This is a challenging kinetic model (highly nonlinear) of polystyrene production. There are several strategies to help find a successful solution including:
- Good initial guess values such as
Mn=1e4 (number average molecular weight) and Mw=1e5 (weight average molecular weight)
- Solve as a steady-state simulation with
IMODE=1 to find good starting points for the dynamic simulation.
- Inspect
infeasibilities.txt to find equations or constraints that may be a problem for the solver.
- Reduce the problem down to only essential variables and equations and gradually add equations until the problem is initialized.
- Use variable bounds to limit the search space such as
lb=0 for concentrations.
Here is a version of the script that can help resolve the infeasible solution:
import numpy as np
import math
from gekko import GEKKO
import matplotlib.pyplot as plt
# define the model
m = GEKKO(remote=False) # create GEKKO model
# Defining the temperature in degrees kelvin.
T=(100+273);
# Pre-defining the kinetic parameters we have:
k_d=1.99e6*math.exp(-14842/T);
k_i=1.712e15*math.exp(-15924/T);
k_iterm=2.019e1*math.exp(-13810/T);
k_p=1.051e7*math.exp(-3577/T);
k_trM=2.31e-6*math.exp(-6377/T);
k_trS=1.8;
k_td=0.99*1.255e9*math.exp(-844/T);
k_tc=1.255e9*math.exp(-844/T);
M0=104.15; # molecular mass of Styrene monomer
f = 0.65; # Initiator efficiency
dm=0.909; # density of styrene monomer
dp=1.000; # density of Polystyrene
e=(dm-dp)/dp;
# create GEKKO variables
X = m.Var(0.1,lb=0,name='x') # monomer conversion
M = m.Var(70.0,lb=0,name='m') # monomer concentration
I = m.Var(30.0,lb=0,name='i') # initiator concentration
l_0 = m.Var(1,name='l_0')
l_1 = m.Var(1,name='l_1')
l_2 = m.Var(1,name='l_2')
u_0 = m.Var(1,name='u_0')
u_1 = m.Var(1,name='u_1')
u_2 = m.Var(1,name='u_2')
p = m.Var(30,lb=0,name='p')
Mn = m.Var(1e4,lb=0,name='mn')
Mw = m.Var(1e5,lb=0,name='mw')
PDI = m.Var(2,lb=0,name='pdi')
Mm=104.15;
# create GEKKO equations
m.Equation(p==(M/M0))
m.Equation((k_tc+k_p*e*p)*l_0==((2*f*k_d*I)))
m.Equation(((k_trM*M+k_tc*l_0+k_p*e*l_0*p)*l_1==\
2*f*k_d + k_p*M*l_0+k_trM*M*l_0))
m.Equation(M.dt()==-k_p*p*l_0*(1+e*p))
m.Equation(I.dt()==-k_d*I-k_p*I*l_0*e*p)
m.Equation((M0+e*M)*X==(M0-M))
m.Equation((u_0+l_0)*Mn==Mm*(u_1+l_1))
m.Equation((u_1+l_1)*Mw==Mm*(u_2+l_2))
m.Equation((k_i*l_0+k_p*e*p*l_0)*l_2==l_1+((2*k_p*M*l_1)))
m.Equation(u_1.dt()==k_trM*M*l_0+k_tc*l_0*l_1-k_p*u_1*l_0*e*p)
m.Equation(u_2.dt()==k_trM*M*l_2+k_tc*l_0*l_2+k_tc*l_1**2-k_p*u_2*l_0*e*p)
m.Equation((l_0**2)-u_0*l_0*e*k_p*p*u_0.dt()==\
k_trM*M*l_0+k_td*l_0**2+k_tc*0.5)
m.Equation(Mn*PDI==Mw)
try:
# solve steady-state first
m.options.IMODE=1
m.options.SOLVER=1
m.solve(disp=True)
# if steady-state is successful, try dynamic simulation
m.options.IMODE=4
#defining the simulation time
tm= np.linspace(0,20,100) # time points
m.time= tm
m.solve(disp=True)
except:
with open(m.path+'/infeasibilities.txt') as f:
print(f.read())
It currently produces the model gk0_model.apm in the run folder (open with m.open_folder()).
Model
Variables
x = 0.0, >= 0
m = 70.0, >= 0
i = 30.0, >= 0
l_0 = 0.01, >= 0
l_1 = 0.01, >= 0
l_2 = 0.01, >= 0
u_0 = 0.01, >= 0
u_1 = 0.01, >= 0
u_2 = 0.01, >= 0
p = 0.01, >= 0
mn = 0.01, >= 0
mw = 0.01, >= 0
pdi = 0.01, >= 0
End Variables
Equations
p=((m)/(104.15))
(((130602226.78465715+((-65.43973468411447)*(p))))*(l_0))=((1.354673906019174e-11)*(i))
((((((8.683403179277924e-14)*(m))+((130602226.78465715)*(l_0)))+((((-65.43973468411447)*(l_0)))*(p))))*(l_1))=((1.354673906019174e-11+((((719.1179635616977)*(m)))*(l_0)))+((((8.683403179277924e-14)*(m)))*(l_0)))
$m=((((((-719.1179635616977)*(p)))*(l_0)))*((1+((-0.09099999999999997)*(p)))))
$i=(((-1.04205685078398e-11)*(i))-((((((((719.1179635616977)*(i)))*(l_0)))*(-0.09099999999999997)))*(p)))
(((104.15+((-0.09099999999999997)*(m))))*(x))=(104.15-m)
(((u_0+l_0))*(mn))=((104.15)*((u_1+l_1)))
(((u_1+l_1))*(mw))=((104.15)*((u_2+l_2)))
(((((0.0004928772821909254)*(l_0))+((((-65.43973468411447)*(p)))*(l_0))))*(l_2))=(l_1+((((1438.2359271233954)*(m)))*(l_1)))
$u_1=((((((8.683403179277924e-14)*(m)))*(l_0))+((((130602226.78465715)*(l_0)))*(l_1)))-((((((((719.1179635616977)*(u_1)))*(l_0)))*(-0.09099999999999997)))*(p)))
$u_2=(((((((8.683403179277924e-14)*(m)))*(l_2))+((((130602226.78465715)*(l_0)))*(l_2)))+((130602226.78465715)*(((l_1)^(2)))))-((((((((719.1179635616977)*(u_2)))*(l_0)))*(-0.09099999999999997)))*(p)))
(((l_0)^(2))-((((((((((u_0)*(l_0)))*(-0.09099999999999997)))*(719.1179635616977)))*(p)))*($u_0)))=((((((8.683403179277924e-14)*(m)))*(l_0))+((129296204.51681055)*(((l_0)^(2)))))+65301113.392328575)
mw=((mn)*(pdi))
End Equations
End Model
The infeasibility report identifies two equations that cannot be solved with the current constraints and equations:
************************************************
***** POSSIBLE INFEASBILE EQUATIONS ************
************************************************
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
11 0.0000E+00 -4.3117E-03 0.0000E+00 4.3117E-03 ss.Eqn(11): 0 = $u_2-((((((((8.683403179277924e-14)*(m)))*(l_2))+((((130602226.78465715)*(l_0)))*(l_2)))+((130602226.78465715)*(((l_1)^(2)))))-((((((((719.1179635616977)*(u_2)))*(l_0)))*(-0.09099999999999997)))*(p))))
Variable Lower Value Upper $Value Name
2 0.0000E+00 0.0000E+00 1.2346E+20 0.0000E+00 ss.m
4 -1.2346E+20 0.0000E+00 1.2346E+20 0.0000E+00 ss.l_0
5 -1.2346E+20 5.7458E-06 1.2346E+20 0.0000E+00 ss.l_1
6 -1.2346E+20 -2.6903E+02 1.2346E+20 0.0000E+00 ss.l_2
9 -1.2346E+20 2.6903E+02 1.2346E+20 0.0000E+00 ss.u_2
10 0.0000E+00 0.0000E+00 1.2346E+20 0.0000E+00 ss.p
9 -1.2346E+20 2.6903E+02 1.2346E+20 0.0000E+00 ss.u_2
____________________________________________________________________________
EQ Number Lower Residual Upper Infeas. Name
12 0.0000E+00 -6.5301E+07 0.0000E+00 6.5301E+07 ss.Eqn(12): 0 = (((l_0)^(2))-((((((((((u_0)*(l_0)))*(-0.09099999999999997)))*(719.1179635616977)))*(p)))*($u_0)))-(((((((8.683403179277924e-14)*(m)))*(l_0))+((129296204.51681057)*(((l_0)^(2)))))+65301113.392328575))
Variable Lower Value Upper $Value Name
2 0.0000E+00 0.0000E+00 1.2346E+20 0.0000E+00 ss.m
4 -1.2346E+20 0.0000E+00 1.2346E+20 0.0000E+00 ss.l_0
7 -1.2346E+20 -4.3536E+02 1.2346E+20 0.0000E+00 ss.u_0
10 0.0000E+00 0.0000E+00 1.2346E+20 0.0000E+00 ss.p
7 -1.2346E+20 -4.3536E+02 1.2346E+20 0.0000E+00 ss.u_0
