Python GEKKO-ODE에서 배열의 값을 어떻게 사용할 수 있습니까?

Nov 13 2020

프로젝트가 있고 정말 도움이 필요합니다.

기본적으로 우리가 시도하는 것은 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,lact,'bx', label = 'Lactate measured')


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


plot1 = plt.figure(4)

plt.plot(tspan,ammonia,'ro', label = 'Ammonia measured')
plt.plot(tspan,glutamine,'bx', 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()

데이터는 Experimental_Data를 사용하여 얻습니다.

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']

여기에 엑셀 파일을 넣을 수 없기 때문에 데이터는 다음과 같습니다.

그리고 miu는이 모듈 (ann_test)을 사용하여 예측됩니다.

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,144,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,cell_br1, 'mo', label ='BR1')
        plt.plot(xp,np.polyval(pbr,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,144,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,deriv1,'m:',linewidth=2, 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, cell_br2, 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,144,99)
xg = np.array([ 0.0 , 23.0 , 47.0 , 71.5 ,\
                95.0 , 119.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,144,198)
xp = np.linspace(0,144,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 (ann_test의) 값을 미분 방정식과 병합 할 수 없다는 것입니다.

이것은 내가 얻은 오류입니다.

TypeError : 'safe'규칙에 따라 dtype ( 'O')에서 dtype ( 'float64')로 배열 데이터를 캐스팅 할 수 없습니다.

누군가 도와 줄 수 있습니까?

답변

1 JohnHedengren Nov 26 2020 at 11:00

문제는 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매개 변수 를 생성합니다 .