Why doesn't the Gekko solver adapt to variations in the system?
The following is related to this question : predictive control model using GEKKO
I am still traying applying MPC to maintain the temperature of a room within a defined range(16,18), as Professor @John explained last time.the gain is normally listed as K=array([[0.93705481,−12.24012156]]). Thus, an increase in β by +1 leads to a decrease in T by -12.24, while an increase in Text by +1 leads to an increase in T by +0.937.
In my case, I tried to implement an external temperature profile for a day (sinusoidal), and I modified the control to vary between 0 and 1, not just 0 or 1, since the equipment can operate within this range. I also adjusted m.bint.MV_STEP_HOR=1 so that the control is updated at each iteration. However, the control still does not react! Even when the external temperature increases or decreases, the control remains unchanged. 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
total_minutes_in_day = 24 * 60
interval = 5 # minutes
num_points = total_minutes_in_day // interval
# Générer les points de temps
time_points = np.linspace(0, total_minutes_in_day, num_points)
temperature = 8.5 + 11.5 * np.sin((2 * np.pi / total_minutes_in_day) * time_points - (np.pi / 2))
temperature = np.clip(temperature, -3, 20)
plt.plot(time_points, temperature)
plt.xlabel('Temps (minutes)')
plt.ylabel('Température (°C)')
plt.title('Variation de Tex sur une journée')
plt.savefig('Tex pour une jrn')
plt.grid(True)
plt.show()
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 = temperature
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]
m.free(m.beta)
m.bint = m.MV(0,lb=0,ub=1)
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=False)
# 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 = 1
m.bint.value = 0
# Controlled variables
m.Tb.STATUS = 1 # drive to set point
m.Tb.FSTATUS = 0 # receive measurement
m.Tb.SPHI = 18 # set point high level
m.Tb.SPLO = 16 # set point low level
m.Tb.WSPHI = 100 # set point high priority
m.Tb.WSPLO = 100 # set point low priority
T_MEAS = 20 # 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=False)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
Solution
The solution that you showed was able to find a solution with a single step to keep within the temperature range (17-19). Tightening the temperature range (16.5-17.5) shows additional movement of the MV over the time period and shows that it is responsive. Unless there is a violation of this temperature range, the controller does not take action if there is even a small DCOST value to move the MV.
The dynamics influence the speed of the response. A better insulated building will have a slower response to changes in the external temperature. One other suggestion is that the bint Manipulated Variable isn't needed when solving for continuous instead of integer solutions.
Here is the complete code with bint removed and the tighter temperature target range:
# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
import matplotlib.pyplot as plt
total_minutes_in_day = 24 * 60
interval = 5 # minutes
num_points = total_minutes_in_day // interval
# Générer les points de temps
time_points = np.linspace(0, total_minutes_in_day, num_points)
temperature = 8.5 + 11.5 * np.sin((2 * np.pi / total_minutes_in_day) * time_points - (np.pi / 2))
temperature = np.clip(temperature, -3, 20)
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])}
T_externel = temperature
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])
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=False)
# 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.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.MV_STEP_HOR = 1
m.beta.value = 0
# Controlled variables
m.Tb.STATUS = 1 # drive to set point
m.Tb.FSTATUS = 0 # receive measurement
m.Tb.SPHI = 17.5 # set point high level
m.Tb.SPLO = 16.5 # set point low level
m.Tb.WSPHI = 100 # set point high priority
m.Tb.WSPLO = 100 # set point low priority
T_MEAS = 20
m.Tb.value = T_MEAS
m.bias.value = T_MEAS - m.T.value[0]
m.options.SOLVER = 3
m.solve(disp=False)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
# 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()
- 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