请求将电池产热计算的Pyomo非线性约束线性化（适配CBC求解器）

首先，咱们得明确核心问题：你的原方程q_batt = R(capacity) * (flow / V)^2里包含了平方项和容量与平方项的乘积项，这两个都是非线性的，而CBC作为线性规划（LP）求解器无法直接处理。你提到的对数法其实在这里并不适用——因为R(capacity)是容量的线性函数，它的对数ln(R(c))是非线性的，取对数后依然无法得到全线性的方程。

下面是更实用的辅助变量+McCormick松弛线性化方案，这是LP/MILP中处理二次非线性项的标准方法：

步骤1：拆解原方程并定义关键参数 #

首先把原方程展开，因为R(capacity)是线性的，我们可以把它写成：
R(c) = a*c + b
其中：

• a是内阻随容量变化的斜率
• b是内阻的截距常量
• c是电池容量变量block.capacity[n,t]
• V是常量n.inflow_voltage
• flow是变量m.flow[n, i[n], t]

原方程因此变为：
q_batt = (a*c + b) * (flow/V)^2 = b*(flow/V)^2 + a*c*(flow/V)^2

这里有两个非线性项需要处理：(flow/V)^2和c*(flow/V)^2。

步骤2：设定变量上下界 #

线性化的前提是给所有变量设定合理的上下界（这会直接影响线性化的精度）：

• 给flow设定上下界：flow_min（比如放电的最大电流，负数）和flow_max（充电的最大电流，正数）
• 给capacity设定上下界：c_min=0（空容量）和c_max（电池额定容量）
• 推导current = flow/V的上下界：current_min = flow_min/V，current_max = flow_max/V
• 推导current_sq = (flow/V)^2的上下界：csq_min = current_min²，csq_max = current_max²

步骤3：引入辅助变量并线性化 #

我们需要引入两个辅助变量来替代非线性项：

1. current_sq[n,t]：替代(flow/V)^2
2. c_current_sq[n,t]：替代c*(flow/V)^2

3.1 线性化current_sq = current * current #

用McCormick松弛给current_sq添加线性约束：

def _current_sq_linear_rule(block, n, t):
    flow_min = n.flow_min
    flow_max = n.flow_max
    V = n.inflow_voltage
    current_min = flow_min / V
    current_max = flow_max / V

    constraints = []
    # 下界约束：取两个线性下界的最大值
    constraints.append(block.current_sq[n,t] >= 2*current_min*block.current[n,t] - current_min**2)
    constraints.append(block.current_sq[n,t] >= 2*current_max*block.current[n,t] - current_max**2)
    # 上界约束：根据current的符号选择不同的线性上界
    if current_min * current_max >= 0:
        # current同号，用线性组合约束
        constraints.append(block.current_sq[n,t] <= current_min*block.current[n,t] + current_max*block.current[n,t] - current_min*current_max)
    else:
        # current可正可负，上界取最大平方值
        constraints.append(block.current_sq[n,t] <= max(current_min**2, current_max**2))
    return constraints

# 绑定约束到模型
self.current_sq_constraint = Constraint(self.STORAGES, reduced_timesteps, rule=_current_sq_linear_rule)

3.2 线性化c_current_sq = c * current_sq #

同样用McCormick松弛处理容量与平方项的乘积：

def _c_current_sq_linear_rule(block, n, t):
    c_min = 0
    c_max = n.capacity_max
    flow_min = n.flow_min
    flow_max = n.flow_max
    V = n.inflow_voltage
    current_min = flow_min / V
    current_max = flow_max / V
    csq_min = current_min**2
    csq_max = current_max**2

    constraints = []
    # 下界约束：取四个线性下界的最大值
    constraints.append(block.c_current_sq[n,t] >= c_min*block.current_sq[n,t] + csq_min*block.capacity[n,t] - c_min*csq_min)
    constraints.append(block.c_current_sq[n,t] >= c_min*block.current_sq[n,t] + csq_max*block.capacity[n,t] - c_min*csq_max)
    constraints.append(block.c_current_sq[n,t] >= c_max*block.current_sq[n,t] + csq_min*block.capacity[n,t] - c_max*csq_min)
    constraints.append(block.c_current_sq[n,t] >= c_max*block.current_sq[n,t] + csq_max*block.capacity[n,t] - c_max*csq_max)
    # 上界约束：取两个线性上界的最小值
    constraints.append(block.c_current_sq[n,t] <= c_min*block.current_sq[n,t] + csq_max*block.capacity[n,t] - c_min*csq_max)
    constraints.append(block.c_current_sq[n,t] <= c_max*block.current_sq[n,t] + csq_min*block.capacity[n,t] - c_max*csq_min)
    return constraints

# 绑定约束到模型
self.c_current_sq_constraint = Constraint(self.STORAGES, reduced_timesteps, rule=_c_current_sq_linear_rule)

步骤4：改写原产热约束为线性形式 #

先从n.inflow_resistance函数中提取线性参数a和b（可以通过两个容量点计算）：

# 示例：从线性内阻函数中提取斜率和截距
c1 = 0
R1 = n.inflow_resistance(c1)
c2 = n.capacity_max
R2 = n.inflow_resistance(c2)
a = (R2 - R1) / (c2 - c1)
b = R1 - a * c1

然后写出线性化后的产热约束：

def _q_batt_linear_rule(block, n, t):
    # 提前提取的线性参数a和b
    a = n.inflow_resistance_slope
    b = n.inflow_resistance_intercept
    # 线性化后的产热方程
    expr = block.q_batt[n,t] - b*block.current_sq[n,t] - a*block.c_current_sq[n,t]
    return expr == 0

# 绑定约束到模型
self.q_batt_linear_constraint = Constraint(self.STORAGES, reduced_timesteps, rule=_q_batt_linear_rule)

步骤5：补充变量定义 #

别忘了在模型中定义新增的辅助变量：

# 定义current变量（flow/V的别名，方便计算）
self.current = Var(self.STORAGES, reduced_timesteps, domain=Reals)
# 定义current平方的辅助变量
self.current_sq = Var(self.STORAGES, reduced_timesteps, domain=NonNegativeReals)
# 定义容量与current平方乘积的辅助变量
self.c_current_sq = Var(self.STORAGES, reduced_timesteps, domain=NonNegativeReals)

# 绑定current的定义约束
def _current_def_rule(block, n, t):
    expr = block.current[n,t] - m.flow[n, i[n], t] / n.inflow_voltage
    return expr == 0
self.current_def_constraint = Constraint(self.STORAGES, reduced_timesteps, rule=_current_def_rule)

关键注意事项 #

1. 上下界的合理性：变量上下界越窄，线性化的精度越高，求解器的性能也越好。
2. 符号处理：如果你的模型只考虑充电（flow正）或只考虑放电（flow负），可以简化current的上下界，减少约束数量。
3. 精度与速度平衡：McCormick松弛会引入更多约束，可能增加求解时间，但这是在LP中处理二次项的精确松弛方法（可行域包含原问题的可行域）。

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