用Gekko求解微分代数方程（DAE）问题

Question

我目前正在尝试求解一个由12个方程组成的系统，其中包括代数和微分表达式，如this article中所示。基本上，我试图在Python中复制文档中显示的结果，这些结果是使用上面链接的补充文件部分中提供的MATLAB代码获得的。我考虑的所有参数、常数和初始条件都取自这些文件。作为系统一部分的方程式如下： DAE equation system

为了尝试用Python实现解决方案，我在Gekko中设置了模型的核心，如下所示（尽可能简单，因为将来需要处理更多的方程和变量）

from gekko import GEKKO
import numpy as np

# Define constants and parameters of the system
# Will be updated in the future to work as variables or changing parameters
T = 298.0                   # [K] Temperature of the cell
Iapp = 1.7                  # [A] Single value for now
tMax  = 3600                # [s] Simulation time

# Physical constants
F = 9.64853399e4            # [sA/mol = C/mol] Faraday constant 
R = 8.314462175             # [J/(K*mol)] Gas constant
NA = 6.0221412927e23        # [1/mol] Avogadro number
e = 1.60217656535e-19       # [C] Electron charge
ne = 4.0                    # [] Electron number per reaction
ns8 = 8.0                   # [] Number of sulfur atoms in the polysulfide
ns4 = 4.0                   # [] Number of sulfur atoms in the polysulfide
ns2 = 2.0                   # [] Number of sulfur atoms in the polysulfide
ns = 1.0                    # [] Number of sulfur atoms in the polysulfide
Mm = 32.0                   # [g/mol] Molecular mass of Sulfur

# Cell parameters
Eh0 = 2.195                 # [V] Reference open circuit potential of high plateau (H)
El0 = 2.35                  # [V] Reference OCP of low plateau (L)
ih0 = 0.96                  # [A/m^2] Exchange current density in H
il0 = ih0/2                 # [A/m^2] Echange current density in L
Ar = 0.96                   # [m^2] Active reaction area of the cell
ms = 2.7                    # [g] Mass of active sulfur on the cell
velyt = 0.0114              # [L] Electrolyte volume in the cell
fh = ns4**2*Mm*velyt/ns8                 # [gL/mol] Dimensionality factor H
fl = ns**2*ns2*Mm**2*velyt**2/ns4        # [[g^2 L^2/mol]] Dimensionality factor L

# Precipitation and shuttle parameters
Vol = 0.0114e-3             # [m^3] Cell volume
RhoS = 2.0e6                # [g/m^3] Density of precipitated sulfur
Ksp = 0.0001                # [g] Saturation mass of sulfur in electrolyte
kp = 100.0                  # [1/s] Precipitation/dissolution rate constant
ks = 0.0002                 # [1/s] Shuttle rate constant

# REAL INITIAL (t = 0s) for discharge
S80 = 2.67300000000000          # [g] S80
S40 = 0.0128002860955374        # [g] S40
S20 = 4.33210229104915e-06      # [g] S20
S0 = 1.63210229104915e-06      # [g] S0
Sp0 = 2.70000000000000e-06      # [g] Sp0
Ih0 = 1.70000000000000          # [A] Ih0
Il0 = 0.00000               # [A] Il0
V0 = 2.40000000000000          # [V] V0
ETAh0 = -0.0102460242980059       # [V] ETAh0
ETAl0 = 0.00000               # [V] Etal0
Eh0 = 2.4102460242980059        # [V] Eh0
El0 = 2.40000000000000          # [V] El0

# Setup Gekko model 
Model0D = GEKKO()            # create GEKKO model
Model0D.time = np.linspace(0,tMax,tMax)

# Define the variables of the problem
S8 = Model0D.Var(value=S80, lb=1e-6, ub=2.7)
S4 = Model0D.Var(value=S40, lb=1e-6, ub=2.7)
S2 = Model0D.Var(value=S20, lb=1e-6, ub=2.7)
S = Model0D.Var(value=S0, lb=1e-6, ub=2.7)
Sp = Model0D.Var(value=Sp0, lb=1e-6, ub=2.7)
Ih = Model0D.Var(value=Ih0)
Il = Model0D.Var(value=Il0)
V = Model0D.Var(value=V0)
ETAh = Model0D.Var(value=ETAh0)
ETAl = Model0D.Var(value=ETAl0)
Eh = Model0D.Var(value=Eh0)
El = Model0D.Var(value=El0)

# Define all the equations of the system (DAE)
Model0D.Equation(S8.dt() == - Ih*(ns8*Mm)/(ne*F) - ks*S8)
Model0D.Equation(S4.dt() == Ih*(ns8*Mm)/(ne*F) + ks*S8 - Il*(ns4*Mm)/(ne*F))
Model0D.Equation(S2.dt() == Il*(ns2*Mm)/(ne*F))
Model0D.Equation(S.dt() == 2.0*Il*(ns*Mm)/(ne*F) - Sp*(kp/(Vol*RhoS))*(S-Ksp))
Model0D.Equation(Sp.dt() == Sp*(kp/(Vol*RhoS))*(S-Ksp))
Model0D.Equation(I == Ih + Il)
Model0D.Equation(Ih == 2.0*ih0*Ar*Model0D.sinh(ne*F*ETAh/(2.0*R*T)))
Model0D.Equation(Il == 2.0*il0*Ar*Model0D.sinh(ne*F*ETAl/(2.0*R*T)))
Model0D.Equation(ETAh == V - Eh)
Model0D.Equation(ETAl == V - El)
Model0D.Equation(Eh == Eh0 + (R*T/(4.0*F))*Model0D.log(fh*S8/S4**2.0))
Model0D.Equation(El == El0 + (R*T/(4.0*F))*Model0D.log(fl*S4/(S2*S**2.0)))

# Set model options
Model0D.options.IMODE = 7               # Set model type (simulation -> 4-simultaneous 7-sequential)
# Solve model
Model0D.solve(disp=True)

看到我向您展示的代码块没有收敛到任何解决方案，我尝试调整一些从文档中可以理解的解算器和模型相关选项（顺序/同步、解算器选择等），但不能说它是否正在改变过程，以收敛到一个合理的解决方案

我确实意识到，对于一组相对复杂的方程，我提出的模型可能过于简单，以至于解算器实际上无法收敛。由于我对Gekko和Python本身非常陌生，而且我在网上没有找到任何有用的信息，因此我正在尝试联系任何对Gekko工具有一定经验的人，他们可以为我提供一些关于如何有效解决DAE问题或问题本身解决方案的指南，以便我可以与我已经进行的失败尝试进行比较

Answer 1

在paper中有关于初始化策略的其他信息：

Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016.

以下是一些可以尝试的东西：

• 从一个非常小的时间步开始，比如m.time = np.linspace(0,1e-3,2)来验证一个成功的解决方案
• 切换到m.options.IMODE=4和m.options.COLDSTART=2。这也会在第一个块中产生不可行的解决方案
• 通过重新排列方程，如x.dt()==1/y到y*x.dt()==1，消除被零除
• 脚本中未定义I的值
• 尝试为所有变量设置一个零的下限，尤其是那些物理上不应该为负的变量
• 尝试删除其他非物理边界（上限为2.7？），看看问题是否可行

from gekko import GEKKO
import numpy as np

# Define constants and parameters of the system
# Will be updated in the future to work as variables or changing parameters
T = 298.0                   # [K] Temperature of the cell
Iapp = 1.7                  # [A] Single value for now
tMax  = 3600                # [s] Simulation time
I = Iapp

# Physical constants
F = 9.64853399e4            # [sA/mol = C/mol] Faraday constant 
R = 8.314462175             # [J/(K*mol)] Gas constant
NA = 6.0221412927e23        # [1/mol] Avogadro number
e = 1.60217656535e-19       # [C] Electron charge
ne = 4.0                    # [] Electron number per reaction
ns8 = 8.0                   # [] Number of sulfur atoms in the polysulfide
ns4 = 4.0                   # [] Number of sulfur atoms in the polysulfide
ns2 = 2.0                   # [] Number of sulfur atoms in the polysulfide
ns = 1.0                    # [] Number of sulfur atoms in the polysulfide
Mm = 32.0                   # [g/mol] Molecular mass of Sulfur

# Cell parameters
Eh0 = 2.195                 # [V] Reference open circuit potential of high plateau (H)
El0 = 2.35                  # [V] Reference OCP of low plateau (L)
ih0 = 0.96                  # [A/m^2] Exchange current density in H
il0 = ih0/2                 # [A/m^2] Echange current density in L
Ar = 0.96                   # [m^2] Active reaction area of the cell
ms = 2.7                    # [g] Mass of active sulfur on the cell
velyt = 0.0114              # [L] Electrolyte volume in the cell
fh = ns4**2*Mm*velyt/ns8                 # [gL/mol] Dimensionality factor H
fl = ns**2*ns2*Mm**2*velyt**2/ns4        # [[g^2 L^2/mol]] Dimensionality factor L

# Precipitation and shuttle parameters
Vol = 0.0114e-3             # [m^3] Cell volume
RhoS = 2.0e6                # [g/m^3] Density of precipitated sulfur
Ksp = 0.0001                # [g] Saturation mass of sulfur in electrolyte
kp = 100.0                  # [1/s] Precipitation/dissolution rate constant
ks = 0.0002                 # [1/s] Shuttle rate constant

# REAL INITIAL (t = 0s) for discharge
S80 = 2.67300000000000          # [g] S80
S40 = 0.0128002860955374        # [g] S40
S20 = 4.33210229104915e-06      # [g] S20
S0 = 1.63210229104915e-06      # [g] S0
Sp0 = 2.70000000000000e-06      # [g] Sp0
Ih0 = 1.70000000000000          # [A] Ih0
Il0 = 0.00000               # [A] Il0
V0 = 2.40000000000000          # [V] V0
ETAh0 = -0.0102460242980059       # [V] ETAh0
ETAl0 = 0.00000               # [V] Etal0
Eh0 = 2.4102460242980059        # [V] Eh0
El0 = 2.40000000000000          # [V] El0

# Setup Gekko model 
Model0D = GEKKO(remote=False)            # create GEKKO model
Model0D.time = np.linspace(0,1e-3,2) #tMax,tMax)

# Define the variables of the problem
S8 = Model0D.Var(value=S80)
S4 = Model0D.Var(value=S40)
S2 = Model0D.Var(value=S20)
S = Model0D.Var(value=S0)
Sp = Model0D.Var(value=Sp0)
Ih = Model0D.Var(value=Ih0)
Il = Model0D.Var(value=Il0)
V = Model0D.Var(value=V0)
ETAh = Model0D.Var(value=ETAh0)
ETAl = Model0D.Var(value=ETAl0)
Eh = Model0D.Var(value=Eh0)
El = Model0D.Var(value=El0)

# Define all the equations of the system (DAE)
Model0D.Equation(S8.dt() == - Ih*(ns8*Mm)/(ne*F) - ks*S8)
Model0D.Equation(S4.dt() == Ih*(ns8*Mm)/(ne*F) + ks*S8 - Il*(ns4*Mm)/(ne*F))
Model0D.Equation(S2.dt() == Il*(ns2*Mm)/(ne*F))
Model0D.Equation(S.dt() == 2.0*Il*(ns*Mm)/(ne*F) - Sp*(kp/(Vol*RhoS))*(S-Ksp))
Model0D.Equation(Sp.dt() == Sp*(kp/(Vol*RhoS))*(S-Ksp))
Model0D.Equation(I == Ih + Il)
Model0D.Equation(Ih == 2.0*ih0*Ar*Model0D.sinh(ne*F*ETAh/(2.0*R*T)))
Model0D.Equation(Il == 2.0*il0*Ar*Model0D.sinh(ne*F*ETAl/(2.0*R*T)))
Model0D.Equation(ETAh == V - Eh)
Model0D.Equation(ETAl == V - El)
Model0D.Equation(Eh == Eh0 + (R*T/(4.0*F))*Model0D.log(fh*S8/S4**2.0))
Model0D.Equation(El == El0 + (R*T/(4.0*F))*Model0D.log(fl*S4/(S2*S**2.0)))

# Set model options
Model0D.options.IMODE = 4 # Set model type (simulation -> 4-simultaneous 7-sequential)
Model0D.options.COLDSTART=2
Model0D.options.SOLVER=1
# Solve model
Model0D.open_folder()
Model0D.solve(disp=True)

这并没有使问题变得可行，但可能有助于搜索。infeasibilities.txt文件也有帮助。当使用运行目录m.path中的remote=False进行求解时，或通过使用m.open_folder()打开运行目录（在m.solve()命令之前）时，此选项可用。如果将变量命名为x = m.Var(name='x')，则可读性更高