带有实时测量的GEKKO MPC求解器
带有实时测量的GEKKO MPC求解器
提问于 2022-11-02 20:42:28
试图用目标函数和实时测量来求解MPC，一次只测量一次。以下几点我有点不知所措：
1-是否有必要将预测视界缩短到n_steps - step + 1，并在每次新测量时重新初始化MVs和CVs？
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()回答 1
Stack Overflow用户
发布于 2022-11-03 04:09:54
每次发出m.solve()命令时，Gekko都会管理时间转移、重新初始化和解决方案。
- 没有必要缩短每一个周期的时间范围。除非是批处理过程缩短了批处理的时间范围，否则时间范围保持不变。这是一个图表，显示了时间视界是如何保持不变的。这两个CVs (顶级地块)有一个由虚空目标区域指示的设定点的预测地平线。
- 预测值是：
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.value[1],
Cw.value[1],
n_pumps.value[1],
Cp.value[1]])- 新预测的状态是：
mpc_x_state[step] = np.array([topoil.MODEL,
hotspot.MODEL,
puload.MODEL])或者，您可以从时间范围内获取任何值，例如初始条件：
mpc_x_state[step] = np.array([topoil.value[0],
hotspot.value[0],
puload.value[0]])温控实验室是实时MPC的一个很好的例子，它与Arduino Leonardo为DAQ一起运行，并具有与Python或Matlab的串行接口。TCLab示例可以使用TCLab()运行，如果TCLab硬件不可用，则可以使用TCLabModel()运行。
对编辑的响应
每个m.Var()、m.SV()和m.CV()都需要一个与m.Equation()对应的方程来确定值。m.Var()的声明创建了一个额外的自由度，而m.Equation()将自由度降低了一个。该模型有三个m.CV()定义，但没有对应的puload、hotspot和topoil方程。需要定义将MVs或其他可调输入与这些输出相关联的方程。然后，优化器选择最佳的MVs或FVs，以最小化结合功率和跟踪成本的目标函数。
一种检查自由度是否正确指定的方便方法是将m.options.COLDSTART=1设置为第一个解。
m.options.COLDSTART = 1
m.solve(disp=True, debug=True)
m.options.COLDSTART = 0
m.solve(disp=True, debug=False)如果没有正确地设置自由度，就会出现以下错误：
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...一旦自由度正确后，另一个建议是避免对履历施加严格的限制。这可能导致不可行。
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))最好使用CV_TYPE=1并设置SPHI和SPLO值，以便能够发生违反这些约束的情况，以保持可行性。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/74295404
复制相关文章
相似问题
腾讯云开发者
Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 
粤公网安备44030502008569号
腾讯云计算（北京）有限责任公司 京ICP证150476号 | 京ICP备11018762号