Python GEKKO-如何在ODE中使用数组中的值

问题描述

我们有一个项目,我们确实需要帮助。

基本上,我们正在尝试使用GEKKO解决一个多方程组。但是,参数(miu)之一是由神经网络预测的。 但是,当我们尝试将预测的数据和方程式放在一起时,会出现多个错误

我有两个程序: 这是第一个,这是主要的:

import numpy as np
from gekko import GEKKO,brain
import pandas as pd
import matplotlib.pyplot as plt
from math import e
m = GEKKO(remote=False)    # create GEKKO model --  optimization and accesses solvers of constrained,unconstrained,continuous,and discrete problems

KdQ = 0.001        #degree of degradation of glutamine (1/h)
mG = 1.1e-12# 1.1e-10   #glucose maintenance coefficient (mmol/cell/hour)
YAQ = 0.1#0.90         #yield of ammonia from glutamine
YLG = 0.1 #2            #yield of lactate from glucose
YXG = 2.2e8    #yield of cells from glucose (cells/mmol)
YXQ = 0.5e9#1.5e9    #yield of cells from glutamine (cells/mmol)
KL = 150           #lactate saturation constant (mM)
KA = 40            #ammonia saturation constant (mM)
Kdmax = 0.01       #maximum death rate (1/h)
mumax = 0.044      #maximum growth rate (1/h)
KG = 30#1             #glucose saturation constant (mM)
KQ = 0.22          #glutamine saturation constant (mM)
mQ = 0             #glutamine maintenance coefficient (mmol/cell/hour)
kmu = 0.01         #intrinsic death rate (1/h)
Klysis = 2e-2  #rate of cell lysis (1/h)
Ci_star = 100      #inhibitor saturation concentration (mM)
qi = 2.5e-10   #specific inhibitor production rate (1/h)

#Flow,volume and concentration
Fo = 0         #Feed-rate (L/h)
Fi = 0        #Feed-rate (L/h)
V = 3              #volume (L)
SG = 653           #glucose concentration in the Feed (mM)
SQ = 58.8          #glutamine concentration in the feced (mM)

#Load experimental data
from Experimental_Data import tspan,glucose,glutamine,glutamate,lact,ammonia,cell_br1,cell_br2
# create GEKKO parameter
t = np.linspace(0,144,99)
m.time = t

XT= m.Var(value=5e8,name='XT')         #total cell density (MMcells/L)
XV = m.Var(value=5e8,lb=0,name='XV')   #viable cell density (MMcells/L)

from test_ann import  b,x
# mu values are given by neural network

mu2 = b.think(x)
mu1 = np.array(mu2)

#mu = m.abs3(mu2)
mu = m.sos1(mu1)
Kd = m.Intermediate(Kdmax*(kmu/(mu+kmu)))    #death rate(1/h)
# create GEEKO equations
m.Equation(XT.dt()== mu*XV )
m.Equation(XV.dt() == ((mu - Kd)*XV ))

# solve ODE
m.options.IMODE  = 4  #Simulation   #2-Regression mode
m.options.soLVER = 1  #Public software version
m.options.NODES  = 3  #Default
m.options.COLDSTART = 2
# objective
m.solve(display=False)

# objective
#m.Obj(sum([ (z[j]-1)**2 + y for j in range(p)]))
#figure,axes = plt.subplots(nrows=5,ncols=1)
plot1 = plt.figure(1)
plt.plot(t,XV.value,label='viable cell')
#axes[0].plot(t,XT.value,label='total cell')


plt.xlabel='Time [hr]' 
plt.ylabel='Concentration [cells/ml]'
plt.legend()

plot1 = plt.figure(2)

plt.xlabel='Time [hr]' 
plt.ylabel='Concentration [mM]'
plt.legend()

plot1 = plt.figure(3)
plt.plot(tspan,'bx',label = 'Lactate measured')


plt.xlabel='Time [hr]' 
plt.ylabel='Concentration [mM]'
plt.legend()


plot1 = plt.figure(4)

plt.plot(tspan,'ro',label = 'Ammonia measured')
plt.plot(tspan,label = 'glutamine measured')

plt.xlabel='Time [hr]' 
plt.ylabel='Concentration [mM]'
plt.legend()

plot1 = plt.figure(5)
plt.plot(m.time,mu,label='\u03BC')
plt.plot(m.time,Kd,label='Kd')

plt.xlabel='Time [hr]' 
plt.ylabel='Miu[1/h]'
plt.legend()




plt.show()

使用实验数据获取数据

import pandas as pd

#Load experimental data
df = pd.read_excel(r'path')
sheet = df[0:9] #we have to include row 235  

tspan = sheet['TIME']

cell_br1= sheet['CELL_BR1']
cell_br2= sheet['CELL_BR2']

由于我无法将excel文件放在此处,因此数据如下:

image

然后使用此模块(ann_test)预测miu

from gekko import GEKKO
from gekko import brain
import numpy as np
import matplotlib.pyplot as plt  
from numpy import diff
from scipy.interpolate import CubicSpline


xm = np.array([ 0.0,23.0,47.0,71.5,95.0,119.0,143.0 ]) # 47.0,deriv1 = 0
from Experimental_Data import  cell_br1,cell_br2
def spline(cell):    
    m = GEKKO()
    m.options.IMODE=2
    c = [m.FV(value=0) for i in range(4)]
    x = m.Param(value=xm)
    cell = np.array(cell)
    y = m.CV(value=cell)
    y.FSTATUS = 1
    # polynomial model
    m.Equation(y==c[0]+c[1]*x+c[2]*x**2+c[3]*x**3)
    c[0].STATUS=1
    m.solve(disp=False)
    c[1].STATUS=1
    m.solve(disp=False)
    c[2].STATUS=1
    c[3].STATUS=1
    m.solve(disp=False)
    pbr = [c[3].value[0],c[2].value[0],\
           c[1].value[0],c[0].value[0]]
   # print(pbr)
    xp = np.linspace(0,100)
    plot1 = plt.figure(1)
    if cell[0] == cell_br2[0]:
        plt.plot(xm,cell_br2,'ko',label ='BR2')
        plt.plot(xp,np.polyval(pbr,xp),'g:',linewidth=2)
    elif cell[0]  == cell_br1[0] :
        plt.plot(xm,'mo',label ='BR1')
        plt.plot(xp,'r:',linewidth=2)

    plt.xlabel('time(hr)')
    plt.ylabel('cells')
    plt.legend()
    dx = diff(xp)
    dy1 = diff(np.polyval(pbr,xp))
    deriv1 = dy1/dx
    time =np.linspace(0,99)
    plot1 = plt.figure(2)
    if cell[0] == cell_br2[0]:
        plt.plot(time,deriv1,'b:',linewidth=2,label ='BR2')
    elif cell[0] == cell_br1[0]:
        plt.plot(time,'m:',label ='BR1')
    plt.xlabel('time(hr)')
    plt.ylabel('miu(1/h)')
    plt.legend()
    #plt.show()
    return(deriv1)

m = GEKKO()



from Experimental_Data import  cell_br1,glucose


b = brain.Brain(remote=True)
b.input_layer(2)
b.layer(linear=5)
b.layer(tanh=3)
b.layer(tanh=5)
b.output_layer(1)

x_s = np.linspace(0,99)
xg = np.array([ 0.0,\
                95.0,144.0 ])
cells_spline = CubicSpline(xm,cell_br1) 
y_cells = cells_spline(x_s)
miu_1 = spline(cell_br1)
miu_2 = spline(cell_br2)
scale = [1.0e6,1.0e4]
x = (x_s,y_cells) #,y_glucose) #Inputs (3)
y1 = (miu_1)    #Output (2)
y2 = (miu_2)    #Output (2)

b.learn(x,y1) # train
b.learn(x,y2) # train
yp = b.think(x) # validate
x_1 = np.linspace(0,198)
xp = np.linspace(0,99)
yyp = np.array(yp)
miu = np.reshape(yyp,(99,))


plot1 = plt.figure(3)
plt.plot(x_s,miu,'r-',label = 'Predicted ')
plt.plot(x_s,miu_1,'.',label = 'Experimental points')
plt.xlabel('Time [hr]')
plt.ylabel('miu [1/h]')
plt.legend()
plt.show()

问题是我无法将miu的值(来自an_test)与微分方程合并。

这是我得到的错误

TypeError:无法根据规则“安全”将数组数据从dtype('O')转换为dtype('float64')

可以请人帮忙吗?

解决方法

问题可能是您正在使用m.sos1()函数为您的微分方程生成mu

mu = m.sos1(mu1)
Kd = m.Intermediate(Kdmax*(kmu/(mu+kmu)))    #death rate(1/h)
# create GEEKO equations
m.Equation(XT.dt()== mu*XV )
m.Equation(XV.dt() == ((mu - Kd)*XV ))

要将参数向量(与m.time相同的长度)转换为微分方程,请使用m.Param()创建mu参数。