Python Gekko 微分方程时域问题中的逐步变量

Question

我正在 Python Gekko 中编写一个 MPC 问题，以便在冬日为建筑物供暖（目前我正在处理一栋只有 2 个区域的简单建筑物）。该建筑用电阻电容 (RC) 模型表示，目标是最大限度地减少最大功率和总能耗的组合。 RC模型通过微分方程来表达。通过使用区域温度作为操纵变量 (MV) 来最小化目标，它已经可以正常工作。

但是为了提高准确性并将其更好地集成到建筑物中，我现在在 MPC 问题中实现了 PID（因为建筑物使用 PID 恒温器，而我使用 MPC 只是为了向 PID 控制器提供更好的设定值）。这意味着 MPC 将更改 PID 设定点以最小化目标函数。在这种情况下，MV 不是区域温度，而是 PID 设定点。区域温度将尝试通过 PID 控制遵循设定值。

到目前为止效果很好，但我需要 PID 设置点每 2 小时而不是 15 分钟（这是我正在使用的时间步长）更改一次。 MPC 必须使用 15 分钟的时间步长，因为 2 小时就太大了。这意味着 MV 应该具有逐步形状，其值只能每 8 个时间步改变一次。我正在尝试使用设置点 (SP) 的 m.Array 变量（m.time 的长度为 97）来尝试每 8 个时间步强加 12 个不同的固定值 (FV) (SP[0:8] = FV_1 ，SP[8:16] = FV_2 等）但我无法解决我的问题。我在 Stack Overflow 上搜索了多个帖子，但没有找到类似的案例。

MPC中带有PID的工作代码是：

m = GEKKO(remote=False) # create GEKKO model
n_days_mpc = 1
starting_datetime_mpc = pd.Timestamp("2015-02-02 00:00", tz="US/Eastern")
finishing_datetime_mpc = starting_datetime_mpc +  timedelta(hours = (24*n_days_mpc))
dates_mpc = pd.date_range(starting_datetime_mpc, finishing_datetime_mpc,freq='15min',tz="US/Eastern")
time_array = np.arange(0,len(dates_mpc)*timestep,timestep)
m.time = time_array

# import input data
df, location, solar_position, Q_solar_df = weather_data (starting_datetime_mpc, finishing_datetime_mpc)
U_test = Power_delivered[starting_datetime_mpc:finishing_datetime_mpc].values
X_test = Temp[starting_datetime_mpc:finishing_datetime_mpc].values
W_cols = ["T_ext", "Tground_1", "Tground_2", 
          "Q_sol_1", "Q_sol_2", "Q_sol_roof",
          "Q_int_1", "Q_int_2", "wind_sq"]  # exogenous_columns
W_test = df[W_cols].values

# reference performances from TRNSYS
power_trnsys = U_test[:,0] + U_test[:,1]
Q_int_trnsys = W_test[:,7] + W_test[:,8]
Q_sol_direct_trnsys = W_test[:,3] + W_test[:,4]

# comfort schedule
df["is_occupied"] = (7. <= df.index.hour) & (df.index.hour < 18) # comfort schedules
heating_SP_unoccupied = 19.0 # Setpoints when occupied/unoccupied
heating_SP_occupied = 21.0
df["heating_SP"] = heating_SP_unoccupied
df.loc[df["is_occupied"], "heating_SP"] = heating_SP_occupied # Update for occupied setpoint
heating_SP = df["heating_SP"].values
T_comfort_21 = m.Param(value=heating_SP)

# PID parameters
Kc = 30000/3600             # controller gain
tauI = 1000                 # controller reset time
tauD = 0                    # derivative constant
Q1_0 = m.Const(value=0.0)   # OP bias Zone 1
Q2_0 = m.Const(value=0.0)   # OP bias Zone 2

# change solver options
m.options.IMODE = 6 # MPC

# parameters from building training
for i in params_list:
    locals()[i] = m.Const(value = locals()[i+'_results'][-1])

# initial conditions
T_wall_int_initial = np.ones(len(dates_mpc))*20
initial_thermal_power = U_test[:,0][0] + U_test[:,1][0]

# state variables (temperature)
T_1 = m.SV(value=X_test[:,0][0]);
T_2 = m.SV(value=X_test[:,1][0]);
T_wall_int = m.SV(value=T_wall_int_typique);

# manipulated variables (PID thermostat setpoints)
SP_1 = m.MV(value=heating_SP[0],fixed_initial=False); SP_1.STATUS = 1   # set point
SP_2 = m.MV(value=heating_SP[0],fixed_initial=False); SP_2.STATUS = 1   # set point
Intgl_1 = m.Var((heating_SP[0]-X_test[:,0][0])*tauI*Kc)    # integral of the error Zone 1
Intgl_2 = m.Var((heating_SP[0]-X_test[:,1][0])*tauI*Kc)    # integral of the error Zone 2

# controlled variables (heating)
Q_1 = m.Var(value=U_test[:,0][0],lb=0);
Q_2= m.Var(value=U_test[:,1][0],lb=0);
thermal_power = m.Var(value=initial_thermal_power)
max_power = m.FV(value=initial_thermal_power); max_power.STATUS = 1

# uncontrolled variables (weather and internal heat gains)
T_ext = m.Param(value=W_test[:,0])
Tground_1 = m.Param(value=W_test[:,1])
Tground_2 = m.Param(value=W_test[:,2])
Q_sol_1 = m.Param(value=W_test[:,3])
Q_sol_2 = m.Param(value=W_test[:,4])
Q_sol_roof = m.Param(value=W_test[:,5])
Q_int_1 = m.Param(value=W_test[:,6])
Q_int_2 = m.Param(value=W_test[:,7])
wind_sq = m.Param(value=W_test[:,8])

# building equations (RC grey box)
m.Equation(C_1 *1e6* T_1.dt() == (Q_1 + lambda_1 * Q_int_1 + sigma_1 * Q_sol_1 + sigma_roof * Q_sol_roof + alpha_1 * wind_sq * (T_ext - T_1) \
                    + Uext_1*(T_ext - T_1) + Uroof_1*(T_ext - T_1) + Ug_1*(Tground_1 - T_1) \
                    + U_wall_int*(T_wall_int - T_1) + U_12*(T_2- T_1)  ))

m.Equation(C_2 *1e6* T_2.dt() == (Q_2 + lambda_2 * Q_int_2 + sigma_2 * Q_sol_2 + sigma_roof * Q_sol_roof + alpha_2 * wind_sq * (T_ext - T_2) \
                    + Uext_2*(T_ext - T_2) + Uroof_2*(T_ext - T_2) + Ug_2*(Tground_2 - T_2) \
                    + U_wall_int*(T_wall_int - T_2) + U_12*(T_1- T_2) ))

m.Equation(C_wall_int *1e6* T_wall_int.dt() == U_wall_int*(T_1 - T_wall_int) + U_wall_int*(T_2 - T_wall_int))

# PID thermostats
err_1 = m.Intermediate(SP_1-T_1) # set point error
err_2 = m.Intermediate(SP_2-T_2) # set point error
m.Equation(Intgl_1.dt()==err_1) # integral of the error
m.Equation(Q_1 == Q1_0 + Kc*err_1 + (Kc/tauI)*Intgl_1 - (Kc*tauD)*T_1.dt())
m.Equation(Intgl_2.dt()==err_2) # integral of the error
m.Equation(Q_2 == Q2_0 + Kc*err_2 + (Kc/tauI)*Intgl_2 - (Kc*tauD)*T_2.dt())

m.Equation(thermal_power == Q_1+Q_2)
m.Equation(max_power >= thermal_power)

m.Equation(SP_1 >= T_comfort_21)
m.Equation(SP_2 >= T_comfort_21)

m.Minimize(14.58*max_power + 5.03/100*thermal_power/ (60/minutes_per_timestep))
m.solve()

我想做的是（我在这里只报告我修改的部分）：

# manipulated variables (PID thermostat setpoints)
SP_1 = m.Array(m.Var, len(m.time))
SP_2 = m.Array(m.Var, len(m.time))
SP_1_values = m.Array(m.FV,12)
SP_2_values = m.Array(m.FV,12)
for SP in SP_1_values:
    SP.STATUS = 1
for SP in SP_2_values:
    SP.STATUS = 1
Intgl_1 = m.Var((heating_SP[0]-X_test[:,0][0])*tauI*Kc)    # integral of the error Zone 1
Intgl_2 = m.Var((heating_SP[0]-X_test[:,1][0])*tauI*Kc)    # integral of the error Zone 2

# PID thermostats
for i in (0,1,2,3,4,5,6,7,8):
    m.Equation(SP_1[i] == SP_1_values[0])
    m.Equation(SP_2[i] == SP_2_values[0])
for j in range(1, 12):
    start_index = 8 * j + 1
    indices = tuple(range(start_index, start_index + 8))
    for i in indices:
        m.Equation(SP_1[i] == SP_1_values[j])
        m.Equation(SP_2[i] == SP_2_values[j])
err_1 = m.Intermediate(SP_1-T_1) # set point error
err_2 = m.Intermediate(SP_2-T_2) # set point error
m.Equation(Intgl_1.dt()==err_1) # integral of the error
m.Equation(Q_1 == Q1_0 + Kc*err_1 + (Kc/tauI)*Intgl_1 - (Kc*tauD)*T_1.dt())
m.Equation(Intgl_2.dt()==err_2) # integral of the error
m.Equation(Q_2 == Q2_0 + Kc*err_2 + (Kc/tauI)*Intgl_2 - (Kc*tauD)*T_2.dt())

m.Equation(thermal_power == Q_1+Q_2)
m.Equation(max_power >= thermal_power)

for i in range(len(m.time)):
    m.Equation(SP_1[i] >= T_comfort_21)
    m.Equation(SP_2[i] >= T_comfort_21)

我首先尝试将 PID 设置点仅定义为 m.Var，但后来我无法使用 SP_1[i]，因为它说“标量变量的索引无效”，所以这就是我使用数组然后尝试插入的原因类似 SP_1[0:9] = 一个可以更改的 FV 参数，然后 SP_1[9:17] = 另一个可以更改的 FV 参数，依此类推，因此 SP_1 是一个逐步的 多变的。但现在 m.Intermediate 行

err_1 = m.Intermediate(SP_1-T_1)

给出了此错误：

@error: Intermediate Definition Error: Intermediate variable with no equality (=) expression

。我想这是因为现在 SP_1 需要在那里建立索引（因为它是一个 m.Array），但是我不能对 T_1 做同样的事情，因为它需要保留为状态变量（SV），因为它受微分方程的影响，然后PID的积分误差也受微分方程的影响。所以我基本上被困住了，不知道该怎么办。

我想问的是：我想实现的逐步变量是否可以实现？我该怎么做？

感谢所有回复这篇文章的人。如果问题不清楚，我可以尝试澄清。另外，如果输入数据来计算代码有用，我可以提供它。

Answer 1

尝试使用

MV

类型的选项，指定操纵变量每个允许的移动之间的时间步数。

SP_1.MV_STEP_HOR=8
SP_2.MV_STEP_HOR=8

在动态优化课程和文档中可以找到有关MPC调整选项的更多信息。

可以为所有 MV 设置全局选项（例如

m.options.MV_STEP_HOR=8

），也可以单独设置它们（例如

SP_1.MV_STEP_HOR=6

、

SP_1.MV_STEP_HOR=10

）。默认情况下，当使用全局选项时，各个 MV 会设置为

0

。