'Solving a differential algebraic equation (DAE) problem with Gekko
I'm currently trying to solve a system of 12 equations that include algebraic and differential expressions as it is presented in this article. Basically I'm trying to replicate in Python the results presented on the document, which are obtained using the MATLAB code provided in the supplementary files section of the link above. All the parameters, constants and initial conditions I'm considering are taken from these files. The equations that are part of the system are the next ones: DAE equation system
To try to reach the solution in Python, I set up the core of the model in Gekko as following (as simple as possible as it will need to handle more equations and variables in the future).
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)
Seeing that the block of the code that I showed you wasn't converging to any solution, I tried to tweak some of the solver and model related options (sequential/simultaneous, solver election, etc.) that I was able to comprehend from the documentation, but can't really say if it is making a change on the process to converge to a reasonable solution.
I do realise that for a relatively complex set of equations, the model I'm proposing may result too simplistic for the solver to actually be able to converge. As I'm pretty new to Gekko and Python itself, and I didn't really found any useful information online, I'm trying to reach anyone with some experience on the Gekko tool that can provide me some guidelines on how to solve DAE problems efficiently or the solution of the problem itself, so I can compare with the unsuccessful trials I already made.
Solution 1:[1]
There is additional information on initialization strategies in the 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.
Here are a couple things to try:
- Start with a very small time step such as
m.time = np.linspace(0,1e-3,2)to verify a successful solution. - Switch to
m.options.IMODE=4andm.options.COLDSTART=2. This also produces an infeasible solution in the first block. - Eliminate divide-by-zero by rearranging the equations such as
x.dt()==1/ytoy*x.dt()==1. - The value of
Iis undefined in the script. - Try setting a lower bound of zero for all variables, especially those that shouldn't physically be negative.
- Try removing other non-physical bounds (upper bounds of 2.7?) to see if the problem becomes feasible.
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)
This didn't make the problem feasible but may help with the search. The infeasibilities.txt file can also help. It is available when solving with remote=False in the run directory m.path or by opening the run directory with m.open_folder() (before the m.solve() command). It is more readable if the variables are named such as x = m.Var(name='x').
Solution 2:[2]
I have tried to use this model as a subject to learn the GEKKO ().
I advise on using I=Model0D.intermediate() instead of;
Model0D.Equation(ETAh == V - Eh)
Model0D.Equation(ETAl == V - El)
Also the time step are very large, creating unstable solution.
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|---|
| Solution 1 | John Hedengren |
| Solution 2 |