在 python 中使用带有 GEKKO 的 ODE 系统中的数组变量
Using array variables in system of ODEs with GEKKO in python
我之前 post 回答过有关在参数估计领域使用 GEKKO 的问题。 这个问题对于理解 GEKKO 的一些工作原理非常有用，但我没有利用的一个方面是数组变量的使用。在重新审视这个问题时，我意识到我的 'measured' 数据是错误的，因为我将初始（和未测量）条件作为我测量变量的一部分作为数组的第一个元素。 （很抱歉，我不确定如何使 python 语法高亮显示）
#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])
在重新审视这方面时，我 运行 遇到了一个问题，我必须向 GEKKO 指定一个积分周期，其中必须包括时间 = 0。模拟时间跨度与测量时间跨度的这种变化由于术语长度不同，数据导致我的 objective 函数无法正常工作：
m.Minimize((A-Am)**2) #Array A has 17 elements, array Am has 16
m.Minimize((P-Pm)**2)
m.Minimize((C-Cm)**2)
为了解决这个问题，我想我需要用 m.Array() 来实现我的变量，这样我就可以按元素对它们进行索引，以期设置表单的 objective 函数:
m.Obj(np.sum([(A[i+1]-Am[i])**2 for i in range(len(Am))]))
完整的代码显示在 post 的末尾，但我的问题似乎在于设置系统的 m.Equations()。目前，我不使用索引来引用变量的各个元素：
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(A.dt() == - r1 - r2 - r4 - r5 - r6 )
AttributeError: 'numpy.ndarray' object has no attribute 'dt'
当我尝试将它们设置为数组时，python 冻结了，因为我认为我一定是生成了太多的方程式或其他一些内存问题。也许它试图分别区分索引中的每个方程，这也是不正确的：
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)) ])
如果有人能解释如何在方程中使用数组变量，特别是 ODE，将不胜感激，因为我找不到任何示例，其中数组变量是 ODE 系统的一部分（我一定错过了许多例子）。此外，如果有另一种可能更简单的方法来编写我的 objective 具有不同长度变量的函数（并且不必像我最初那样将它们指定为数组变量），那也很好知道。
完整脚本如下：
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()
使用数组变量没有问题，但这需要 collocation equations 的自定义实现以将导数与变量相关联。为缺失的测量值放置一个占位符并创建一个新参数 mobj 以忽略某些测量值（例如初始测量值）要容易得多。
# ignore first measurement
tobj = np.ones(len(t))
tobj[0] = 0
mobj = m.Param(tobj)
数组对于压缩代码很有用，例如变量的定义。
A,B,C,P = m.Array(m.Var,4,lb=0,fixed_initial=False)
解决方法如下：
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()
我之前 post 回答过有关在参数估计领域使用 GEKKO 的问题。
#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])
在重新审视这方面时，我 运行 遇到了一个问题，我必须向 GEKKO 指定一个积分周期，其中必须包括时间 = 0。模拟时间跨度与测量时间跨度的这种变化由于术语长度不同，数据导致我的 objective 函数无法正常工作：
m.Minimize((A-Am)**2) #Array A has 17 elements, array Am has 16
m.Minimize((P-Pm)**2)
m.Minimize((C-Cm)**2)
为了解决这个问题，我想我需要用 m.Array() 来实现我的变量，这样我就可以按元素对它们进行索引，以期设置表单的 objective 函数:
m.Obj(np.sum([(A[i+1]-Am[i])**2 for i in range(len(Am))]))
完整的代码显示在 post 的末尾，但我的问题似乎在于设置系统的 m.Equations()。目前，我不使用索引来引用变量的各个元素：
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(A.dt() == - r1 - r2 - r4 - r5 - r6 )
AttributeError: 'numpy.ndarray' object has no attribute 'dt'
当我尝试将它们设置为数组时，python 冻结了，因为我认为我一定是生成了太多的方程式或其他一些内存问题。也许它试图分别区分索引中的每个方程，这也是不正确的：
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)) ])
如果有人能解释如何在方程中使用数组变量，特别是 ODE，将不胜感激，因为我找不到任何示例，其中数组变量是 ODE 系统的一部分（我一定错过了许多例子）。此外，如果有另一种可能更简单的方法来编写我的 objective 具有不同长度变量的函数（并且不必像我最初那样将它们指定为数组变量），那也很好知道。
完整脚本如下：
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()
使用数组变量没有问题，但这需要 collocation equations 的自定义实现以将导数与变量相关联。为缺失的测量值放置一个占位符并创建一个新参数 mobj 以忽略某些测量值（例如初始测量值）要容易得多。
# ignore first measurement
tobj = np.ones(len(t))
tobj[0] = 0
mobj = m.Param(tobj)
数组对于压缩代码很有用，例如变量的定义。
A,B,C,P = m.Array(m.Var,4,lb=0,fixed_initial=False)
解决方法如下：
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()