predictive control model using GEKKO
I am modeling an MPC to maintain the temperature in a building within a given interval while minimizing the energy consumption. I am using GEKKO to model my algorithm.
I wrote the following code. First, I identified my model using data with the input (the disturbance: external temperature and the control), and the output y , which is the temperature. Then, I built an ARX model (using the arx function in GEKKO. This is my code :
# Import library
import numpy as np
import pandas as pd
import time
# Initialize Model
ts = 300
t = np.arange(0,len(data_1)*ts, ts)
u_id = data_1[['T_ext','beta']]
y_id = data_1[['T_int']]
# system identification
#meas : the time-series next step is predicted from prior measurements as in ARX
na=5; nb=5 # ARX coefficients
print('Identify model')
start = time.time()
yp,p,K = m.sysid(t,u_id,y_id,na,nb,objf=100,scale=False,diaglevel=0,pred='meas')
print('temps de prediction :'+str(time.time()-start)+'s')
#%% create control ARX model
T_externel = np.array([5.450257,5.448852,5.447447,5.446042,5.444637,5.443232,5.441826,5.440421,5.439016,
5.440421,5.437610,5.436205,5.434799,5.433394,5.431988,5.430583,5.429177,5.427771,
5.426365, 5.424959, 5.423553 ])
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)
# lower,heigh bound for MV
TL = m.Param(values = 16)
TH = m.Param(values = 18)
# steady state initialization
m.options.IMODE = 1
m.solve(disp=False)
# set up MPC
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 = 2.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.options.CV_TYPE=2 # the objective is an l2-norm (squared error)
m.T.SP = 17 # set point
TL.values = np.ones(len(T_externel))*16
TH.values = np.ones(len(T_externel))*18
m.T.value = 17 # Temprature starts at 17
for i in range(len(T_externel)):
m.d = T_externel[i]
m.solve(disp = False)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
I get the following error :
File c:\algoritmempc\gekko_mpc.py:84
m.arx(p,m.y,m.u)
ValueError: operands could not be broadcast together with shapes (2,) (0,)
I also have another question. Since one of my inputs is a disturbance, I'm not sure if the way I declared my variables is correct or not (I want to provide the disturbances myself).
Solution
When creating a question, please include example data that reproduces the error. It appears that the code is correct with this randomly generated data:
# Example data for data_1
np.random.seed(42) # For reproducibility
n = 100 # Number of data points
# Generating example data
T_ext = np.random.uniform(low=5, high=25, size=n)
beta = np.random.uniform(low=0, high=1, size=n)
T_int = T_ext * 0.5 + beta * 10 + np.random.normal(loc=0, scale=2, size=n)
# Create DataFrame
data_1 = pd.DataFrame({
'T_ext': T_ext,
'beta': beta,
'T_int': T_int
})
I've also made minor corrections so that T_externel is correctly defined for the steady-state initialization (one value) and for the model predictive control application.
# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
# Example data for data_1
np.random.seed(42) # For reproducibility
n = 100 # Number of data points
# Generating example data
T_ext = np.random.uniform(low=5, high=25, size=n)
beta = np.random.uniform(low=0, high=1, size=n)
T_int = T_ext * 0.5 + beta * 10 + np.random.normal(loc=0, scale=2, size=n)
# Create DataFrame
data_1 = pd.DataFrame({
'T_ext': T_ext,
'beta': beta,
'T_int': T_int
})
# Initialize Model
ts = 300
t = np.arange(0,len(data_1)*ts, ts)
u_id = data_1[['T_ext','beta']]
y_id = data_1[['T_int']]
# system identification
m = GEKKO()
#meas : the time-series next step is predicted from prior measurements as in ARX
na=5; nb=5 # ARX coefficients
print('Identify model')
start = time.time()
yp,p,K = m.sysid(t,u_id,y_id,na,nb,objf=100,scale=False,diaglevel=0,pred='meas')
print('temps de prediction :'+str(time.time()-start)+'s')
#%% create control ARX model
T_externel = np.array([5.450257,5.448852,5.447447,5.446042,5.444637,5.443232,5.441826,5.440421,5.439016,
5.440421,5.437610,5.436205,5.434799,5.433394,5.431988,5.430583,5.429177,5.427771,
5.426365, 5.424959, 5.423553 ])
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 = 2.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.d = T_externel[i]
m.solve(disp = False)
if m.options.APPSTATUS == 1:
# Retrieve new values
beta = m.beta.NEWVAL
else:
# Solution failed
beta = 0.0
I recommend the Temperature Control Lab (TCLab) for learning about model predictive control with ARX models. This is a heat transfer application, similar to the home application where there is an external temperature (T2) that can be used as a disturbance, and internal temperature (T1) that can be used as the internal temperature. This code uses TCLabModel() (digital twin) that can be switched to TCLab() if you have the device. It demonstrates MIMO MPC, but you could switch T2 to the disturbance.
import numpy as np
import time
import matplotlib.pyplot as plt
import pandas as pd
import json
# get gekko package with:
# pip install gekko
from gekko import GEKKO
# get tclab package with:
# pip install tclab
from tclab import TCLab
# Connect to Arduino
a = TCLabModel()
# Make an MP4 animation?
make_mp4 = False
if make_mp4:
import imageio # required to make animation
import os
try:
os.mkdir('./figures')
except:
pass
# Final time
tf = 10 # min
# number of data points (every 2 seconds)
n = tf * 30 + 1
# Percent Heater (0-100%)
Q1s = np.zeros(n)
Q2s = np.zeros(n)
# Temperatures (degC)
T1m = a.T1 * np.ones(n)
T2m = a.T2 * np.ones(n)
# Temperature setpoints
T1sp = T1m[0] * np.ones(n)
T2sp = T2m[0] * np.ones(n)
# Heater set point steps about every 150 sec
T1sp[3:] = 50.0
T2sp[40:] = 35.0
T1sp[80:] = 30.0
T2sp[120:] = 50.0
T1sp[160:] = 45.0
T2sp[200:] = 35.0
T1sp[240:] = 60.0
#########################################################
# Initialize Model
#########################################################
# load data (20 min, dt=2 sec) and parse into columns
url = 'http://apmonitor.com/do/uploads/Main/tclab_2sec.txt'
data = pd.read_csv(url)
t = data['Time']
u = data[['H1','H2']]
y = data[['T1','T2']]
# generate time-series model
m = GEKKO()
##################################################################
# system identification
na = 2 # output coefficients
nb = 2 # input coefficients
print('Identify model')
yp,p,K = m.sysid(t,u,y,na,nb,objf=10000,scale=False,diaglevel=1)
##################################################################
# plot sysid results
plt.figure()
plt.subplot(2,1,1)
plt.plot(t,u)
plt.legend([r'$H_1$',r'$H_2$'])
plt.ylabel('MVs')
plt.subplot(2,1,2)
plt.plot(t,y)
plt.plot(t,yp)
plt.legend([r'$T_{1meas}$',r'$T_{2meas}$',\
r'$T_{1pred}$',r'$T_{2pred}$'])
plt.ylabel('CVs')
plt.xlabel('Time')
plt.savefig('sysid.png')
plt.show()
##################################################################
# create control ARX model
y = m.Array(m.CV,2)
u = m.Array(m.MV,2)
m.arx(p,y,u)
# rename CVs
TC1 = y[0]
TC2 = y[1]
# rename MVs
Q1 = u[0]
Q2 = u[1]
# steady state initialization
m.options.IMODE = 1
m.solve(disp=False)
# set up MPC
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # Objective type
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.time=np.linspace(0,120,61)
# Manipulated variables
Q1.STATUS = 1 # manipulated
Q1.FSTATUS = 0 # not measured
Q1.DMAX = 50.0
Q1.DCOST = 0.1
Q1.UPPER = 100.0
Q1.LOWER = 0.0
Q2.STATUS = 1 # manipulated
Q2.FSTATUS = 0 # not measured
Q2.DMAX = 50.0
Q2.DCOST = 0.1
Q2.UPPER = 100.0
Q2.LOWER = 0.0
# Controlled variables
TC1.STATUS = 1 # drive to set point
TC1.FSTATUS = 1 # receive measurement
TC1.TAU = 20 # response speed (time constant)
TC1.TR_INIT = 2 # reference trajectory
TC1.TR_OPEN = 0
TC2.STATUS = 1 # drive to set point
TC2.FSTATUS = 1 # receive measurement
TC2.TAU = 20 # response speed (time constant)
TC2.TR_INIT = 2 # dead-band
TC2.TR_OPEN = 1
##################################################################
# Create plot
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
# Main Loop
start_time = time.time()
prev_time = start_time
tm = np.zeros(n)
try:
for i in range(1,n-1):
# Sleep time
sleep_max = 2.0
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep-0.01)
else:
time.sleep(0.01)
# Record time and change in time
t = time.time()
dt = t - prev_time
prev_time = t
tm[i] = t - start_time
# Read temperatures in Celsius
T1m[i] = a.T1
T2m[i] = a.T2
# Insert measurements
TC1.MEAS = T1m[i]
TC2.MEAS = T2m[i]
# Adjust setpoints
db1 = 1.0 # dead-band
TC1.SPHI = T1sp[i] + db1
TC1.SPLO = T1sp[i] - db1
db2 = 0.2
TC2.SPHI = T2sp[i] + db2
TC2.SPLO = T2sp[i] - db2
# Adjust heaters with MPC
m.solve()
if m.options.APPSTATUS == 1:
# Retrieve new values
Q1s[i+1] = Q1.NEWVAL
Q2s[i+1] = Q2.NEWVAL
# get additional solution information
with open(m.path+'//results.json') as f:
results = json.load(f)
else:
# Solution failed
Q1s[i+1] = 0.0
Q2s[i+1] = 0.0
# Write new heater values (0-100)
a.Q1(Q1s[i])
a.Q2(Q2s[i])
# Plot
plt.clf()
ax=plt.subplot(3,1,1)
ax.grid()
plt.plot(tm[0:i+1],T1sp[0:i+1]+db1,'k-',\
label=r'$T_1$ target',lw=3)
plt.plot(tm[0:i+1],T1sp[0:i+1]-db1,'k-',\
label=None,lw=3)
plt.plot(tm[0:i+1],T1m[0:i+1],'r.',label=r'$T_1$ measured')
plt.plot(tm[i]+m.time,results['v1.bcv'],'r-',\
label=r'$T_1$ predicted',lw=3)
plt.plot(tm[i]+m.time,results['v1.tr_hi'],'k--',\
label=r'$T_1$ trajectory')
plt.plot(tm[i]+m.time,results['v1.tr_lo'],'k--')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(3,1,2)
ax.grid()
plt.plot(tm[0:i+1],T2sp[0:i+1]+db2,'k-',\
label=r'$T_2$ target',lw=3)
plt.plot(tm[0:i+1],T2sp[0:i+1]-db2,'k-',\
label=None,lw=3)
plt.plot(tm[0:i+1],T2m[0:i+1],'b.',label=r'$T_2$ measured')
plt.plot(tm[i]+m.time,results['v2.bcv'],'b-',\
label=r'$T_2$ predict',lw=3)
plt.plot(tm[i]+m.time,results['v2.tr_hi'],'k--',\
label=r'$T_2$ range')
plt.plot(tm[i]+m.time,results['v2.tr_lo'],'k--')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(3,1,3)
ax.grid()
plt.plot([tm[i],tm[i]],[0,100],'k-',\
label='Current Time',lw=1)
plt.plot(tm[0:i+1],Q1s[0:i+1],'r.-',\
label=r'$Q_1$ history',lw=2)
plt.plot(tm[i]+m.time,Q1.value,'r-',\
label=r'$Q_1$ plan',lw=3)
plt.plot(tm[0:i+1],Q2s[0:i+1],'b.-',\
label=r'$Q_2$ history',lw=2)
plt.plot(tm[i]+m.time,Q2.value,'b-',
label=r'$Q_2$ plan',lw=3)
plt.plot(tm[i]+m.time[1],Q1.value[1],color='red',\
marker='.',markersize=15)
plt.plot(tm[i]+m.time[1],Q2.value[1],color='blue',\
marker='X',markersize=8)
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc=2)
plt.draw()
plt.pause(0.05)
if make_mp4:
filename='./figures/plot_'+str(i+10000)+'.png'
plt.savefig(filename)
# Turn off heaters and close connection
a.Q1(0)
a.Q2(0)
a.close()
# Save figure
plt.savefig('tclab_mpc.png')
# generate mp4 from png figures in batches of 350
if make_mp4:
images = []
iset = 0
for i in range(1,n-1):
filename='./figures/plot_'+str(i+10000)+'.png'
images.append(imageio.imread(filename))
if ((i+1)%350)==0:
imageio.mimsave('results_'+str(iset)+'.mp4', images)
iset += 1
images = []
if images!=[]:
imageio.mimsave('results_'+str(iset)+'.mp4', images)
# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
# Turn off heaters and close connection
a.Q1(0)
a.Q2(0)
a.close()
print('Shutting down')
plt.savefig('tclab_mpc.png')
# Make sure serial connection still closes when there's an error
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
a.close()
print('Error: Shutting down')
plt.savefig('tclab_mpc.png')
raise
- how to average in a specific dimension with numpy.mean?
- Removed Realm, but still getting this error: module importing failed: invalid token (rlm_lldb.py, line 37) File "temp.py", line 1,
- Bool object does not support item assignment
- Conda env vs venv / pyenv / virtualenv / poetry / docker, etc
- How to import values from excel with pandas into tkinter faster?
- Speed up concatenation of excel files with Pandas
- Faster way to read Excel files to pandas dataframe
- Create Azure TimerTrigger Durable Function in python
- Read same file from src and test in PyCharm
- How to use Cython to compile Python 3 into C
- Outputting variable from one thread to other
- How to POST a JSON having a single body parameter in FastAPI?
- touchscreen raspberry pi phantom touch (reset fixes)
- How to determine the encoding of text
- Specify extras_require with pip install -e
- Plotting a second scaled y axis in matplotlib from one set of data
- How can I call scikit-learn classifiers from Java?
- How does Pandas.Series.nbytes work for strings? Results don't seem to match expectations
- How to install ursina on pydroid?
- auditlog with Django and DRF
- Python: How to RECURSIVELY remove None values from a NESTED data structure (lists and dictionaries)?
- How to Mark Repeated Entries as True Starting from the Second Occurrence Using NumPy?
- Python3.5 send object over socket (Pygame cam image)
- Value based partial slicing with non-existing keys is now deprecated
- Camelot: DeprecationError: PdfFileReader is deprecated
- How can I partially fill a numpy array with a value given a range?
- Conda wont install pdfplumber
- tomli-w, adding arrays from a python script
- Override value in pydantic model with environment variable
- Polars-Python: Compare list columns