GEKKO MPC 求解器 生成@error：不等式定义 无效不等式：z > x < y 最小化<生成器对象<..>

Question

我正在学习使用GEKKO MPC求解器，并将上面的代码作为测试编写。经过多次尝试，我仍然无法使其运行，并不断出现以下异常。

APMonitor, Version 1.0.1
APMonitor Optimization Suite
----------------------------------------------------------------

@error: Inequality Definition
invalid inequalities: z > x < y
minimize<generatorobject<genexpr>at0x7f8cb1a1c950>
STOPPING . . .
---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-11-e92d32444662> in <module>
122 # Solver
123 m.options.IMODE     = 6 # control
--> 124 m.solve(disp=True, debug=True)
125 
126 # Plot the results
/usr/local/lib/python3.7/dist-packages/gekko/gekko.py in solve(self, disp, debug, GUI, **kwargs)
2183             #print APM error message and die
2184             if (debug >= 1) and ('@error' in response):
-> 2185                 raise Exception(response)
2186 
2187             #load results
Exception:  @error: Inequality Definition
invalid inequalities: z > x < y
minimize<generatorobject<genexpr>at0x7f8cb1a1c950>
STOPPING . . .

任何帮助都将不胜感激。这是我写的代码：

import numpy as np
import pandas as pd
m = GEKKO()           
ALPHA               = 0.5
NUM_FANS            = 8         # MAX Number of fans
NUM_PUMPS           = 1         # MAX number of pumps
n_steps             = 25
m.time              = np.linspace(0, n_steps - 1, n_steps)
DELTA_TOP           = 5         # 5 degC
DELTA_HOT           = 5         # 5 degC
DELTA_PU            = 0.05      # 0.05 p.u
fan_powers          = np.array([145, 130, 120, 100, 500, 460, 430, 370, 860, 800, 720, 610, 1500, 1350, 1230, 1030])   #  kW
pump_powers         = np.array([0.43, 1.07, 2.95, 6.92, 8.83])                                                         #  kW
C_base              = NUM_FANS * np.max(fan_powers) + NUM_PUMPS * np.max(pump_powers)                                  #  kW
x_state             = np.array([61.29027692, 70.15582365,  0.86972331])
u_state             = np.array([5, 1.00, 1500.00, 0.43])  # np.array(pd.DataFrame([[5, 1.00, 1500.00, 0.43]], columns=['nfans', 'npumps', 'fpower', 'ppower'] )) #np.array([1.00, 1500.00, 0.43])  
ref_state   = pd.DataFrame([
[0,   '2022-08-30T19:33:07.637217',  50.949829,  56.055570, 0.70],
[1,   '2022-08-30T19:33:07.719390',  46.113708,  48.741882, 0.60],
[2,   '2022-08-30T19:33:07.754899',  43.921465,  49.425708, 0.60],
[3,   '2022-08-30T19:33:07.782037',  44.792515,  49.895490, 0.60],
[4,   '2022-08-30T19:33:07.831646',  45.814439,  51.055404, 0.60],
[5,   '2022-08-30T19:33:07.910940',  46.677830,  51.900248, 0.60],
[6,   '2022-08-30T19:33:07.951684',  47.500278,  52.609172, 0.60],
[7,   '2022-08-30T19:33:08.050460',  48.187270,  53.240813, 0.60],
[8,   '2022-08-30T19:33:08.126050',  48.866124,  53.806335, 0.60],
[9,   '2022-08-30T19:33:08.205533',  49.395292,  54.303250, 0.60],
[10,  '2022-08-30T19:33:08.237825',  49.908234,  54.732465, 0.60],
[11,  '2022-08-30T19:33:08.261200',  50.315668,  55.112417, 0.60],
[12,  '2022-08-30T19:33:08.303079',  50.750658,  55.793464, 0.70],
[13,  '2022-08-30T19:33:08.370096',  51.341523,  57.619243, 0.70],
[14,  '2022-08-30T19:33:08.463300',  51.666736,  58.602764, 0.80],
[15,  '2022-08-30T19:33:08.524749',  52.738678,  60.785766, 0.80],
[16,  '2022-08-30T19:33:08.552913',  53.460458,  62.226178, 0.90],
[17,  '2022-08-30T19:33:08.589561',  55.055422,  64.867184, 0.90],
[18,  '2022-08-30T19:33:08.633709',  56.231096,  66.286857, 0.90],
[19,  '2022-08-30T19:33:08.671211',  57.876352,  67.644409, 0.90],
[20,  '2022-08-30T19:33:08.708004',  59.015503,  68.931404, 0.90],
[21,  '2022-08-30T19:33:08.729763',  60.586943,  70.227752, 0.90],
[22,  '2022-08-30T19:33:08.753146',  61.809524,  71.492779, 0.90],
[23,  '2022-08-30T19:33:08.779459',  63.232974,  72.710304, 0.90],
[24,  '2022-08-30T19:33:08.808419',  64.324357, 74.556550, 1.05]], columns=['id', 'sampdate', 'optopoil', 'ophotspot', 'opload'])
# Initial State at i = 0
puload_0    = x_state[2]
hotspot_0   = x_state[1]
topoil_0    = x_state[0]
# Initial Controls at i =0
n_fans_0  = u_state[0]
n_pumps_0 = u_state[1]
Cw_0 = u_state[0] * u_state[2]
Cp_0 = u_state[1] * u_state[3]
# References
ref_puload  = np.array(ref_state['opload'])
ref_hotspot = np.array(ref_state['ophotspot'])
ref_topoil  = np.array(ref_state['optopoil'])
# Controlled variables
tophigh = m.Param(value = ref_topoil)
toplow = m.Param(value =  ref_topoil  - DELTA_TOP)
hothigh = m.Param(value=ref_hotspot)
hotlow = m.Param(value=ref_topoil  - DELTA_TOP)
pulow = m.Param(value=ref_puload)
puhigh= m.Param(value=ref_puload  - DELTA_PU)
puload     = m.CV (value = np.array([puload_0]*n_steps),  lb = ref_puload  - DELTA_PU,   ub = ref_puload)
hotspot    = m.CV (value = np.array([hotspot_0]*n_steps), lb = ref_hotspot - DELTA_HOT,  ub = ref_hotspot)
topoil     = m.CV (value = np.array([topoil_0]*n_steps),  lb = ref_topoil  - DELTA_TOP,  ub = ref_topoil)
m.Equations([topoil >= toplow, topoil <= tophigh])
m.Equations([hotspot >= hotlow, hotspot <= hothigh])
m.Equations([puload >= pulow, puload <= puhigh])
puload.STATUS   = 1
hotspot.STATUS  = 1
topoil.STATUS   = 1
# Manipulated variables
fan_low   = m.Param(value = np.zeros(n_steps))
fan_high  = m.Param(value =  np.array([NUM_FANS]*n_steps))
pump_low  = m.Param(value = np.ones(n_steps))
pump_high = m.Param(value =  np.array([NUM_PUMPS]*n_steps))
Cw_low    = m.Param(value = np.zeros(n_steps))
Cw_high   = m.Param(value= np.array([NUM_FANS * np.max(fan_powers)]*n_steps))
Cp_low    = m.Param(value = np.ones(n_steps))
Cp_high   = m.Param(value = np.array([NUM_PUMPS * np.max(pump_powers)]*n_steps))
n_fans      = m.MV (value = np.array([n_fans_0]*n_steps),  lb = np.zeros(n_steps),  ub = np.array([NUM_FANS]*n_steps))
n_pumps     = m.MV (value = np.array([n_pumps_0]*n_steps), lb = np.ones(n_steps),   ub = np.array([NUM_PUMPS]*n_steps))
Cw          = m.MV (value = np.array([Cw_0]*n_steps),      lb = np.zeros(n_steps),  ub = np.array([NUM_FANS * np.max(fan_powers)]*n_steps))
Cp          = m.MV (value = np.array([Cp_0]*n_steps),      lb = np.ones(n_steps),   ub = np.array([NUM_PUMPS * np.max(pump_powers)]*n_steps))
m.Equation([n_fans >= fan_low, n_fans <= fan_high])
m.Equation([n_pumps >= pump_low, n_pumps <= pump_high])
m.Equation([Cw >= Cw_low, Cw <= Cw_high])
m.Equation([Cp>= Cp_low, Cp <= Cp_high])

n_fans.STATUS = 1
n_pumps.STATUS = 1
Cw.STATUS = 1
Cp.STATUS = 1
# Objective Function
Fuv = m.Var(value=0)
m.Minimize(Fuv = sum(ALPHA * (Cw[i] + Cp[i]) / C_base
+ (1 - ALPHA) * (
pow ((ref_puload[i] - puload[i]) / DELTA_PU, 2)
+  pow ((ref_hotspot[i]- hotspot[i]) / DELTA_HOT, 2)
+  pow ((ref_topoil[i] - topoil[i]) / DELTA_TOP, 2)
)) for i in range(n_steps))
# Solver
m.options.IMODE     = 6 # control
m.solve(disp=True, debug=True)

Answer 1

Gekko进行时间索引。以下是成功解决问题的代码的修改版本。解决方案是通过首先求解而不存在如下不等式：

#%% Solve without inequalities
m.options.IMODE = 6 # control
m.options.SOLVER=1
m.solve(disp=True, debug=True)

----------------------------------------------
Dynamic Control with APOPT Solver
----------------------------------------------

Iter    Objective  Convergence
0  6.51001E+07  0.00000E+00
1  1.19052E+16  0.00000E+00
2  7.77312E+01  0.00000E+00
3  7.49548E+00  0.00000E+00
4  7.49548E+00  0.00000E+00
6  7.49548E+00  0.00000E+00
Successful solution

---------------------------------------------------
Solver         :  APOPT (v1.0)
Solution time  :  0.043000000000000003 sec
Objective      :  7.4954845726019945
Successful solution
---------------------------------------------------

然后使用该解决方案初始化具有不等式的情况：

#%% Turn on DOF and include inequalities
n_fans.STATUS = 1
n_pumps.STATUS = 1
Cw.STATUS = 1
Cp.STATUS = 1
m.Equations([topoil >= toplow, topoil <= tophigh])
m.Equations([hotspot >= hotlow, hotspot <= hothigh])
m.Equations([puload >= pulow, puload <= puhigh])
m.Equations([n_fans >= fan_low, n_fans <= fan_high])
m.Equations([n_pumps >= pump_low, n_pumps <= pump_high])
m.Equations([Cw >= Cw_low, Cw <= Cw_high])
m.Equations([Cp>= Cp_low, Cp <= Cp_high])
m.options.TIME_SHIFT=0
m.solve(disp=True, debug=True)

----------------------------------------------
Dynamic Control with APOPT Solver
----------------------------------------------

Iter    Objective  Convergence
0  7.91421E-02  7.49900E+03
1  7.70285E-02  1.30000E-09
2  7.70285E-02  2.00000E-10
3  7.70285E-02  2.00000E-10
Successful solution

---------------------------------------------------
Solver         :  APOPT (v1.0)
Solution time  :  0.033999999999999996 sec
Objective      :  0.07702852753968847
Successful solution
---------------------------------------------------

以下是完整的脚本：

import numpy as np
import pandas as pd
from gekko import GEKKO
m = GEKKO(remote=False)           
ALPHA               = 0.5
NUM_FANS            = 8         # MAX Number of fans
NUM_PUMPS           = 1         # MAX number of pumps
n_steps             = 25
m.time              = np.linspace(0, n_steps - 1, n_steps)
DELTA_TOP           = 5         # 5 degC
DELTA_HOT           = 5         # 5 degC
DELTA_PU            = 0.05      # 0.05 p.u
fan_powers          = np.array([145, 130, 120, 100, 500, 460, 430, 370, 860, 800, 720, 610, 1500, 1350, 1230, 1030])   #  kW
pump_powers         = np.array([0.43, 1.07, 2.95, 6.92, 8.83])                                                         #  kW
C_base              = NUM_FANS * np.max(fan_powers) + NUM_PUMPS * np.max(pump_powers)                                  #  kW
x_state             = np.array([61.29027692, 70.15582365,  0.86972331])
u_state             = np.array([5, 1.00, 1500.00, 0.43]) 
ref_state   = pd.DataFrame([
[0,   '2022-08-30T19:33:07.637217',  50.949829,  56.055570, 0.70],
[1,   '2022-08-30T19:33:07.719390',  46.113708,  48.741882, 0.60],
[2,   '2022-08-30T19:33:07.754899',  43.921465,  49.425708, 0.60],
[3,   '2022-08-30T19:33:07.782037',  44.792515,  49.895490, 0.60],
[4,   '2022-08-30T19:33:07.831646',  45.814439,  51.055404, 0.60],
[5,   '2022-08-30T19:33:07.910940',  46.677830,  51.900248, 0.60],
[6,   '2022-08-30T19:33:07.951684',  47.500278,  52.609172, 0.60],
[7,   '2022-08-30T19:33:08.050460',  48.187270,  53.240813, 0.60],
[8,   '2022-08-30T19:33:08.126050',  48.866124,  53.806335, 0.60],
[9,   '2022-08-30T19:33:08.205533',  49.395292,  54.303250, 0.60],
[10,  '2022-08-30T19:33:08.237825',  49.908234,  54.732465, 0.60],
[11,  '2022-08-30T19:33:08.261200',  50.315668,  55.112417, 0.60],
[12,  '2022-08-30T19:33:08.303079',  50.750658,  55.793464, 0.70],
[13,  '2022-08-30T19:33:08.370096',  51.341523,  57.619243, 0.70],
[14,  '2022-08-30T19:33:08.463300',  51.666736,  58.602764, 0.80],
[15,  '2022-08-30T19:33:08.524749',  52.738678,  60.785766, 0.80],
[16,  '2022-08-30T19:33:08.552913',  53.460458,  62.226178, 0.90],
[17,  '2022-08-30T19:33:08.589561',  55.055422,  64.867184, 0.90],
[18,  '2022-08-30T19:33:08.633709',  56.231096,  66.286857, 0.90],
[19,  '2022-08-30T19:33:08.671211',  57.876352,  67.644409, 0.90],
[20,  '2022-08-30T19:33:08.708004',  59.015503,  68.931404, 0.90],
[21,  '2022-08-30T19:33:08.729763',  60.586943,  70.227752, 0.90],
[22,  '2022-08-30T19:33:08.753146',  61.809524,  71.492779, 0.90],
[23,  '2022-08-30T19:33:08.779459',  63.232974,  72.710304, 0.90],
[24,  '2022-08-30T19:33:08.808419',  64.324357, 74.556550, 1.05]], 
columns=['id', 'sampdate', 'optopoil', 'ophotspot', 'opload'])
# Initial State at i = 0
puload_0    = x_state[2]
hotspot_0   = x_state[1]
topoil_0    = x_state[0]
# Initial Controls at i =0
n_fans_0  = u_state[0]
n_pumps_0 = u_state[1]
Cw_0 = u_state[0] * u_state[2]
Cp_0 = u_state[1] * u_state[3]
# References
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']))
# Controlled variables
tophigh = m.Param(value = ref_topoil)
toplow = m.Param(value =  ref_topoil  - DELTA_TOP)
hothigh = m.Param(value=ref_hotspot)
hotlow = m.Param(value=ref_topoil  - DELTA_TOP)
pulow = m.Param(value=ref_puload)
puhigh= m.Param(value=ref_puload  - DELTA_PU)
puload     = m.Var(value = puload_0)
hotspot    = m.Var(value = hotspot_0)
topoil     = m.Var(value = topoil_0)
# Manipulated variables
fan_low   = m.Param(value = 0)
fan_high  = m.Param(value = NUM_FANS)
pump_low  = m.Param(value = 1)
pump_high = m.Param(value = NUM_PUMPS)
Cw_low    = m.Param(value = 0)
Cw_high   = m.Param(value= NUM_FANS * np.max(fan_powers))
Cp_low    = m.Param(value = 1)
Cp_high   = m.Param(value = NUM_PUMPS * np.max(pump_powers))
n_fans      = m.MV(value = n_fans_0,  lb = 0,  ub = NUM_FANS)
n_pumps     = m.MV(value = n_pumps_0, lb = 1,  ub = NUM_PUMPS)
Cw          = m.MV(value = Cw_0,      lb = 0,  ub = NUM_FANS * np.max(fan_powers))
Cp          = m.MV(value = Cp_0,      lb = 1,  ub = NUM_PUMPS * np.max(pump_powers))
# Objective Function
Fuv = m.Intermediate(ALPHA * (Cw + Cp) / C_base
+ (1 - ALPHA) * ((ref_puload - puload) / DELTA_PU)**2
+  ((ref_hotspot- hotspot) / DELTA_HOT)**2
+  ((ref_topoil - topoil) / DELTA_TOP)**2)
m.Minimize(Fuv)

#%% Solve without inequalities
m.options.IMODE = 6 # control
m.options.SOLVER=1
m.solve(disp=True, debug=True)

#%% Turn on DOF and include inequalities
n_fans.STATUS = 1
n_pumps.STATUS = 1
Cw.STATUS = 1
Cp.STATUS = 1
m.Equations([topoil >= toplow, topoil <= tophigh])
m.Equations([hotspot >= hotlow, hotspot <= hothigh])
m.Equations([puload >= pulow, puload <= puhigh])
m.Equations([n_fans >= fan_low, n_fans <= fan_high])
m.Equations([n_pumps >= pump_low, n_pumps <= pump_high])
m.Equations([Cw >= Cw_low, Cw <= Cw_high])
m.Equations([Cp>= Cp_low, Cp <= Cp_high])
m.options.TIME_SHIFT=0
m.solve(disp=True, debug=True)