使用 Gekko 的 ARX 模型的 MPC
MPC with ARX Model Using Gekko
我正在对 MPC 进行建模以控制冰箱并将温度保持在给定的时间间隔内，同时最大限度地降低成本。我正在使用 GEKKO 为我的算法建模。
我写了下面的代码。首先，我使用系统中的传感器数据识别了我的模型（我使用了 GEKKO 的函数 sysif）。然后我构建了一个 ARX 模型（使用 GEKKO 中的 arx 函数），它成为 sysid() 作为输入的结果。
我正在尝试编写一个“虚拟”算法以在将其实现到 Pi 之前在本地进行测试。
我收到以下错误：
KeyError Traceback (most recent call last)
<ipython-input-13-108148376700> in <module>
107 #Solve the optimization problem.
108
--> 109 m.solve()
~/opt/anaconda3/lib/python3.8/site-packages/gekko/gekko.py in solve(self, disp, debug, GUI, **kwargs)
2214 if timing == True:
2215 t = time.time()
-> 2216 self.load_JSON()
2217 if timing == True:
2218 print('load JSON', time.time() - t)
~/opt/anaconda3/lib/python3.8/site-packages/gekko/gk_post_solve.py in load_JSON(self)
48 vp.__dict__[o] = dpred
49 else: #everything besides value, dpred and pred
---> 50 vp.__dict__[o] = data[vp.name][o]
51 for vp in self._variables:
52 if vp.type != None: #(FV/MV/SV/CV) not Param or Var
KeyError: 'int_p6'
这是我的代码
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO(remote = True)
#initialize variables
#Room Temprature:
T_external = [23,23,23,23,23.5,23.5,23.4,23.5,23.9,23.7,\
23,23.9,23.9,23.4,23.9,24,23.6,23.7,23.8,\
23,23,23,23,23]
# Temprature Lower Limit:
temp_low = 10*np.ones(24)
# Temprature Upper Limit:
temp_upper = 12*np.ones(24)
#Hourly Energy prices:
TOU_v = [39.09,34.93,38.39,40.46,40.57,43.93,25,11,9,24,51.28,45.22,45.72,\
36,35.03,10,12,13,32.81,42.55,8,29.58,29.52,29.52]
###########################################
#System Identification:
#Time
t = np.linspace(0,10,117)
#State of the Fridge
ud = np.append(np.zeros(78) ,np.ones(39),0)
#Temprature Data
y = [14.600000000000001,14.600000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.8,14.8,14.8,14.8,14.8,14.8,14.8,14.8,\
14.8,14.8,14.9,14.9,14.9,14.9,14.9,14.9,14.9,15,15,15,15,15,15,15,15,15,15,15,15,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,
15,15,15,15,15,15,15,15,15,15,14.9,14.9,14.9,14.9,14.8,14.9,14.8,14.8,14.8,14.8,14.8,14.8,\
14.8,14.700000000000001,14.8,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.600000000000001,14.600000000000001,14.600000000000001,\
14.600000000000001,14.600000000000001,14.60]
na = 1 # output coefficients
nb = 1 # input coefficients
print('Identification')
yp,p,K = m.sysid(t,ud,y,na,nb,objf=10000,scale=False,diaglevel=1)
#create control ARX model:
y = m.Array(m.CV,1)
uc = m.Array(m.MV,1)
m.arx(p,y,uc)
# rename CVs
T= y[0]
# rename MVs
uc = uc[0]
# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)
###########################################
#Parameter
P = m.Param(value =100) #power
TL = m.Param(value=temp_low)
TH = m.Param(value=temp_upper)
c = m.Param(value=TOU_v)
# Manipilated variable:
u = m.MV(lb=0, ub=1, integer=True)
u.STATUS = 1 # allow optimizer to change the variable to attein the optimum.
# Controlled Variable (Affected with changes in the manipulated variable)
T = m.CV(value=11) # Temprature will start at 11.
# Soft constraints on temprature.
eH = m.CV(value=0)
eL = m.CV(value=0)
eH.SPHI=0 #Set point high for linear error model.
eH.WSPHI=100 #Objective function weight on upper set point for linear error model.
eH.WSPLO=0 # Objective function weight on lower set point for linear error model
eH.STATUS =1 # eH : Error is considered in the objective function.
eL.SPLO=0
eL.WSPHI=0
eL.WSPLO=100
eL.STATUS = 1
#Linear error (Deviation from the limits)
m.Equations([eH==T-TH,eL==T-TL])
#Objective : minimize the costs.
m.Minimize(c*P*u)
#Optimizer Options.
m.options.IMODE = 6 # MPC mode in Gekko.
m.options.NODES = 2 # Collocation nodes.
m.options.SOLVER = 1 # APOT solver for mixed integer linear programming.
m.time = np.linspace(0,23,24)
#Solve the optimization problem.
m.solve()
问题在于：
T = m.CV(value=11) # Temperature will start at 11.
您正在重新定义 T 变量，但它在内部存储了这两个变量。如果你需要 re-initialize 到 11 然后使用 T.value=11。另外，我在稳态初始化之前添加了 eH 和 eL 变量。这是一个成功运行的完整脚本。
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO(remote = True)
#initialize variables
#Room Temprature:
T_external = [23,23,23,23,23.5,23.5,23.4,23.5,23.9,23.7,\
23,23.9,23.9,23.4,23.9,24,23.6,23.7,23.8,\
23,23,23,23,23]
# Temprature Lower Limit:
temp_low = 10*np.ones(24)
# Temprature Upper Limit:
temp_upper = 12*np.ones(24)
#Hourly Energy prices:
TOU_v = [39.09,34.93,38.39,40.46,40.57,43.93,25,11,9,24,51.28,45.22,45.72,\
36,35.03,10,12,13,32.81,42.55,8,29.58,29.52,29.52]
###########################################
#System Identification:
#Time
t = np.linspace(0,10,117)
#State of the Fridge
ud = np.append(np.zeros(78) ,np.ones(39),0)
#Temprature Data
y = [14.600000000000001,14.600000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.8,14.8,14.8,14.8,14.8,14.8,14.8,14.8,\
14.8,14.8,14.9,14.9,14.9,14.9,14.9,14.9,14.9,15,15,15,15,15,15,15,15,15,15,15,15,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,
15,15,15,15,15,15,15,15,15,15,14.9,14.9,14.9,14.9,14.8,14.9,14.8,14.8,14.8,14.8,14.8,14.8,\
14.8,14.700000000000001,14.8,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.600000000000001,14.600000000000001,14.600000000000001,\
14.600000000000001,14.600000000000001,14.60]
na = 1 # output coefficients
nb = 1 # input coefficients
print('Identification')
yp,p,K = m.sysid(t,ud,y,na,nb,objf=10000,scale=False,diaglevel=1)
#create control ARX model:
y = m.Array(m.CV,1)
uc = m.Array(m.MV,1)
m.arx(p,y,uc)
# rename CVs
T= y[0]
# rename MVs
uc = uc[0]
###########################################
#Parameter
P = m.Param(value =100) #power
TL = m.Param(value=temp_low[0])
TH = m.Param(value=temp_upper[0])
c = m.Param(value=TOU_v[0])
# Manipilated variable:
u = m.MV(lb=0, ub=1, integer=True)
u.STATUS = 1 # allow optimizer to change the variable to attein the optimum.
# Controlled Variable (Affected with changes in the manipulated variable)
# Soft constraints on temprature.
eH = m.CV(value=0)
eL = m.CV(value=0)
eH.SPHI=0 #Set point high for linear error model.
eH.WSPHI=100 #Objective function weight on upper set point for linear error model.
eH.WSPLO=0 # Objective function weight on lower set point for linear error model
eH.STATUS =1 # eH : Error is considered in the objective function.
eL.SPLO=0
eL.WSPHI=0
eL.WSPLO=100
eL.STATUS = 1
#Linear error (Deviation from the limits)
m.Equations([eH==T-TH,eL==T-TL])
#Objective : minimize the costs.
m.Minimize(c*P*u)
#Optimizer Options.
# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)
TL.value = temp_low
TH.value = temp_upper
c.value = TOU_v
T.value = 11 # Temprature starts at 11
m.options.IMODE = 6 # MPC mode in Gekko.
m.options.NODES = 2 # Collocation nodes.
m.options.SOLVER = 1 # APOT solver for mixed integer linear programming.
m.time = np.linspace(0,23,24)
#Solve the optimization problem.
m.solve()
控制器输出如下：
--------- APM Model Size ------------
Each time step contains
Objects : 1
Constants : 0
Variables : 9
Intermediates: 0
Connections : 2
Equations : 3
Residuals : 3
Number of state variables: 1035
Number of total equations: - 1012
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 23
----------------------------------------------
Dynamic Control with APOPT Solver
----------------------------------------------
Iter: 1 I: 0 Tm: 0.07 NLPi: 3 Dpth: 0 Lvs: 0 Obj: 6.76E+03 Gap: 0.00E+00
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 8.319999999366701E-002 sec
Objective : 6763.77971670735
Successful solution
---------------------------------------------------
我正在对 MPC 进行建模以控制冰箱并将温度保持在给定的时间间隔内，同时最大限度地降低成本。我正在使用 GEKKO 为我的算法建模。
我写了下面的代码。首先，我使用系统中的传感器数据识别了我的模型（我使用了 GEKKO 的函数 sysif）。然后我构建了一个 ARX 模型（使用 GEKKO 中的 arx 函数），它成为 sysid() 作为输入的结果。
我正在尝试编写一个“虚拟”算法以在将其实现到 Pi 之前在本地进行测试。
我收到以下错误：
KeyError Traceback (most recent call last)
<ipython-input-13-108148376700> in <module>
107 #Solve the optimization problem.
108
--> 109 m.solve()
~/opt/anaconda3/lib/python3.8/site-packages/gekko/gekko.py in solve(self, disp, debug, GUI, **kwargs)
2214 if timing == True:
2215 t = time.time()
-> 2216 self.load_JSON()
2217 if timing == True:
2218 print('load JSON', time.time() - t)
~/opt/anaconda3/lib/python3.8/site-packages/gekko/gk_post_solve.py in load_JSON(self)
48 vp.__dict__[o] = dpred
49 else: #everything besides value, dpred and pred
---> 50 vp.__dict__[o] = data[vp.name][o]
51 for vp in self._variables:
52 if vp.type != None: #(FV/MV/SV/CV) not Param or Var
KeyError: 'int_p6'
这是我的代码
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO(remote = True)
#initialize variables
#Room Temprature:
T_external = [23,23,23,23,23.5,23.5,23.4,23.5,23.9,23.7,\
23,23.9,23.9,23.4,23.9,24,23.6,23.7,23.8,\
23,23,23,23,23]
# Temprature Lower Limit:
temp_low = 10*np.ones(24)
# Temprature Upper Limit:
temp_upper = 12*np.ones(24)
#Hourly Energy prices:
TOU_v = [39.09,34.93,38.39,40.46,40.57,43.93,25,11,9,24,51.28,45.22,45.72,\
36,35.03,10,12,13,32.81,42.55,8,29.58,29.52,29.52]
###########################################
#System Identification:
#Time
t = np.linspace(0,10,117)
#State of the Fridge
ud = np.append(np.zeros(78) ,np.ones(39),0)
#Temprature Data
y = [14.600000000000001,14.600000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.8,14.8,14.8,14.8,14.8,14.8,14.8,14.8,\
14.8,14.8,14.9,14.9,14.9,14.9,14.9,14.9,14.9,15,15,15,15,15,15,15,15,15,15,15,15,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,
15,15,15,15,15,15,15,15,15,15,14.9,14.9,14.9,14.9,14.8,14.9,14.8,14.8,14.8,14.8,14.8,14.8,\
14.8,14.700000000000001,14.8,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.600000000000001,14.600000000000001,14.600000000000001,\
14.600000000000001,14.600000000000001,14.60]
na = 1 # output coefficients
nb = 1 # input coefficients
print('Identification')
yp,p,K = m.sysid(t,ud,y,na,nb,objf=10000,scale=False,diaglevel=1)
#create control ARX model:
y = m.Array(m.CV,1)
uc = m.Array(m.MV,1)
m.arx(p,y,uc)
# rename CVs
T= y[0]
# rename MVs
uc = uc[0]
# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)
###########################################
#Parameter
P = m.Param(value =100) #power
TL = m.Param(value=temp_low)
TH = m.Param(value=temp_upper)
c = m.Param(value=TOU_v)
# Manipilated variable:
u = m.MV(lb=0, ub=1, integer=True)
u.STATUS = 1 # allow optimizer to change the variable to attein the optimum.
# Controlled Variable (Affected with changes in the manipulated variable)
T = m.CV(value=11) # Temprature will start at 11.
# Soft constraints on temprature.
eH = m.CV(value=0)
eL = m.CV(value=0)
eH.SPHI=0 #Set point high for linear error model.
eH.WSPHI=100 #Objective function weight on upper set point for linear error model.
eH.WSPLO=0 # Objective function weight on lower set point for linear error model
eH.STATUS =1 # eH : Error is considered in the objective function.
eL.SPLO=0
eL.WSPHI=0
eL.WSPLO=100
eL.STATUS = 1
#Linear error (Deviation from the limits)
m.Equations([eH==T-TH,eL==T-TL])
#Objective : minimize the costs.
m.Minimize(c*P*u)
#Optimizer Options.
m.options.IMODE = 6 # MPC mode in Gekko.
m.options.NODES = 2 # Collocation nodes.
m.options.SOLVER = 1 # APOT solver for mixed integer linear programming.
m.time = np.linspace(0,23,24)
#Solve the optimization problem.
m.solve()
问题在于：
T = m.CV(value=11) # Temperature will start at 11.
您正在重新定义 T 变量，但它在内部存储了这两个变量。如果你需要 re-initialize 到 11 然后使用 T.value=11。另外，我在稳态初始化之前添加了 eH 和 eL 变量。这是一个成功运行的完整脚本。
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO(remote = True)
#initialize variables
#Room Temprature:
T_external = [23,23,23,23,23.5,23.5,23.4,23.5,23.9,23.7,\
23,23.9,23.9,23.4,23.9,24,23.6,23.7,23.8,\
23,23,23,23,23]
# Temprature Lower Limit:
temp_low = 10*np.ones(24)
# Temprature Upper Limit:
temp_upper = 12*np.ones(24)
#Hourly Energy prices:
TOU_v = [39.09,34.93,38.39,40.46,40.57,43.93,25,11,9,24,51.28,45.22,45.72,\
36,35.03,10,12,13,32.81,42.55,8,29.58,29.52,29.52]
###########################################
#System Identification:
#Time
t = np.linspace(0,10,117)
#State of the Fridge
ud = np.append(np.zeros(78) ,np.ones(39),0)
#Temprature Data
y = [14.600000000000001,14.600000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.8,14.8,14.8,14.8,14.8,14.8,14.8,14.8,\
14.8,14.8,14.9,14.9,14.9,14.9,14.9,14.9,14.9,15,15,15,15,15,15,15,15,15,15,15,15,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,\
15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,15.100000000000001,
15,15,15,15,15,15,15,15,15,15,14.9,14.9,14.9,14.9,14.8,14.9,14.8,14.8,14.8,14.8,14.8,14.8,\
14.8,14.700000000000001,14.8,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.700000000000001,14.700000000000001,14.700000000000001,\
14.700000000000001,14.600000000000001,14.600000000000001,14.600000000000001,\
14.600000000000001,14.600000000000001,14.60]
na = 1 # output coefficients
nb = 1 # input coefficients
print('Identification')
yp,p,K = m.sysid(t,ud,y,na,nb,objf=10000,scale=False,diaglevel=1)
#create control ARX model:
y = m.Array(m.CV,1)
uc = m.Array(m.MV,1)
m.arx(p,y,uc)
# rename CVs
T= y[0]
# rename MVs
uc = uc[0]
###########################################
#Parameter
P = m.Param(value =100) #power
TL = m.Param(value=temp_low[0])
TH = m.Param(value=temp_upper[0])
c = m.Param(value=TOU_v[0])
# Manipilated variable:
u = m.MV(lb=0, ub=1, integer=True)
u.STATUS = 1 # allow optimizer to change the variable to attein the optimum.
# Controlled Variable (Affected with changes in the manipulated variable)
# Soft constraints on temprature.
eH = m.CV(value=0)
eL = m.CV(value=0)
eH.SPHI=0 #Set point high for linear error model.
eH.WSPHI=100 #Objective function weight on upper set point for linear error model.
eH.WSPLO=0 # Objective function weight on lower set point for linear error model
eH.STATUS =1 # eH : Error is considered in the objective function.
eL.SPLO=0
eL.WSPHI=0
eL.WSPLO=100
eL.STATUS = 1
#Linear error (Deviation from the limits)
m.Equations([eH==T-TH,eL==T-TL])
#Objective : minimize the costs.
m.Minimize(c*P*u)
#Optimizer Options.
# steady state initialization
m.options.IMODE = 1
m.solve(disp=True)
TL.value = temp_low
TH.value = temp_upper
c.value = TOU_v
T.value = 11 # Temprature starts at 11
m.options.IMODE = 6 # MPC mode in Gekko.
m.options.NODES = 2 # Collocation nodes.
m.options.SOLVER = 1 # APOT solver for mixed integer linear programming.
m.time = np.linspace(0,23,24)
#Solve the optimization problem.
m.solve()
控制器输出如下：
--------- APM Model Size ------------
Each time step contains
Objects : 1
Constants : 0
Variables : 9
Intermediates: 0
Connections : 2
Equations : 3
Residuals : 3
Number of state variables: 1035
Number of total equations: - 1012
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 23
----------------------------------------------
Dynamic Control with APOPT Solver
----------------------------------------------
Iter: 1 I: 0 Tm: 0.07 NLPi: 3 Dpth: 0 Lvs: 0 Obj: 6.76E+03 Gap: 0.00E+00
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 8.319999999366701E-002 sec
Objective : 6763.77971670735
Successful solution
---------------------------------------------------