Gekko: Does Not Respond To Constraints Nor Solves The Obj Function

The following is related to this question :MPC with ARX Model Using Gekko. I am trying to identify my system with data of 15 minutes. And I Am trying to update my MPC MV every hour

Solution 1:

You need to define u = m.MV() and T=m.CV() before calling the m.arx() model so that these values will be used as inputs and outputs. I also increased the WSPHI value so that the cost objective does not cause the temperature limit to be ignored. The current refrigeration system appears to be insufficient to cool to this level. It needs a system that is about 3 times more powerful to maintain the temperature limit. I set the upper bound to the refrigeration system to be 4 so that it could maintain the temperature in the limits. It gives up on temperature control at the end because it finds that the energy savings is more valuable than meeting the temperature limit so a small period of time. You can enforce the limit by either increasing WSPHI and WSPLO or else with TH.UPPER = 0 as a hard constraint. The hard constraint may lead to an infeasible solution if the refrigeration system can't meet that constraint.

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 for 10 min 
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 coefficientsprint('Identification')
yp,p,K = m.sysid(t,ud,y,na,nb,objf=10000,scale=False,diaglevel=1)
#create control ARX model:# Controlled variable:
T = m.CV()
# Manipulated variable:
u = m.MV(value=0,lb=0, ub=4, integer=True)
# Create ARX Model
m.arx(p,T,u)

############################################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])

u.STATUS = 1# allow optimizer to change the variable to attein the optimum.# Controlled Variable (Affected with changes in the manipulated variable)#T = m.CV()# 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=100000#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=100000 
eL.STATUS = 1#Linear error (Deviation from the limits)
m.Equations([eH==T-TH,eL==T-TL])

#Objective: minimize 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#Set Up MPC
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()
m.solve() 

#Calculate the costs.
c= 0
cost_list = []
for i inrange(0,len(u)):
    c = c + TOU_v[i]*u[i]
    cost_list.append(c)
print('The daily energy cost is' ,c/100, 'Euro') 

plt.subplot(4,1,1)
plt.plot(m.time,temp_low,'k--', label='Lower limit')
plt.plot(m.time,temp_upper,'k--',label='Upper limit')
plt.plot(m.time,T.value,'r-')
plt.ylabel('Temperature')
plt.legend()
plt.subplot(4,1,2)
plt.step(m.time,u.value,'b:',label='u')
plt.ylabel('Fridge State')
#plt.grid()
plt.legend()
plt.subplot(4,1,3)
plt.plot(m.time, eH.value, 'k--', label='Upper Temperatue Limit Error')
plt.plot(m.time, eL.value, 'b--', label='Lower Temperature Limit Error')
plt.ylabel('Cumulative Linear Error')
plt.legend()
plt.subplot(4,1,4)
plt.plot(m.time, cost_list, 'r-')
plt.ylabel('Costs in cent')

plt.show()