How do I change the MV status of a GEKKO variable after solving for steady state?

Question

I am trying to find the value for an inlet flow rate that keeps a gravity drained tank at a required level. After the solution is found using IMODE 3, I want to simulate the system at this state for a period of time using IMODE 4. Below is the code I am using:

from gekko import GEKKO
import numpy as np
m = GEKKO(remote = False)

F_in = m.MV(value = 1, name= 'F_in')
level = m.CV(value = 1, name = 'level')
F_out = m.Var(value = 1, name = 'F_out')

level.STATUS = 1
level.SPHI = 6
level.SPLO = 4

m.Equation(level.dt() == (F_in - F_out)/5)
m.Equation(F_out == 2*(level)**0.5)

# Find the steady state.

F_in.STATUS = 1

m.options.IMODE = 3
m.options.SOLVER = 3
m.solve(disp = False)
print(F_in.value)

# Run at steady state for 5 time units.

F_in.STATUS = 0 # Turn STATUS off for 0 DOF.
m.time = np.linspace(0,5,10)
m.options.IMODE = 4
m.solve(disp = False)
print(level.value)

However, the code produces the error stating that IMODE 4 requires 0 DOF. Does changing the F_in status to 0 not result in 0 DOF? Below is the full error message:

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In[4], line 27
     25 m.time = np.linspace(0,5,10)
     26 m.options.IMODE = 4
---> 27 m.solve(disp = False)
     28 print(level.value)

File ~\VRMEng\VRMEng\Lib\site-packages\gekko\gekko.py:2140, in GEKKO.solve(self, disp, debug, GUI, **kwargs)
   2138         print("Error:", errs)
   2139     if (debug >= 1) and record_error:
-> 2140         raise Exception(apm_error)
   2142 else: #solve on APM server
   2143     def send_if_exists(extension):

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

Any help will be appreciated greatly!

Answer 1

Set m.options.SPECS = 0 to ignore the calculated / fixed specifications when changing modes (see documentation). Here is a complete script that solves successfully.

from gekko import GEKKO
import numpy as np
m = GEKKO(remote = False)
F_in = m.MV(value = 1, name= 'F_in')
level = m.CV(value = 1, name = 'level')
F_out = m.Var(value = 1, name = 'F_out')
level.STATUS = 1
level.SPHI = 6
level.SPLO = 4
m.Equation(level.dt() == (F_in - F_out)/5)
m.Equation(F_out == 2*(level)**0.5)
# Find the steady state.
F_in.STATUS = 1
m.options.IMODE = 3
m.options.SOLVER = 3
m.solve(disp = False)
print(level.value)
print(F_in.value)
print(F_out.value)

# Run at steady state for 5 time units.
m.options.SPECS = 0
F_in.value = F_in.value.value
F_in.STATUS = 0 # Turn STATUS off for 0 DOF.
m.time = np.linspace(0,5,10)
m.options.IMODE = 4
m.solve(disp = True)
print(level.value)
print(F_in.value)
print(F_out.value)

A few more details: The t0 file in the run directory m._path stores whether a variable should be calculated or fixed. The rto.t0 file is generated with IMODE=3 and stores the steady-state values and the fixed/calculated information. Setting SPECS=0 ignores the specifications from the t0 file and uses the default specifications. For IMODE=4 the STATUS for each of the FVs and MVs are over-written to turn off. Set F_in.value = F_in.value.value to retain the F_in value from the steady state simulation.