I am trying to use GEKKO to solve a dynamic parameter identification problem. For context it is a CSTR producing biodiesel in which I am trying to find the average reaction rate.
x[3] is the biodiesel fraction on the reactor output that I am trying to fit to ymeas that come from lab measurements. u is a matrix with the measurements of four load variables.
Given the following formaulation:
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
t = np.arange(0, 36000, 3600)
ymeas = np.array([0.429, 0.429, 0.427, 0.429, 0.428, 0.429, 0.429, 0.428, 0.428, 0.429])
u = np.array([
[2996.9, 2989.4, 2997.0, 3005.6, 3001.7, 2988.8, 2986.6, 3005.5, 2993.9, 2994.0,],
[656.9, 656.8, 657.7, 658.8, 657.3, 660.2, 659.5, 662.0, 659.7, 659.0,],
[333.1, 334.1, 332.6, 333.6, 331.0, 330.9, 331.0, 330.8, 331.7, 331.3,],
[323.2, 324.2, 324.9, 325.9, 326.9, 327.11, 327.6, 327.6, 327.4, 326.8],
])
par0=0.046
m = GEKKO(remote=False) # create GEKKO model
m.time = t # time points
M = np.array([0.853, 0.032, 0.092, 0.286,])
ro = np.array([954.0, 757.0, 1340.0, 844.0,])
cp = np.array([2110.0, 2785.0, 2556.0, 2146.0,])
xo = np.array([1, 0, 0, 0,])
xm = np.array([0, 1, 0, 0,])
vmol = M / ro
cpmol = cp * M
# create GEKKO constants
nc = 4
Mo = m.Const(value=0.853)
Mm = m.Const(value=0.032)
cpmolo = m.Const(value=cpmol[0])
cpmolm = m.Const(value=cpmol[1])
VR = m.Const(20)
dHr = m.Const(-6309)
# create GEKKO parameter
r = m.FV(value=par0)
u0 = m.Param(value=u[0, :])
u1 = m.Param(value=u[1, :])
u2 = m.Param(value=u[2, :])
u3 = m.Param(value=u[3, :])
# create GEKKO variables
x = m.Array(m.Var, 4)
x[0].value = 0.0031
x[1].value = 0.4235
x[2].value = 0.1432
x[3] = m.CV(value=ymeas)
x[3].FSTATUS = 1
T = m.Var(333.5500)
rx = [-r, -3*r, r, 3*r,]
# create GEKKO equations
To = u2
Tm = u3
No = m.Intermediate(u0/Mo/3600)
Nm = m.Intermediate(u1/Mm/3600)
cpmolR = m.Intermediate(np.sum(x * cpmol))
nR = m.Intermediate(VR / np.sum(vmol * x))
m.Equation(nR*x[i].dt() == Nm*(xm[i] - x[i]) + No*(xo[i] - x[i]) + rx[i]*VR for i in range(nc))
m.Equation(nR*cpmolR*T.dt() == Nm*cpmolm*(Tm - T) + No*cpmolo*(To - T) + VR*(-dHr)*r)
# solve ODE
r.STATUS = 1
m.options.IMODE = 5 # dynamic estimation
m.options.NODES = 3 # collocation nodes
m.solve(disp=False)
print('Parameter identification: done!')
I get the following error:
Exception: @error: Inequality Definition
invalid inequalities: z > x < y
<generatorobject<genexpr>at0x7f88ecb7cf90>
STOPPING . . .
I can't think of a way to fix this, because I did not define inequalities.
Anyone has ideas to fix this and get the solver to find r?
I tried create x[3] as an independent y and append to the x array but got a shape error in Intermediate equation cpmolR = m.Intermediate(np.sum(x * cpmol)).
Use m.open_folder() and open the gk_model0.apm model file with a text editor. It shows that a generator is used for the 2nd-to-last equation.
Equations
<generator object <genexpr> at 0x7f9adc07bc10>
((((i9)*(i8)))*($v6))=((((((i7)*(i3)))*((p5-v6)))
+((((i6)*(i2)))*((p4-v6))))+((((i4)*((-i5))))*(p1)))
End Equations
Use square brackets to make it a list of equations:
m.Equation([nR*x[i].dt() == Nm*(xm[i] - x[i]) \
+ No*(xo[i] - x[i]) \
+ rx[i]*VR \
for i in range(nc)])
There were a few other issues with the syntax that are resolved with this script:
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
t = np.arange(0, 36000, 3600)
ymeas = np.array([0.429, 0.429, 0.427, 0.429, 0.428, 0.429, 0.429, 0.428, 0.428, 0.429])
u = np.array([
[2996.9, 2989.4, 2997.0, 3005.6, 3001.7, 2988.8, 2986.6, 3005.5, 2993.9, 2994.0,],
[656.9, 656.8, 657.7, 658.8, 657.3, 660.2, 659.5, 662.0, 659.7, 659.0,],
[333.1, 334.1, 332.6, 333.6, 331.0, 330.9, 331.0, 330.8, 331.7, 331.3,],
[323.2, 324.2, 324.9, 325.9, 326.9, 327.11, 327.6, 327.6, 327.4, 326.8]])
par0=0.046
m = GEKKO(remote=False) # create GEKKO model
m.time = t # time points
M = np.array([0.853, 0.032, 0.092, 0.286])
ro = np.array([954.0, 757.0, 1340.0, 844.0])
cp = np.array([2110.0, 2785.0, 2556.0, 2146.0])
xo = np.array([1, 0, 0, 0])
xm = np.array([0, 1, 0, 0])
vmol = M / ro
cpmol = cp * M
# create GEKKO constants
nc = 4
Mo = 0.853
Mm = 0.032
cpmolo = cpmol[0]
cpmolm = cpmol[1]
VR = 20
dHr = -6309
# create GEKKO parameter
r = m.FV(value=par0)
u0 = m.Param(value=u[0, :])
u1 = m.Param(value=u[1, :])
u2 = m.Param(value=u[2, :])
u3 = m.Param(value=u[3, :])
# create GEKKO variables
x = np.array([None]*4)
x[0] = m.Var(0.0031)
x[1] = m.Var(0.4235)
x[2] = m.Var(0.1432)
x[3] = m.Var(0.429)
y = m.Param(ymeas)
T = m.Var(333.5500)
rx = [-r, -3*r, r, 3*r,]
# create GEKKO equations
To = u2
Tm = u3
No = m.Intermediate(u0/Mo/3600)
Nm = m.Intermediate(u1/Mm/3600)
cpmolR = m.Intermediate(np.sum(x * cpmol))
nR = m.Intermediate(VR / np.sum(vmol * x))
m.Equations([nR*x[i].dt() == Nm*(xm[i] - x[i]) + No*(xo[i] - x[i]) + rx[i]*VR for i in range(nc)])
m.Equation(nR*cpmolR*T.dt() == Nm*cpmolm*(Tm - T) + No*cpmolo*(To - T) + VR*(-dHr)*r)
# minimize sum of squared errors
m.Minimize((x[3]-y)**2)
# solve ODE
r.STATUS = 1
m.options.IMODE = 5 # dynamic estimation
m.options.NODES = 3 # collocation nodes
#m.open_folder()
m.solve(disp=True)
print('Parameter identification: done!')
This produces a successful solution:
----------------------------------------------
Dynamic Estimation with APOPT Solver
----------------------------------------------
Iter Objective Convergence
0 1.65358E-05 2.02852E-03
1 1.19327E-05 1.09210E-07
2 1.18844E-05 3.12011E-10
3 1.18844E-05 2.64239E-12
4 1.18844E-05 2.64239E-12
Successful solution
---------------------------------------------------
Solver : IPOPT (v3.12)
Solution time : 1.279999999678694E-002 sec
Objective : 1.188444178924983E-005
Successful solution
---------------------------------------------------