python - 使用 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。这也会在第一个 block 中产生不可行的解决方案。
• 通过重新排列方程式(例如 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') 则更具可读性。