pythongekko

Consistent initialisation in simple DAE system using GEKKO

I'm messing around with GEKKO for dynamic simulations to see how it compares to commercial simulation software e.g. gPROMS.

As a simple test I'm attempting to model the blowdown of a pressurised vessel. However I cannot seem to get the solver to consistently initialise the simulation.

I have modelled the system using the following system of DAE's

dM/dt = -outlet_flowrate (1)
outlet_flowrate = f(Pressure,Outlet_Pressure) (2)
Pressure = f(M,T) (3)

Where Outlet_Pressure is fixed at a value, T is fixed at a value. I expect that I can specify the initial mass (M), the initial pressure could then be calculated using equation (3) , the initial outlet flowrate can then be calculated using equation (2).

I have used arbitrary function for f(Pressure,Outlet_pressure) and f(M,T) to simplify the model as a minimum reproducible example:

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

m = GEKKO(remote=False)
m.time = np.linspace(0,20,100)

mass = m.Var(value=2,lb=0,ub=100)
Outlet_flowrate = m.Var(lb=0,ub=10000000)
pressure = m.Var(lb=0,ub=1000000)

m.Equation(mass.dt() == -Outlet_flowrate)
m.Equation(Outlet_flowrate == (pressure-2))
m.Equation(pressure == mass*2)

m.free_initial(pressure)
m.free_initial(Outlet_flowrate)

m.options.IMODE=4
m.options.RTOL=1e-15
m.solve()

plt.plot(m.time,pressure.value)
print(pressure.value)
plt.show()

However when I attempt to solve this model I get the following output from the solver:

 --------- APM Model Size ------------
 Each time step contains
   Objects      :  0
   Constants    :  0
   Variables    :  3
   Intermediates:  0
   Connections  :  2
   Equations    :  3
   Residuals    :  3

 Number of state variables:    398
 Number of total equations: -  396
 Number of slack variables: -  0
 ---------------------------------------
 Degrees of freedom       :    2

 @error: Degrees of Freedom
 * Error: DOF must be zero for this mode
 STOPPING...

I have tried solving the model using IMODE=7 which does solve the model however this initialises the pressure at 0. I think the initial mass should be sufficient to specify the initial pressure at 4 not zero.

I am aware that I can trivially reduce the model to remove pressure and flowrate as variabes, however in the future I would like to apply this methodology to other systems which aren't as trivial to reduce. I'm also aware that I could specify the initial pressure by manually computing the initial pressure from the initial mass.

Is manually computing the initial condition something that most people do or have I missed a method in GEKKO to calculate the initial values of these variables ?

Thanks

Solution

  • The first node of a dynamic simulation is deactivated in Gekko because initial conditions are typically fixed. Gekko doesn't require consistent initial condition and can solve higher-index DAEs without Pantelides algorithm (differentiate higher-order algebraic equations to return to ODE or index-1 DAE form). There is no easy way to undo this default behavior. As an alternative, a small time step (e.g. 1e-10) can be placed at the beginning to calculate the correct initial condition of variables solved with algebraic equations.

    results

    from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

m = GEKKO(remote=False)
m.time = np.linspace(0,20,100)
m.time = np.insert(m.time,1,1e-10)

mass = m.Var(value=2,lb=0,ub=100)
Outlet_flowrate = m.Var(lb=0,ub=10000000)
pressure = m.Var(lb=0,ub=1000000)

m.Equation(mass.dt() == -Outlet_flowrate)
m.Equation(Outlet_flowrate == (pressure-2))
m.Equation(pressure == mass*2)

m.options.IMODE=4
m.options.SOLVER=1
m.solve()

plt.plot(m.time[1:],pressure.value[1:])
print(pressure.value)
plt.show()

    The first time point is removed from the results with:

    plt.plot(m.time[1:],pressure.value[1:])

    The initial conditions are removed and the results are numerically the same if the initial time step is small enough.