python - Gekko 优化包和 numpy 反函数

Question

我正在使用 Gekko 为一组 react 动力学选择 A 最优实验。目标函数是最小化 trace(inv(Z'Z))，其中 Z 是通过围绕其参数线性化 ODE 计算的尺度灵敏度矩阵。如您所见，目标函数涉及 Z'Z 的倒数。我使用了 numpy(甚至是 scipy)反函数，但遇到了以下错误:“没有找到与指定签名和转换匹配的循环，用于 ufunc inv”

我真的不知道怎么了。没有反函数，优化器工作正常。请，请帮助我。我被这个问题困扰了两个多星期。

代码如下:

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

m = GEKKO()

# setting up the paraemter values
k = [2.11, 8.229, 0.0034, 0.000000000001, 4.1233e-05, 3.798835e-05, 0.0000000001]

k1f = k[0] # Var(bounds=(0,1))
Keq = k[1] # Var(bounds=(0,20))
k2 = k[2] # Var(bounds=(0,0.1)) 
k3 = k[3] # Var(bounds=(0,0.01)) 
k4 = k[4] # Var(bounds=(0,0.0001)) 
k5 = k[5] # Var(bounds=(0,0.0001)) 
k6 = k[6] # Var(bounds=(0,0.0001)) 

# Seting up the reaction time for 6hr
m.time = np.linspace(0,21600)
t = m.Param(value=m.time)

# The decision variable is the inital condition for te component SM1
SM1_0 = m.MV(value=1.0,lb=0.1,ub=100.0)
SM1_0.STATUS = 1

# Setting up the State variables
SM1 = m.Var(value=SM1_0)
D = m.Var(value=0.05)
SM2 = m.Var(value= 1.15)
P = m.Var(value=0.0)
SM1D = m.Var(value=0.0)
H2O = m.Var(value=0.3*595.7/(18.0*100.0))
I1 = m.Var(value=0.0)
I2 = m.Var(value=0.0)
I3 = m.Var(value=0.0)
I4 = m.Var(value=0.0)

# Defining the ODES
m.Equation(SM1.dt() ==  - k1f * SM1 * D + (k1f/Keq) * SM1D - k4 * SM1 * H2O)
m.Equation(D.dt() ==  -k1f * SM1 * D + (k1f/Keq) * SM1D + k2 * SM2 * SM1D - k5 * D * H2O)
m.Equation(SM2.dt() == - k2 * SM2 * SM1D - k3 * SM2 * P)
m.Equation(P.dt() == k2 * SM2 * SM1D - k3 * SM2 * P - k6 * P)
m.Equation(SM1D.dt() == k1f * SM1* D - (k1f/Keq)  * SM1D - k2 * SM2 * SM1D)
m.Equation(H2O.dt() == -k4 * SM1 * H2O - k5 * D * H2O)
m.Equation(I1.dt() == k3 * SM2 * P)
m.Equation(I2.dt() == k4 * SM1 * H2O)
m.Equation(I3.dt() == k5 * D * H2O)
m.Equation(I4.dt() == k6 * P)

# Defining the parameteric sensitivites
S_SM1k1f = m.Var(value=0.0)
S_Dk1f = m.Var(value=0.0)
S_SM2k1f = m.Var(value=0.0)
S_Pk1f = m.Var(value=0.0)
#    
S_SM1Keq = m.Var(value=0.0)
S_DKeq = m.Var(value=0.0)
S_SM2Keq = m.Var(value=0.0)
S_PKeq = m.Var(value=0.0)
#    
#        
S_SM1k2 = m.Var(value=0.0)
S_Dk2 = m.Var(value=0.0)
S_SM2k2 = m.Var(value=0.0)
S_Pk2 = m.Var(value=0.0)

S_SM1k3 = m.Var(value=0.0)
S_Dk3 = m.Var(value=0.0)
S_SM2k3 = m.Var(value=0.0)
S_Pk3 = m.Var(value=0.0)

S_SM1k4 = m.Var(value=0.0)
S_Dk4 = m.Var(value=0.0)
S_SM2k4 = m.Var(value=0.0)
S_Pk4 = m.Var(value=0.0)

S_SM1k5 = m.Var(value=0.0)
S_Dk5 = m.Var(value=0.0)
S_SM2k5 = m.Var(value=0.0)
S_Pk5 = m.Var(value=0.0)

S_SM1k6 = m.Var(value=0.0)
S_Dk6 = m.Var(value=0.0)
S_SM2k6 = m.Var(value=0.0)
S_Pk6 = m.Var(value=0.0)


m.Equation(S_SM1k1f.dt() == -(k1f*D + k4*H2O)*S_SM1k1f - k1f*SM1 * S_Dk1f-SM1*D + SM1D/Keq)

m.Equation(S_Dk1f.dt()==-k1f*D *S_SM1k1f - (k1f*SM1+k5*H2O) * S_Dk1f+ k2*SM1D*S_SM2k1f-SM1*D + SM1D/Keq)

m.Equation(S_SM2k1f.dt() == -(k2*SM1D+ k3*P)*S_SM2k1f - k3*SM2*S_Pk1f)

m.Equation(S_Pk1f.dt() == (k2*SM1D- k3*P)*S_SM2k1f - k3*SM2*S_Pk1f- k6* S_Pk1f)


m.Equation(S_SM1Keq.dt() == -(k1f*D + k4*H2O)*S_SM1Keq - k1f*SM1 * S_DKeq - SM1D/(Keq)**2)

m.Equation(S_DKeq.dt() == -k1f*D *S_SM1Keq - (k1f*SM1+k5*H2O) * S_DKeq+ k2*SM1D*S_SM2Keq- SM1D/Keq**2)

m.Equation(S_SM2Keq.dt() == -(k2*SM1D+ k3*P)*S_SM2Keq - k3*SM2*S_PKeq)

m.Equation(S_PKeq.dt() == (k2*SM1D- k3*P)*S_SM2Keq - k3*SM2*S_PKeq- k6* S_PKeq)


m.Equation(S_SM1k2.dt() == -(k1f*D + k4*H2O)*S_SM1k2 - k1f*SM1 * S_Dk2 )

m.Equation(S_Dk2.dt() == -k1f*D *S_SM1k2 - (k1f*SM1+k5*H2O) * S_Dk2+ k2*SM1D*S_SM2k2 + SM2* SM1D)

m.Equation(S_SM2k2.dt() == -(k2*SM1D+ k3*P)*S_SM2k2 - k3*SM2*S_Pk2-SM2*SM1D)

m.Equation(S_Pk2.dt() == (k2*SM1D- k3*P)*S_SM2k2 - k3*SM2*S_Pk2- k6* S_Pk2+ SM2*SM1D)


m.Equation(S_SM1k3.dt() == -(k1f*D + k4*H2O)*S_SM1k3 - k1f*SM1 * S_Dk3 )

m.Equation(S_Dk3.dt() == -k1f*D *S_SM1k3 - (k1f*SM1+k5*H2O) * S_Dk3+ k2*SM1D*S_SM2k3 )

m.Equation(S_SM2k3.dt() == -(k2*SM1D+ k3*P)*S_SM2k3 - k3*SM2*S_Pk3-SM2*P)

m.Equation(S_Pk3.dt() == (k2*SM1D- k3*P)*S_SM2k3 - k3*SM2*S_Pk3- k6* S_Pk3- SM2*P)


m.Equation(S_SM1k4.dt() == -(k1f*D + k4*H2O)*S_SM1k4 - k1f*SM1 * S_Dk4 -SM1*H2O)

m.Equation(S_Dk4.dt() == -k1f*D *S_SM1k4 - (k1f*SM1+k5*H2O) * S_Dk4+ k2*SM1D*S_SM2k4)

m.Equation(S_SM2k4.dt() == -(k2*SM1D+ k3*P)*S_SM2k4 - k3*SM2*S_Pk4)

m.Equation(S_Pk4.dt() == (k2*SM1D- k3*P)*S_SM2k4 - k3*SM2*S_Pk4- k6* S_Pk4)


m.Equation(S_SM1k5.dt() == -(k1f*D + k4*H2O)*S_SM1k5 - k1f*SM1 * S_Dk5 )

m.Equation(S_Dk5.dt() == -k1f*D *S_SM1k5 - (k1f*SM1+k5*H2O) * S_Dk5+ k2*SM1D*S_SM2k5-D*H2O)

m.Equation(S_SM2k5.dt() == -(k2*SM1D+ k3*P)*S_SM2k5 - k3*SM2*S_Pk5)

m.Equation(S_Pk5.dt() == (k2*SM1D- k3*P)*S_SM2k5 - k3*SM2*S_Pk5- k6* S_Pk5)


m.Equation(S_SM1k6.dt() == -(k1f*D + k4*H2O)*S_SM1k6 - k1f*SM1 * S_Dk6 )

m.Equation(S_Dk6.dt() == -k1f*D *S_SM1k6 - (k1f*SM1+k5*H2O) * S_Dk6+ k2*SM1D*S_SM2k6)

m.Equation(S_SM2k6.dt() == -(k2*SM1D+ k3*P)*S_SM2k6 - k3*SM2*S_Pk6)

m.Equation(S_Pk6.dt() == (k2*SM1D- k3*P)*S_SM2k6 - k3*SM2*S_Pk6- k6* S_Pk6-P)


# defining the Z matrix(Scaled senisitivity matrix)
sk = 5.0*np.array(k)/2.0 # scaling factor for parameters
sy2 = np.array([0.048,0.00021,0.052,0.05]) # scaling factors for measurments

dSM1dk = np.matrix(np.transpose([S_SM1k1f*sk[0], S_SM1Keq*sk[1], S_SM1k2*sk[2], S_SM1k3*sk[3], S_SM1k4*sk[4], S_SM1k5*sk[5], S_SM1k6*sk[6]]))/np.sqrt(sy2[0])
dDdk   = np.matrix(np.transpose([S_Dk1f*sk[0],S_DKeq*sk[1],S_Dk2*sk[2],S_Dk3*sk[3],S_Dk4*sk[4],S_Dk5*sk[5],S_Dk6*sk[6]]))/np.sqrt(sy2[1])
dSM2dk = np.matrix(np.transpose([S_SM2k1f*sk[0], S_SM2Keq*sk[1], S_SM2k2*sk[2], S_SM2k3*sk[3], S_SM2k4*sk[4], S_SM2k5*sk[5], S_SM2k6*sk[6]]))/np.sqrt(sy2[2])
dPdk   = np.matrix(np.transpose([S_Pk1f*sk[0],   S_PKeq*sk[1],   S_Pk2*sk[2],   S_Pk3*sk[3],   S_Pk4*sk[4],   S_Pk5*sk[5],   S_Pk6*sk[6]]))/np.sqrt(sy2[3])

Z = np.vstack((dSM1dk,dDdk,dSM2dk,dPdk))

# definging the objecitve function as trace(inv(Z'Z))
m.Obj(np.trace((np.linalg.inv(np.dot(np.transpose(Z),Z)))))


m.options.IMODE = 6
m.solve()

Answer 1

阿里，

看起来您将无法在您的模型中使用 np.inv。这是因为求解器需要具有模型的第一个和第二个梯度才能求解，而这些对于外部函数来说并不容易获得。这就是为什么许多常用函数(如 np.sqrt)被 Gekko 特定方法(如 GEKKO 中的 m.sqrt)替换的原因。

我建议尝试找到一种方法，用其他东西替换模型中的 np.inv。如果这是应该添加到 GEKKO 的方法，请考虑在 GitHub repo 上提出功能请求.