LKV_Nemawashi_v02 - NemaBosch
JOM
20/12/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/NemaBosch/")
dat = read.csv(file="Nemawashi 4 KPI - 3 Value Streams.csv",sep=",",
header=TRUE,stringsAsFactors=FALSE)
place = dat[1,]
signf = dat[2,]
dat=dat[-c(1,2),]
units = dat[1,]
dat = dat[-1,]Ahora normalizamos a escala consistente los datos:
#
normaliza_ternario=function(x,minv,maxv){
y=(as.numeric(gsub(',','.',x))-minv)/(maxv-minv)
ty=sum(y)
return(y/ty)
}
LotVmod = function (Time, State, Pars) {
with(as.list(c(State, Pars)), {
dx = x*(r1 -a11*x -a12*y -a13*z)
dy = y*(r2 -a21*x -a22*y -a23*z)
dz = z*(r3 -a31*x -a32*y -a33*z)
return(list(c(dx, dy, dz)))
})
}
LKVcost=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])
State = c(x = dat[1,1], y = dat[1,2], z=dat[1,3])
Time = seq(0, nrow(dat), by = 1)
lres = as.data.frame(ode(func = LotVmod, y = State, parms = Pars, times = Time,
method = "ode45", atol = 1e-14, rtol = 1e-14))
err = 0.
for (i in 1:3) {
err = err + sum(abs(dat[,i]-lres[(1:nrow(dat)),(i+1)]))
}
return(-err)
}
LKVtime=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])
State = c(x = dat[1,1], y = dat[1,2], z=dat[1,3])
Time = seq(0, nrow(dat), by = 1)
lres = as.data.frame(ode(func = LotVmod, y = State, parms = Pars, times = Time))
return(lres)
}
#
# 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)
#
dd=t(apply(dat[,c(5:8,13:16,21:24,26)],1,as.numeric))
dda=t(apply(dat[,c(1:4,9:12,17:20)],1,as.numeric))
oug= as.numeric(dat[,26])-as.numeric(min(dat[,26]))/
(as.numeric(max(dat[,26]))-as.numeric(min(dat[,26])))
#
pca=list()
pca[[1]]=prcomp(as.matrix(dd[,1:4]))
pca[[2]]=prcomp(as.matrix(dd[,5:8]))
pca[[3]]=prcomp(as.matrix(dd[,9:12]))
pca[[4]]=prcomp(as.matrix(dd[,c(1:4,5:8,9:13)]))
#
evl=function(x){
return(cumsum(summary(x)$importance[2,])[3]*100)
}
lapply(pca,evl)## [[1]]
## PC3
## 99.212
##
## [[2]]
## PC3
## 99.042
##
## [[3]]
## PC3
## 97.288
##
## [[4]]
## PC3
## 92.43lapply(pca,print)## Standard deviations:
## [1] 0.39975377 0.21913359 0.15168185 0.04281096
##
## Rotation:
## PC1 PC2 PC3 PC4
## [1,] -0.4681851 0.0377319 0.8685747 -0.157977914
## [2,] 0.3795438 0.9102706 0.1653795 0.001860239
## [3,] 0.6821892 -0.3692521 0.4619605 0.429957296
## [4,] -0.4139650 0.1834027 -0.0694274 0.888918580
## Standard deviations:
## [1] 0.34522807 0.13656164 0.06456002 0.03704541
##
## Rotation:
## PC1 PC2 PC3 PC4
## [1,] -0.58426076 0.05762939 -0.68446612 -0.43223182
## [2,] 0.05834482 0.99010238 0.09044107 -0.09007536
## [3,] -0.51830179 0.11466062 -0.08066429 0.84362875
## [4,] -0.62177046 -0.05682491 0.71890181 -0.30553660
## Standard deviations:
## [1] 0.2701962 0.1137074 0.0544836 0.0497862
##
## Rotation:
## PC1 PC2 PC3 PC4
## [1,] -0.6544042 0.2275442 -0.12446424 0.71027283
## [2,] -0.3675108 -0.7860189 -0.46758967 -0.16873054
## [3,] -0.6006670 0.4152808 0.02043567 -0.68287874
## [4,] -0.2754817 -0.3974163 0.87490053 0.02681676
## Standard deviations:
## [1] 0.62348410 0.24824563 0.18344142 0.13782297 0.10010829 0.05765288
## [7] 0.05256413 0.03676924 0.03290518 0.02922468 0.02720799 0.01881352
## [13] 0.01235712
##
## Rotation:
## PC1 PC2 PC3 PC4 PC5
## [1,] -0.31220861 0.07169975 -0.577763666 0.07467899 -0.453203295
## [2,] 0.19043872 -0.86443420 -0.221352339 -0.25105891 0.168233310
## [3,] 0.43166124 0.14717952 -0.524232691 0.32331554 0.307313916
## [4,] -0.26830593 -0.08229789 0.174986507 0.10939675 -0.126589939
## [5,] -0.30643929 0.06424780 -0.148265565 0.34018684 0.536982350
## [6,] 0.03309702 -0.34848664 0.015869504 0.68313578 -0.383217139
## [7,] -0.28026906 -0.01832456 0.030299118 0.22282164 0.298588121
## [8,] -0.34458109 -0.06665805 -0.116848937 -0.10564743 0.199029368
## [9,] -0.28243589 -0.12224415 0.058870492 -0.16878880 -0.005678861
## [10,] -0.14169594 0.07684302 -0.463280646 -0.25431092 -0.223606214
## [11,] -0.24898585 -0.25550245 0.009386156 0.18306700 0.046670682
## [12,] -0.11062084 0.01407082 -0.234723984 -0.21129658 0.197166427
## [13,] -0.37626235 -0.05000119 0.040837817 -0.01886710 0.027924110
## PC6 PC7 PC8 PC9 PC10
## [1,] 0.41006136 -0.30328986 0.179129681 -0.13621160 -0.066634767
## [2,] 0.19217417 -0.03361464 0.118800874 0.06042700 -0.068402471
## [3,] 0.05627389 0.12984645 -0.170313455 0.24836688 0.177166487
## [4,] 0.20436874 0.14527847 0.254573093 0.06710746 0.501039884
## [5,] 0.03447866 0.04956922 0.449158173 0.06008248 -0.007406178
## [6,] -0.31658786 -0.11493292 -0.235768526 0.03339346 -0.065134445
## [7,] 0.21811733 0.10472194 -0.324720638 -0.30821391 -0.549191947
## [8,] 0.09064198 -0.15278563 -0.654787791 -0.01524677 0.520932489
## [9,] 0.08581609 0.23411242 -0.170426449 0.06364050 -0.173895883
## [10,] -0.33469985 0.60289072 -0.064227295 0.06766262 -0.098888699
## [11,] -0.30420383 0.33617502 0.162960545 -0.30641246 0.234705290
## [12,] -0.58714749 -0.50308316 0.093256121 -0.28632580 0.004337884
## [13,] -0.19001871 -0.18587738 0.006419243 0.79192976 -0.189863206
## PC11 PC12 PC13
## [1,] 0.07986004 -0.14956378 -0.09295514
## [2,] -0.09372545 0.02091098 0.08018833
## [3,] -0.01586296 0.05646726 -0.41552893
## [4,] -0.44556132 0.47602596 -0.23297250
## [5,] 0.25217018 0.13348870 0.43063283
## [6,] 0.08919602 0.22430144 0.17615620
## [7,] -0.44614394 0.03815710 -0.14247013
## [8,] 0.05430553 -0.09490552 0.26177680
## [9,] 0.65115273 0.40716377 -0.39750267
## [10,] -0.20494131 0.14452186 0.29638420
## [11,] 0.08784281 -0.57011129 -0.35087975
## [12,] -0.10091245 0.31167217 -0.23041951
## [13,] -0.16917934 -0.24489879 -0.17500430## [[1]]
## Standard deviations:
## [1] 0.39975377 0.21913359 0.15168185 0.04281096
##
## Rotation:
## PC1 PC2 PC3 PC4
## [1,] -0.4681851 0.0377319 0.8685747 -0.157977914
## [2,] 0.3795438 0.9102706 0.1653795 0.001860239
## [3,] 0.6821892 -0.3692521 0.4619605 0.429957296
## [4,] -0.4139650 0.1834027 -0.0694274 0.888918580
##
## [[2]]
## Standard deviations:
## [1] 0.34522807 0.13656164 0.06456002 0.03704541
##
## Rotation:
## PC1 PC2 PC3 PC4
## [1,] -0.58426076 0.05762939 -0.68446612 -0.43223182
## [2,] 0.05834482 0.99010238 0.09044107 -0.09007536
## [3,] -0.51830179 0.11466062 -0.08066429 0.84362875
## [4,] -0.62177046 -0.05682491 0.71890181 -0.30553660
##
## [[3]]
## Standard deviations:
## [1] 0.2701962 0.1137074 0.0544836 0.0497862
##
## Rotation:
## PC1 PC2 PC3 PC4
## [1,] -0.6544042 0.2275442 -0.12446424 0.71027283
## [2,] -0.3675108 -0.7860189 -0.46758967 -0.16873054
## [3,] -0.6006670 0.4152808 0.02043567 -0.68287874
## [4,] -0.2754817 -0.3974163 0.87490053 0.02681676
##
## [[4]]
## Standard deviations:
## [1] 0.62348410 0.24824563 0.18344142 0.13782297 0.10010829 0.05765288
## [7] 0.05256413 0.03676924 0.03290518 0.02922468 0.02720799 0.01881352
## [13] 0.01235712
##
## Rotation:
## PC1 PC2 PC3 PC4 PC5
## [1,] -0.31220861 0.07169975 -0.577763666 0.07467899 -0.453203295
## [2,] 0.19043872 -0.86443420 -0.221352339 -0.25105891 0.168233310
## [3,] 0.43166124 0.14717952 -0.524232691 0.32331554 0.307313916
## [4,] -0.26830593 -0.08229789 0.174986507 0.10939675 -0.126589939
## [5,] -0.30643929 0.06424780 -0.148265565 0.34018684 0.536982350
## [6,] 0.03309702 -0.34848664 0.015869504 0.68313578 -0.383217139
## [7,] -0.28026906 -0.01832456 0.030299118 0.22282164 0.298588121
## [8,] -0.34458109 -0.06665805 -0.116848937 -0.10564743 0.199029368
## [9,] -0.28243589 -0.12224415 0.058870492 -0.16878880 -0.005678861
## [10,] -0.14169594 0.07684302 -0.463280646 -0.25431092 -0.223606214
## [11,] -0.24898585 -0.25550245 0.009386156 0.18306700 0.046670682
## [12,] -0.11062084 0.01407082 -0.234723984 -0.21129658 0.197166427
## [13,] -0.37626235 -0.05000119 0.040837817 -0.01886710 0.027924110
## PC6 PC7 PC8 PC9 PC10
## [1,] 0.41006136 -0.30328986 0.179129681 -0.13621160 -0.066634767
## [2,] 0.19217417 -0.03361464 0.118800874 0.06042700 -0.068402471
## [3,] 0.05627389 0.12984645 -0.170313455 0.24836688 0.177166487
## [4,] 0.20436874 0.14527847 0.254573093 0.06710746 0.501039884
## [5,] 0.03447866 0.04956922 0.449158173 0.06008248 -0.007406178
## [6,] -0.31658786 -0.11493292 -0.235768526 0.03339346 -0.065134445
## [7,] 0.21811733 0.10472194 -0.324720638 -0.30821391 -0.549191947
## [8,] 0.09064198 -0.15278563 -0.654787791 -0.01524677 0.520932489
## [9,] 0.08581609 0.23411242 -0.170426449 0.06364050 -0.173895883
## [10,] -0.33469985 0.60289072 -0.064227295 0.06766262 -0.098888699
## [11,] -0.30420383 0.33617502 0.162960545 -0.30641246 0.234705290
## [12,] -0.58714749 -0.50308316 0.093256121 -0.28632580 0.004337884
## [13,] -0.19001871 -0.18587738 0.006419243 0.79192976 -0.189863206
## PC11 PC12 PC13
## [1,] 0.07986004 -0.14956378 -0.09295514
## [2,] -0.09372545 0.02091098 0.08018833
## [3,] -0.01586296 0.05646726 -0.41552893
## [4,] -0.44556132 0.47602596 -0.23297250
## [5,] 0.25217018 0.13348870 0.43063283
## [6,] 0.08919602 0.22430144 0.17615620
## [7,] -0.44614394 0.03815710 -0.14247013
## [8,] 0.05430553 -0.09490552 0.26177680
## [9,] 0.65115273 0.40716377 -0.39750267
## [10,] -0.20494131 0.14452186 0.29638420
## [11,] 0.08784281 -0.57011129 -0.35087975
## [12,] -0.10091245 0.31167217 -0.23041951
## [13,] -0.16917934 -0.24489879 -0.17500430lapply(pca,summary)## [[1]]
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 0.3998 0.2191 0.15168 0.04281
## Proportion of Variance 0.6868 0.2064 0.09889 0.00788
## Cumulative Proportion 0.6868 0.8932 0.99212 1.00000
##
## [[2]]
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 0.3452 0.1366 0.06456 0.03705
## Proportion of Variance 0.8313 0.1301 0.02907 0.00957
## Cumulative Proportion 0.8313 0.9614 0.99043 1.00000
##
## [[3]]
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 0.2702 0.1137 0.05448 0.04979
## Proportion of Variance 0.7989 0.1415 0.03248 0.02712
## Cumulative Proportion 0.7989 0.9404 0.97288 1.00000
##
## [[4]]
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 0.6235 0.2482 0.18344 0.13782 0.10011 0.05765
## Proportion of Variance 0.7423 0.1177 0.06426 0.03627 0.01914 0.00635
## Cumulative Proportion 0.7423 0.8600 0.92430 0.96058 0.97971 0.98606
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.05256 0.03677 0.03291 0.02922 0.02721 0.01881
## Proportion of Variance 0.00528 0.00258 0.00207 0.00163 0.00141 0.00068
## Cumulative Proportion 0.99134 0.99392 0.99599 0.99762 0.99903 0.99971
## PC13
## Standard deviation 0.01236
## Proportion of Variance 0.00029
## Cumulative Proportion 1.00000dpca=lapply(pca,function(x){x$x[,1:3]})
#Ahora tendríamos los valores del sistema ternario normalizado para 22 pasos de tiempo (semanas).
#
ipca = 1
ddn=t(apply(dpca[[ipca]],1,normaliza_ternario,
apply(dpca[[ipca]],2,min,na.rm=TRUE),
apply(dpca[[ipca]],2,max,na.rm=TRUE)))
sdat=data.frame(ddn,oug)
panel.hist <- function(x, ...) {
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y,
col = "cyan", ...)
}
## put (absolute) correlations on the upper panels,
## with size proportional to the correlations.
panel.cor <- function(x, y, digits = 2,
prefix = "", cex.cor, ...){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(sdat,panel=panel.smooth,
upper.panel=panel.cor, diag.panel=panel.hist)
dat=ddn
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,1], y = dat[1,2], z=dat[1,3])
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';}
#Caso 1 de los 4 PCAs
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.8084221 0.7170726 0.6668084 0.6371229 0.6190097 0.6079033
## [8] 0.6012612 0.5975745 0.5959084 0.5956669 0.5964628 0.5980419 0.6002364
## [15] 0.6029352 0.6060646 0.6095754 0.6134339 0.6176160 0.6221027 0.6268773
## [22] 0.6319229como 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(1:3)], type = "l", xlab = "Time", ylab = "Relative value per unit")
legend("topright", c("PcaKPI1", "PcaKPI2","PcaKPI3"), 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.5619450 0.7082886 0.4262639 0.2377784 0.4491938 0.4632352 0.5602191
## a23 r3 a31 a32 a33
## 0.4717402 0.4155801 0.1021375 0.4607434 0.9853695La 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)## 0.0223632269132467 0.0170586298420778 3.49479276632197e-16
## 1.0000000 0.7817870 0.7085941
## 0.195323416815065 0.311330202563056 0.353024591056636
## 0.6892486 0.6972868 0.7216161
## 0.439702511037452 0.810417161453791 0.468119791702424
## 0.7559483 0.7956517 0.8368067
## 0.337015471837637 0.34335173137218 0.365376379936351
## 0.8761388 0.9112486 0.9407884
## 0.397030245403745 0.390224630395006 0.453844057870419
## 0.9644109 0.9825164 0.9959268
## 0.509395010798617 0.574307254827581 0.572944829183242
## 1.0056017 1.0124510 1.0172443
## 0.527326949425464 0.669499930665261 0.525628446830254
## 1.0205843 1.0229191 1.0245692
## 0.467823164851314
## 1.0257572La conclusión es obvia, en este caso el error es mejor y el ajuste es mucho más satisfactorio.
# Gráfico ternario
r 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),"PcaKPI-3", pos=3) text(0,0,"PcaKPI-1", pos=1) text(1,0,"PcaKPI-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,22,4),sep=""),text.col=rainbow(22)[seq(1,22,4)],pt.bg=rainbow(22)[seq(1,22,4)],pch=rep(21,25), inset=c(0.1,0.0))
r # save.image(file=paste("~/git/JVD/NemaBosch/all_",ipca,"_GA_2015.RData",sep="")) #
Ahora tendríamos los valores del sistema ternario normalizado para 22 pasos de tiempo (semanas).
#
ipca = 2
ddn=t(apply(dpca[[ipca]],1,normaliza_ternario,
apply(dpca[[ipca]],2,min,na.rm=TRUE),
apply(dpca[[ipca]],2,max,na.rm=TRUE)))
sdat=data.frame(ddn,oug)
panel.hist <- function(x, ...) {
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y,
col = "cyan", ...)
}
## put (absolute) correlations on the upper panels,
## with size proportional to the correlations.
panel.cor <- function(x, y, digits = 2,
prefix = "", cex.cor, ...){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(sdat,panel=panel.smooth,
upper.panel=panel.cor, diag.panel=panel.hist)
dat=ddn
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,1], y = dat[1,2], z=dat[1,3])
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';}
#Caso 2 de los 4 PCAs
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.7862229 0.6691112 0.5963803 0.5476432 0.5132994 0.4882315
## [8] 0.4694605 0.4551340 0.4440395 0.4353510 0.4284866 0.4230256 0.4186571
## [15] 0.4151471 0.4123168 0.4100282 0.4081732 0.4066670 0.4054420 0.4044446
## [22] 0.4036316como 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(1:3)], type = "l", xlab = "Time", ylab = "Relative value per unit")
legend("topright", c("PcaKPI1", "PcaKPI2","PcaKPI3"), 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.5331793 0.5394517 0.5791342 0.5410247 0.2668811 0.4690450 0.9979227
## a23 r3 a31 a32 a33
## 0.4888278 0.5588858 0.4919228 0.4916377 0.4507969La 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)## -7.48733836237751e-16 0.111196964926142 0.164772086334952
## 1.0000000 0.6091553 0.4688357
## 0.246915075346115 0.350362484914679 0.344136381126356
## 0.3985360 0.3574881 0.3313539
## 0.534984328489872 0.709678701349658 0.525725975986138
## 0.3137875 0.3015461 0.2927989
## 0.500255379693831 0.423805536804522 0.461056830560283
## 0.2864361 0.2817479 0.2782601
## 0.412837659626307 0.45740761937979 0.403492537444766
## 0.2756471 0.2736790 0.2721909
## 0.441821591128845 0.484224505937143 0.649132856439195
## 0.2710623 0.2702042 0.2695509
## 0.552489826318026 0.541807477864439 0.532222303232553
## 0.2690527 0.2686724 0.2683819
## 0.533915225521476
## 0.2681599La conclusión es obvia, en este caso el error es mejor y el ajuste es mucho más satisfactorio.
# Gráfico ternario
r 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),"PcaKPI-3", pos=3) text(0,0,"PcaKPI-1", pos=1) text(1,0,"PcaKPI-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,22,4),sep=""),text.col=rainbow(22)[seq(1,22,4)],pt.bg=rainbow(22)[seq(1,22,4)],pch=rep(21,25), inset=c(0.1,0.0))
r # save.image(file=paste("~/git/JVD/NemaBosch/all_",ipca,"_GA_2015.RData",sep="")) #
Ahora tendríamos los valores del sistema ternario normalizado para 22 pasos de tiempo (semanas).
#
ipca = 3
ddn=t(apply(dpca[[ipca]],1,normaliza_ternario,
apply(dpca[[ipca]],2,min,na.rm=TRUE),
apply(dpca[[ipca]],2,max,na.rm=TRUE)))
sdat=data.frame(ddn,oug)
panel.hist <- function(x, ...) {
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y,
col = "cyan", ...)
}
## put (absolute) correlations on the upper panels,
## with size proportional to the correlations.
panel.cor <- function(x, y, digits = 2,
prefix = "", cex.cor, ...){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(sdat,panel=panel.smooth,
upper.panel=panel.cor, diag.panel=panel.hist)
dat=ddn
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,1], y = dat[1,2], z=dat[1,3])
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';}
#Caso 3 de los 4 PCAs
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.8130963 0.7158665 0.6631373 0.6371209 0.6294602 0.6355493
## [8] 0.6522549 0.6768699 0.7066732 0.7388872 0.7708936 0.8005379 0.8263579
## [15] 0.8476413 0.8643183 0.8767646 0.8855947 0.8914971 0.8951291 0.8970616
## [22] 0.8977612como 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(1:3)], type = "l", xlab = "Time", ylab = "Relative value per unit")
legend("topright", c("PcaKPI1", "PcaKPI2","PcaKPI3"), 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
## 0.52437123 0.98114180 0.05572541 0.47365783 0.30757553 0.05203012
## a22 a23 r3 a31 a32 a33
## 0.80762160 0.51559053 0.49844031 0.44529521 0.47072948 0.47263139La 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)## 0.0468053919177978 -2.90225301393846e-16 0.120024189147987
## 1.0000000 0.7517233 0.6679435
## 0.220677856678109 0.227465493092571 0.263588930044604
## 0.6502745 0.6679113 0.7040360
## 0.360435066874517 0.462715623050949 0.585928111518639
## 0.7458157 0.7840091 0.8139410
## 0.52677788478258 0.527979922312596 0.655510499301836
## 0.8348916 0.8483354 0.8563575
## 0.539279410664284 0.505805999179368 0.430597909107428
## 0.8608192 0.8630941 0.8640984
## 0.40886428197891 0.405758875114255 0.444926916974781
## 0.8644050 0.8643555 0.8641449
## 0.480564798223731 0.50972973992706 0.514250972276234
## 0.8638798 0.8636151 0.8633756
## 0.440398807654523
## 0.8631707La conclusión es obvia, en este caso el error es mejor y el ajuste es mucho más satisfactorio.
# Gráfico ternario
r 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),"PcaKPI-3", pos=3) text(0,0,"PcaKPI-1", pos=1) text(1,0,"PcaKPI-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,22,4),sep=""),text.col=rainbow(22)[seq(1,22,4)],pt.bg=rainbow(22)[seq(1,22,4)],pch=rep(21,25), inset=c(0.1,0.0))
r # save.image(file=paste("~/git/JVD/NemaBosch/all_",ipca,"_GA_2015.RData",sep="")) #
Ahora tendríamos los valores del sistema ternario normalizado para 22 pasos de tiempo (semanas).
#
ipca = 4
ddn=t(apply(dpca[[ipca]],1,normaliza_ternario,
apply(dpca[[ipca]],2,min,na.rm=TRUE),
apply(dpca[[ipca]],2,max,na.rm=TRUE)))
sdat=data.frame(ddn,oug)
panel.hist <- function(x, ...) {
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y,
col = "cyan", ...)
}
## put (absolute) correlations on the upper panels,
## with size proportional to the correlations.
panel.cor <- function(x, y, digits = 2,
prefix = "", cex.cor, ...){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(sdat,panel=panel.smooth,
upper.panel=panel.cor, diag.panel=panel.hist)
dat=ddn
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,1], y = dat[1,2], z=dat[1,3])
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';}
#Caso 4 de los 4 PCAs
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.7867793 0.6700218 0.5976156 0.5492345 0.5153131 0.4907633
## [8] 0.4726368 0.4591170 0.4490342 0.4416141 0.4363376 0.4328590 0.4309560
## [15] 0.4304974 0.4314198 0.4337102 0.4373887 0.4424908 0.4490465 0.4570569
## [22] 0.4664688como 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(1:3)], type = "l", xlab = "Time", ylab = "Relative value per unit")
legend("topright", c("PcaKPI1", "PcaKPI2","PcaKPI3"), 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.4296417 0.4647824 0.4339652 0.4306910 0.2399285 0.4561576 0.9816484
## a23 r3 a31 a32 a33
## 0.1401166 0.4484962 0.5073900 0.2019588 0.9718218La 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)## 4.19202927313733e-16 0.0111703149498767 0.127213791394191
## 1.0000000 0.6041920 0.4617283
## 0.345624188928509 0.39668464793221 0.38454557095388
## 0.3906118 0.3497923 0.3250231
## 0.299655337992101 0.233273124428604 0.301303798832656
## 0.3102583 0.3027383 0.3013031
## 0.429340694691289 0.608273326458586 0.661551004813404
## 0.3056924 0.3161559 0.3331130
## 0.657187699770924 0.687130744198574 0.557616544667311
## 0.3567433 0.3865161 0.4208349
## 0.522075346622979 0.451352933842529 0.456419573040862
## 0.4570757 0.4921607 0.5234102
## 0.484649798141984 0.399936723982161 0.478589776141675
## 0.5492091 0.5691412 0.5836956
## 0.521507544587344
## 0.5938223La conclusión es obvia, en este caso el error es mejor y el ajuste es mucho más satisfactorio.
# Gráfico ternario
r 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),"PcaKPI-3", pos=3) text(0,0,"PcaKPI-1", pos=1) text(1,0,"PcaKPI-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,22,4),sep=""),text.col=rainbow(22)[seq(1,22,4)],pt.bg=rainbow(22)[seq(1,22,4)],pch=rep(21,25), inset=c(0.1,0.0))
r # save.image(file=paste("~/git/JVD/NemaBosch/all_",ipca,"_GA_2015.RData",sep="")) #