Question: 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

Answer: 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.  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.