pythonoptimizationgekkompc

Why doesn't the Gekko solver adapt to variations in the system?

The following is related to this question : predictive control model using GEKKO

I am still traying applying MPC to maintain the temperature of a room within a defined range(16,18), as Professor @John explained last time.the gain is normally listed as K=array([[0.93705481,−12.24012156]]). Thus, an increase in β by +1 leads to a decrease in T by -12.24, while an increase in Text by +1 leads to an increase in T by +0.937.

In my case, I tried to implement an external temperature profile for a day (sinusoidal), and I modified the control to vary between 0 and 1, not just 0 or 1, since the equipment can operate within this range. I also adjusted m.bint.MV_STEP_HOR=1 so that the control is updated at each iteration. However, the control still does not react! Even when the external temperature increases or decreases, the control remains unchanged. 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

total_minutes_in_day = 24 * 60
interval = 5  # minutes
num_points = total_minutes_in_day // interval

# Générer les points de temps
time_points = np.linspace(0, total_minutes_in_day, num_points)

temperature = 8.5 + 11.5 * np.sin((2 * np.pi / total_minutes_in_day) * time_points - (np.pi / 2))
temperature = np.clip(temperature, -3, 20)


plt.plot(time_points, temperature)
plt.xlabel('Temps (minutes)')
plt.ylabel('Température (°C)')
plt.title('Variation de Tex sur une journée')
plt.savefig('Tex pour une jrn')
plt.grid(True)
plt.show()


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 = temperature
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]
m.free(m.beta)
m.bint = m.MV(0,lb=0,ub=1)
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=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)*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 = 1
m.bint.value = 0

# Controlled variables
m.Tb.STATUS = 1        # drive to set point
m.Tb.FSTATUS = 0       # receive measurement
m.Tb.SPHI = 18         # set point high level
m.Tb.SPLO = 16         # set point low level
m.Tb.WSPHI = 100       # set point high priority
m.Tb.WSPLO = 100       # set point low priority

T_MEAS = 20 # 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=False)

if m.options.APPSTATUS == 1:
   # Retrieve new values
   beta  = m.beta.NEWVAL
else:
   # Solution failed
   beta  = 0.0

And is the output and the optimal control obtained by GEKKO:

Solution

  • The solution that you showed was able to find a solution with a single step to keep within the temperature range (17-19). Tightening the temperature range (16.5-17.5) shows additional movement of the MV over the time period and shows that it is responsive. Unless there is a violation of this temperature range, the controller does not take action if there is even a small DCOST value to move the MV.

    results

    The dynamics influence the speed of the response. A better insulated building will have a slower response to changes in the external temperature. One other suggestion is that the bint Manipulated Variable isn't needed when solving for continuous instead of integer solutions.

    Here is the complete code with bint removed and the tighter temperature target range:

    # 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

# Générer les points de temps
time_points = np.linspace(0, total_minutes_in_day, num_points)

temperature = 8.5 + 11.5 * np.sin((2 * np.pi / total_minutes_in_day) * time_points - (np.pi / 2))
temperature = np.clip(temperature, -3, 20)

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])}

T_externel = temperature
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)*300,300) # 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
else:
   # Solution failed
   beta  = 0.0

# 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()