Pyomo与EnergyPlus协同优化：迭代更新EnergyPlus数值的方法咨询

为啥现在只调用一次？ #

你当前代码里把室内温度定义成了Expression，Pyomo在构建这个表达式的时候，会用power_w的初始值（5000）调用一次Eplus算出温度，之后优化迭代时变量变了，Pyomo不会自动重新跑Eplus——因为Eplus的温度计算是外部黑箱，不是Pyomo能解析的代数运算，它不知道这个温度值会跟着变量变化。

怎么改才能让它每次迭代都更新温度？ #

核心是把Eplus的温度计算包装成Pyomo求解器能识别的外部黑箱函数，让求解器每次需要约束值的时候，都主动调用Eplus重新计算温度。下面是修改后的完整代码：

import pyomo.environ as pyo
import numpy as np
import matplotlib.pyplot as plt
from EPlus.model import TemperatureModel
from pyomo.opt import SolverFactory

# 初始化EnergyPlus模型
EPLUS_model = TemperatureModel("C:\\EnergyPlusV9-4-0")

# 创建Pyomo模型
model = pyo.ConcreteModel()

# 时间集合
model.time = pyo.Set(initialize=(range(30, 31 + 1, 30)))

# 优化变量：输入功率
model.power_w = pyo.Var(model.time, within=pyo.NonNegativeReals, bounds=(0, 40000), initialize=5000)

# 定义温度计算函数：输入功率，返回对应室内温度
def get_indoor_temp(power):
    # 调用EnergyPlus模拟，取第一个温度结果
    return EPLUS_model.T(np.zeros(96) + power)[0]

# 把这个函数注册成Pyomo的外部函数，让求解器能识别
model.calculate_temp = pyo.ExternalFunction(get_indoor_temp)

# 温度下限约束：直接用当前功率调用外部函数算温度，确保每次都更新
def temp_limit_rule(model, t):
    return model.calculate_temp(model.power_w[t]) >= 13
model.temp_constraint = pyo.Constraint(model.time, rule=temp_limit_rule)

# 目标函数：最小化总功耗
def minimize_power(model):
    return sum(model.power_w[t] for t in model.time)
model.total_power = pyo.Objective(rule=minimize_power, sense=pyo.minimize)

# 用IPOPT求解，开tee=True能看日志，确认Eplus是不是多次调用了
solver = SolverFactory('ipopt')
result = solver.solve(model, tee=True)

# 提取结果
solved_p = {t: model.power_w[t].value for t in model.time}
solved_t = {t: model.calculate_temp(model.power_w[t]()).value for t in model.time}

print("求解后的功率值：", solved_p)
print("对应的室内温度：", solved_t)

关键改动说明 #

• 换掉Expression，用ExternalFunction：把温度计算逻辑包成外部函数后，IPOPT这类求解器在迭代时，每次需要计算约束值，都会调用这个函数，自然就会重新跑Eplus。
• 直接在约束里调用外部函数：不再通过中间的Expression传递温度，避免Pyomo提前计算后就不再更新。
• 开启tee看日志：可以在求解过程里看到Eplus是不是被多次调用，方便调试。

额外提醒 #

如果Eplus模拟很慢，频繁调用会拖慢优化速度，可以试试这两个办法：

• 加个缓存：只有当功率变化超过某个阈值（比如100W）时，才重新跑Eplus，否则用之前的结果
• 先搞个代理模型：用响应面、神经网络之类的快速模型替代Eplus做优化，收敛后再用Eplus验证最终结果
• 梯度问题：IPOPT默认需要梯度信息，如果你的Eplus没法直接输出梯度，可以让IPOPT用有限差分近似梯度，加一行solver.options['finite_diff_rel_step'] = 1e-3就行

内容的提问来源于stack exchange，提问作者Brown