MPC With ARX Model Using Gekko

I am modelling an MPC to control a fridge and keep the temperature within a given interval while minimizing the cost. I am using GEKKO to model my algorithm. I wrote the following

Solution 1:

The problem is with:

T = m.CV(value=11) # Temperature will start at 11.

You are redefining the T variable but it stores both internally. If you need to re-initialize to 11 then use T.value=11. Also, I added the eH and eL variables before the steady state initialization. Here is a complete script that runs successfully.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO(remote = True)


#initialize variables

#Room Temprature:
T_external = [23,23,23,23,23.5,23.5,23.4,23.5,23.9,23.7,\
              23,23.9,23.9,23.4,23.9,24,23.6,23.7,23.8,\
              23,23,23,23,23]

# Temprature Lower Limit:
temp_low = 10*np.ones(24)

# Temprature Upper Limit:
temp_upper = 12*np.ones(24)

#Hourly Energy prices:
TOU_v = [39.09,34.93,38.39,40.46,40.57,43.93,25,11,9,24,51.28,45.22,45.72,\
            36,35.03,10,12,13,32.81,42.55,8,29.58,29.52,29.52]

###########################################
#System Identification:

#Time 
t = np.linspace(0,10,117)
#State of the Fridge
ud = np.append(np.zeros(78) ,np.ones(39),0)
#Temprature Data
y = [14.600000000000001,14.600000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
     14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
     14.700000000000001,14.700000000000001,14.700000000000001,14.8,14.8,14.8,14.8,14.8,14.8,14.8,14.8,\
    14.8,14.8,14.9,14.9,14.9,14.9,14.9,14.9,14.9,15,15,15,15,15,15,15,15,15,15,15,15,15.100000000000001,\
    15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
    15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
    15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
    15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
    15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
    15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
    15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,
    15,15,15,15,15,15,15,15,15,15,14.9,14.9,14.9,14.9,14.8,14.9,14.8,14.8,14.8,14.8,14.8,14.8,\
    14.8,14.700000000000001,14.8,14.700000000000001,14.700000000000001,14.700000000000001,\
    14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
    14.700000000000001,14.600000000000001,14.600000000000001,14.600000000000001,\
    14.600000000000001,14.600000000000001,14.60]

na = 1 # output coefficients
nb = 1 # input coefficients
print('Identification')
yp,p,K = m.sysid(t,ud,y,na,nb,objf=10000,scale=False,diaglevel=1)
#create control ARX model:

y = m.Array(m.CV,1)
uc = m.Array(m.MV,1)
m.arx(p,y,uc)
# rename CVs
T= y[0]

# rename MVs
uc = uc[0]


###########################################

#Parameter
P = m.Param(value =100) #power
TL = m.Param(value=temp_low[0]) 
TH = m.Param(value=temp_upper[0])
c = m.Param(value=TOU_v[0])
# Manipilated variable:

u = m.MV(lb=0, ub=1, integer=True)
u.STATUS = 1  # allow optimizer to change the variable to attein the optimum.

# Controlled Variable (Affected with changes in the manipulated variable)

# Soft constraints on temprature.

eH = m.CV(value=0)
eL = m.CV(value=0)

eH.SPHI=0       #Set point high for linear error model.
eH.WSPHI=100    #Objective function weight on upper set point for linear error model.
eH.WSPLO=0      # Objective function weight on lower set point for linear error model
eH.STATUS =1    # eH : Error is considered in the objective function.
eL.SPLO=0
eL.WSPHI=0
eL.WSPLO=100 
eL.STATUS = 1   
#Linear error (Deviation from the limits)
m.Equations([eH==T-TH,eL==T-TL])

#Objective : minimize the costs.

m.Minimize(c*P*u)

#Optimizer Options.

# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)

TL.value = temp_low
TH.value = temp_upper
c.value  = TOU_v
T.value = 11 # Temprature starts at 11

m.options.IMODE = 6 # MPC mode in Gekko.
m.options.NODES = 2  # Collocation nodes.
m.options.SOLVER = 1 # APOT solver for mixed integer linear programming.
m.time = np.linspace(0,23,24)

#Solve the optimization problem.

m.solve() 

Here is the controller output:

 --------- APM Model Size ------------
 Each time step contains
   Objects      :            1
   Constants    :            0
   Variables    :            9
   Intermediates:            0
   Connections  :            2
   Equations    :            3
   Residuals    :            3
 
 Number of state variables:           1035
 Number of total equations: -         1012
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :             23
 
 ----------------------------------------------
 Dynamic Control with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.07 NLPi:    3 Dpth:    0 Lvs:    0 Obj:  6.76E+03 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   8.319999999366701E-002 sec
 Objective      :    6763.77971670735     
 Successful solution
 ---------------------------------------------------