Python GEKKO - Comment puis-je utiliser les valeurs d'un tableau dans mes ODE

Question

Nous avons un projet et nous avons vraiment besoin d'aide.

Fondamentalement, ce que nous essayons de faire est de résoudre un système d'équations multiples en utilisant GEKKO. Cependant, l'un des paramètres (miu) est prédit par les réseaux de neurones. Cependant, lorsque nous essayons de mettre ensemble les données prédites et les équations, nous obtenons de multiples erreurs.

J'ai deux programmes: C'est le premier, qui est le principal:

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()

Les données sont obtenues à l'aide de 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']

Comme je ne peux pas mettre le fichier Excel ici, les données sont les suivantes:

Et le miu est prédit en utilisant ce module (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()

Le problème est que je ne peux pas fusionner les valeurs de miu (de ann_test) avec les équations différentielles.

Voici l'erreur que j'ai obtenue:

TypeError: Impossible de convertir les données du tableau de dtype ('O') en dtype ('float64') selon la règle 'safe'

Quelqu'un peut-il aider?

Answer 1

Le problème peut être que vous utilisez la m.sos1()fonction pour générer muvos équations différentielles:

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 ))

Pour obtenir un vecteur de paramètres (de même longueur que m.time) dans une équation différentielle, utilisez m.Param()pour créer le muparamètre.