LKV2_v01
JOM
23/4/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 Life Cycle 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(7,15,65)
maxKPI=c(15,50,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
## CW31 CW32 CW33 CW34 CW35 CW36 CW37 CW38 CW39 CW40 CW41 CW42 CW43 CW44 CW45
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## CW46 CW47 CW48 CW49 CW50 CW51 CW52
## 1 1 1 1 1 1 1#Ahora tendríamos los valores del sistema ternario normalizado para 52 pasos de tiempo (semanas).
Pars <- c(r1 = .4, a11 = .5, a12 = .2, a13 = .6,
r2 = .2, a21 = .2, a22 = .5, a23 = .3,
r3 = .4, a31 = .4, a32 = .4, a33 = .8)
State <- c(x = dat[1,5], y = dat[1,6], z=dat[1,7])
Time <- seq(1, nrow(dat), by = 1)
A=matrix(NA,nrow=3,ncol = 3)
A[1,1]=Pars["a11"]
A[1,2]=Pars["a12"]
A[1,3]=Pars["a13"]
A[2,1]=Pars["a21"]
A[2,2]=Pars["a22"]
A[2,3]=Pars["a23"]
A[3,1]=Pars["a31"]
A[3,2]=Pars["a32"]
A[3,3]=Pars["a33"]
m11 = as.numeric(Pars["a22"]*Pars["a33"] - Pars["a23"]*Pars["a32"])
m22 = as.numeric(Pars["a11"]*Pars["a33"] - Pars["a13"]*Pars["a31"])
m33 = as.numeric(Pars["a11"]*Pars["a22"] - Pars["a12"]*Pars["a21"])
da = det(A)
if (m11 > 0) { res1='se cumple';} else { res1='no se cumple';}
if (m22 > 0) { res2='se cumple';} else { res2='no se cumple';}
if (m33 > 0) { res3='se cumple';} else { res3='no se cumple';}
if (da > 0) { res4='se cumple ya que es positivo';} else
{ res4='no se cumple, ya que es negativo';}
vizda=A[1,3]*A[2,1]*A[3,2]+A[1,2]*A[2,3]*A[3,1]-2*(A[1,1]*A[2,2]*A[3,3])
vdcha=2*(sqrt(A[1,1]*m11*m22)+sqrt(A[2,2]*A[3,3]*m22*m33)+sqrt(A[1,1]*A[3,3]*m11*m33))
if (vizda < vdcha) {res5='se cumple';} else { res5='ńo se cumple';}
#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:
Esto sitúa el problema en su verdadera magnitud, es decir el sistema modelado NO SIEMPRE ES WL, ya que
## [1] 1.0000000 0.9518622 0.9224710 0.9041662 0.8925872 0.8851203 0.8801567
## [8] 0.8766940 0.8741104 0.8720250 0.8702101 0.8685332 0.8669225 0.8653419
## [15] 0.8637757 0.8622202 0.8606774 0.8591516 0.8576479 0.8561709 0.8547246
## [22] 0.8533117 0.8519344 0.8505936 0.8492900 0.8480235 0.8467936 0.8455996
## [29] 0.8444405 0.8433152 0.8422224 0.8411611 0.8401299 0.8391276 0.8381530
## [36] 0.8372051 0.8362826 0.8353847 0.8345103 0.8336586 0.8328286 0.8320195
## [43] 0.8312307 0.8304614 0.8297110 0.8289789 0.8282644 0.8275671 0.8268865
## [50] 0.8262220 0.8255733 0.8249399como 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
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)
En este caso de Life Cycle Value Stream el comportamiento de los KPIs SI se parece al de tres presas compitiendo por alimento común limitado.
Una vez resuelto el algoritmo genético, tendremos que la mejor solución encontrada proviene del objeto bestSolution de la última iteración. Vamos a sacar los parámetros y valorar las soluciones:
LKVpars=function(x) {
Pars = c(r1 = x[10], a11 = x[1], a12 = x[2], a13 = x[3],
r2 = x[11], a21 = x[4], a22 = x[5], a23 = x[6],
r3 = x[12], a31 = x[7], a32 = x[8], a33 = x[9])
return(Pars)
}
LKVpars(GA1@bestSol[[GA1@maxiter]])## r1 a11 a12 a13 r2 a21 a22
## 0.4720175 0.5498038 0.4797541 0.4597939 0.4839662 0.4667363 0.5450558
## a23 r3 a31 a32 a33
## 0.4882042 0.5286260 0.4715993 0.3962419 0.5666388La 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:
matplot(LKVtime(GA1@bestSol[[GA1@maxiter]])[,-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)
Y la invariancia a la disipación:
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
## 1.0000000 0.9898091 0.9848088 0.9829326 0.9829572 0.9841390 0.9860134
## CW8 CW9 CW10 CW11 CW12 CW13 CW14
## 0.9882867 0.9907693 0.9933313 0.9958857 0.9983699 1.0007393 1.0029619
## CW15 CW16 CW17 CW18 CW19 CW20 CW21
## 1.0050153 1.0068846 1.0085607 1.0100405 1.0113246 1.0124178 1.0133273
## CW22 CW23 CW24 CW25 CW26 CW27 CW28
## 1.0140626 1.0146347 1.0150562 1.0153398 1.0154986 1.0155456 1.0154934
## CW29 CW30 CW31 CW32 CW33 CW34 CW35
## 1.0153538 1.0151383 1.0148572 1.0145202 1.0141360 1.0137126 1.0132569
## CW36 CW37 CW38 CW39 CW40 CW41 CW42
## 1.0127753 1.0122734 1.0117560 1.0112275 1.0106916 1.0101516 1.0096101
## CW43 CW44 CW45 CW46 CW47 CW48 CW49
## 1.0090696 1.0085321 1.0079993 1.0074726 1.0069533 1.0064423 1.0059404
## CW50 CW51 CW52
## 1.0054482 1.0049663 1.0044950La conclusión es obvia, en este caso el error es mejor y el ajuste es mucho más satisfactorio.