Question: 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).

Answer: 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