LKV_Nemawashi_v03_hospitales

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/HOSP/")
dat   = read.csv2(file="Hospitals_for_HKF.csv",sep=";",dec=",",
             header=TRUE,stringsAsFactors=FALSE)
#
normaliza_ternario=function(x,lims){
  y=(as.numeric(gsub(',','.',x))-lims[1])/(lims[2]-lims[1])
  return(y)
}
#
units = dat[1,]
dat   = dat[-1,]
dat[,3] = as.numeric(gsub(",",".",dat[,3]))
dat[,4] = as.numeric(gsub(",",".",dat[,4]))
dat[,5] = as.numeric(gsub(",",".",dat[,5]))
idx=c(1,13,25,nrow(dat)+1)
rg = list()
dd = list()
for (i in 1:(length(idx)-1)) {
  dd[[i]] = dat[idx[i]:(idx[i]-1+diff(idx)[i]),3:5]  
  rownames(dd[[i]]) = dat[idx[i]:(idx[i]-1+diff(idx)[i]),2]
  rg[[i]] = apply(dd[[i]],2,range)
}
ddn= dd
ddx=ddn
for (i in 1:(length(idx)-1)) {
  for (j in 1:ncol(dd[[i]])) {
    ddn[[i]][,j] = normaliza_ternario(dd[[i]][,j],rg[[i]][,j])
  }
  ddx[[i]] = t(apply(ddn[[i]],1,function(x){return(x/sum(x))}))
}
#
print(xtable(dat),type="html")

	Hospital	X	PLT	OSR	R	Number.of.Nodes.connected.through..CPD.nA	Number.of..CPD.nA.active.in.the.V.t.	Number.of.new..CPD.nA	Number.of.new..CPD.nA.within.existing.POs	Number.of.new..CPD.nA.with.new.POs	Number.of.re.wired..CPD.nA	m.t..Number.of.new.activated.nodes.in.V.t.	p.t.	q.t.	X1.p.t..q.t.	LAMBDA	Diameter.of.HKT..Avg..Distance.between.nodes.in.the..CPD.nA.Structural.Network	Condition.for.Scale.Free.d.ln.ln.N…if.2.LAMBDA.3	Condition.for.SWN.2.λ.3
2	H1	Q1	26.00	9300.00	19.70	5	6	6	6	0	0	5.00	1.00	0.00	0.00	3.00	2.00	0.48	1.61
3		Q2	24.00	8733.00	20.60	32	69	52	48	4	6	27.00	0.69	0.09	0.22	2.80	4.20	1.24	3.47
4		Q3	23.00	7600.00	17.60	62	109	41	40	1	29	30.00	0.37	0.27	0.36	2.50	9.40	1.42	4.13
5		Q4	19.00	5000.00	16.70	81	189	76	72	4	65	19.00	0.38	0.34	0.27	2.30	8.20	1.48	4.39
6		Q5	18.00	3150.00	17.20	123	232	31	30	1	31	42.00	0.13	0.13	0.74	2.80	7.30	1.57	4.81
7		Q6	24.20	4300.00	19.80	134	310	51	49	2	49	11.00	0.16	0.16	0.69	2.70	6.10	1.59	4.90
8		Q7	14.00	2100.00	19.00	145	395	76	72	4	41	11.00	0.18	0.10	0.72	2.90	5.40	1.60	4.98
9		Q8	15.20	2700.00	23.80	152	431	41	39	2	32	7.00	0.09	0.07	0.84	3.00	5.01	1.61	5.02
10		Q9	13.00	2500.00	27.00	183	521	31	30	1	87	31.00	0.06	0.17	0.78	2.70	4.20	1.65	5.21
11		Q10	12.00	2400.00	25.00	211	621	29	25	4	101	28.00	0.04	0.16	0.80	2.70	3.70	1.68	5.35
12		Q11	10.50	2100.00	31.00	317	714	34	32	2	86	106.00	0.05	0.12	0.83	2.80	2.90	1.75	5.76
13		Q12	9.20	1900.00	33.00	348	813	87	84	3	6	178.00	0.10	0.01	0.89	3.00	2.09	1.77	5.85
14	H2	Q1	18.30	8700.00	18.30	6	7	7	7	0	0	6.00	0.95	0.00	0.05	3.00	2.00	0.58	1.79
15		Q2	22.10	7600.00	19.20	37	72	31	30	1	18	31.00	0.41	0.25	0.34	2.50	4.70	1.28	3.61
16		Q3	21.00	7900.00	28.00	79	111	39	37	2	23	51.00	0.33	0.21	0.46	2.60	8.90	1.47	4.37
17		Q4	13.00	5300.00	26.30	83	124	65	62	3	45	4.00	0.50	0.36	0.14	2.30	7.10	1.49	4.42
18		Q5	14.20	6700.00	29.40	91	139	71	65	6	34	8.00	0.47	0.24	0.29	2.50	6.30	1.51	4.51
19		Q6	16.20	9200.00	31.00	102	189	81	58	23	65	11.00	0.30	0.34	0.35	2.30	5.90	1.53	4.62
20		Q7	15.00	11000.00	29.00	113	198	61	57	4	96	11.00	0.29	0.48	0.23	2.10	4.20	1.55	4.73
21		Q8	11.00	9000.00	30.00	119	207	67	64	3	27	6.00	0.31	0.13	0.56	2.80	4.08	1.56	4.78
22		Q9	10.30	8700.00	33.00	153	257	61	56	5	75	34.00	0.22	0.29	0.49	2.40	4.30	1.62	5.03
23		Q10	10.00	8100.00	35.00	187	249	19	18	1	69	34.00	0.07	0.28	0.65	2.50	3.90	1.65	5.23
24		Q11	10.20	7300.00	34.00	208	312	87	80	7	102	21.00	0.26	0.33	0.42	2.40	2.50	1.67	5.34
25		Q12	11.00	7400.00	33.00	308	459	56	52	4	152	100.00	0.11	0.33	0.56	2.30	1.80	1.75	5.73
26	H3	Q1	32.20	19300.00	27.30	13	24	24	24	0	0	13.00	1.00	0.00	0.00	3.00	2.00	0.94	2.56
27		Q2	29.20	17200.00	25.40	35	104	32	29	3	11	22.00	0.28	0.11	0.61	2.80	4.80	1.27	3.56
28		Q3	24.20	18300.00	31.10	91	176	43	37	6	76	56.00	0.21	0.43	0.36	2.10	8.30	1.51	4.51
29		Q4	21.50	16200.00	33.80	112	239	38	35	3	109	21.00	0.14	0.46	0.40	2.10	7.20	1.55	4.72
30		Q5	19.20	12300.00	36.00	178	319	65	63	2	86	66.00	0.20	0.27	0.53	2.50	6.30	1.65	5.18
31		Q6	16.20	7600.00	39.00	210	329	41	37	4	107	32.00	0.11	0.33	0.56	2.40	5.90	1.68	5.35
32		Q7	18.00	6600.00	41.00	258	430	59	57	2	38	48.00	0.13	0.09	0.78	2.80	5.70	1.71	5.55
33		Q8	17.30	5500.00	39.00	284	512	31	30	1	102	26.00	0.06	0.20	0.74	2.60	5.53	1.73	5.65
34		Q9	16.80	4400.00	40.00	311	581	48	46	2	56	27.00	0.08	0.10	0.82	2.80	4.20	1.75	5.74
35		Q10	14.30	3900.00	39.00	389	619	76	65	11	6	78.00	0.10	0.01	0.89	3.00	3.80	1.79	5.96
36		Q11	14.20	3700.00	38.00	412	654	89	71	18	6	23.00	0.11	0.01	0.88	3.00	2.40	1.80	6.02
37		Q12	14.30	3500.00	39.00	509	987	89	41	48	6	97.00	0.04	0.01	0.95	3.00	1.90	1.83	6.23

print(xtable(ddn[[1]]),type="html")

	PLT	OSR	R
Q1	1.00	1.00	0.18
Q2	0.88	0.92	0.24
Q3	0.82	0.77	0.06
Q4	0.58	0.42	0.00
Q5	0.52	0.17	0.03
Q6	0.89	0.32	0.19
Q7	0.29	0.03	0.14
Q8	0.36	0.11	0.44
Q9	0.23	0.08	0.63
Q10	0.17	0.07	0.51
Q11	0.08	0.03	0.88
Q12	0.00	0.00	1.00

print(xtable(ddn[[2]]),type="html")

	PLT	OSR	R
Q1	0.69	0.60	0.00
Q2	1.00	0.40	0.05
Q3	0.91	0.46	0.58
Q4	0.25	0.00	0.48
Q5	0.35	0.25	0.66
Q6	0.51	0.68	0.76
Q7	0.41	1.00	0.64
Q8	0.08	0.65	0.70
Q9	0.02	0.60	0.88
Q10	0.00	0.49	1.00
Q11	0.02	0.35	0.94
Q12	0.08	0.37	0.88

print(xtable(ddn[[3]]),type="html")

	PLT	OSR	R
Q1	1.00	1.00	0.12
Q2	0.83	0.87	0.00
Q3	0.56	0.94	0.37
Q4	0.41	0.80	0.54
Q5	0.28	0.56	0.68
Q6	0.11	0.26	0.87
Q7	0.21	0.20	1.00
Q8	0.17	0.13	0.87
Q9	0.14	0.06	0.94
Q10	0.01	0.03	0.87
Q11	0.00	0.01	0.81
Q12	0.01	0.00	0.87

print(xtable(ddx[[1]]),type="html")

	PLT	OSR	R
Q1	0.46	0.46	0.08
Q2	0.43	0.45	0.12
Q3	0.50	0.47	0.03
Q4	0.58	0.42	0.00
Q5	0.72	0.23	0.04
Q6	0.63	0.23	0.14
Q7	0.63	0.06	0.31
Q8	0.40	0.12	0.48
Q9	0.24	0.09	0.67
Q10	0.22	0.09	0.68
Q11	0.08	0.03	0.89
Q12	0.00	0.00	1.00

print(xtable(ddx[[2]]),type="html")

	PLT	OSR	R
Q1	0.53	0.47	0.00
Q2	0.69	0.28	0.04
Q3	0.47	0.23	0.30
Q4	0.34	0.00	0.66
Q5	0.28	0.20	0.53
Q6	0.26	0.35	0.39
Q7	0.20	0.49	0.31
Q8	0.06	0.45	0.49
Q9	0.02	0.40	0.59
Q10	0.00	0.33	0.67
Q11	0.01	0.27	0.72
Q12	0.06	0.28	0.66

print(xtable(ddx[[3]]),type="html")

	PLT	OSR	R
Q1	0.47	0.47	0.06
Q2	0.49	0.51	0.00
Q3	0.30	0.50	0.20
Q4	0.23	0.46	0.31
Q5	0.18	0.37	0.45
Q6	0.09	0.21	0.70
Q7	0.15	0.14	0.71
Q8	0.15	0.11	0.74
Q9	0.13	0.05	0.82
Q10	0.01	0.03	0.97
Q11	0.00	0.02	0.98
Q12	0.01	0.00	0.99

Hospital number 1

Ahora tendríamos los valores del sistema ternario normalizado para 36 pasos de tiempo (semanas).

dat=ddx[[1]]
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';}
#

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.9562908 0.9301776 0.9142734 0.9043940 0.8980634 0.8937890
##  [8] 0.8906714 0.8881770 0.8859979 0.8839647 0.8819908

como 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,2,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.45635646 0.45870100 0.09265862 0.96479635 0.35141091 0.56588408 
##        a22        a23         r3        a31        a32        a33 
## 0.34899121 0.73431869 0.63287436 0.21002559 0.67178928 0.45660073

La 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:

GA1=GAs[[1]]
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.457865168539326 0.431079822239339 0.498768477334589 0.582022471910112 
##         1.0000000         1.0093798         1.0123102         1.0101518 
## 0.724090602735764 0.634417379522762 0.629540709812109 0.396458105993873 
##         1.0046866         0.9982813         0.9941681         0.9966916 
## 0.240839952163669 0.224184075533418 0.078822731042803                 0 
##         1.0110200         1.0415910         1.0892145         1.1488020

La conclusión es obvia, en este caso el error es mejor y el ajuste es mucho más satisfactorio.

Gráfico ternario

par(xpd=NA)
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,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))

Hospital number 2

Ahora tendríamos los valores del sistema ternario normalizado para 12 pasos de tiempo (semanas).

dat=ddx[[2]]
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';}
#

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.9662921 0.9464612 0.9344762 0.9269347 0.9218704 0.9181477
##  [8] 0.9151232 0.9124458 0.9099357 0.9075100 0.9051384

como 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,2,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.04899708 0.01809333 0.65177118 0.48383965 0.35704583 0.59545457 
##        a22        a23         r3        a31        a32        a33 
## 0.74197257 0.44554756 0.53246611 0.56834349 0.59554881 0.50659437

La 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:

GA1=GAs[[1]]
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.53487846240814  0.686152958985079  0.467142036894357 
##          1.0000000          0.7953060          0.6794101 
##  0.341048332198775  0.276053451617558  0.261816100826468 
##          0.6050081          0.5538187          0.5172776 
##  0.201185427910563 0.0576979727215089 0.0165121477645546 
##          0.4908300          0.4718098          0.4585126 
##                  0 0.0126413755048947 0.0620778976859207 
##          0.4497467          0.4446001          0.4423181

La conclusión es obvia, en este caso el error es mejor y el ajuste es mucho más satisfactorio.

Gráfico ternario

par(xpd=NA)
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,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))

Hospital number 2

Ahora tendríamos los valores del sistema ternario normalizado para 12 pasos de tiempo (semanas).

dat=ddx[[2]]
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';}
#

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.9662921 0.9464612 0.9344762 0.9269347 0.9218704 0.9181477
##  [8] 0.9151232 0.9124458 0.9099357 0.9075100 0.9051384

como 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,2,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.04899708 0.01809333 0.65177118 0.48383965 0.35704583 0.59545457 
##        a22        a23         r3        a31        a32        a33 
## 0.74197257 0.44554756 0.53246611 0.56834349 0.59554881 0.50659437

La 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:

GA1=GAs[[1]]
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.53487846240814  0.686152958985079  0.467142036894357 
##          1.0000000          0.7953060          0.6794101 
##  0.341048332198775  0.276053451617558  0.261816100826468 
##          0.6050081          0.5538187          0.5172776 
##  0.201185427910563 0.0576979727215089 0.0165121477645546 
##          0.4908300          0.4718098          0.4585126 
##                  0 0.0126413755048947 0.0620778976859207 
##          0.4497467          0.4446001          0.4423181

La conclusión es obvia, en este caso el error es mejor y el ajuste es mucho más satisfactorio.

Gráfico ternario

par(xpd=NA)
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,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))

Hospital Number 3

Ahora tendríamos los valores del sistema ternario normalizado para 12 pasos de tiempo (semanas).

dat=ddx[[3]]
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';}
#

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.9585230 0.9341785 0.9196326 0.9107526 0.9051173 0.9012908
##  [8] 0.8984283 0.8960452 0.8938752 0.8917822 0.8897052

como 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,2,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.2981745 0.5640629 0.5235605 0.5518292 0.3309207 0.1163930 0.3377623 
##       a23        r3       a31       a32       a33 
## 0.9632983 0.7931337 0.1149552 0.1761578 0.7799840

La 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:

GA1=GAs[[1]]
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.471299093655589   0.490074441687345   0.299063787655975 
##           1.0000000           0.9690857           0.9673741 
##   0.232035803046106   0.183445270077165  0.0894326642574128 
##           0.9822297           1.0041280           1.0247808 
##   0.150009994003598    0.14712311168824   0.127006021632214 
##           1.0388473           1.0452659           1.0459131 
## 0.00615460216819392                   0 0.00633219678519249 
##           1.0433268           1.0394443           1.0353935

La conclusión es obvia, en este caso el error es mejor y el ajuste es mucho más satisfactorio.

Gráfico ternario

par(xpd=NA)
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,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))