GEKKO RTO vs MPC MODES

Question

This is a question derived from this one. After posting my question I found a solution (more like a patch to force the optimizer to optimize). There is something that baffles me. John Hedengren correctly points out that a b=1.0 in the ODE leads to an infeasible solution with IMODE=6. However in my patchy work around with IMODE=3 I do get a solution.

I'm trying to understand what's happening here reading GEKKO's documentation for IMODE=3 and 6 but it is not clear to me

IMODE=3

RTO Real-Time Optimization (RTO) is a steady-state mode that allows decision variables (FV or MV types with STATUS=1) or additional variables in excess of the number of equations. An objective function guides the selection of the additional variables to select the optimal feasible solution. RTO is the default mode for Gekko if m.options.IMODE is not specified.

IMODE=6

MPC Model Predictive Control (MPC) is implemented with IMODE=6 as a simultaneous solution or with IMODE=9 as a sequential shooting method.

Why b=1. works in one mode but not in the other?

This is my patchy work around with IMODE=3 and b=1.0:

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

m = GEKKO(remote=False)
m.time = np.linspace(0,23,24)

#initialize variables
T_e = [50.,50.,50.,50.,45.,45.,45.,60.,60.,63.,\
       64.,45.,45.,50.,52.,53.,53.,54.,54.,53.,52.,51.,50.,45.]
temp_low = [55.,55.,55.,55.,55.,55.,55.,68.,68.,68.,68.,55.,55.,68.,\
            68.,68.,68.,55.,55.,55.,55.,55.,55.,55.]
temp_upper = [75.,75.,75.,75.,75.,75.,75.,70.,70.,70.,70.,75.,75.,\
              70.,70.,70.,70.,75.,75.,75.,75.,75.,75.,75.]
TOU = [0.05,0.05,0.05,0.05,0.05,0.05,0.05,200.,200.,\
       200.,200.,200.,200.,200.,200.,200.,200.,200.,\
       200.,200.,200.,0.05,0.05,0.05]

b = m.Param(value=1.)
k = m.Param(value=0.05)

u = [m.MV(0.,lb=0.,ub=1.) for i in range(24)]

# Controlled Variable
T = [m.SV(60.,lb=temp_low[i],ub=temp_upper[i]) for i in range(24)]

for i in range(24):
    u[i].STATUS = 1

for i in range(23):
    m.Equation( T[i+1]-T[i]-k*(T_e[i]-T[i])-b*u[i]==0.0 )

m.Obj(np.dot(TOU,u))

m.options.IMODE = 3
m.solve(debug=True)
myu =[u[0:][i][0] for i in range(24)]
print myu
myt =[T[0:][i][0] for i in range(24)]
plt.plot(myt)
plt.plot(temp_low)
plt.plot(temp_upper)
plt.show()
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax1.plot(myu,color='b')
ax2.plot(TOU,color='k')
plt.show()

Results:

Answer 1

The difference between the infeasible IMODE=6 and the feasible IMODE=3 is that the IMODE=3 case allows the temperature initial condition to be adjusted by the optimizer. The optimizer recognizes that the initial condition can be changed and so it modifies it to 75 to both stay feasible and also minimize the future energy consumption.

from gekko import GEKKO
import numpy as np

m = GEKKO(remote=False)
m.time = np.linspace(0,23,24)

#initialize variables
T_external = [50.,50.,50.,50.,45.,45.,45.,60.,60.,63.,\
              64.,45.,45.,50.,52.,53.,53.,54.,54.,\
              53.,52.,51.,50.,45.]
temp_low = [55.,55.,55.,55.,55.,55.,55.,68.,68.,68.,68.,\
            55.,55.,68.,68.,68.,68.,55.,55.,55.,55.,55.,55.,55.]
temp_upper = [75.,75.,75.,75.,75.,75.,75.,70.,70.,70.,70.,75.,\
              75.,70.,70.,70.,70.,75.,75.,75.,75.,75.,75.,75.]
TOU_v = [0.05,0.05,0.05,0.05,0.05,0.05,0.05,200.,200.,200.,200.,\
         200.,200.,200.,200.,200.,200.,200.,200.,200.,200.,0.05,\
         0.05,0.05]

b = m.Param(value=1.)
k = m.Param(value=0.05)
T_e = m.Param(value=T_external)
TL = m.Param(value=temp_low)
TH = m.Param(value=temp_upper)
TOU = m.Param(value=TOU_v)

u = m.MV(lb=0, ub=1)
u.STATUS = 1  # allow optimizer to change

# Controlled Variable
T = m.SV(value=75)

m.Equations([T>=TL,T<=TH])
m.Equation(T.dt() == k*(T_e-T) + b*u)

m.Minimize(TOU*u)

m.options.IMODE = 6
m.solve(disp=True,debug=True)

import matplotlib.pyplot as plt
plt.subplot(2,1,1)
plt.plot(m.time,temp_low,'k--')
plt.plot(m.time,temp_upper,'k--')
plt.plot(m.time,T.value,'r-')
plt.ylabel('Temperature')
plt.subplot(2,1,2)
plt.step(m.time,u.value,'b:')
plt.ylabel('Heater')
plt.xlabel('Time (hr)')
plt.show()

If you went for another day (48 hours), you'd probably see that the problem would eventually be infeasible because the smaller heater b=1 wouldn't be able to meet the temperature lower constraint.

One of the advantages of using IMODE=6 is that you can write the differential equation instead of doing the discretization yourself. With IMODE=3, you use an Euler's method for the differential equation. The default discretization for IMODE>=4 is NODES=2, equivalent to your Euler finite difference method. Setting NODES=3-6 increases the accuracy with orthogonal collocation on finite elements.