Question: The following is related to this question: Why doesn't the Gekko solver adapt to variations in the system? What I still don't understand is why, sometimes, when the outside temperature rises, the inside temperature remains constant. Normally, it should also increase since beta remains constant. Here are examples of this:   In the two images, we can see that on the second day, the outside temperature increases and β remains constant, but the inside temperature decreases. I checked my data and found that it had been taken every 2 seconds, which represents 43,200 measurements per day. However, during the system identification with the ARX model, I had set the sampling time to 300 sec. I ran the corrected MPC code with 2 sec sampling from my previous question. # 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 # import data tx = pd.read_csv(r'C:\Users\mpc\tx_b.csv') tz = pd.read_csv(r'C:\Users\mpc\tz.csv') data = pd.concat([tx,tz],axis=1) data_1 = data[0:8596800] data_2 = data[8596800:] #%% Initialize Model m = GEKKO(remote=False) # system identification ts = 2 t = np.arange(0,len(data_1)*ts, ts) u_id = data_1[['Tx_i','beta_i']] y_id = data_1[['Tz_i']] #meas : the time-series next step is predicted from prior measurements as in ARX na=10; nb=10 # 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') #%% parametres K = array([[ 0.93688819, -12.22410568]]) p = {'a': array([[ 1.08945931], [-0.00243571], [-0.00247112], [-0.00273341], [-0.00296342], [-0.00319516], [-0.00343794], [-0.00366398], [-0.00394255], [-0.06661506]]), 'b': array([[[-0.05134201, -0.01035174], [ 0.00170311, -0.01551259], [ 0.00172715, -0.01178932], [ 0.00178147, -0.01051817], [ 0.00184694, -0.00821511], [ 0.00192371, -0.00570574], [ 0.00201409, -0.00344425], [ 0.00210016, -0.0014708 ], [ 0.00222189, 0.00021622], [ 0.03789636, 0.04235503]]]), 'c': array([0.0266222])} #%% I used the last day's external temperature data as a disturbance. T_externel = data_2[["Tx_i"]].values 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)*2,2) # 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 for i in range(43200): print(i) else: # Solution failed beta = 0.0 # Plot the results plt.figure(figsize=(8,3.5)) plt.subplot(3,1,1) plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{int}$') 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.legend(loc=1); plt.grid() plt.ylabel('Tin (°C)') plt.subplot(3,1,2) plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$') plt.ylabel('Tex (°C)') plt.subplot(3,1,3) 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('results7.png',dpi=300) plt.show()  For the external temperature over one day, I took one day's worth of data (data_2) and concatenated it with itself to create a profile for two consecutive days.  I haven't solved the problem yet.

Answer: There is a small error where m.d is redefined and the connection is broken with the m.arx() model. Here is a correction: # distrubance and parametres #m.d = m.Param(T_externel[0]) m.d.value = T_externel[0]  There is also a dynamics issue. The response to beta (and likely external temperature) is very slow as shown in the step response here:  This could be from using a different sampling time interval for the SysID (System Identification) and the controller. This is especially the case if you used a short time interval for system identification and use a longer time interval for the MPC application. Because it is a discrete ARX model, the dynamics will take longer for the inputs to affect the outputs. Response with Additional Information Thanks for confirming that the identification was performed with 2 second data and the MPC application with 300 second time intervals. This is the source of one of the issues. One method is to use 2 sec for both (as you have shown) or else use 300 sec for both. I recommend going to 300 sec intervals for both because of the slow dynamics and need to calculate over a long time horizon. There is a way to reduce the system identification data to 300 second intervals with a rolling average and sampling every 150 data points with davg.iloc[::150]. Here are the first 500 points of the identification.  # 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 # import data tx = pd.read_csv('tx_b.csv') tz = pd.read_csv('tz.csv') data = pd.concat([tx,tz],axis=1) # reduce from 2 sec to 5 min intervals # calculate rolling average davg = data.rolling(window=150).mean() # sample every 150 cycles (2 sec to 5 min) dnew = davg.iloc[::150] dnew.dropna(inplace=True) slice = int(len(dnew)*0.99) data_1 = dnew.iloc[:slice].copy() # 99% for training data_2 = dnew.iloc[slice:].copy() # 1% for testing #%% Initialize Model m = GEKKO(remote=False) # system identification ts = 300 t = np.arange(0,len(data_1)*ts, ts) u_id = data_1[['Tx_i','beta_i']] y_id = data_1[['Tz_i']] #meas : the time-series next step is predicted from prior measurements as in ARX na=10; nb=10 # 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') # new parameters print(f'Parameters: {p}') # new gains print(f'Gains: {K}') # Plot the results nview = 500 plt.figure(figsize=(8,3.5)) plt.subplot(3,1,1) plt.plot(t[0:nview],data_1[['Tz_i']].iloc[0:nview],'r-',label=r'$T_{meas}$') plt.plot(t[0:nview],yp[0:nview],'b--',label=r'$T_{pred}$') plt.legend(loc=1); plt.grid() plt.ylabel('Tin (°C)') plt.subplot(3,1,2) plt.plot(t[0:nview],data_1['Tx_i'].iloc[0:nview],'g:',label=r'$T_{ext}$') plt.ylabel('Tex (°C)') plt.legend(loc=1); plt.grid() plt.subplot(3,1,3) plt.step(t[0:nview],data_1['beta_i'].iloc[0:nview],'b--',label=r'$\beta$') plt.ylabel('optimal control') plt.xlabel('Time (sec)') plt.legend(loc=1); plt.grid() plt.savefig('results_id.png',dpi=300) plt.show()  The new gains look good: Gains: [[ 0.97529892 -14.28674666]]  The new parameters are: Parameters: {'a': array([[ 0.65059382], [-0.03027903], [ 0.06080624], [-0.02233488], [ 0.06877417], [ 0.37139369], [-0.22535397], [-0.01737772], [-0.03430238], [ 0.02420281]]), 'b': array([[[ 3.71568256e-01, -8.09407185e+00], [-7.62375332e-02, 4.97505860e+00], [-1.34657379e-01, -1.56110911e+00], [ 7.32255512e-02, 1.89452386e+00], [-2.50583697e-02, -1.68565102e+00], [-2.06554465e-01, 2.03970469e+00], [ 1.50625153e-01, 8.79190526e-01], [-9.20974935e-03, -1.34505749e-01], [ 1.58409954e-03, -6.32745151e-01], [ 4.79076587e-03, 1.21199721e-01]]]), 'c': array([2.08698176])}  Here is a step response with T_externel.  # 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 #%% parameters p = {'a': array([[ 0.65059382], [-0.03027903], [ 0.06080624], [-0.02233488], [ 0.06877417], [ 0.37139369], [-0.22535397], [-0.01737772], [-0.03430238], [ 0.02420281]]), 'b': array([[[ 3.71568256e-01, -8.09407185e+00], [-7.62375332e-02, 4.97505860e+00], [-1.34657379e-01, -1.56110911e+00], [ 7.32255512e-02, 1.89452386e+00], [-2.50583697e-02, -1.68565102e+00], [-2.06554465e-01, 2.03970469e+00], [ 1.50625153e-01, 8.79190526e-01], [-9.20974935e-03, -1.34505749e-01], [ 1.58409954e-03, -6.32745151e-01], [ 4.79076587e-03, 1.21199721e-01]]]), 'c': array([2.08698176])} 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.value = 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 n = 120 # 10 hours at 5 min sampling rate T_externel = np.zeros(n) T_externel[10:] = 15 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.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=True) 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(3,1,1) plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{int}$') plt.legend(loc=1); plt.grid() plt.ylabel('Tin (°C)') plt.subplot(3,1,2) plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$') plt.ylabel('Tex (°C)') plt.legend(loc=1); plt.grid() plt.subplot(3,1,3) 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.tight_layout() plt.savefig('results_step.png',dpi=300) plt.show()  Here is the final controller with a sample sinusoidal profile for the external temperature.  # 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 time_points = np.linspace(0, total_minutes_in_day, num_points) # Générer les points de temps T_externel = 8.5 + 11.5 * np.sin((2 * np.pi / total_minutes_in_day) * time_points - (np.pi / 2)) #%% parameters p = {'a': array([[ 0.65059382], [-0.03027903], [ 0.06080624], [-0.02233488], [ 0.06877417], [ 0.37139369], [-0.22535397], [-0.01737772], [-0.03430238], [ 0.02420281]]), 'b': array([[[ 3.71568256e-01, -8.09407185e+00], [-7.62375332e-02, 4.97505860e+00], [-1.34657379e-01, -1.56110911e+00], [ 7.32255512e-02, 1.89452386e+00], [-2.50583697e-02, -1.68565102e+00], [-2.06554465e-01, 2.03970469e+00], [ 1.50625153e-01, 8.79190526e-01], [-9.20974935e-03, -1.34505749e-01], [ 1.58409954e-03, -6.32745151e-01], [ 4.79076587e-03, 1.21199721e-01]]]), 'c': array([2.08698176])} 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.value = 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) 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 = time_points # step time = 300s # Manipulated variables m.beta.STATUS = 1 # calculated by the optimizer m.beta.FSTATUS = 0 # use measured value m.beta.DCOST = 1e-4 # 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 = 19 # set point high level m.Tb.SPLO = 17 # 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=True) 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(3,1,1) plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{int}$') 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.legend(loc=1); plt.grid() plt.ylabel('Tin (°C)') plt.subplot(3,1,2) plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$') plt.ylabel('Tex (°C)') plt.legend(loc=1); plt.grid() plt.subplot(3,1,3) 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.tight_layout() plt.savefig('results_mpc.png',dpi=300) plt.show()