GEKKO MPC Solver with real-time measurements

Question

Trying to solve MPC with an objective function and real-time measurements, one measurement getting in at a time. I am a bit at a loss on the followings:

1 - Is it necessary to shorten the prediction horizon to n_steps - step + 1 and reinitialize the MVs and CVs at every time interval when new measurement comes in?

2 - Not sure how to collect the next step predicted actuation inputs/ states values after the model is solved.

Should that the predicted actuation inputs be:

self.mpc_u_state[step]  = np.array([n_fans.NEWVAL, 
                                    Cw.NEWVAL, 
                                    n_pumps.NEWVAL, 
                                    Cp.NEWVAL])

or

self.mpc_u_state[step] = np.array([n_fans[step], 
                                   Cw [step], 
                                   n_pumps[step],
                                   Cp[step]])

3 - How about the newly predicted state? Should that be:

mpc_x_state[step]   = np.array([topoil.VALUE[step], 
                                hotspot.VALUE[step],
                                puload.VALUE[step]])

Here is my real-time MPC code. Any help would be much appreciated.

#!/usr/bin/python

from datetime import datetime
import numpy as np
import pandas as pd
import csv as csv
from gekko import GEKKO
import numpy as np
import matplotlib
import matplotlib.pyplot as plt  


ALPHA               = 0.5
DELTA_TOP           = 5       # 5 degC
DELTA_HOT           = 5       # 5 degC
DELTA_PU            = 0.05    # 0.05 p.u

NUM_FANS            = 8         # MAX Number of fans
NUM_PUMPS           = 3         # MAX number of pumps
FAN_POWERS          = [145, 130, 120, 100, 500, 460, 430, 370, 860, 800, 720, 610, 1500, 1350, 1230, 1030]
PUMP_POWERS         = [430.0, 1070.0, 2950.0, 6920.0, 8830.0] #  [0.43, 1.07, 2.95, 6.92, 8.83]


# set up matplotlib
is_ipython = 'inline' in matplotlib.get_backend()
if is_ipython:
    from IPython import display


class MPCooController:
    
    def __init__(self):

        self.ref_state      = pd.DataFrame([
            [0  , '2022-11-11T15:12:17.476577',  67.78,   77.94,   0.6],  
            [1  , '2022-11-11T15:12:17.535194',  64.31,   73.03,   0.6],  
            [2  , '2022-11-11T15:12:17.566615',  61.44,   69.90,   0.6], 
            [3  , '2022-11-11T15:12:17.613887',  58.41,   67.16,   0.6], 
            [4  , '2022-11-11T15:12:17.653718',  55.98,   64.62,   0.6],  
            [5  , '2022-11-11T15:12:17.696774',  53.47,   62.41,   0.6], 
            [6  , '2022-11-11T15:12:17.726733',  51.41,   60.38,   0.6], 
            [7  , '2022-11-11T15:12:17.765546',  49.37,   58.57,   0.6], 
            [8  , '2022-11-11T15:12:17.809288',  47.63,   56.93,   0.6], 
            [9  , '2022-11-11T15:12:17.841497',  46.04,   55.50,   0.6], 
            [10 , '2022-11-11T15:12:17.878795',  44.61,   54.22,   0.6],  
            [11 , '2022-11-11T15:12:17.921976',  43.46,   53.14,   0.6],  
            [12 , '2022-11-11T15:12:17.964345',  42.32,   52.75,   0.7],  
            [13 , '2022-11-11T15:12:17.997516',  42.10,   54.73,   0.7], 
            [14 , '2022-11-11T15:12:18.037895',  41.82,   55.56,   0.8],  
            [15 , '2022-11-11T15:12:18.076159',  42.63,   58.60,   0.8],  
            [16 , '2022-11-11T15:12:18.119739',  43.19,   60.29,   0.9], 
            [17 , '2022-11-11T15:12:18.153816',  44.96,   64.24,   0.9], 
            [18 , '2022-11-11T15:12:18.185398',  46.34,   66.69,   1.0],  
            [19 , '2022-11-11T15:12:18.219051',  49.00,   71.43,   1.0],  
            [20 , '2022-11-11T15:12:18.249319',  51.10,   73.73,   1.0],  
            [21 , '2022-11-11T15:12:18.278797',  53.67,   75.80,   1.0],
            [22 , '2022-11-11T15:12:18.311761',  55.53,   77.71,   1.0],
            [23 , '2022-11-11T15:12:18.339181',  57.86,   79.58,   1.0],
            [24 , '2022-11-11T15:12:18.386485',  59.56,   81.72,   1.05],
            [25 , '2022-11-11T15:12:18.421970',  62.10,   85.07,   1.05],
            [26 , '2022-11-11T15:12:18.451925',  64.14,   87.55,   1.1],
            [27 , '2022-11-11T15:12:18.502646',  66.91,   91.12,   1.1], 
            [28 , '2022-11-11T15:12:18.529126',  69.22,   93.78,   1.15],
            [29 , '2022-11-11T15:12:18.557800',  72.11,   97.48,   1.15], 
            [30 , '2022-11-11T15:12:18.591488',  74.60,   100.25,  1.2],
            [31 , '2022-11-11T15:12:18.620894',  77.50,   103.99,  1.2],
            [32 , '2022-11-11T15:12:18.652168',  80.04,   105.84,  1.15],
            [33 , '2022-11-11T15:12:18.692116',  81.82,   106.17,  1.15],
            [34 , '2022-11-11T15:12:18.739722',  83.28,   106.96,  1.1],
            [35 , '2022-11-11T15:12:18.786310',  83.99,   106.39,  1.1], 
            [36 , '2022-11-11T15:12:18.839116',  84.62,   106.82,  1.1],
            [37 , '2022-11-11T15:12:18.872161',  84.91,   107.12,  1.1], 
            [38 , '2022-11-11T15:12:18.908019',  85.34,   107.36,  1.1],
            [39 , '2022-11-11T15:12:18.938229',  85.30,   107.40,  1.1],
            [40 , '2022-11-11T15:12:18.967031',  85.46,   106.54,  1.0],
            [41 , '2022-11-11T15:12:19.001552',  84.21,   103.19,  1.0], 
            [42 , '2022-11-11T15:12:19.035265',  83.19,   101.22,  0.9], 
            [43 , '2022-11-11T15:12:19.069475',  80.95,   97.04,   0.9], 
            [44 , '2022-11-11T15:12:19.094408',  79.11,   94.33,   0.8], 
            [45 , '2022-11-11T15:12:19.123621',  76.21,   89.62,   0.8],
            [46 , '2022-11-11T15:12:19.158660',  73.81,   86.42,   0.7],
            [47 , '2022-11-11T15:12:19.192915',  70.51,   81.42,   0.7], 
            [48 , '2022-11-11T15:12:19.231802',  67.78,   77.94,   0.6]], columns=['id', 'sampdate', 'optopoil', 'ophotspot', 'opload'])


        self.puload              = np.zeros(len(self.ref_state))
        self.hot_noise           = np.zeros(len(self.ref_state))
        self.top_noise           = np.zeros(len(self.ref_state))
        
        self.ref_puload         = []
        self.ref_hotspot        = []
        self.ref_topoil         = []

        self.mpc_play_time      = []
        self.mpc_ref_state      = []
        self.mpc_x_state        = []
        self.mpc_u_state        = []


    # This function simulates observations

    def get_observation(self, step, u_state):

        # Slee 5 seconds to pretend to actuate something with (u_state) and get the resulting state back
        # here the resulting state is simulated with the reference curve affected by a random noise
       # time.sleep(5)

        optopoil     = float(self.ref_state['optopoil'][step])  +  self.top_noise[step]               # Top oil temperature
        ophotspot    = float(self.ref_state['ophotspot'][step]) +  self.hot_noise[step]               # Winding X temperature                                        # Water activity
        opuload      = float(self.ref_state['opload'][step])  + self.puload[step]                     # pu load current X Winding
        return  np.array([optopoil, ophotspot, opuload]) 
    
    def mpc_free_resources(self):
        n_steps                  = len(self.ref_state)
        self.mpc_play_time       = list(np.empty(n_steps))
        self.mpc_x_state         = list(np.empty(n_steps))
        self.mpc_u_state         = list(np.empty(n_steps))
        self.mpc_x_meas          = list(np.empty(n_steps))

        self.pu_noise            = np.random.normal(0, .05, len(self.ref_state))
        self.hot_noise           = np.random.normal(0, 5, len(self.ref_state))
        self.top_noise           = np.random.normal(0, 5, len(self.ref_state))
    
    def mpc_real_mpc(self):   
                
        m                       = GEKKO(remote=False) 
        n_steps                 = len(self.ref_state )
        m.time                  = np.linspace(0, n_steps -1 , n_steps)
        self.mpc_ref_state      = self.ref_state
        mpc_play_time           = list(np.empty(n_steps))
        mpc_x_state             = list(np.empty(n_steps))
        mpc_u_state             = list(np.empty(n_steps))
        mpc_x_meas              = list(np.empty(n_steps))

        alpha                   = m.Const(value = ALPHA)
        delta_top               = m.Const(value = DELTA_TOP)
        delta_hot               = m.Const(value = DELTA_HOT)
        delta_pu                = m.Const(value = DELTA_PU)
        C_base                  = m.Const(value = NUM_FANS * np.max(FAN_POWERS) + NUM_PUMPS * np.max(PUMP_POWERS))   #  kW

        # Reference parameters

        ref_puload              = m.Param(np.array(self.ref_state['opload']))
        ref_hotspot             = m.Param(np.array(self.ref_state['ophotspot']))
        ref_topoil              = m.Param(np.array(self.ref_state['optopoil']))


        # Reference curves lower and higher bounds

        tophigh             = m.Param(value = ref_topoil.VALUE)
        toplow              = m.Param(value = ref_topoil.VALUE  - delta_top.VALUE)

        hothigh             = m.Param(value = ref_hotspot.VALUE)
        hotlow              = m.Param(value = ref_hotspot.VALUE - delta_hot.VALUE)

        puhigh              = m.Param(value = ref_puload.VALUE)
        pulow               = m.Param(value = ref_puload.VALUE  - delta_pu.VALUE)

        # Controlled Variables 

        puload              = m.CV(lb = np.min(pulow.VALUE),  ub = np.max(puhigh.VALUE))
        hotspot             = m.CV(lb = np.min(hotlow.VALUE), ub = np.max(hothigh.VALUE))
        topoil              = m.CV(lb = np.min(toplow.VALUE), ub = np.max(tophigh.VALUE))


        # Manipulated variables 

        n_fans              = m.MV(value = 0,                   lb = 0,                     ub = NUM_FANS, integer=True)
        n_pumps             = m.MV(value = 1,                   lb = 1,                     ub = NUM_PUMPS, integer=True)
        Cw                  = m.MV(value = np.min(FAN_POWERS),  lb = np.min(FAN_POWERS),    ub = np.max(FAN_POWERS))
        Cp                  = m.MV(value = np.min(PUMP_POWERS), lb = np.min(PUMP_POWERS),   ub = np.max(PUMP_POWERS))


        # CVs Status (both measured and calculated)
        
        puload.FSTATUS      = 1
        hotspot.FSTATUS     = 1
        topoil.FSTATUS      = 1

        puload.STATUS       = 1
        hotspot.STATUS      = 1
        topoil.STATUS       = 1

        # Action status
        n_fans.STATUS       = 1
        n_pumps.STATUS      = 1
        Cw.STATUS           = 1
        Cp.STATUS           = 1

        # Not measured
        n_fans.FSTATUS       = 0
        n_pumps.FSTATUS      = 0
        Cw.FSTATUS           = 0
        Cp.FSTATUS           = 0
        
        # The Objective Function (Fuv) cumulating overtime
        power_cost      = m.Intermediate((((n_fans * Cw + n_pumps * Cp) - C_base) / C_base)**2)
 
        tracking_cost   = m.Intermediate (((ref_puload - puload) / ref_puload)**2 
                                        + ((ref_hotspot - hotspot) / ref_hotspot)**2 
                                        + ((ref_topoil  - topoil) / ref_topoil)**2)
                         
        Fuv = m.Intermediate(alpha * power_cost + (1 - alpha) * tracking_cost) 
        
        # Initial solution
        step  = 0
        u_state                  = np.array([0, np.min(FAN_POWERS), 1, np.min(PUMP_POWERS)])  
        x_state                  = self.get_observation(step, u_state)
        topoil.MEAS              = x_state[0] 
        hotspot.MEAS             = x_state[1] 
        puload.MEAS              = x_state[2] 


        m.options.TIME_SHIFT = 1
        m.options.CV_TYPE = 2
        m.Obj(Fuv)
        m.options.IMODE  = 6    
        m.options.SOLVER = 1
        m.solve(disp=True, debug=False)


        mpc_x_state[0]     = np.array([topoil.MODEL, hotspot.MODEL, puload.MODEL])
        mpc_u_state[0]     = np.array([n_fans.NEWVAL, Cw.NEWVAL, n_pumps.NEWVAL, Cp.NEWVAL])
        mpc_x_meas[0]      = np.array([topoil.MEAS, hotspot.MEAS, puload.MEAS]) 
        u_state            = mpc_u_state[0] 
        mpc_play_time[0]   = 0

        # Actuation Input at time step = 0

        while(True):

            for step in range(1, n_steps):
                x_state                  = self.get_observation(step, u_state)
                topoil.MEAS              = x_state[0] 
                hotspot.MEAS             = x_state[1] 
                puload.MEAS              = x_state[2] 

                topoil.SP                = tophigh[step]
                hotspot.SP               = hothigh[step]
                puload.SP                = puhigh[step]
     
                m.solve(disp=True, debug=False)
                        
                mpc_x_state[step]   = np.array([topoil.MODEL, hotspot.MODEL, puload.MODEL]) 
                mpc_x_meas[step]    = np.array([topoil.MEAS, hotspot.MEAS, puload.MEAS]) 
                mpc_u_state[step]   = np.array([n_fans.NEWVAL, Cw.NEWVAL, n_pumps.NEWVAL, Cp.NEWVAL]) 
                
                # New actuation inputs 
                u_state             = mpc_u_state[step] 
                mpc_play_time[step] = step
            
            self.mpc_x_state    = mpc_x_state
            self.mpc_x_meas     = mpc_x_meas
            self.mpc_u_state    = mpc_u_state
            self.mpc_play_time  = mpc_play_time
            self.plot_ctl_mpc()
            self.mpc_free_resources()


    def plot_ctl_mpc(self):
        print("\n\n\n\n===== mpc_u_state ========\n", self.mpc_u_state)
        print("\n\n===== mpc_x_state ========\n", self.mpc_x_state)
                  
        self.mpc_x_state    = pd.DataFrame(self.mpc_x_state, columns=['optopoil','ophotspot','opload'])
        self.mpc_x_meas     = pd.DataFrame(self.mpc_x_meas, columns=['optopoil','ophotspot','opload'])
        self.mpc_u_state    = pd.DataFrame(self.mpc_u_state, columns=['nfans', 'fpower', 'npumps', 'ppower'])
        
        print("\n\n===== mpc_u_state ========\n", self.mpc_u_state)
        print("\n\n===== mpc_x_state ========\n", self.mpc_x_state)
        print("\n\n===== mpc_x_meas ========\n",  self.mpc_x_meas)

        # Results Collection over play time           
        fig1, ax        = plt.subplots()
        ref_lns_hot,    = ax.plot(self.mpc_play_time, self.mpc_ref_state['ophotspot'], 'r', label="ref-hot spot") 
        mpc_lns_hot,    = ax.plot(self.mpc_play_time, self.mpc_x_state['ophotspot'], 'r--', label="mpc-hot spot") 
        # mpc_hot_meas,   = ax.plot(self.mpc_play_time, self.mpc_x_meas['ophotspot'], 'r+-', label="mpc_hot_meas") 

        ref_lns_top,    = ax.plot(self.mpc_play_time, self.mpc_ref_state['optopoil'], 'y', label="ref-top oil")
        mpc_lns_top,    = ax.plot(self.mpc_play_time, self.mpc_x_state['optopoil'], 'y--', label="mpc-top oil")
        # mpc_top_meas,   = ax.plot(self.mpc_play_time, self.mpc_x_meas['optopoil'], 'y+-', label="mpc_top_meas") 

        ax2 = ax.twinx()
        ref_lns_load, = ax2.plot(self.mpc_play_time, self.mpc_ref_state['opload'], 'k', drawstyle='steps-post', label='ref-pu-load') 
        mpc_lns_load, = ax2.plot(self.mpc_play_time, self.mpc_x_state['opload'], 'k--', drawstyle='steps-post', label="mpc-pu-load") 
        # mpc_load_meas, = ax2.plot(self.mpc_play_time, self.mpc_x_meas['opload'], 'k+-', drawstyle='steps-post', label="meas-pu-load") 
        
        ax2.set_ylabel('Load[p.u]')
        ax.set_xlabel('Time [min]')
        ax.set_ylabel('Temperatures[degC]')
        ax.set_title('Thermal and loads stimuli distribution')

        # ax2.legend(handles=[ref_lns_hot, mpc_lns_hot, rl_lns_hot, ref_lns_top, mpc_lns_top, rl_lns_top, ref_lns_load, mpc_lns_load, rl_lns_load], loc='best')
        fig2, ax3        = plt.subplots()
        ax3.plot(self.mpc_play_time, self.mpc_u_state['fpower'] * self.mpc_u_state['nfans'], drawstyle='steps-post', label="Fans Power") 
        ax3.plot(self.mpc_play_time, self.mpc_u_state['ppower'] * self.mpc_u_state['npumps'], drawstyle='steps-post', label="Pumps Power")
        plt.show() 

          
if __name__ == '__main__':

    mpco_controller                         = MPCooController()
    mpco_controller.mpc_real_mpc()

John Hedengren · Answer 1 · 2022-11-20T14:04:13.980

Every time the m.solve() command is issued, Gekko manages the time shifting, re-initialization, and solution.

It is not necessary to shorten the time horizon with every cycle. The time horizon remains constant unless it is a batch process that shortens the horizon as the batch proceeds. Here is a graphic that shows how the time horizon remains constant. The two CVs (top plots) have a prediction horizon with a setpoint indicated by the dashed target region.

The predicted value is:

self.mpc_u_state[step]  = np.array([n_fans.NEWVAL, 
                                    Cw.NEWVAL, 
                                    n_pumps.NEWVAL, 
                                    Cp.NEWVAL])

this is equivalent to:

self.mpc_u_state[step]  = np.array([n_fans.value[1], 
                                    Cw.value[1], 
                                    n_pumps.value[1], 
                                    Cp.value[1]])

The newly predicted state is:

mpc_x_state[step]   = np.array([topoil.MODEL, 
                                hotspot.MODEL,
                                puload.MODEL])

or you can take any value from the time horizon such as the initial condition:

mpc_x_state[step]   = np.array([topoil.value[0], 
                                hotspot.value[0],
                                puload.value[0]])

The Temperature Control Lab is a good example of real-time MPC that runs with an Arduino Leonardo for DAQ and has a serial interface to Python or Matlab for the plots. The TCLab examples can be run with TCLab() or with TCLabModel() if the TCLab hardware is not available.

Response to Edit

Each m.Var(), m.SV(), and m.CV() needs a corresponding equation with m.Equation() to determine the value. The declaration of an m.Var() creates an additional degree of freedom and m.Equation() reduces the degree of freedom by one. The model has three m.CV() definitions but no corresponding equations for puload, hotspot, and topoil. Equations need to be defined that relate the MVs or other adjustable inputs to these outputs. The optimizer then selects the best MVs or FVs to minimize the objective function that combines power and tracking costs.

A convenient way to check that the degrees of freedom are specified correctly is to set m.options.COLDSTART=1 for the first solve.

m.options.COLDSTART = 1
m.solve(disp=True, debug=True)

m.options.COLDSTART = 0
m.solve(disp=True, debug=False)

If the degrees of freedom are not set properly, there is an error:

 Number of state variables:           1104
 Number of total equations: -          960
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :            144
 
 @error: Degrees of Freedom
 * Error: DOF must be zero for this mode
 STOPPING...

Once the degrees of freedom are correct, another suggestion is to avoid hard constraints on the CVs. This can lead to an infeasibility.

puload = m.CV() #lb = np.min(pulow.VALUE),  ub = np.max(puhigh.VALUE))
hotspot = m.CV() #lb = np.min(hotlow.VALUE), ub = np.max(hothigh.VALUE))
topoil = m.CV() #lb = np.min(toplow.VALUE), ub = np.max(tophigh.VALUE))

It is better to use CV_TYPE=1 and set SPHI and SPLO values so that violations of these constraints can occur to maintain feasibility.

The curious part is that the actuation manipulated variables (MVs) values never change and remains static at every time step, despite large fluctuations in the Control variables (CVs) value. What could possibly be causing such behaviors? — info2metlab, Nov 10 '22 at 14:55
It is likely one of three issues: (1) The objective function (check `m.options.OBJFCNVAL`) is not defined and needs to be activated with `cv.STATUS=1` for each controlled variable, (2) The manipulated variable is not turned on - check that `mv.STATUS=1` for each manipulated variable, (3) There is no mathematical relationship between the MV and CV. It looks like #1 and #2 are there. It's hard to troubleshoot without a script that can be run. Could you replace `get_reference()` with something else? — John Hedengren, Nov 11 '22 at 00:05
The code is now updated with the static reference curve. and dummy observations to simulate the states value. — info2metlab, Nov 12 '22 at 03:08

GEKKO MPC Solver with real-time measurements

1 Answers1