GEKKO: Does not respond to constraints nor solves the obj function

Question

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 during one day. Does this affect my controller?

I run the corrected code from my previous question but it does not seem to maintain the constraints or to change the MV over the day.

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 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
u = 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)
#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=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 costs.

m.Obj(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() 

#Calculate the costs.
c= 0
cost_list = []
for i in range(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(5,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(5,1,2)
plt.step(m.time,u.value,'b:')
plt.ylabel('Fridge State')
plt.legend()
plt.subplot(5,1,3)
plt.plot(m.time, eH.value, 'k--', label='Upper Tempratue Limit Error')
plt.plot(m.time, eL.value, 'b--', label='Lower Temprature Limit Error')
plt.ylabel('Cumulative Linar Error')
plt.legend()
plt.subplot(5,1,4)
plt.plot(m.time, cost_list, 'r-')
plt.ylabel('Costs in cent')

plt.show()

The results look like this :

I will appreciate any kind of help :)

score 2 · Answer 1 · answered Sep 05 '20 at 21:34

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 coefficients
print('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 in range(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()

The problem is I can only decide if the fridge is on or off, using u = {0,1,2,3,4} is not solving my problem. Plus, the changes in temperature do not seem logical to me. It basically says that if I Keep my fridge for about 6 hours off the temperature will not rise inside of it. Maybe I need to identify my model for the MPC mit 24 hours data ? — tt40kiwi, Sep 06 '20 at 16:45
Yes, it sounds like it is a model identification issue. You may want to start with a first order differential equation model such as tau dT/dt = -(T_ambient - T) + K * u where tau is the time constant (perhaps 30 min) and K is the compressor gain (perhaps -40). You can implement that equation in Gekko with m.Equation(30 * T.dt()==-(30-T)-40*u). You can also do step test simulations on your ARX model to make sure the identification results are reasonable. — John Hedengren, Sep 07 '20 at 12:59

GEKKO: Does not respond to constraints nor solves the obj function

1 Answers1

Linked