在 GEKKO 中使用 MPC 进行自由度警告
Degrees of freedom warning with MPC in GEKKO
我正在尝试使用 GEKKO 和 MPC 预测控制器在给定温度范围内的行为。我是控制系统的新手,以前从未使用过 MPC 或 GEKKO。
我想尽量减少通风系统中风扇的功率,同时将温度保持在所需范围内,特别是在图中显示的范围内:
我使用 PySINDy 获得了一个方程式,该方程式描述了通风系统的风扇功率与房间温度之间的关系。给出为:
dx/dt = 0.205+0.029y-0.003xy-0.001x^2
其中 x 是室温,y 是风扇功率。
我的代码如下:
# Initialize GEKKO model
m = GEKKO()
m.time = np.linspace(0,144,144)
# Manipulated variable
y = m.MV(value=0,lb=0,ub=100) # *** This is the fan power y ***
y.STATUS = 1
y.DCOST = 0.1
y.DMAX = 20
# Controlled variable
x = m.CV(value=18) #Initial value of control variable. *** This is temperature x ***
x.STATUS = 1
m.option.CV_TYPE = 1 #L1 norm
x.RT_INIT = 1
x.TAU = 20
# Add contraints
temp_u = np.concatenate((19*np.ones(32),22.5*np.ones(80),19*np.ones(32))) # Upper temperature limit
temp_l = np.concatenate((17.5*np.ones(38),21.0*np.ones(68),17.5*np.ones(38))) # Lower temperature limit
T_low = m.Param(value=temp_l)
T_high = m.Param(value=temp_u)
m.Equations([x>=T_low, x<=T_high])
# Equation
m.Equation(x.dt() == 0.205+0.029*y-0.003*x*y-0.001*(x**2))
#MPC Mode
m.options.IMODE = 6
#Solve
m.solve()
我得到以下输出:
apm 77.241.105.43_gk_model4 <br><pre> ----------------------------------------------------------------
APMonitor, Version 1.0.1
APMonitor Optimization Suite
----------------------------------------------------------------
--------- APM Model Size ------------
Each time step contains
Objects : 0
Constants : 0
Variables : 6
Intermediates: 0
Connections : 0
Equations : 3
Residuals : 3
Number of state variables: 2860
Number of total equations: - 2717
Number of slack variables: - 286
---------------------------------------
Degrees of freedom : -143
* Warning: DOF <= 0
**********************************************
Dynamic Control with Interior Point Solver
...
Creating file: infeasibilities.txt
Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
@error: Solution Not Found
---------------------------------------------------------------------------
Exception Traceback (most recent call last)
~\mpc.ipynb Cell 11' in <cell line: 1>()
----> 1 m.solve()
File ~\AppData\Local\Programs\Python\Python39\lib\site-packages\gekko\gekko.py:2185, in GEKKO.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)
2187 #load results
2188 def byte2str(byte):
Exception: @error: Solution Not Found
对于那些处于类似位置的人,用这段代码解决了我自己的问题:
m = GEKKO()
m.time = np.linspace(0,144,145)
# Upper temperature limit
temp_u = np.concatenate((18*np.ones(32),22*np.ones(81),18*np.ones(32)))
# Lower temperature limit
temp_l = np.concatenate((17.0*np.ones(42),21.0*np.ones(61),17.0*np.ones(42)))
T_low = m.Param(value=temp_l)
T_high = m.Param(value=temp_u)
#Manipulated variable
u = m.MV(value=0,lb=0,ub=100)
u.STATUS = 1
#Controlled variable
T = m.SV(value=18)
#Soft constraints
eH = m.CV(value=0)
eL = m.CV(value=0)
eH.SPHI=0; eH.WSPHI=10; eH.WSPLO=0 ; eH.STATUS = 1
eL.SPLO=0; eL.WSPHI=0 ; eL.WSPLO=10; eL.STATUS = 1
m.Equations([eH==T-T_high,eL==T-T_low])
#m.Equations([eH==T-T_high])
m.Equation(T.dt() == 02.05+0.029*u-0.003*T*u-0.001*(T**2))
m.Minimize(u)
m.options.IMODE = 6
m.solve(disp=True,debug=True)
我正在尝试使用 GEKKO 和 MPC 预测控制器在给定温度范围内的行为。我是控制系统的新手,以前从未使用过 MPC 或 GEKKO。
我想尽量减少通风系统中风扇的功率,同时将温度保持在所需范围内,特别是在图中显示的范围内:
我使用 PySINDy 获得了一个方程式,该方程式描述了通风系统的风扇功率与房间温度之间的关系。给出为:
dx/dt = 0.205+0.029y-0.003xy-0.001x^2
其中 x 是室温,y 是风扇功率。
我的代码如下:
# Initialize GEKKO model
m = GEKKO()
m.time = np.linspace(0,144,144)
# Manipulated variable
y = m.MV(value=0,lb=0,ub=100) # *** This is the fan power y ***
y.STATUS = 1
y.DCOST = 0.1
y.DMAX = 20
# Controlled variable
x = m.CV(value=18) #Initial value of control variable. *** This is temperature x ***
x.STATUS = 1
m.option.CV_TYPE = 1 #L1 norm
x.RT_INIT = 1
x.TAU = 20
# Add contraints
temp_u = np.concatenate((19*np.ones(32),22.5*np.ones(80),19*np.ones(32))) # Upper temperature limit
temp_l = np.concatenate((17.5*np.ones(38),21.0*np.ones(68),17.5*np.ones(38))) # Lower temperature limit
T_low = m.Param(value=temp_l)
T_high = m.Param(value=temp_u)
m.Equations([x>=T_low, x<=T_high])
# Equation
m.Equation(x.dt() == 0.205+0.029*y-0.003*x*y-0.001*(x**2))
#MPC Mode
m.options.IMODE = 6
#Solve
m.solve()
我得到以下输出:
apm 77.241.105.43_gk_model4 <br><pre> ----------------------------------------------------------------
APMonitor, Version 1.0.1
APMonitor Optimization Suite
----------------------------------------------------------------
--------- APM Model Size ------------
Each time step contains
Objects : 0
Constants : 0
Variables : 6
Intermediates: 0
Connections : 0
Equations : 3
Residuals : 3
Number of state variables: 2860
Number of total equations: - 2717
Number of slack variables: - 286
---------------------------------------
Degrees of freedom : -143
* Warning: DOF <= 0
**********************************************
Dynamic Control with Interior Point Solver
...
Creating file: infeasibilities.txt
Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
@error: Solution Not Found
---------------------------------------------------------------------------
Exception Traceback (most recent call last)
~\mpc.ipynb Cell 11' in <cell line: 1>()
----> 1 m.solve()
File ~\AppData\Local\Programs\Python\Python39\lib\site-packages\gekko\gekko.py:2185, in GEKKO.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)
2187 #load results
2188 def byte2str(byte):
Exception: @error: Solution Not Found
对于那些处于类似位置的人,用这段代码解决了我自己的问题:
m = GEKKO()
m.time = np.linspace(0,144,145)
# Upper temperature limit
temp_u = np.concatenate((18*np.ones(32),22*np.ones(81),18*np.ones(32)))
# Lower temperature limit
temp_l = np.concatenate((17.0*np.ones(42),21.0*np.ones(61),17.0*np.ones(42)))
T_low = m.Param(value=temp_l)
T_high = m.Param(value=temp_u)
#Manipulated variable
u = m.MV(value=0,lb=0,ub=100)
u.STATUS = 1
#Controlled variable
T = m.SV(value=18)
#Soft constraints
eH = m.CV(value=0)
eL = m.CV(value=0)
eH.SPHI=0; eH.WSPHI=10; eH.WSPLO=0 ; eH.STATUS = 1
eL.SPLO=0; eL.WSPHI=0 ; eL.WSPLO=10; eL.STATUS = 1
m.Equations([eH==T-T_high,eL==T-T_low])
#m.Equations([eH==T-T_high])
m.Equation(T.dt() == 02.05+0.029*u-0.003*T*u-0.001*(T**2))
m.Minimize(u)
m.options.IMODE = 6
m.solve(disp=True,debug=True)