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
- RuntimeError: Trying to backward through the graph a second time on loss tensor
- ERROR: post() got an unexpected keyword argument 'proxies'
- How do I automatically fix an invalid JSON string?
- Python generating a list of dates between two dates
- ModuleNotFoundError: No module named 'json'
- Check if string is None, empty or has spaces only
- invalid command name ".!canvas"
- SqlAlchemy: How to make a LONGBLOB column in mysql?
- Validate if input number is prime
- OpenAI remove key access while using AAD authentication
- Converting a Text File to a String in Python
- Problems with repeated draws of Sobol sequence in scipy
- Render NumPy array in FastAPI
- Any faster and memory-efficient alternative of torch.autograd.functional.jacobian(model.decoder, latent_l)?
- Convert integer to string in Python
- starting container process caused: exec: "fastapi": executable file not found in $PATH: unknown
- How to change version of Python in PyCharm
- convert image from CV_64F to CV_8U
- Pyspark new column when otherwise results in "should be a column" error
- FastAPI as a Windows service
- Regex to match a whole number not ending in some digits
- How can I measure the width of a string rendering via tkFont without creating a window first?
- The tilde operator in Python
- 2D stacking and optimizing
- Compute percentage of positive rows in a group_by polars DataFrame
- Extract domain name from URL in Python
- How to use AssertRaisesMessage() in Django tests
- Langchain ChatOllama always produces invalid format: expected "json" error
- Python client library for Language Server Protocol (LSP)
- How can I convert a python urandom to a string?