Using array variables in system of ODEs with GEKKO in python
I have previously posted questions regarding the use of GEKKO in the area of parameter estimation. The question here was very useful in understanding some of the workings of GEKKO, but one aspect I did not leverage was the use of array variables. In revisiting this problem, I realized my 'measured' data was wrong, because I included the initial (and unmeasured) conditions as part of my measured variables as the first element of the arrays. (I apologize as I am not sure how to make the python syntax highlighting happen)
#Experimental data
times = np.array([0.0, 0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250,
0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500,
0.934375, 1.006250, 1.078125, 1.150000])
A_obs = np.array([1.0, 0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096,
0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698,
0.001944, 0.001116, 0.000732, 0.000426])
C_obs = np.array([0.0, 0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264,
0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211,
0.385798, 0.358132, 0.380497, 0.383051])
P_obs = np.array([0.0, 0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637,
0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645,
0.303198, 0.277822, 0.284194, 0.301471])
In revisitng this aspect, I ran into a problem where I had to specify to GEKKO an integration period, which had to include time = 0. This change in the time span for the simulation vs the measured data leads to issues with my objective functions not working due to different lengths of the terms:
m.Minimize((A-Am)**2) #Array A has 17 elements, array Am has 16
m.Minimize((P-Pm)**2)
m.Minimize((C-Cm)**2)
To try to get around this, I think I need to implement my variables with m.Array() so I can index them by element in the hopes of setting up objective functions of the form:
m.Obj(np.sum([(A[i+1]-Am[i])**2 for i in range(len(Am))]))
The complete code is shown at the end of the post, but my issue seems to be in setting up the m.Equations() of my system. Currently, I do not use indexes to refer to the individual elements of the variables:
m.Equation(r1 == k1 * A )
m.Equation(r2 == k2 * A * B )
m.Equation(r3 == k3 * C * B )
m.Equation(r4 == k4 * A )
m.Equation(r5 == k5 * A )
m.Equation(r6 == k6 * A * B )
#mass balance diff eqs, function calls rxn function
m.Equation(A.dt() == - r1 - r2 - r4 - r5 - r6 )
m.Equation(B.dt() == r1 - r2 - r3 - r6 )
m.Equation(C.dt() == r2 - r3 + r4)
m.Equation(P.dt() == r3 + r5 + r6)
The program give an error when using equations of this form:
m.Equation(A.dt() == - r1 - r2 - r4 - r5 - r6 )
AttributeError: 'numpy.ndarray' object has no attribute 'dt'
When I try to set them up as an array, python freezes, as I think I must be generating too many equations or some other memory issue. Perhaps it is trying to differentiate each equation in the index individually, which is also not correct:
m.Equation([r1[k] == k1 * A[k] for k in range(len(r1))])
m.Equation([r2[k] == k2 * A[k] * B[k] for k in range(len(r2))])
m.Equation([r3[k] == k3 * C[k] * B[k] for k in range(len(r3))])
m.Equation([r4[k] == k4 * A[k] for k in range(len(r4))])
m.Equation([r5[k] == k5 * A[k] for k in range(len(r5))])
m.Equation([r6[k] == k6 * A[k] * B[k] for k in range(len(r6))])
#mass balance diff eqs, function calls rxn function
m.Equation([A[j].dt() == - r1[j] - r2[j] - r4[j] - r5[j] - r6[j] for j in range(len(A))] )
m.Equation([B[j].dt() == r1[j] - r2[j] - r3[j] - r6[j] for j in range(len(B))] )
m.Equation([C[j].dt() == r2[j] - r3[j] + r4[j] for j in range(len(C))])
m.Equation([P[j].dt() == r3[j] + r5[j] + r6[j] for j in range(len(P)) ])
If anybody could explain how to use array variables within equations, particulalry ODEs, it would be greatly appreciated, as I could not find any examples where array variables are part of a system of ODEs (I must have missed it in the many examples). Additionally, If there is another, perhaps simpler way, to write my objective functions with variables of different lengths (and not have to specify them as array variables as I had originally), that would be good to know as well.
Full script below:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from gekko import GEKKO
data = {'times':[0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250,
0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500,
0.934375, 1.006250, 1.078125, 1.150000],
'A_obs':[0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096,
0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698,
0.001944, 0.001116, 0.000732, 0.000426],
'C_obs':[0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264,
0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211,
0.385798, 0.358132, 0.380497, 0.383051],
'P_obs':[0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637,
0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645,
0.303198, 0.277822, 0.284194, 0.301471]}
df = pd.DataFrame(data)
df.set_index('times',inplace = True)
m = GEKKO(remote=False)
time_obs = df.index
t = m.time = np.append(0,time_obs) #adding a time = 0 start value to time discretization for integration
Am= m.Array(m.Param,(len(time_obs)))
Cm= m.Array(m.Param,(len(time_obs)))
Pm= m.Array(m.Param,(len(time_obs)))
for i in range(len(time_obs)):
Am[i].value = df.iloc[i]['A_obs']
Cm[i].value = df.iloc[i]['C_obs']
Pm[i].value = df.iloc[i]['P_obs']
A = m.Array(m.Var,(len(t)),lb=0,value=0.0)
A[0].value=1.0
B = m.Array(m.Var,(len(t)),lb=0,value=0.0)
C = m.Array(m.Var,(len(t)),lb=0,value=0.0)
P = m.Array(m.Var,(len(t)),lb=0,value=0.0)
k = m.Array(m.FV,6,value=1,lb=0)
for ki in k:
ki.STATUS = 1
k1,k2,k3,k4,k5,k6 = k
r1 = m.Array(m.Var,(len(t)),lb=0,value=0.0)
r2 = m.Array(m.Var,(len(t)),lb=0,value=0.0)
r3 = m.Array(m.Var,(len(t)),lb=0,value=0.0)
r4 = m.Array(m.Var,(len(t)),lb=0,value=0.0)
r5 = m.Array(m.Var,(len(t)),lb=0,value=0.0)
r6 = m.Array(m.Var,(len(t)),lb=0,value=0.0)
m.Equation(r1 == k1 * A )
m.Equation(r2 == k2 * A * B )
m.Equation(r3 == k3 * C * B )
m.Equation(r4 == k4 * A )
m.Equation(r5 == k5 * A )
m.Equation(r6 == k6 * A * B )
#mass balance diff eqs, function calls rxn function
m.Equation(A.dt() == - r1 - r2 - r4 - r5 - r6 )
m.Equation(B.dt() == r1 - r2 - r3 - r6 )
m.Equation(C.dt() == r2 - r3 + r4)
m.Equation(P.dt() == r3 + r5 + r6)
'''
m.Equation([r1[k] == k1 * A[k] for k in range(len(r1))])
m.Equation([r2[k] == k2 * A[k] * B[k] for k in range(len(r2))])
m.Equation([r3[k] == k3 * C[k] * B[k] for k in range(len(r3))])
m.Equation([r4[k] == k4 * A[k] for k in range(len(r4))])
m.Equation([r5[k] == k5 * A[k] for k in range(len(r5))])
m.Equation([r6[k] == k6 * A[k] * B[k] for k in range(len(r6))])
#mass balance diff eqs, function calls rxn function
m.Equation([A[j].dt() == - r1[j] - r2[j] - r4[j] - r5[j] - r6[j] for j in range(len(A))] )
m.Equation([B[j].dt() == r1[j] - r2[j] - r3[j] - r6[j] for j in range(len(B))] )
m.Equation([C[j].dt() == r2[j] - r3[j] + r4[j] for j in range(len(C))])
m.Equation([P[j].dt() == r3[j] + r5[j] + r6[j] for j in range(len(P)) ])
'''
m.Obj(np.sum([(A[i+1]-Am[i])**2 for i in range(len(Am))]))
m.Obj(np.sum([(C[i+1]-Cm[i])**2 for i in range(len(Cm))]))
m.Obj(np.sum([(P[i+1]-Pm[i])**2 for i in range(len(Pm))]))
m.options.IMODE = 5
m.options.SOLVER = 3 #IPOPT optimizer
m.options.RTOL = 1E-6
m.options.OTOL = 1E-6
m.options.NODES = 6
m.options.DIAGLEVEL = 10
#m.open_folder()
m.solve()
It is no problem to use an Array variable, but this would require a customized implementation of the collocation equations to relate derivatives to the variables. It is much easier to put in a placeholder for the missing measurement and create a new parameter mobj to ignore certain measurements such as the initial measurement.
# ignore first measurement
tobj = np.ones(len(t))
tobj[0] = 0
mobj = m.Param(tobj)
Arrays are useful to condense the code such as with the definition of variables.
A,B,C,P = m.Array(m.Var,4,lb=0,fixed_initial=False)
Here is the solution:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from gekko import GEKKO
data = {'times':[0,0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250,
0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500,
0.934375, 1.006250, 1.078125, 1.150000],
'A_obs':[0,0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096,
0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698,
0.001944, 0.001116, 0.000732, 0.000426],
'C_obs':[0,0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264,
0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211,
0.385798, 0.358132, 0.380497, 0.383051],
'P_obs':[0,0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637,
0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645,
0.303198, 0.277822, 0.284194, 0.301471]}
df = pd.DataFrame(data)
df.set_index('times',inplace = True)
m = GEKKO(remote=False)
time_obs = df.index
t = m.time = time_obs
# ignore first measurement
tobj = np.ones(len(t))
tobj[0] = 0
mobj = m.Param(tobj)
Am= m.Param(df['A_obs'].values)
Cm= m.Param(df['C_obs'].values)
Pm= m.Param(df['P_obs'].values)
A,B,C,P = m.Array(m.Var,4,lb=0,fixed_initial=False)
k = m.Array(m.FV,6,value=1,lb=0)
for ki in k:
ki.STATUS = 1
k1,k2,k3,k4,k5,k6 = k
r = m.Array(m.Var,6,lb=0,value=0.0)
r1,r2,r3,r4,r5,r6 = r
m.Equation(r1 == k1 * A )
m.Equation(r2 == k2 * A * B )
m.Equation(r3 == k3 * C * B )
m.Equation(r4 == k4 * A )
m.Equation(r5 == k5 * A )
m.Equation(r6 == k6 * A * B )
#mass balance diff eqs, function calls rxn function
m.Equation(A.dt() == - r1 - r2 - r4 - r5 - r6 )
m.Equation(B.dt() == r1 - r2 - r3 - r6 )
m.Equation(C.dt() == r2 - r3 + r4)
m.Equation(P.dt() == r3 + r5 + r6)
m.Minimize(mobj*(A-Am)**2)
m.Minimize(mobj*(C-Cm)**2)
m.Minimize(mobj*(P-Pm)**2)
m.options.IMODE = 5
m.options.SOLVER = 3 #IPOPT optimizer
m.options.RTOL = 1E-6
m.options.OTOL = 1E-6
m.options.NODES = 3
m.options.DIAGLEVEL = 0
#m.open_folder()
m.solve()
import matplotlib.pyplot as plt
plt.plot(m.time,Am.value,'ro')
plt.plot(m.time,A.value,'r-',label='A')
plt.plot(m.time,Cm.value,'bo')
plt.plot(m.time,C.value,'b-',label='C')
plt.plot(m.time,Pm.value,'go')
plt.plot(m.time,P.value,'g-',label='P')
plt.xlabel('Time')
plt.legend()
plt.show()
- Change value in Ini to empty using Python configparser
- Script works with error line but won't run corectly when bad line is removed
- FFMPEG not saving logs when converting to audio format
- Numerically obtaining response of a damped driven oscillator gives peak at wrong frequency
- What is the deal with the pony in Python community?
- 'Area' object is not callable
- ValueError: Attribute Users.request is required
- How to place class in its own file when it appears to be inheriting from an instance?
- Python 3.6 with pony 0.7 gives error on commit to oracle db
- Why pyqt tablewigdet is only displaying row number and no data?
- How to check if a value exists in a dictionary?
- OpenCV: draw a rectangle around a region
- Python, Tkinter, trying to pull random numbers from a list based off user input for number and have results open in mew window
- How to remove repeated elements in a vector, similar to 'set' in Python
- Understanding descriptor protocol for 'wrapper-descriptor' itself
- Passing a NumPy 3d array to a C function with a triple pointer as an argument
- How do I create variable variables?
- Turn a tf.data.Dataset to a jax.numpy iterator
- Django 5.1 + Postgresql (debian server)
- Passing in arguments to Dependency function makes it recognized as a Query Parameter
- How to pass parameters to an endpoint using `add_route()` in FastAPI?
- How to make Depends optional in FastAPI?
- Streaming multiple videos through FastAPI to Web Browser causes HTTP requests to stall
- FastAPI error when handling file together with form-data defined in a Pydantic model
- Use serial port in Python without installing external packages
- How can i get the output to print 10 circles instead of 9 in python turtle module
- Unable to install TA-Lib on Ubuntu
- Unable to install Python without sudo access
- Size of an open file object
- two if or more with one ELSE