LKV_PdPsPs_v01
JOM
02/5/2015
Introducción
Se pretende ajustar un modelo de LKV a los datos de un sistema compuesto de varios POs que interaccionan entre si.
Lo primero será leer los datos y después definir la estrategia.
setwd("~/git/JVD_HEALTH/")
dat=read.csv(file="Case Study Production Value Stream-Table 1.csv",sep=",",
header=TRUE,stringsAsFactors=FALSE)
units=dat[1,]
dat=dat[-1,]Ahora normalizamos a escala consistente los datos:
#
minKPI=c(0.2,50,15)
maxKPI=c(2,100,100)
#
normaliza_ternario=function(x,minv,maxv){
y=(as.numeric(gsub(',','.',x))-minv)/(maxv-minv)
ty=sum(y)
return(y/ty)
}
#
dd= t(apply(dat[,2:4],1,normaliza_ternario,minKPI,maxKPI))
dat[,5:7]=dd
rownames(dat)=dat[,1]
apply(dat[,5:7],1,sum)## CW1 CW2 CW3 CW4 CW5 CW6 CW7 CW8 CW9 CW10 CW11 CW12 CW13 CW14 CW15
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## CW16 CW17 CW18 CW19 CW20 CW21 CW22 CW23 CW24 CW25 CW26 CW27 CW28 CW29 CW30
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1#Ahora tendríamos los valores del sistema ternario normalizado para 30 pasos de tiempo (semanas).
Las condiciones de estabilidad para m11 se cumple [0.28]; para m22 se cumple [0.16] y para m33 se cumple [0.21]. La otra condición exigida, el determinante de A se cumple ya que es positivo [0.06]. Adicionalmente la evaluación de la última condición indica que ésta se cumple.
Procedemos a resolver el sistema:
apply(out[,-1],1,sum)## [1] 1.0000000 0.8293398 0.7278250 0.6679559 0.6335942 0.6145552 0.6044225
## [8] 0.5992584 0.5967248 0.5954864 0.5948134 0.5943292 0.5938560 0.5933224
## [15] 0.5927112 0.5920298 0.5912954 0.5905267 0.5897403 0.5889503 0.5881677
## [22] 0.5874003 0.5866541 0.5859326 0.5852382 0.5845723 0.5839353 0.5833273
## [29] 0.5827477 0.5821958Como se ve, no siempre existe reparto entre los KPIs, sino que al ser todos “presas” crecen hasta el agotamiento de los rescursos disponibles, caracterizados por los términos -aii*var^2. Es decir el caracter disipativo dependerá de los coeficientes.
El impacto es que los términos de reparto entre los KPIi no siempre podrán ser presentados en un modelo ternario, dependiendo de los valores de esos parámetros.
Un criterio indirecto de calidad para el genético es observar si la solución encontrada es conservativa, que será un requisito de mínimos.
Búsqueda de Parámetros
Una vez que tenemos el sistema de solución de las ecuaciones operativo vamos a establecer el procedimiento para ese ajuste. Se empleará una técnica de algoritmo genético donde la función de coste tenga que ver, para un vector de coeficientes (12) dado, una vez resueltas las ecuaciones de LKV, se evalúa la distancia en cada paso de tiempo entre el objetivo y la solcuión actual. Esa distancia será sumada para cada uno de los KPI’s. En el algoritmo genético se codificarán como números reales cada uno de los parámetros de los que depende el sistema a ser resuelto mediante LKV. En caso de que alguno sea negativo la función de coste adoptará el valor de -100000
SPCmod <- function(t, x, parms, input) {
with(as.list(c(parms, x)), {
import <- input(t)
dS <- import - b*S*P + g*C # substrate
dP <- c*S*P - d*C*P # producer
dC <- e*P*C - f*C # consumer
res <- c(dS, dP, dC)
list(res)
})
}
signal = data.frame(times = Time,
import = rep(0.2, length(Time)))
signal$import = sample(c(0.3,0.4,0.35,0.325,0.375,0.25,0.45),length(Time),replace=TRUE)
sigimp <- approxfun(signal$times, signal$import, rule = 2)
LKVcost=function(x) {
xstart <- c(S = dat[1,5], P = dat[1,6], C = dat[1,7])
parms <- c(b = x[1], c = x[2], d = x[3], e = x[4], f = x[5], g = x[6])
lres = as.data.frame(ode(func = SPCmod, y = xstart, parms = parms, times = Time,
input=sigimp))
if (nrow(lres) != length(Time) || is.na(lres)) {
return (-100000)
}
err = 0.
for (i in 1:3) {
err = err + sum(abs(dat[,(i+4)]-lres[1:nrow(dat),(i+1)]))
}
return(-err)
}
LKVtime=function(x){
Pars <- c(b = x[1], c = x[2], d = x[3], e = x[4], f = x[5], g = x[6])
xstart <- c(S = dat[1,5], P = dat[1,6], C = dat[1,7])
lres = as.data.frame(ode(func = SPCmod, y = xstart, parms = Pars, times = Time,
input=sigimp))
return(lres)
}
#
NITER = 5
tt=list()
GAs = list()
if ( file.exists("~/git/JVD_HEALTH/GA_2015.RData")) {
j=1
while(j <= NITER) {
load(file=paste("~/git/JVD_HEALTH/GA_2015_LKV_",sprintf("%02d",j),".RData",sep=""))
tt[[j]]=t_tot
GAs[[j]]=GA1
signal$import = sample(c(0.3,0.4,0.35,0.325,0.375,0.25,0.45),length(Time),replace=TRUE)
sigimp <- approxfun(signal$times, signal$import, rule = 2)
j=j+1
}
} else {
j=1
while(j <= NITER) {
cat(paste("Starting Iteration:",j,"<br>",sep=""))
flush.console()
tt[[j]]= system.time({GA1 = ga(type="real-valued", fitness = LKVcost,
min=rep(0,6),max=rep(1.,6),
monitor = FALSE,parallel=3,maxiter=400,
popSize=8000,keepBest=TRUE)})
t_tot=tt[[j]]
GAs[[j]]=GA1
save(GA1,t_tot,j,file=paste("~/git/JVD_HEALTH/GA_2015_LKV_",
sprintf("%02d",j),".RData",sep=""))
j=j+1
}
save(GAs,tt,NITER,signal,sigimp,file="~/git/JVD_HEALTH/GA_2015.RData")
}Valoración de la solución
El comportamiento original de los KPIs normalziados es
matplot(dat[,c(5,6,7)], type = "l", xlab = "Time", ylab = "Relative value per unit")
legend("topright", c("KPI1", "KPI2","KPI3"), lty = c(1,2,3), col = c(1,2,3), box.lwd = 0)
La mejor solución proviene del objeto bestSolution de la última iteración. Vamos a sacar los parámetros y valorar las soluciones:
LKVpars=function(x) {
Pars <- c(b = x[1], c = x[2], d = x[3], e = x[4], f = x[5], g = x[6])
return(Pars)
}
LKVpars(GA1@bestSol[[GA1@maxiter]])## b c d e f g
## 0.999375287 0.435201102 0.760433457 0.731432921 0.410940666 0.002933935La evolución del error en la población de parámetros puede observarse en la figura siguiente:
matplot(GA1@summary[,c(1,3)], type = "l", xlab = "Iterations", ylab = "Error",ylim=c(-30,0))
legend("bottomright", c("min", "mean"), lty = c(1,2), col = c(1,2), box.lwd = 0)
Para los parámetros identificados la solución del LKV es:
for (i in 1:length(GAs)) {
cat(paste("Pintando resultados de la iteración: ",i,"\n",sep=""))
print( matplot(LKVtime(GAs[[i]]@bestSol[[GAs[[i]]@maxiter]])[1:nrow(dat),-1],
type = "l", xlab = "Time", ylab = "LKV relative value per unit"))
legend("topright", c("KPI1", "KPI2","KPI3"), lty = c(1,2,3),
col = c(1,2,3), box.lwd = 0)
}## Pintando resultados de la iteración: 1
## NULL
## Pintando resultados de la iteración: 2
## NULL
## Pintando resultados de la iteración: 3
## NULL
## Pintando resultados de la iteración: 4
## NULL
## Pintando resultados de la iteración: 5
## NULLY la invariancia a la disipación:
LKVpars(GA1@bestSol[[GA1@maxiter]])## b c d e f g
## 0.999375287 0.435201102 0.760433457 0.731432921 0.410940666 0.002933935#
dd=LKVtime(GA1@bestSol[[GA1@maxiter]])[1:nrow(dat),-1]
rownames(dd)=dat[,1]
apply(dd,1,sum)## CW1 CW2 CW3 CW4 CW5 CW6 CW7 CW8
## 1.000000 1.132636 1.287840 1.429939 1.534880 1.599203 1.655173 1.621646
## CW9 CW10 CW11 CW12 CW13 CW14 CW15 CW16
## 1.582429 1.603526 1.645202 1.696711 1.662402 1.576890 1.550026 1.554994
## CW17 CW18 CW19 CW20 CW21 CW22 CW23 CW24
## 1.600267 1.588379 1.496361 1.519136 1.573080 1.558388 1.535863 1.529972
## CW25 CW26 CW27 CW28 CW29 CW30
## 1.546265 1.511155 1.435078 1.379117 1.360558 1.391621Gráfico ternario
plot(NA,NA,xlim=c(0,1),ylim=c(0,sqrt(3)/2),asp=1,bty="n",axes=F,xlab="",ylab="")
segments(0,0,0.5,sqrt(3)/2)
segments(0.5,sqrt(3)/2,1,0)
segments(1,0,0,0)
text(0.5,(sqrt(3)/2),"KPI-3", pos=3)
text(0,0,"KPI-1", pos=1)
text(1,0,"KPI-2", pos=1)
#
tern2cart <- function(coord){
coord[1]->x
coord[2]->y
coord[3]->z
x+y+z -> tot
x/tot -> x # First normalize the values of x, y and z
y/tot -> y
z/tot -> z
(2*y + z)/(2*(x+y+z)) -> x1 # Then transform into cartesian coordinates
sqrt(3)*z/(2*(x+y+z)) -> y1
return(c(x1,y1))
}
# Apply this equation to each set of coordinates
t(apply(dd,1,tern2cart)) -> tern
points(tern,pch=19,col=rainbow(nrow(dd)),cex=1:nrow(dd)/20)
legend("topright",cex=0.8,bty="n",legend=paste("CW",seq(1,52,4),sep=""),text.col=rainbow(52)[seq(1,52,4)],pt.bg=rainbow(52)[seq(1,52,4)],pch=rep(21,25), inset=c(0.1,0.0))