Question: The following is related to this question: predictive control model using GEKKO I am trying to apply the MPC to maintain the temperature of a room within a defined range, but GEKKO gives me null commands even if the output diverges. 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 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 = np.linspace(9.51,9.78,200) 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 = 0.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.solve(disp = False) if m.options.APPSTATUS == 1: # Retrieve new values beta = m.beta.NEWVAL else: # Solution failed beta = 0.0 import matplotlib.pyplot as plt # Plot the results plt.figure(figsize=(12,6)) plt.subplot(2,1,1) plt.plot(m.time,m.T.value,'r-',label=r'$T_{int}$') plt.plot(m.time,TL.value,'k--',label='Lower Bound') plt.plot(m.time,TH.value,'k--',label='Upper Bound') plt.ylabel('Temperature (°C)') plt.legend() plt.subplot(2,1,2) plt.plot(m.time,m.beta.value,'b--',label=r'$\beta$') plt.ylabel('optimal control') plt.xlabel('Time (sec)') plt.legend() plt.show()  And is the output and the optimal control obtained by GEKKO:    I need to to maintain the temperature of a room within a defined range.

Answer: The gain is listed as K = array([[ 0.93705481, -12.24012156]]) so an increase in beta by +1 leads to a decrease in the T by -12.24. An increase in T_ext by +1 leads to an increase in T by +0.937. The steady-state value of T is 13.32 so to get a starting value of 17, a common practice is to create an additive (or multiplicative bias value that adjusts the starting value. An additional variable Tb is added, although Gekko does this internally. When investigating a controller response, a good first step is to confirm the model gain between the MV and CV with a step response. Here is a step response on beta.  # Import library import numpy as np import pandas as pd import time from gekko import GEKKO from numpy import array 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 = np.linspace(9.51,9.78,200) m = GEKKO(remote=True) 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=True) # 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 = 0 # calculated by the optimizer m.beta.FSTATUS = 0 # use measured value m.beta.DMAX = 0.2 # Delta MV maximum step per horizon interval 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.MEAS = 0.0 # Measured value # step test m.beta.value = np.zeros_like(m.time) m.beta.value[20:] = 1 # Controlled variables m.Tb.STATUS = 1 # drive to set point m.Tb.FSTATUS = 1 # receive measurement m.Tb.SPHI = 21 # set point high level m.Tb.SPLO = 19 # set point low level T_MEAS = 17 # Temperature starts at 17 m.Tb.value = T_MEAS m.bias.value = T_MEAS - m.T.value[0] m.solve(disp=True) if m.options.APPSTATUS == 1: # Retrieve new values beta = m.beta.NEWVAL else: # Solution failed beta = 0.0 import matplotlib.pyplot as plt # 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,'b:',label=r'T_{external}') 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(); 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(); plt.grid() plt.savefig('results.png',dpi=300) plt.show()  The next step to investigate the controller response is to add the disturbance and set the appropriate set point range. The controller is switched to m.options.CV_TYPE=1 to use the l1-norm objective that gives m.Tb.SPHI and m.Tb.SPLO options to specify the range.  # Import library import numpy as np import pandas as pd import time from gekko import GEKKO from numpy import array 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 = np.linspace(9.51,9.78,200) m = GEKKO(remote=True) 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=True) # 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.DMAX = 1.0 # Delta MV maximum step per horizon interval 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 # Controlled variables m.Tb.STATUS = 1 # drive to set point m.Tb.FSTATUS = 0 # receive measurement m.Tb.SPHI = 15 # set point high level m.Tb.SPLO = 13 # set point low level m.Tb.WSPHI = 100 # set point high priority m.Tb.WSPLO = 100 # set point low priority T_MEAS = 17 # Temperature starts at 17 m.Tb.value = T_MEAS m.bias.value = T_MEAS - m.T.value[0] m.solve(disp=True) if m.options.APPSTATUS == 1: # Retrieve new values beta = m.beta.NEWVAL else: # Solution failed beta = 0.0 import matplotlib.pyplot as plt # 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()  The controller keeps the Tb value within the setpoint range by adjusting beta. It is also possible to make the beta value an integer if the cooling system only has an on/off state instead of continuous values between 0 and 1.  # Import library import numpy as np import pandas as pd import time from gekko import GEKKO from numpy import array 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 = np.linspace(9.51,9.78,200) m = GEKKO(remote=True) 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,integer=True) 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=True) # 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 = 10 m.bint.value = 0 # Controlled variables m.Tb.STATUS = 1 # drive to set point m.Tb.FSTATUS = 0 # receive measurement m.Tb.SPHI = 15 # set point high level m.Tb.SPLO = 13 # set point low level m.Tb.WSPHI = 100 # set point high priority m.Tb.WSPLO = 100 # set point low priority T_MEAS = 17 # 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=True) if m.options.APPSTATUS == 1: # Retrieve new values beta = m.beta.NEWVAL else: # Solution failed beta = 0.0 import matplotlib.pyplot as plt # 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.bint.value,'b--',label=r'$\beta_{int}$') plt.ylabel('optimal control') plt.xlabel('Time (sec)') plt.legend(loc=1); plt.grid() plt.savefig('results.png',dpi=300) plt.show()  Solution time for binary variables is significantly longer so many approximate a 0.37 as 37% ON with some post-processing such as off 3.7 minutes and on 6.3 minutes in 10-minute interval blocks. Iter: 71 I: 0 Tm: 16.99 NLPi: 10 Dpth: 4 Lvs: 87 Obj: 3.04E+02 Gap: 6.24E-01 Iter: 72 I: 0 Tm: 11.37 NLPi: 3 Dpth: 4 Lvs: 88 Obj: 3.04E+02 Gap: 6.24E-01 --Integer Solution: 3.04E+02 Lowest Leaf: 3.04E+02 Gap: 9.61E-09 Iter: 73 I: 0 Tm: 13.62 NLPi: 3 Dpth: 13 Lvs: 88 Obj: 3.04E+02 Gap: 9.61E-09 Successful solution --------------------------------------------------- Solver : APOPT (v1.0) Solution time : 1064.45290000000 sec Objective : 304.259454634496 Successful solution ---------------------------------------------------