【发布时间】：2022-11-12 11:11:47
【问题描述】：
尝试使用目标函数和实时测量来解决 MPC，一次测量一个测量。我对以下几点有点不知所措：
1 - 是否有必要将预测范围缩短到n_steps - step + 1 并在新测量进入时的每个时间间隔重新初始化 MV 和 CV？
2 - 不确定如何在模型求解后收集下一步预测的驱动输入/状态值。
如果预测的驱动输入是：
self.mpc_u_state[step] = np.array([n_fans.NEWVAL,
Cw.NEWVAL,
n_pumps.NEWVAL,
Cp.NEWVAL])
或者
self.mpc_u_state[step] = np.array([n_fans[step],
Cw [step],
n_pumps[step],
Cp[step]])
3 - 新预测的状态如何？那应该是：
mpc_x_state[step] = np.array([topoil.VALUE[step],
hotspot.VALUE[step],
puload.VALUE[step]])
这是我的实时 MPC 代码。任何帮助将非常感激。
#!/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()
【问题讨论】：