Я использую Гекко для выбора А-оптимальных экспериментов для набора кинетики реакции.Задача состоит в том, чтобы минимизировать трассу (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()