Why does GEKKO not provide optimal commands even though the output does not match the reference?
The following is related to this question: predictive control model using GEKKO
I am trying to apply the MPC to maintain the temperature of a room within a defined range, but GEKKO gives me null commands even if the output diverges. I run the corrected code from my previous question:
# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
K = array([[ 0.93705481, -12.24012156]])
p = {'a': array([[ 1.08945247],
[-0.00242145],
[-0.00245978],
[-0.00272713],
[-0.00295845],
[-0.00319119],
[-0.00343511],
[-0.00366243],
[-0.00394247],
[-0.06665054]]),
'b': array([[[-0.05160235, -0.01039767],
[ 0.00173511, -0.01552485],
[ 0.00174602, -0.01179519],
[ 0.00180031, -0.01052658],
[ 0.00186416, -0.00822121],
[ 0.00193947, -0.00570905],
[ 0.00202877, -0.00344507],
[ 0.00211395, -0.00146947],
[ 0.00223514, 0.00021945],
[ 0.03800987, 0.04243736]]]),
'c': array([0.0265903])}
# i have used only 200 mes of T_externel
T_externel = np.linspace(9.51,9.78,200)
m = GEKKO(remote=False)
m.y = m.Array(m.CV,1)
m.u = m.Array(m.MV,2)
m.arx(p,m.y,m.u)
# rename CVs
m.T = m.y[0]
# rename MVs
m.beta = m.u[1]
# distrubance
m.d = m.u[0]
# distrubance and parametres
m.d = m.Param(T_externel[0])
# lower,heigh bound for MV
TL = m.Param(value = 16)
TH = m.Param(value = 18)
# steady state initialization
m.options.IMODE = 1
m.solve(disp=False)
# set up MPC
m.d.value = T_externel
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 2 # the objective is an l2-norm (squared error)
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 1 # APOPT
m.time = np.arange(0,len(T_externel)*300,300) # step time = 300s
# Manipulated variables
m.beta.STATUS = 1 # calculated by the optimizer
m.beta.FSTATUS = 1 # use measured value
m.beta.DMAX = 1.0 # Delta MV maximum step per horizon interval
m.beta.DCOST = 0.0 # Delta cost penalty for MV movement
m.beta.UPPER = 1.0 # Lower bound
m.beta.LOWER = 0.0
m.beta.MEAS = 0 # set u=0
# Controlled variables
m.T.STATUS = 1 # drive to set point
m.T.FSTATUS = 1 # receive measurement
m.T.SP = 17 # set point
TL.value = np.ones(len(T_externel))*16
TH.value = np.ones(len(T_externel))*18
m.T.value = 17 # Temprature starts at 17
for i in range(len(T_externel)):
m.solve(disp = False)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
import matplotlib.pyplot as plt
# Plot the results
plt.figure(figsize=(12,6))
plt.subplot(2,1,1)
plt.plot(m.time,m.T.value,'r-',label=r'$T_{int}$')
plt.plot(m.time,TL.value,'k--',label='Lower Bound')
plt.plot(m.time,TH.value,'k--',label='Upper Bound')
plt.ylabel('Temperature (°C)')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m.time,m.beta.value,'b--',label=r'$\beta$')
plt.ylabel('optimal control')
plt.xlabel('Time (sec)')
plt.legend()
plt.show()
And is the output and the optimal control obtained by GEKKO:
I need to to maintain the temperature of a room within a defined range.
Solution
The gain is listed as K = array([[ 0.93705481, -12.24012156]]) so an increase in beta by +1 leads to a decrease in the T by -12.24. An increase in T_ext by +1 leads to an increase in T by +0.937. The steady-state value of T is 13.32 so to get a starting value of 17, a common practice is to create an additive (or multiplicative bias value that adjusts the starting value. An additional variable Tb is added, although Gekko does this internally.
When investigating a controller response, a good first step is to confirm the model gain between the MV and CV with a step response. Here is a step response on beta.
# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
K = array([[ 0.93705481, -12.24012156]])
p = {'a': array([[ 1.08945247],
[-0.00242145],
[-0.00245978],
[-0.00272713],
[-0.00295845],
[-0.00319119],
[-0.00343511],
[-0.00366243],
[-0.00394247],
[-0.06665054]]),
'b': array([[[-0.05160235, -0.01039767],
[ 0.00173511, -0.01552485],
[ 0.00174602, -0.01179519],
[ 0.00180031, -0.01052658],
[ 0.00186416, -0.00822121],
[ 0.00193947, -0.00570905],
[ 0.00202877, -0.00344507],
[ 0.00211395, -0.00146947],
[ 0.00223514, 0.00021945],
[ 0.03800987, 0.04243736]]]),
'c': array([0.0265903])}
# i have used only 200 mes of T_externel
T_externel = np.linspace(9.51,9.78,200)
m = GEKKO(remote=True)
m.y = m.Array(m.CV,1)
m.u = m.Array(m.MV,2)
m.arx(p,m.y,m.u)
# rename CVs
m.T = m.y[0]
# rename MVs
m.beta = m.u[1]
# distrubance
m.d = m.u[0]
# distrubance and parametres
m.d = m.Param(T_externel[0])
m.bias = m.Param(0)
m.Tb = m.CV()
m.Equation(m.Tb==m.T+m.bias)
# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)
# set up MPC
#m.d.value = T_externel
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # the objective is an l1-norm (region)
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.time = np.arange(0,len(T_externel)*300,300) # step time = 300s
# Manipulated variables
m.beta.STATUS = 0 # calculated by the optimizer
m.beta.FSTATUS = 0 # use measured value
m.beta.DMAX = 0.2 # Delta MV maximum step per horizon interval
m.beta.DCOST = 0.0 # Delta cost penalty for MV movement
m.beta.UPPER = 1.0 # Upper bound
m.beta.LOWER = 0.0 # Lower bound
m.beta.MEAS = 0.0 # Measured value
# step test
m.beta.value = np.zeros_like(m.time)
m.beta.value[20:] = 1
# Controlled variables
m.Tb.STATUS = 1 # drive to set point
m.Tb.FSTATUS = 1 # receive measurement
m.Tb.SPHI = 21 # set point high level
m.Tb.SPLO = 19 # set point low level
T_MEAS = 17 # Temperature starts at 17
m.Tb.value = T_MEAS
m.bias.value = T_MEAS - m.T.value[0]
m.solve(disp=True)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
import matplotlib.pyplot as plt
# Plot the results
plt.figure(figsize=(8,3.5))
plt.subplot(2,1,1)
plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{biased}$')
plt.plot(m.time,m.T.value,'r--',label=r'$T_{unbiased}$')
plt.plot(m.time,m.d.value,'b:',label=r'T_{external}')
plt.plot([0,m.time[-1]],[m.Tb.SPHI,m.Tb.SPHI],'k--',label='Upper Bound')
plt.plot([0,m.time[-1]],[m.Tb.SPLO,m.Tb.SPLO],'k--',label='Lower Bound')
plt.ylabel('Temperature (°C)')
plt.legend(); plt.grid()
plt.subplot(2,1,2)
plt.step(m.time,m.beta.value,'b--',label=r'$\beta$')
plt.ylabel('optimal control')
plt.xlabel('Time (sec)')
plt.legend(); plt.grid()
plt.savefig('results.png',dpi=300)
plt.show()
The next step to investigate the controller response is to add the disturbance and set the appropriate set point range. The controller is switched to m.options.CV_TYPE=1 to use the l1-norm objective that gives m.Tb.SPHI and m.Tb.SPLO options to specify the range.
# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
K = array([[ 0.93705481, -12.24012156]])
p = {'a': array([[ 1.08945247],
[-0.00242145],
[-0.00245978],
[-0.00272713],
[-0.00295845],
[-0.00319119],
[-0.00343511],
[-0.00366243],
[-0.00394247],
[-0.06665054]]),
'b': array([[[-0.05160235, -0.01039767],
[ 0.00173511, -0.01552485],
[ 0.00174602, -0.01179519],
[ 0.00180031, -0.01052658],
[ 0.00186416, -0.00822121],
[ 0.00193947, -0.00570905],
[ 0.00202877, -0.00344507],
[ 0.00211395, -0.00146947],
[ 0.00223514, 0.00021945],
[ 0.03800987, 0.04243736]]]),
'c': array([0.0265903])}
# i have used only 200 mes of T_externel
T_externel = np.linspace(9.51,9.78,200)
m = GEKKO(remote=True)
m.y = m.Array(m.CV,1)
m.u = m.Array(m.MV,2)
m.arx(p,m.y,m.u)
# rename CVs
m.T = m.y[0]
# rename MVs
m.beta = m.u[1]
# distrubance
m.d = m.u[0]
# distrubance and parametres
m.d = m.Param(T_externel[0])
m.bias = m.Param(0)
m.Tb = m.CV()
m.Equation(m.Tb==m.T+m.bias)
# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)
# set up MPC
m.d.value = T_externel
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # the objective is an l1-norm (region)
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.time = np.arange(0,len(T_externel)*300,300) # step time = 300s
# Manipulated variables
m.beta.STATUS = 1 # calculated by the optimizer
m.beta.FSTATUS = 0 # use measured value
m.beta.DMAX = 1.0 # Delta MV maximum step per horizon interval
m.beta.DCOST = 0.0 # Delta cost penalty for MV movement
m.beta.UPPER = 1.0 # Upper bound
m.beta.LOWER = 0.0 # Lower bound
# Controlled variables
m.Tb.STATUS = 1 # drive to set point
m.Tb.FSTATUS = 0 # receive measurement
m.Tb.SPHI = 15 # set point high level
m.Tb.SPLO = 13 # set point low level
m.Tb.WSPHI = 100 # set point high priority
m.Tb.WSPLO = 100 # set point low priority
T_MEAS = 17 # Temperature starts at 17
m.Tb.value = T_MEAS
m.bias.value = T_MEAS - m.T.value[0]
m.solve(disp=True)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
import matplotlib.pyplot as plt
# Plot the results
plt.figure(figsize=(8,3.5))
plt.subplot(2,1,1)
plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{biased}$')
plt.plot(m.time,m.T.value,'r--',label=r'$T_{unbiased}$')
plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$')
plt.plot([0,m.time[-1]],[m.Tb.SPHI,m.Tb.SPHI],'k--',label='Upper Bound')
plt.plot([0,m.time[-1]],[m.Tb.SPLO,m.Tb.SPLO],'k--',label='Lower Bound')
plt.ylabel('Temperature (°C)')
plt.legend(loc=1); plt.grid()
plt.subplot(2,1,2)
plt.step(m.time,m.beta.value,'b--',label=r'$\beta$')
plt.ylabel('optimal control')
plt.xlabel('Time (sec)')
plt.legend(loc=1); plt.grid()
plt.savefig('results.png',dpi=300)
plt.show()
The controller keeps the Tb value within the setpoint range by adjusting beta.
It is also possible to make the beta value an integer if the cooling system only has an on/off state instead of continuous values between 0 and 1.
# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
K = array([[ 0.93705481, -12.24012156]])
p = {'a': array([[ 1.08945247],
[-0.00242145],
[-0.00245978],
[-0.00272713],
[-0.00295845],
[-0.00319119],
[-0.00343511],
[-0.00366243],
[-0.00394247],
[-0.06665054]]),
'b': array([[[-0.05160235, -0.01039767],
[ 0.00173511, -0.01552485],
[ 0.00174602, -0.01179519],
[ 0.00180031, -0.01052658],
[ 0.00186416, -0.00822121],
[ 0.00193947, -0.00570905],
[ 0.00202877, -0.00344507],
[ 0.00211395, -0.00146947],
[ 0.00223514, 0.00021945],
[ 0.03800987, 0.04243736]]]),
'c': array([0.0265903])}
# i have used only 200 mes of T_externel
T_externel = np.linspace(9.51,9.78,200)
m = GEKKO(remote=True)
m.y = m.Array(m.CV,1)
m.u = m.Array(m.MV,2)
m.arx(p,m.y,m.u)
# rename CVs
m.T = m.y[0]
# rename MVs
m.beta = m.u[1]
m.free(m.beta)
m.bint = m.MV(0,lb=0,ub=1,integer=True)
m.Equation(m.beta==m.bint)
# distrubance
m.d = m.u[0]
# distrubance and parametres
m.d = m.Param(T_externel[0])
m.bias = m.Param(0)
m.Tb = m.CV()
m.Equation(m.Tb==m.T+m.bias)
# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)
# set up MPC
m.d.value = T_externel
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # the objective is an l1-norm (region)
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.time = np.arange(0,len(T_externel)*300,300) # step time = 300s
# Manipulated variables
m.bint.STATUS = 1 # calculated by the optimizer
m.bint.FSTATUS = 0 # use measured value
m.bint.DCOST = 0.0 # Delta cost penalty for MV movement
m.bint.UPPER = 1.0 # Upper bound
m.bint.LOWER = 0.0 # Lower bound
m.bint.MV_STEP_HOR = 10
m.bint.value = 0
# Controlled variables
m.Tb.STATUS = 1 # drive to set point
m.Tb.FSTATUS = 0 # receive measurement
m.Tb.SPHI = 15 # set point high level
m.Tb.SPLO = 13 # set point low level
m.Tb.WSPHI = 100 # set point high priority
m.Tb.WSPLO = 100 # set point low priority
T_MEAS = 17 # Temperature starts at 17
m.Tb.value = T_MEAS
m.bias.value = T_MEAS - m.T.value[0]
m.options.SOLVER = 1
m.solve(disp=True)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
import matplotlib.pyplot as plt
# Plot the results
plt.figure(figsize=(8,3.5))
plt.subplot(2,1,1)
plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{biased}$')
plt.plot(m.time,m.T.value,'r--',label=r'$T_{unbiased}$')
plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$')
plt.plot([0,m.time[-1]],[m.Tb.SPHI,m.Tb.SPHI],'k--',label='Upper Bound')
plt.plot([0,m.time[-1]],[m.Tb.SPLO,m.Tb.SPLO],'k--',label='Lower Bound')
plt.ylabel('Temperature (°C)')
plt.legend(loc=1); plt.grid()
plt.subplot(2,1,2)
plt.step(m.time,m.bint.value,'b--',label=r'$\beta_{int}$')
plt.ylabel('optimal control')
plt.xlabel('Time (sec)')
plt.legend(loc=1); plt.grid()
plt.savefig('results.png',dpi=300)
plt.show()
Solution time for binary variables is significantly longer so many approximate a 0.37 as 37% ON with some post-processing such as off 3.7 minutes and on 6.3 minutes in 10-minute interval blocks.
Iter: 71 I: 0 Tm: 16.99 NLPi: 10 Dpth: 4 Lvs: 87 Obj: 3.04E+02 Gap: 6.24E-01
Iter: 72 I: 0 Tm: 11.37 NLPi: 3 Dpth: 4 Lvs: 88 Obj: 3.04E+02 Gap: 6.24E-01
--Integer Solution: 3.04E+02 Lowest Leaf: 3.04E+02 Gap: 9.61E-09
Iter: 73 I: 0 Tm: 13.62 NLPi: 3 Dpth: 13 Lvs: 88 Obj: 3.04E+02 Gap: 9.61E-09
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 1064.45290000000 sec
Objective : 304.259454634496
Successful solution
---------------------------------------------------
- AttributeError: install_layout when attempting to install a package in a virtual environment
- Python list comprehension - want to avoid repeated evaluation
- Hash algorithm for dynamic growing/streaming data?
- matplotlib - making labels for violin plots
- Python How to I check if last element has been reached in iterator tool chain?
- Polars and the Lazy API: How to drop columns that contain only null values?
- Why are my Mean, Var, and Std outputs from NumPy different from what the online grader expects?
- Correlation dataframe convertion from results from pl.corr
- Polars DataFrame transformation
- Discord rate limiting while only sending 1 request per minute
- Check if column contains (/,-,_, *or~) and split in another column - Pandas
- How to draw a rectangle at (x,y) in a PyQt GraphicsView?
- how to calculate correlation between ten columns with polars
- How to set class attribute with await in __init__
- Detect hindi encoding, response received from Facebook API in Python
- Is it possible to write a horizontal if statement with a multi-line body?
- Max length of items in list
- Cannot subclass multiprocessing Queue in Python 3.5
- How can I get notified of updates to Python packages in a unified way?
- Using python AST to traverse code and extract return statements
- merge groups of columns in a polars dataframe to single columns
- Group Pandas DataFrame by Continuous Date Ranges
- Flask login @login_required not working
- Odoo: one2many and many2one? KeyError:'___'
- merge some columns in a Polars dataframe and duplicate the others
- Python: Create table from string mixed with separators using FOR loops
- How do I type hint a method with the type of the enclosing class?
- How can I verify an emails DKIM signature in Python?
- Writing a class that accepts a callback in Python?
- Python Paramiko channel.exec_command not returning output intermittently