I'm new to Gekko. Following @John Hedengren documentation and answers to this questions from @jmiguel, I'm trying to estimate two parameters of a biological reaction described as follows:
- R+L<-->Y
- R+I<-->X
This result in a system of five coupled differential equations:
- dX/dt=k1*R(t)L(t)-k2Y(t)
- dY/dt=k3*R(t)I(t)-k4X(t)
- dR/dt=k2Y(t)-k1R(t)L(t)+k4X(t)-k3*R(t)*I(t)
- dL/dt=k2Y(t)-k1R(t)*L(t)
- dI/dt=k4X(t)-k3R(t)*I(t) where k1 and k2 are knowed constants and k3,k4 are the two parameters that I want to estimate.
To estimate I have 4 sets of 'Y' data that correspond to four different initial concentrations of 'I' measured at the same times. I want to use all this data to estimate k3 and k4, but I don't know how to implement the four data sets together.
Here is my try:
#Modules
import matplotlib.pyplot as plt
import numpy as np
from gekko import GEKKO
import math as math
import pandas as pd
#4 sets of data for Y with same tame and differente concentrations
texp=[0,0.5,0.75,1,1.25,2,2.77,3.55,4.32,5.1,5.87,6.65,7.42,8.2,8.97,13.92,18.92,23.92,28.92,33.92,38.92,43.92,48.92,53.92,58.92,63.92,68.92,83.9,98.9,113.9,128.9,143.9,158.9,173.9,188.9,203.9,218.9,233.9,248.9]
yexp=[[0,134.0972365,101.6778993,95.40580822,77.46417317,97.38426357,94.95124669,92.01743212,103.1524256,101.6335581,100.8064012,97.36957339,98.54218566,105.1649371,100.7135276,91.76116765,98.50303784,99.09282871,101.2884698,98.16682444,98.95862295,100.1781904,97.67407619,98.05420149,101.7811854,99.65343296,99.35574195,101.5667478,97.8827634,99.21928788,100.2506635,103.3401759,97.77635857,95.59188964,97.63152972,98.68365575,97.94605393,94.40686956,101.4638065],
[0,52.8368317,60.63226219,57.1596205,77.40224428,76.59558061,83.08583996,92.71540378,92.74825352,97.77466238,87.90927062,91.55955552,93.40910581,97.38680721,98.68892712,96.28587494,97.24139983,99.14913748,97.05786666,99.52031603,98.25378711,98.45507123,98.48665412,96.58596306,98.49718892,97.03260859,97.78876552,96.96130531,96.16927889,99.05385817,99.81984318,98.34275219,97.06707695,98.13554762,96.40027439,97.02992383,97.92157304,97.45153492,100.7723151],
[0,6.651270535,-0.894974357,28.70354922,41.78345531,28.39710692,38.44803531,52.61165239,31.33900251,47.7098898,46.13716223,60.80324347,63.06296709,61.14672153,62.56302041,73.57634262,80.68134402,84.09161715,86.42168805,83.62885021,82.69730894,87.38430639,92.26394519,87.78013332,85.96624579,87.84265572,85.32430667,87.74945544,87.06258236,88.05326643,86.29714124,90.465153,86.36689116,81.69960378,87.69867171,82.08550027,85.6811316,88.07994935,87.69384792],
[0,21.00084301,-54.20967226,-12.0118567,-25.27198718,-1.764831016,10.29814076,-5.340599221,6.988265971,9.56252586,-3.705303123,1.063813346,12.32611118,7.763248428,9.074028389,20.60003402,22.1001936,23.13229101,27.31536018,25.00455108,31.70315201,35.10288809,38.0816535,35.30253723,36.81655545,36.11171691,41.57221204,42.47852506,46.28315167,42.66070948,44.73318881,37.36241544,39.69557981,38.71667563,37.49757832,42.35943236,41.68017195,44.91883581,47.80088108]]
#Initial condition
Y0=[0,0,0,0]
X0=[0,0,0,0]
R0=[0.5,0.5,0.5,0.5]
L0=[30e-9,30e-9,30e-9,30e-9]
I0=[1000e-9,100e-9,10e-9,1e-9]#One initial condicion for I, one per dataset
#Combine the data in the same dataframe
#···Define each dataset in a dataframe
df1=pd.DataFrame({'time':texp,'y1':yexp[0]})
df1.set_index('time',inplace=True)
df2=pd.DataFrame({'time':texp,'y2':yexp[1]})
df2.set_index('time',inplace=True)
df3=pd.DataFrame({'time':texp,'y3':yexp[2]})
df3.set_index('time',inplace=True)
df4=pd.DataFrame({'time':texp,'y4':yexp[3]})
df4.set_index('time',inplace=True)
#···Merge dataframes
data=df1.join(df2,how='outer').join(df3,how='outer').join(df4,how='outer')
#Model
m=GEKKO(remote=False)
t=m.time=data.index
#···Arrays to store observed values of Y
Ym=m.Array(m.Param,4)
Ym[0].value=data['y1'].values
Ym[1].value=data['y2'].values
Ym[2].value=data['y3'].values
Ym[3].value=data['y4'].values
#···Variable for differential equations
Y=m.Array(m.Var,4,lb=0)
X=m.Array(m.Var,4,lb=0)
R=m.Array(m.Var,4,lb=0)
L=m.Array(m.Var,4,lb=0)
I=m.Array(m.Var,4,lb=0)
#···Initial values
for j in range(4):
Y[j].value=Y0[j]
X[j].value=X0[j]
R[j].value=R0[j]
L[j].value=L0[j]
I[j].value=I0[j]
#···Parameters
#····Fixed
k1=m.Const(value=3.58e6)
k2=m.Const(value=1.25e-1)
#····Variable
k3=m.FV(value=1e8,lb=10e1,ub=10e10)
k3.STATUS=1
k4=m.FV(value=1e-2,lb=10e-5,ub=10e3)
k4.STATUS=1
#···Define reaction rate variables
r1=m.Array(m.Var,4,value=1,lb=0)
r2=m.Array(m.Var,4,value=1,lb=0)
r3=m.Array(m.Var,4,value=1,lb=0)
r4=m.Array(m.Var,4,value=1,lb=0)
#···Model equations
for k in range(4):
#····Rate equations
m.Equation(r1[k] == k1*R[k]*L[k])
m.Equation(r2[k] == k2*Y[k])
m.Equation(r3[k] == k3*R[k]*I[k])
m.Equation(r4[k] == k4*X[k])
#····Differential species balances
m.Equation(X[k].dt() == -r4[k]+r3[k])
m.Equation(Y[k].dt() == -r2[k]+r1[k])
m.Equation(R[k].dt() == r2[k]-r1[k]+r4[k]-r3[k])
m.Equation(L[k].dt() == r2[k]-r1[k])
m.Equation(I[k].dt() == r4[k]-r3[k])
#····Minimization objective functions
m.Minimize((Y[k]-Ym[k])**2)
#···Solver
m.options.IMODE=2
m.options.SOLVER=3
m.options.MAX_ITER=10000
m.options.NODES=3
#m.options. COLDSTART=2
m.solve()
print(k3.value[0], k4.value[0])
labels=['X0','X1','X2','X3','Y0','Y1','Y2','Y3','R0','R1','R2','R3','L0','L1','L2','L3','I0','I1','I2','I3','yexp0','yexp1','yexp2','yexp3']
datas=[X[0].value,X[1].value,X[2].value,X[3].value,Y[0].value,Y[1].value,Y[2].value,Y[3].value,R[0].value,R[1].value,
R[2].value,R[3].value,L[0].value,L[1].value,L[2].value,L[3].value,I[0].value,I[1].value,I[2].value,I[3].value,
yexp[0],yexp[1],yexp[2],yexp[3]]
plt.figure(1,figsize=(24,10))
for f in range(24):
plt.subplot(6,4,f+1)
plt.plot(t,datas[f],label=labels[f])
plt.legend()
plt.show()
And here is the answer:
Warning: model time only used for dynamic modes (IMODE>3) ---------------------------------------------------------------- APMonitor, Version 1.0.0 APMonitor Optimization Suite ----------------------------------------------------------------
--------- APM Model Size ------------ Each time step contains Objects : 0 Constants : 2 Variables : 42 Intermediates: 0 Connections : 0 Equations : 40 Residuals : 40 Number of state variables: 1406 Number of total equations: - 1404 Number of slack variables: - 0 --------------------------------------- Degrees of freedom : 2 ********************************************** Model Parameter Estimation with Interior Point Solver ********************************************** Info: Exact Hessian ****************************************************************************** This program contains Ipopt, a library for large-scale nonlinear optimization. Ipopt is released as open source code under the Eclipse Public License (EPL). For more information visit http://projects.coin-or.org/Ipopt ****************************************************************************** This is Ipopt version 3.10.2, running with linear solver mumps. Number of nonzeros in equality constraint Jacobian...: 3744 Number of nonzeros in inequality constraint Jacobian.: 0 Number of nonzeros in Lagrangian Hessian.............: 936 MUMPS returned INFO(1) = -9 and requires more memory, reallocating. Attempt 1 Increasing icntl[13] from 1000 to 2000. MUMPS returned INFO(1) = -9 and requires more memory, reallocating. Attempt 2 Increasing icntl[13] from 2000 to 4000. Total number of variables............................: 1406 variables with only lower bounds: 1404 variables with lower and upper bounds: 2 variables with only upper bounds: 0 Total number of equality constraints.................: 1404 Total number of inequality constraints...............: 0 inequality constraints with only lower bounds: 0 inequality constraints with lower and upper bounds: 0 inequality constraints with only upper bounds: 0 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 0 9.4857890e+005 5.00e+005 1.01e+002 0.0 0.00e+000 - 0.00e+000 0.00e+000 0 1 9.4878851e+005 1.24e+005 2.38e+000 -5.1 1.00e+000 - 2.40e-001 9.84e-001h 1 2 9.4879369e+005 3.17e+004 6.49e+001 -1.3 2.66e-001 2.0 9.90e-001 9.90e-001h 1 3 9.4879397e+005 4.07e+002 2.23e+001 -7.0 2.93e-001 1.5 7.05e-001 7.04e-001h 1 4 9.4879401e+005 7.00e+000 1.48e+001 -7.5 7.83e-002 - 3.59e-001 3.55e-001h 1 ... Number of Iterations....: 107 (scaled) (unscaled) Objective...............: 1.4048811068946009e+003 3.7678134809125418e+003 Dual infeasibility......: 8.5533874233035452e-007 2.2939712323577221e-006 Constraint violation....: 4.2404425129693523e-009 7.5660587128822954e-007 Complementarity.........: 7.3513361195303513e-010 1.9715877164233078e-009 Overall NLP error.......: 8.5533874233035452e-007 2.2939712323577221e-006 Number of objective function evaluations = 108 Number of objective gradient evaluations = 108 Number of equality constraint evaluations = 108 Number of inequality constraint evaluations = 0 Number of equality constraint Jacobian evaluations = 108 Number of inequality constraint Jacobian evaluations = 0 Number of Lagrangian Hessian evaluations = 107 Total CPU secs in IPOPT (w/o function evaluations) = 4.992 Total CPU secs in NLP function evaluations = 1.235 EXIT: Optimal Solution Found. The solution was found. The final value of the objective function is 3767.8134809125418 --------------------------------------------------- Solver : IPOPT (v3.12) Solution time : 6.2372000000000005 sec Objective : 3767.8134828332963 Successful solution --------------------------------------------------- k3=4563.4117885 k4=1732.3167918
To check that it is working, I have solved this using an simplified equivalent system of equations, solving it symbolically and doing a non-linear regression using curve_fit from scipy and the results are very different.
I suspect that I am not getting the data entered correctly in the model. Moreover, maybe the model is not well built at all.
Thank you for your time reading my question and for any help/ideas you may have!
