OverflowError: (34, «Численный результат вне диапазона») при решении и оптимизации системы ODE с использованием Pyomo? - PullRequest
Новая
       18

OverflowError: (34, «Численный результат вне диапазона») при решении и оптимизации системы ODE с использованием Pyomo?

1 голос
/ 25 апреля 2019

В моей попытке оптимизировать процесс, моделируемый системой ODE, всякий раз, когда я пытаюсь добавить температурную зависимость от параметра, я продолжаю получать ошибки.Как я могу исправить этот код для правильной работы и где я иду не так?

Впервые мои коллеги и я попытались добавить функцию в модель ODE.Я попытался взять зависимые от температуры (T) переменные и просто установить их как функции, не определяя их как Constraint-s.Я также попытался не определять их явно, а просто вставить формулы для переменных, зависящих от температуры, в оду.

#importing all of the packages that we will need for optimization
from pyomo.environ import *
from pyomo.dae import *
import math as m
import numpy as np

#Initialize Constants

#rate constants 

#constants for energy balance
rho = 1040        #g/L
C_p = 4.016       #J/grams/K
DH_FG = -91.2     #KJ/gmol
DH_FM = -226.3    #kJ/gmol
DH_FN = -361.3    #kJ/gmol
Tc = 5 + 273.15   #K
R = 8.314         #J/mol/K
#We are not exactly sure what these are. We think they the gains for the Temperature Controls
p1 = 1.6  
p2 = 0.188
p3 = 0.1
p = (p1*p2)/p3
e = m.exp(1)

 #creating the final time parameter
t_f = 300  #final time

#Creating initial conditions for each of the equations

#Initial values...
u_x = .5            #h^-1         
X_int = 125    #biomass 
G_int = 70     #glucose
M_int = 220    #Maltose
N_int = 40     #Maltotriose
T_int = 20     #Temperature of system
u = 400        #Control Variable
uG_int = m.exp(35.77)
KprimeG_int = 2
KM_int = m.exp(-19.5)
KG_int = m.exp(-121.3)
u1_int = uG_int*G_int/(KG_int +G_int)
u2_int = 1
u3_int = 1


#Activation Energies
E_uG = 22.6*4184 # J/mol
E_KG = -68.6*4184 #J/mol
g = ConcreteModel()

#Growth Model

#Establishing the time variable
g.t = ContinuousSet(bounds=(0, t_f))    #creating time parameter with bounds from 0 to tf

#Establishing the dependent variables on time

g.X = Var(g.t)  #Amount of Biomass in the system as a function of t
g.G = Var(g.t)  #amount of Glucose in the system as a function of t
g.M = Var(g.t)  #amount of Maltose in the system as a function of t
g.N = Var(g.t)  #Amount of Maltotriose in the system is a function of t 
g.T = Var(g.t)

g.KG = Var(g.t)
g.uG = Var(g.t)
g.u1 = Var(g.t)

#Creating Derivatives of each of the depend variables
g.dXdt = DerivativeVar(g.X, wrt=g.t)
g.dGdt = DerivativeVar(g.G, wrt=g.t)
g.dMdt = DerivativeVar(g.M, wrt=g.t)
g.dNdt = DerivativeVar(g.N, wrt=g.t)
g.dTdt = DerivativeVar(g.T, wrt=g.t)

#Temperature dependent variables
g.uGeq  = Constraint(g.t, rule=lambda g, t: g.uG[t]   == uG_int*e**(-E_uG/R/g.T[t]) )
g.KGeq  = Constraint(g.t, rule=lambda g, t: g.KG[t]   == KG_int*e**(-E_KG/R/g.T[t]) )

#creating the system of ode's
g.u1eq  = Constraint(g.t, rule=lambda g, t: g.u1[t]   == g.uG[t]*g.G[t]/(g.KG[t] + g.G[t]) )
g.X_ode = Constraint(g.t, rule=lambda g, t: g.dXdt[t] == u_x*g.X[t])
g.G_ode = Constraint(g.t, rule=lambda g, t: g.dGdt[t] == g.u1[t]*g.X[t])
g.M_ode = Constraint(g.t, rule=lambda g, t: g.dMdt[t] == -u2_int*g.X[t])
g.N_ode = Constraint(g.t, rule=lambda g, t: g.dNdt[t] == -u3_int*g.X[t])
g.T_ode = Constraint(g.t, rule=lambda g, t: g.dTdt[t] == (g.u1[t]*DH_FG + u2_int*DH_FM + u3_int*DH_FN)*g.X[t] + u*(Tc - g.T[t])/(rho*C_p))

g.KG[0].fix(KG_int)
g.uG[0].fix(uG_int)
g.X[0].fix(X_int)
g.G[0].fix(G_int)
g.M[0].fix(M_int)
g.N[0].fix(N_int)
g.T[0].fix(T_int)


def solve(g):
    TransformationFactory('dae.finite_difference').apply_to(g, nfe=200, method='backwards')
    SolverFactory('ipopt',  executable=ipopt_executable).solve(g)

    #creating lists to store the different values as a function of time
    tsim  = [t for t in g.t]
    Xsim  = [g.X[t]() for t in g.t]
    Gsim  = [g.G[t]() for t in g.t]
    Msim  = [g.M[t]() for t in g.t]
    Nsim  = [g.N[t]() for t in g.t]
    Tsim  = [g.T[t]() for t in g.t]
    uGsim = [g.uG[t]() for t in g.t]
    KGsim = [g.KG[t]() for t in g.t]
    u1sim = [g.u1[t]() for t in g.t]

#run the solver and see what happens    
solve(g)

Мы ожидаем, что различные переменные изменятся со временем, и отобразим результаты с помощью mathplotlib

1 Ответ

0 голосов
/ 25 апреля 2019

Я думаю, вы видите эту ошибку из-за вашего ограничения g.KGeq. Если вы инициализируете g.T в T_int=20, показатель степени в этом уравнении будет примерно 1711, что приведет к ошибке переполнения Python. У вас может быть опечатка с отрицательным знаком.

Еще одна вещь, которую я вижу, это то, что вы делите на g.T и g.KG + g.G в некоторых ваших ограничениях, но вы никогда не инициализируете или не связываете эти переменные. Инициализация по умолчанию в IPOPT - установить их все в 0, что приведет к ошибкам оценки.

Существует также ошибка в том, как вы применяете преобразование дискретизации. Должно быть: 

TransformationFactory('dae.finite_difference').apply_to(g, nfe=200, scheme='BACKWARD')
...