Python GEKKO: ¿Cómo puedo usar valores de la matriz en mis ODE?

Question

Tenemos un proyecto y realmente necesitamos ayuda.

Básicamente, lo que estamos tratando de hacer es resolver un sistema de ecuaciones múltiples usando GEKKO. Sin embargo, las redes neuronales predicen uno de los parámetros (miu). Sin embargo, cuando intentamos juntar los datos predichos y las ecuaciones, obtenemos múltiples errores.

Tengo dos programas: este es el primero, que es el 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()

Los datos se obtienen utilizando 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']

Como no puedo poner el archivo de Excel aquí, los datos son los siguientes:

Y el miu se predice usando este módulo (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()

El problema es que no puedo fusionar los valores de miu (de ann_test) con las ecuaciones diferenciales.

Este es el error que obtuve:

TypeError: no se pueden convertir datos de matriz de dtype ('O') a dtype ('float64') de acuerdo con la regla 'safe'

¿Puede alguien ayudarme?

Answer 1

El problema puede ser que está utilizando la m.sos1()función para generar musus ecuaciones diferenciales:

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

Para obtener un vector de parámetro (de la misma longitud que m.time) en una ecuación diferencial, utilice m.Param()para crear el muparámetro.