GEKKO MPC Solver avec mesures en temps réel

Question

J'essaie de résoudre un MPC avec une fonction objective et des mesures en temps réel, une mesure à la fois. Je suis un peu perdu sur les points suivants :

1 - Est-il nécessaire de raccourcir l'horizon de prédiction pour n_steps - step + 1 et réinitialiser les MVs et CVs à chaque intervalle de temps lorsqu'une nouvelle mesure arrive ?

2 - Je ne suis pas sûr de savoir comment collecter les valeurs prédites des entrées/états d'actionnement de l'étape suivante après la résolution du modèle.

Il faudrait que les entrées d'actionnement prédites soient :

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

o

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

3 - Qu'en est-il de l'état nouvellement prédit ? Cela devrait être :

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

Voici mon code MPC en temps réel. Toute aide serait très appréciée.

import numpy as np
from gekko import GEKKO

m = GEKKO()

# get the reference states curve to follow
ref_state = get_reference() 
n_steps = len(ref_state )
m.time = np.linspace(0, n_steps -1, n_steps) 

mpc_ref_state = np.array([0, np.min(FAN_POWERS), 
                          1, np.min(PUMP_POWERS)])
mpc_play_time = list(np.empty(n_steps))
mpc_x_state = list(np.empty(n_steps))
mpc_u_state = 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.Param(value = NUM_FANS * np.max(FAN_POWERS) 
                       + NUM_PUMPS * np.max(PUMP_POWERS)) 

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

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

# Variables status

puload.FSTATUS = 1
hotspot.FSTATUS = 1
topoil.FSTATUS = 1

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

# Actuation (plant inputs) status
n_fans.STATUS = 1
n_pumps.STATUS = 1
Cw.STATUS = 1
Cp.STATUS = 1

# Objective Function (to be continuously minimize over
#   the time period of measurements being taken)
Fuv = m.Intermediate((alpha * (n_fans * Cw + n_pumps * Cp) / C_base)
 + (1 - alpha) * ((ref_puload - puload) / delta_pu)**2
 + ((ref_hotspot - hotspot) / delta_hot)**2
 + ((ref_topoil - topoil) / delta_top)**2) 

# Objective function and solver settings
m.Obj(Fuv)
m.Minimize(Fuv) 
m.options.IMODE = 6 
m.options.SOLVER = 1

m.solve(disp=True, debug=False)

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

# Starts acquiring the real-time state measurements x_meas, 
#  and solves to predict the next states and actuation inputs 

while (True):
   # get new measurement
   x_meas = get_observation(mpc_u_state[step])
   topoil.MEAS = x_meas[0] 
   hotspot.MEAS = x_meas[1] 
   puload.MEAS = x_meas[2] 

   topoil.SPHI = tophigh.VALUE[step]
   topoil.SPLO = toplow.VALUE[step]

   hotspot.SPHI = hothigh.VALUE[step]
   hotspot.SPLO = hotlow.VALUE[step]

   puload.SPHI = puhigh.VALUE[step]
   puload.SPLO = pulow.VALUE[step]

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

   if (step < n_steps): 
      step = step + 1
   else: # reinitialize to the next cycle
      step = 0

   # Predicted state values for the next time step
   mpc_x_state[step] = np.array([topoil.VALUE[step], 
                                 hotspot.VALUE[step], 
                                 puload.VALUE[step]])

   # Next actuation inputs
   mpc_u_state[step] = np.array([n_fans.NEWVAL, 
                                 Cw.NEWVAL, 
                                 n_pumps.NEWVAL, 
                                 Cp.NEWVAL]) 
   mpc_play_time[step] = step

Answer 1

Chaque fois que le m.solve() est émise, Gekko gère le décalage temporel, la réinitialisation et la solution.

1. Il est no nécessaire de raccourcir l'horizon temporel à chaque cycle. L'horizon temporel reste constant, sauf s'il s'agit d'un processus par lots qui raccourcit l'horizon au fur et à mesure que le lot avance. Voici un graphique qui montre comment l'horizon temporel reste constant. Les deux CV (graphiques du haut) ont un horizon de prédiction avec un point de consigne indiqué par la région cible en pointillés.

2. La valeur prédite est :

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

ce qui est équivalent à :

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

3. Le nouvel état prédit est :

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

ou vous pouvez prendre n'importe quelle valeur de l'horizon temporel, comme la condition initiale :

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

En Laboratoire de contrôle de la température est un bon exemple de MPC en temps réel qui fonctionne avec un Arduino Leonardo pour le DAQ et possède une interface série vers Python ou Matlab pour les graphiques. Les exemples TCLab peuvent être exécutés à l'aide de TCLab() ou de TCLabModel() si l'interface TCLab est utilisée. Matériel TCLab n'est pas disponible.