Estimation des paramètres et minimisation des moindres carrés - Python et Gekko

Question

J'essaie d'estimer 4 paramètres (kdoc1, kdoc2, S1, S2), en utilisant GEKKO pour résoudre un système d'équations algébriques ordinaires, minimisant la somme des carrés de l'erreur résiduelle entre les valeurs observées et prédites. Je n'obtiens pas un bon ajustement entre les données mesurées et les données prédites, avez-vous une suggestion sur ce que je fais mal dans le code ?
En outre, j'essaie également d'ajouter un autre ensemble de données de mesure. Je veux que le modèle fasse l'estimation des paramètres en utilisant deux ensembles de données mesurées, je dois minimiser la somme des résidus en utilisant les deux ensembles de données. Avez-vous une idée sur la façon de le faire dans le code ?

J'ai déjà ajouté des points supplémentaires pour aider le solveur à converger en utilisant la suggestion suivante python - Mesures peu fréquentes dans Gekko avec des points de simulation supplémentaires - Stack Overflow

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt 
import math as math
import pandas as pd

#TOC=5.93 mg/L
#M2 C1 5mg/L
#time: [0,0.08333,0.5,1,4,22.61667]; C2 measured: [0,9.33e-5,8.55e-5,8.55e-5,7.77e-5,7.00e-5]
#M4 C1 20mg/l
#time: [0,0.08333,0.5,1,4,22.61667]; C2 measured: [0,2.92e-4,2.72e-4,2.92e-4,2.72e-4,2.14e-4]

#measurements
t_data = [0.08333,0.5,1,4,22.61667]
x_data = [9.33e-5,8.55e-5,8.55e-5,7.77e-5,7.00e-5]
#x_data2 = [2.92e-4,2.72e-4,2.92e-4,2.72e-4,2.14e-4]
df1 = pd.DataFrame({'time':t_data,'x':x_data})
df1.set_index('time', inplace=True)
# simulation time points
df2 = pd.DataFrame({'time':np.linspace(0,25,51)}) #(0,25,51)
df2.set_index('time', inplace=True)
# merge dataframes
df = df2.join(df1,how='outer')
# get True (1) or False (0) for measurement
df['meas'] = (df['x'].values==df['x'].values).astype(int)
# replace NaN with zeros
df0 = df.fillna(value=0)

#Estimator Model
m = GEKKO()
m.time = df.index.values

# adjustable parameters
kdoc1 = m.FV() #m-1h-1
kdoc2 = m.FV() #m-1h-1
S1 = m.FV()
S2 = m.FV()

#State variables
# M (mol/L)
C0 = m.Var(value=8.00e-5)
C1 = m.Var(value=1.36e-3)
C2 = m.Var(value=x_data[0])
C2m = m.Param(df0['x'].values) 
meas = m.Param(df0['meas'].values)
m.Minimize(meas*(C2-C2m)**2)
C3 = m.Var(value=0)
C4 = m.Var(value=3.16e-8)
C5 = m.Var(value=3.83e-7)
C6 = m.Var(value=0)
C7 = m.Var(value=0)
C8 = m.Var(value=0)
C9 = m.Var(value=0)
C10 = m.Var(value=0)
C11 = m.Var(value=0)
C12 = m.Var(value=3.21e-3)
TOC = m.Var(value=5.93)#mg-C/L 
cC1 = m.Var(value=0)
cC2 = m.Var(value=0)

#temperature ºC
T = 20 

#Rate constants
k1 = 2040000000*math.exp(-1887/(273.15+T))*3600
k2 = 138000000*math.exp(-8800/(273.15+T))*3600
k3 = 300000*math.exp(-2010/(273.15+T))*3600
k4 = 0.00000065*3600
k6 = 26700*3600
k7 = 167*3600
k8 = 27700*3600
k9 = 8300*3600
k10 = 0
k5 = 37800000000*math.exp(-2169/(273.15+T))*C4+2.52e25*math.exp(-16860/(273.15+T))*C10+0.87*math.exp(-503/(273.15+T))*C9

r1 = k1 * C0 * C1
r2 = k2 * C2
r3 = k3 * C0 * C2
r4 = k4 * C3
r5 = k5 * C2 * C2
r6 = k6 * C3 * C1* C4
r7 = k7 * C3 * C5
r8 = k8 * C6 * C3
r9 = k9 * C6 * C2
r10 = k10 * C2 * C3
r11 = kdoc1*S1*TOC*C2/12000
r12 = kdoc2*S2*TOC*C0/12000

t = m.Param(value=m.time)
m.Equation(C0.dt()== -r1 + r2 - r3 + r4 + r8 - r12)
m.Equation(C1.dt()== -r1 + r2 + r5 - r6 + r11)
m.Equation(C2.dt()== r1 - r2 - r3 + r4 - r5 + r6 - r9 - r10 - r11)
m.Equation(C3.dt()== r3 - r4 + r5 - r6 - r7 - r8 - r10)
m.Equation(C4.dt()== 0)
m.Equation(C5 == 1e-14/C4)
m.Equation(C6.dt()== r7 - r8 - r9)
m.Equation(C7 == (3.16e-8*C0)/C4)
m.Equation(C1 == (5.01e-10*C8)/C4)
m.Equation(C11 == (5.01e-11*C10)/C4)
m.Equation(C9 == (5.01e-7*C9)/C4)
m.Equation(C10 == C12 - 2*C11 - C5 + C4)
m.Equation(C12.dt()== 0)
m.Equation(TOC.dt()== 0)
m.Equation(cC1 == 17000*C1)
m.Equation(cC2 == 51500*C2)

#Application options
m.options.SOLVER = 1 #APOPT solver
m.options.IMODE = 5 #Dynamic Simultaneous - estimation
m.options.EV_TYPE = 2 #absolute error
m.options.NODES = 5 #collocation nodes (2,5)

if True:
    kdoc1.STATUS=1
    kdoc2.STATUS=1
    S1.STATUS=1
    S2.STATUS=1

#Solve
m.solve(disp=False)

#Results
data={'Time (h)':t,'C0 (M)':C0.value,'C1 (M)':C1.value, 'C2 (M)':C2.value,'C3 (M)':C3.value,'C4 (M)':C4.value,'C5 (M)':C5.value,'C6 (M)':C6.value,'C7 (M)':C7.value,'C8 (M)':C8.value,'C9 (M)':C9.value,'C10 (M)':C10.value,'C11 (M)':C11.value,'c12 (M)':C12.value}
dfr = pd.DataFrame(data, columns=['Time (h)','C0 (M)','C1 (M)','C2 (M)','C3 (M)','C4 (M)','C5 (M)','C6 (M)','C7 (M)','C8 (M)','C9 (M)','C10 (M)','C11 (M)','C12 (M)'])
#dfr.to_csv(r'C:\dfestimation.csv', index = False)

print('Solution')
print('kdoc1 = ' + str(kdoc1.value[0]))
print('kdoc2 = ' + str(kdoc2.value[0]))
print('S1 = ' + str(S1.value[0]))
print('S2 = ' + str(S2.value[0]))

plt.figure(1,figsize=(12,8))
plt.subplot(3,1,1)
plt.plot(m.time,C2.value,'bo',label='Predicted')
plt.plot(m.time,df['x'].values,'rs',label='Meas')
plt.plot(m.time,C2.value,label ='C2')
plt.legend(loc='best')
plt.subplot(3,1,2)
plt.plot(m.time,cC2.value,label ='C2')
plt.plot(m.time,cC1.value,label ='C1')
plt.legend(loc='best')
plt.subplot(3,1,3)
plt.plot(m.time,C0.value,label ='C0')
plt.plot(m.time,C1.value,label ='C1')
plt.plot(m.time,C3.value,label ='C3')
plt.plot(m.time,C4.value,label ='C4')
plt.plot(m.time,C5.value,label ='C5')
plt.plot(m.time,C6.value,label ='C6')
plt.plot(m.time,C7.value,label ='C7-')
plt.plot(m.time,C8.value,label ='C8')
plt.plot(m.time,C9.value,label ='C9')
plt.plot(m.time,C10.value,label ='C10')
plt.plot(m.time,C11.value,label ='C11')
plt.plot(m.time,C12.value,label ='C12')

plt.xlabel('time (h)')
plt.ylabel('Concentration (mol/L)')
plt.legend(loc='best')

Answer 1

Voici de nouveaux résultats qui sont plus proches.

J'ai ajouté des contraintes aux paramètres et leur ai donné des valeurs de départ non nulles (1).

kdoc1 = m.FV(1,lb=0.01,ub=10) #m-1h-1
kdoc2 = m.FV(1,lb=0.01,ub=10) #m-1h-1
S1 = m.FV(1,lb=0.01,ub=10)
S2 = m.FV(1,lb=0.01,ub=10)

Calculé les conditions initiales pour C0 y C1 parce que ceux-ci créaient un grand saut instantané dans C2 . Si ces concentrations initiales sont connues, vous devrez peut-être revenir à une condition initiale fixe.

m.free_initial(C0)
m.free_initial(C1)

Mise à l'échelle de la fonction objectif car 1e-4**2 est déjà un petit nombre et peut déjà répondre aux critères de convergence de l'objectif. La valeur de la fonction objectif du problème original était très petite. J'ai rendu chaque valeur carrée approximativement égale à 1 comme point de départ.

m.Minimize(meas*(1e4*(C2-C2m))**2)

Cela fournit une solution raisonnable avec laquelle vous pouvez commencer à obtenir les résultats de régression dont vous avez besoin.

kdoc1 = 3.7210285146
kdoc2 = 2.5743961631
S1 = 3.7210285146
S2 = 2.5743961631

Les valeurs des paramètres ne se terminent pas par une contrainte, ce qui implique que cette solution est optimale. Même si les contraintes des paramètres ne sont pas actives à la fin, les contraintes aident à guider le solveur sur la zone de recherche appropriée.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt 
import math as math
import pandas as pd

#measurements
t_data = [0.08333,0.5,1,4,22.61667]
x_data = [9.33e-5,8.55e-5,8.55e-5,7.77e-5,7.00e-5]
#x_data2 = [2.92e-4,2.72e-4,2.92e-4,2.72e-4,2.14e-4]
df1 = pd.DataFrame({'time':t_data,'x':x_data})
df1.set_index('time', inplace=True)
# simulation time points
df2 = pd.DataFrame({'time':np.linspace(0,25,51)})
df2.set_index('time', inplace=True)
# merge dataframes
df = df2.join(df1,how='outer')
# get True (1) or False (0) for measurement
df['meas'] = (df['x'].values==df['x'].values).astype(int)
# replace NaN with zeros
df0 = df.fillna(value=0)
#Estimator Model
m = GEKKO()
m.time = df.index.values

# adjustable parameters
kdoc1 = m.FV(1,lb=0.01,ub=10) #m-1h-1
kdoc2 = m.FV(1,lb=0.01,ub=10) #m-1h-1
S1 = m.FV(1,lb=0.01,ub=10)
S2 = m.FV(1,lb=0.01,ub=10)

#State variables
# M (mol/L)
C0 = m.Var(value=8.00e-5,lb=0,ub=1e-3)
m.free_initial(C0)
C1 = m.Var(value=1.36e-3,lb=0,ub=2e-3)
m.free_initial(C1)
C2 = m.Var(value=x_data[0])
C2m = m.Param(df0['x'].values,lb=0)
meas = m.Param(df0['meas'].values)
m.Minimize(meas*(1e4*(C2-C2m))**2)
C3 = m.Var(value=0)
C4 = m.Var(value=3.16e-8)
C5 = m.Var(value=3.83e-7)
C6 = m.Var(value=0)
C7 = m.Var(value=0)
C8 = m.Var(value=0)
C9 = m.Var(value=0)
C10 = m.Var(value=0)
C11 = m.Var(value=0)
C12 = m.Var(value=3.21e-3)
TOC = m.Var(value=5.93)#mg-C/L 
cC1 = m.Var(value=0)
cC2 = m.Var(value=0)

#temperature ºC
T = 20 

#Rate constants
k1 = 2040000000*math.exp(-1887/(273.15+T))*3600
k2 = 138000000*math.exp(-8800/(273.15+T))*3600
k3 = 300000*math.exp(-2010/(273.15+T))*3600
k4 = 0.00000065*3600
k6 = 26700*3600
k7 = 167*3600
k8 = 27700*3600
k9 = 8300*3600
k10 = 0
k5 = 37800000000*math.exp(-2169/(273.15+T))*C4\
     +2.52e25*math.exp(-16860/(273.15+T))*C10\
     +0.87*math.exp(-503/(273.15+T))*C9

r1 = k1 * C0 * C1
r2 = k2 * C2
r3 = k3 * C0 * C2
r4 = k4 * C3
r5 = k5 * C2 * C2
r6 = k6 * C3 * C1* C4
r7 = k7 * C3 * C5
r8 = k8 * C6 * C3
r9 = k9 * C6 * C2
r10 = k10 * C2 * C3
r11 = kdoc1*S1*TOC*C2/12000
r12 = kdoc2*S2*TOC*C0/12000

t = m.Param(value=m.time)
m.Equation(C0.dt()== -r1 + r2 - r3 + r4 + r8 - r12)
m.Equation(C1.dt()== -r1 + r2 + r5 - r6 + r11)
m.Equation(C2.dt()== r1 - r2 - r3 + r4 - r5 + r6 - r9 - r10 - r11)
m.Equation(C3.dt()== r3 - r4 + r5 - r6 - r7 - r8 - r10)
m.Equation(C4.dt()== 0)
m.Equation(C5 == 1e-14/C4)
m.Equation(C6.dt()== r7 - r8 - r9)
m.Equation(C7 == (3.16e-8*C0)/C4)
m.Equation(C1 == (5.01e-10*C8)/C4)
m.Equation(C11 == (5.01e-11*C10)/C4)
m.Equation(C9 == (5.01e-7*C9)/C4)
m.Equation(C10 == C12 - 2*C11 - C5 + C4)
m.Equation(C12.dt()== 0)
m.Equation(TOC.dt()== 0)
m.Equation(cC1 == 17000*C1)
m.Equation(cC2 == 51500*C2)

#Application options
m.options.SOLVER = 3 #APOPT solver
m.options.IMODE = 5 #Dynamic Simultaneous - estimation
m.options.EV_TYPE = 2 #absolute error
m.options.NODES = 3 #collocation nodes (2,5)

if True:
    kdoc1.STATUS=1
    kdoc2.STATUS=1
    S1.STATUS=1
    S2.STATUS=1

#Solve
m.solve(disp=True)

print('Solution')
print('kdoc1 = ' + str(kdoc1.value[0]))
print('kdoc2 = ' + str(kdoc2.value[0]))
print('S1 = ' + str(S1.value[0]))
print('S2 = ' + str(S2.value[0]))

plt.figure(1,figsize=(8,5))
plt.subplot(2,1,1)
plt.plot(m.time,C2.value,'bo',label='Predicted')
plt.plot(m.time,df['x'].values,'rs',label='Meas')
plt.plot(m.time,C2.value,label ='C2')
plt.legend(loc='best')
plt.subplot(2,1,2)
plt.plot(m.time,cC2.value,label ='C2')
plt.plot(m.time,cC1.value,label ='C1')
plt.legend(loc='best')
plt.xlabel('time (h)')

plt.figure(2,figsize=(12,8))
plt.subplot(4,3,1)
plt.plot(m.time,C0.value,label ='C0')
plt.legend()

plt.subplot(4,3,2)
plt.plot(m.time,C1.value,label ='C1')
plt.legend()

plt.subplot(4,3,3)
plt.plot(m.time,C3.value,label ='C3')
plt.legend()

plt.subplot(4,3,4)
plt.plot(m.time,C4.value,label ='C4')
plt.legend()

plt.subplot(4,3,5)
plt.plot(m.time,C5.value,label ='C5')
plt.legend()

plt.subplot(4,3,6)
plt.plot(m.time,C6.value,label ='C6')
plt.legend()

plt.subplot(4,3,7)
plt.plot(m.time,C7.value,label ='C7-')
plt.legend()

plt.subplot(4,3,8)
plt.plot(m.time,C8.value,label ='C8')
plt.legend()

plt.subplot(4,3,9)
plt.plot(m.time,C9.value,label ='C9')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,10)
plt.plot(m.time,C10.value,label ='C10')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,11)
plt.plot(m.time,C11.value,label ='C11')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,12)
plt.plot(m.time,C12.value,label ='C12')
plt.legend()
plt.xlabel('time (h)')

plt.show()

Answer 2

@John Hedengren, merci pour votre soutien. J'ai modifié mon code en suivant vos idées et ajouté un autre ensemble de données en suivant vos explications ici :
Comment configurer GEKKO pour l'estimation des paramètres à partir de plusieurs ensembles indépendants de données ?
Cependant, je n'arrive pas à trouver une bonne solution, notamment pour la mesure 4. Avez-vous des idées sur la façon dont je peux mieux l'ajuster ? Je vous remercie de votre temps et de votre aide.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt 
import math as math
import pandas as pd

#data set 1
t_data1 = [0.08333,0.5,1,4,22.61667]
x_data1 = [7.77e-5,7e-5,7e-5,6.22e-5,3.89e-5]
x_data1mgl = [4,3.6,3.6,3.2,2]

#data set 2
t_data4 = [0.08333,0.5,1,4,22.61667]
x_data4 = [2.92e-4,2.72e-4,2.92e-4,2.72e-4,2.14e-4]
x_data4mgl = [15,14,15,14,11]

#combine data and time in the same dataframe
data1 = pd.DataFrame({'time':t_data1,'x1':x_data1})
data4 = pd.DataFrame({'time':t_data4,'x4':x_data4})
data1.set_index('time', inplace=True)
data4.set_index('time', inplace=True)

# merge dataframes
data = data1.join(data4, how='outer')

# simulation time points
dftp = pd.DataFrame({'time':np.linspace(0,25,51)})
dftp.set_index('time', inplace=True)

# merge dataframes
df = data.join(dftp,how='outer')

# get True (1) or False (0) for measurement
#df['meas'] = (df['x'].values==df['x'].values).astype(int)
z1 = (df['x1']==df['x1']).astype(int)
z4 = (df['x4']==df['x4']).astype(int)

# replace NaN with zeros
df0 = df.fillna(value=0)

#Estimator Model
m = GEKKO(remote=True)
#m = GEKKO(remote=False) # remote=False to produce local folder with results

m.time = df.index.values

# measurements - Gekko arrays to store measured data
xm = m.Array(m.Param,2)
xm[0].value = df0['x1'].values
xm[1].value = df0['x4'].values

cl=m.Array(m.Var,2)
cl[0].value= 4e-5
cl[1].value= 1.33e-4

TOC=m.Array(m.Var,2)
TOC[0].value= 11.85
TOC[1].value= 11.85

C0=m.Array(m.Var,2)
C1=m.Array(m.Var,2)
cC1=m.Array(m.Var,2)
C2=m.Array(m.Var,2)
C2m=m.Array(m.Var,2)
cC2=m.Array(m.Var,2)
C3=m.Array(m.Var,2)
C4=m.Array(m.Var,2)
C5=m.Array(m.Var,2)
C6=m.Array(m.Var,2)
C7=m.Array(m.Var,2)
C8=m.Array(m.Var,2)
C9=m.Array(m.Var,2)
C10=m.Array(m.Var,2)
C11=m.Array(m.Var,2)
C12=m.Array(m.Var,2)

r1=m.Array(m.Var,2)
r2=m.Array(m.Var,2)
r3=m.Array(m.Var,2)
r4=m.Array(m.Var,2)
r5=m.Array(m.Var,2)
r6=m.Array(m.Var,2)
r7=m.Array(m.Var,2)
r8=m.Array(m.Var,2)
r9=m.Array(m.Var,2)
r10=m.Array(m.Var,2)
r11=m.Array(m.Var,2)
r12=m.Array(m.Var,2)

k5=m.Array(m.Var,2)

#Define GEKKO variables that determine if time point contains data to be used in regression
#index for objective (0=not measured, 1=measured)
zm = m.Array(m.Param,2)
zm[0].value=z1
zm[1].value=z4

# fit to measurement
x=m.Array(m.Var,2)                         
x[0].value=x_data1[0]
x[1].value=x_data4[0]

#adjustable parameters
kdoc1 = 0 #m-1h-1
kdoc2 = m.FV(1,lb=0.01,ub=10) #m-1h-1
S1 = 0
S2 = m.FV(1,lb=0.01,ub=10)

#Define Variables
for i in range(2):

#State variables
# M (mol/L)
    #C0 = m.Var(cl[i])
    C0[i] = m.Var(value=cl[i],lb=0,ub=1e-3)
    C1[i] = m.Var(value=6.9e-6,lb=0,ub=2e-3)
    C2[i] = m.Var(value=x[i],lb=0, ub=1e-3)
    C2m[i] = m.Param(xm[i],lb=0)
    meas = m.Param(zm[i])
    m.Minimize(meas*((C2[i]-C2m[i]))**2)
    C3[i] = m.Var(value=0)
    C4[i] = m.Var(value=3.16e-8)
    C5[i] = m.Var(value=3.83e-7)
    C6[i] = m.Var(value=0)
    C7[i]= m.Var(value=0)
    C8[i] = m.Var(value=0)
    C9[i] = m.Var(value=0)
    C10[i] = m.Var(value=0)
    C11[i] = m.Var(value=0)
    C12[i] = m.Var(value=3.21e-3)
    TOC[i] = m.Var(value=5.93)#mg-C/L 
    cC1[i] = m.Var(value=0)
    cC2[i] = m.Var(value=0)

#temperature ºC
    T = 20 

#Rate constants
    k1 = 2040000000*math.exp(-1887/(273.15+T))*3600
    k2 = 138000000*math.exp(-8800/(273.15+T))*3600
    k3 = 300000*math.exp(-2010/(273.15+T))*3600
    k4 = 0.00000065*3600
    k6 = 26700*3600
    k7 = 167*3600
    k8 = 27700*3600
    k9 = 8300*3600
    k10 = 0
    k5[i] = 37800000000*math.exp(-2169/(273.15+T))*C4[i]\
         +2.52e25*math.exp(-16860/(273.15+T))*C10[i]\
         +0.87*math.exp(-503/(273.15+T))*C9[i]

    r1[i] = k1 * C0[i] * C1[i]
    r2[i] = k2 * C2[i]
    r3[i] = k3 * C0[i] * C2[i]
    r4[i] = k4 * C3[i]
    r5[i] = k5[i] * C2[i] * C2[i]
    r6[i] = k6 * C3[i] * C1[i]* C4[i]
    r7[i] = k7 * C3[i] * C5[i]
    r8[i] = k8 * C6[i] * C3[i]
    r9[i] = k9 * C6[i] * C2[i]
    r10[i] = k10 * C2[i] * C3[i]
    r11[i] = kdoc1*S1*TOC[i]*C2[i]/12000
    r12[i] = kdoc2*S2*TOC[i]*C0[i]/12000

    t = m.Param(value=m.time)
    m.Equation(C0[i].dt()== -r1[i] + r2[i] - r3[i] + r4[i] + r8[i] - r12[i])
    m.Equation(C1[i].dt()== -r1[i] + r2[i] + r5[i] - r6[i] + r11[i])
    m.Equation(C2[i].dt()== r1[i] - r2[i] - r3[i] + r4[i] - r5[i] + r6[i] - r9[i] - r10[i] - r11[i])
    m.Equation(C3[i].dt()== r3[i] - r4[i] + r5[i] - r6[i] - r7[i] - r8[i] - r10[i])
    m.Equation(C4[i].dt()== 0)
    m.Equation(C5[i] == 1e-14/C4[i])
    m.Equation(C6[i].dt()== r7[i] - r8[i] - r9[i])
    m.Equation(C7[i] == (3.16e-8*C0[i])/C4[i])
    m.Equation(C1[i] == (5.01e-10*C8[i])/C4[i])
    m.Equation(C11[i] == (5.01e-11*C10[i])/C4[i])
    m.Equation(C9[i] == (5.01e-7*C9[i])/C4[i])
    m.Equation(C10[i] == C12[i] - 2*C11[i] - C5[i] + C4[i])
    m.Equation(C12[i].dt()== 0)
    m.Equation(TOC[i].dt()== 0)
    m.Equation(cC1[i] == 17000*C1[i])
    m.Equation(cC2[i] == 51500*C2[i])
    m.Equation(S1+S2<1)

#Application options
m.options.SOLVER = 3 #IPOPT solver
m.options.IMODE = 5 #Dynamic Simultaneous - estimation = MHE
m.options.EV_TYPE = 2 #absolute error
m.options.NODES = 3 #collocation nodes (2,5)

if True:
    #kdoc1.STATUS=1
    kdoc2.STATUS=1
    #S1.STATUS=1
    S2.STATUS=1

#Solve
#m.open_folder() # open folder if remote=False to see infeasibilities.txt
m.solve(disp=True)

print('Final SSE Objective: ' + str(m.options.objfcnval))

print('Solution')
#print('kdoc1 = ' + str(kdoc1.value[0]))
print('kdoc2 = ' + str(kdoc2.value[0]))
#print('S1 = ' + str(S1.value[0]))
print('S2 = ' + str(S2.value[0]))

plt.figure(1,figsize=(8,5))
plt.subplot(2,1,1)
plt.plot(m.time,C2[0],'b.--',label='Predicted 1')
plt.plot(m.time,C2[1],'r.--',label='Predicted 4')

plt.plot(t_data1,x_data1,'bx',label='Meas 1')
plt.plot(t_data4,x_data4,'rx',label='Meas 4')
plt.legend(loc='best')

plt.subplot(2,1,2)
plt.plot(m.time,cC2[0].value,'b',label ='C2_M1')
plt.plot(m.time,cC2[1].value,'r',label ='C2_M4')
plt.plot(t_data1,x_data1mgl,'bx',label='Meas 1')
plt.plot(t_data4,x_data4mgl,'rx',label='Meas 4')

#plt.plot(m.time,cC1.value,label ='C1')
plt.legend(loc='best')
plt.xlabel('time (h)')

plt.figure(2,figsize=(12,8))
plt.subplot(4,3,1)
plt.plot(m.time,C0[i].value,label ='C0')
plt.legend()

plt.subplot(4,3,2)
plt.plot(m.time,C1[i].value,label ='C1')
plt.legend()

plt.subplot(4,3,3)
plt.plot(m.time,C3[i].value,label ='C3')
plt.legend()

plt.subplot(4,3,4)
plt.plot(m.time,C4[i].value,label ='C4')
plt.legend()

plt.subplot(4,3,5)
plt.plot(m.time,C5[i].value,label ='C5')
plt.legend()

plt.subplot(4,3,6)
plt.plot(m.time,C6[i].value,label ='C6')
plt.legend()

plt.subplot(4,3,7)
plt.plot(m.time,C7[i].value,label ='C7')
plt.legend()

plt.subplot(4,3,8)
plt.plot(m.time,C8[i].value,label ='C8')
plt.legend()

plt.subplot(4,3,9)
plt.plot(m.time,C9[i].value,label ='C9')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,10)
plt.plot(m.time,C10[i].value,label ='C10')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,11)
plt.plot(m.time,C11[i].value,label ='C11')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,12)
plt.plot(m.time,C12[i].value,label ='C12')
plt.legend()
plt.xlabel('time (h)')

plt.show()