Pyomo中带负截距电池充电性能分段线性实现问询

问题背景 #

在能源优化场景中，尝试用Pyomo的Piecewise库实现电池充电性能的分段线性表达式，充电性能函数为：
y = 0.6*x - 0.2*x² - 0.01
其中y为充电输出功率，x为充电输入功率。该函数存在负截距，低x值时y为负，导致电池无需充电时x被强制设为约0.0125才能使y=0（y(0.0125)=0）。

现有两个问题：

1. 能否在Pyomo.Piecewise中实现这类带负截距的方程，同时使用无二进制变量的简化线性表示？
2. 能否仅在y>0的区间定义该方程，同时支持充电功率设为0的非运行场景？

原实现忽略截距时可正常运行，但带截距的场景未找到可行方案，原代码如下：

import pyomo.environ as pyo
import random

## Generate some random data for PV and Load
pv = [random.randint(0, 5) for _ in range(48)]
pv_dict = (dict(enumerate(pv,1)))

load_el = [random.randint(0, 8) for _ in range(48)]
load_el_dict = (dict(enumerate(load_el,1)))

# Define model
model = pyo.ConcreteModel()

# Define timeperiod set
model.T = pyo.RangeSet(len(pv_dict))

# Define model parameters
model.pv = pyo.Param(model.T, initialize=pv_dict)
model.load_el= pyo.Param(model.T, initialize=load_el_dict)
model.grid_cost_buy = pyo.Param(model.T, initialize=0.4)

model.battery_eoCH = pyo.Param(initialize=1.0)
model.battery_eoDCH = pyo.Param(initialize=0.1)
model.battery_capacity = pyo.Param(initialize=5)
model.battery_power_max = pyo.Param(initialize=100)

# Define the variables             
model.battery_soc = pyo.Var(model.T, bounds=(model.battery_eoDCH, model.battery_eoCH))     # battery soc with end of ch/DCH levels 
model.grid_power_import = pyo.Var(model.T, domain=pyo.NonNegativeReals)        # grid import power 
model.grid_power_export = pyo.Var(model.T, domain=pyo.NonNegativeReals)        # grid export power 
model.battery_power_DCH = pyo.Var(model.T, domain=pyo.NonNegativeReals)        # battery discharging power
# PWA variables
model.battery_power_CH_in = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.battery_power_CH_out = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.battery_power_CH_in_norm = pyo.Var(model.T, domain=pyo.NonNegativeReals, bounds=(0,1))
model.battery_power_CH_out_norm = pyo.Var(model.T, domain=pyo.NonNegativeReals)

# Linearization of charge efficiency

# Define function for PWA
def f(model,t,x):
    # Normalized charge performance function y=P charge out and x=P charge in
    y = (0.6*x - 0.2*x**2 - 0.01)
    return y
# Define breakpoints
breakpoints = [0.0, 0.3, 0.6, 1.0]
# Create breakpoints dict with same index as variables
PW_PTS = {}
for idx in model.battery_power_CH_in_norm.index_set():
    PW_PTS[idx] = breakpoints
    
# Define the PWA function
model.battery_power_CH_PWA_func = pyo.Piecewise(model.T, 
                                                model.battery_power_CH_out_norm, #y
                                                model.battery_power_CH_in_norm, #x
                                                pw_pts=PW_PTS, 
                                                f_rule=f, 
                                                pw_constr_type='UB', 
                                                pw_repn='CC',
                                                force_pw=False)

# Change normalized values to absolute values
def battery_power_ch_in_rule(m,t):
    return (m.battery_power_CH_in[t] == model.battery_power_CH_in_norm[t] * m.battery_power_max)
model.battery_power_ch_in_rule_c = pyo.Constraint(model.T, rule=battery_power_ch_in_rule)

def battery_power_ch_out_rule(m,t):
    return (m.battery_power_CH_out[t] == model.battery_power_CH_out_norm[t] * m.battery_power_max)
model.battery_power_ch_out_rule_c = pyo.Constraint(model.T, rule=battery_power_ch_out_rule)

# Further battery constraints
# Battery SoC constraint
def battery_soc_rule(m, t):
    if t == m.T.first():
        return m.battery_soc[t] == ((m.battery_power_CH_out[t] - m.battery_power_DCH[t]) / model.battery_capacity)
    
    return m.battery_soc[t] == m.battery_soc[t-1] + ((m.battery_power_CH_out[t] - m.battery_power_DCH[t]) / model.battery_capacity)
model.battery_soc_c = pyo.Constraint(model.T, rule=battery_soc_rule)

# Define balanced electricity bus rule
def balanced_bus_rule(m, t):
    return (0 == (m.pv[t] - m.load_el[t] 
                 + m.battery_power_DCH[t] - m.battery_power_CH_in[t]
                 + m.grid_power_import[t] - m.grid_power_export[t]))
model.bus_c = pyo.Constraint(model.T, rule=balanced_bus_rule)

## Define the cost function
def obj_rule(m):
    return sum(m.grid_power_import[t]*m.grid_cost_buy[t] for t in m.T)
model.obj = pyo.Objective(rule=obj_rule, sense=1)

## Solve the problem
solver = pyo.SolverFactory('gurobi')
results = solver.solve(model)
print('Total operation costs:',pyo.value(model.obj))

问题解答 #

问题1：带负截距的无二进制分段线性实现 #

可以实现，核心是调整分段点和约束逻辑，避免y为负的冲突：

1. 修正分段点：加入y=0对应的x值（约0.0125），确保分段区间覆盖函数的有效区域（y≥0）。
2. 调整约束类型：使用pw_constr_type='EQ'保证分段线性表达式严格匹配原函数的有效部分，结合battery_power_CH_out_norm的非负约束，自动排除y为负的解。
3. 保留无二进制表示：继续使用pw_repn='CC'（凸组合表示），无需引入二进制变量，适配该函数在有效区间的特性。

问题2：仅在y>0区间定义方程并支持x=0 #

核心是建立x和y的逻辑关联：当x=0时y必须为0；当x>0时，x≥0.0125且y遵循原函数的分段线性近似。具体实现：

• 修改分段函数，在x∈[0,0.0125]时强制y=0，x∈[0.0125,1]时使用原函数计算y值。
• 无需额外二进制变量，通过分段点和约束的组合即可实现逻辑控制。

修改后的代码 #

import pyomo.environ as pyo
import random

## Generate some random data for PV and Load
pv = [random.randint(0, 5) for _ in range(48)]
pv_dict = dict(enumerate(pv, 1))

load_el = [random.randint(0, 8) for _ in range(48)]
load_el_dict = dict(enumerate(load_el, 1))

# Define model
model = pyo.ConcreteModel()

# Define timeperiod set
model.T = pyo.RangeSet(len(pv_dict))

# Define model parameters
model.pv = pyo.Param(model.T, initialize=pv_dict)
model.load_el = pyo.Param(model.T, initialize=load_el_dict)
model.grid_cost_buy = pyo.Param(model.T, initialize=0.4)

model.battery_eoCH = pyo.Param(initialize=1.0)
model.battery_eoDCH = pyo.Param(initialize=0.1)
model.battery_capacity = pyo.Param(initialize=5)
model.battery_power_max = pyo.Param(initialize=100)
# 计算y=0对应的临界x值
model.x_critical = pyo.Param(initialize=(0.6 - (0.6**2 - 4*(-0.2)*(-0.01))**0.5)/(2*(-0.2)))  # ≈0.0125

# Define the variables             
model.battery_soc = pyo.Var(model.T, bounds=(model.battery_eoDCH, model.battery_eoCH))
model.grid_power_import = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.grid_power_export = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.battery_power_DCH = pyo.Var(model.T, domain=pyo.NonNegativeReals)
# PWA变量
model.battery_power_CH_in = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.battery_power_CH_out = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.battery_power_CH_in_norm = pyo.Var(model.T, domain=pyo.NonNegativeReals, bounds=(0, 1))
model.battery_power_CH_out_norm = pyo.Var(model.T, domain=pyo.NonNegativeReals)

# 充电效率分段线性化

# 调整后的分段函数
def f(model, t, x):
    if x <= model.x_critical:
        return 0.0
    else:
        return 0.6*x - 0.2*x**2 - 0.01

# 定义包含临界x值的分段点
breakpoints = [0.0, model.x_critical, 0.3, 0.6, 1.0]
PW_PTS = {t: breakpoints for t in model.T}

# 定义分段线性函数：等式约束+凸组合表示（无二进制变量）
model.battery_power_CH_PWA_func = pyo.Piecewise(
    model.T,
    model.battery_power_CH_out_norm,  # y: 归一化充电输出
    model.battery_power_CH_in_norm,   # x: 归一化充电输入
    pw_pts=PW_PTS,
    f_rule=f,
    pw_constr_type='EQ',
    pw_repn='CC',
    force_pw=False
)

# 归一化值与绝对值的转换约束
def battery_power_ch_in_rule(m, t):
    return m.battery_power_CH_in[t] == m.battery_power_CH_in_norm[t] * m.battery_power_max
model.battery_power_ch_in_rule_c = pyo.Constraint(model.T, rule=battery_power_ch_in_rule)

def battery_power_ch_out_rule(m, t):
    return m.battery_power_CH_out[t] == m.battery_power_CH_out_norm[t] * m.battery_power_max
model.battery_power_ch_out_rule_c = pyo.Constraint(model.T, rule=battery_power_ch_out_rule)

# 电池SOC约束
def battery_soc_rule(m, t):
    if t == m.T.first():
        return m.battery_soc[t] == (m.battery_power_CH_out[t] - m.battery_power_DCH[t]) / m.battery_capacity
    return m.battery_soc[t] == m.battery_soc[t-1] + (m.battery_power_CH_out[t] - m.battery_power_DCH[t]) / m.battery_capacity
model.battery_soc_c = pyo.Constraint(model.T, rule=battery_soc_rule)

# 电力平衡约束
def balanced_bus_rule(m, t):
    return 0 == (m.pv[t] - m.load_el[t] 
                 + m.battery_power_DCH[t] - m.battery_power_CH_in[t]
                 + m.grid_power_import[t] - m.grid_power_export[t])
model.bus_c = pyo.Constraint(model.T, rule=balanced_bus_rule)

# 目标函数：最小化购电成本
def obj_rule(m):
    return sum(m.grid_power_import[t] * m.grid_cost_buy[t] for t in m.T)
model.obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize)

# 求解模型
solver = pyo.SolverFactory('gurobi')
results = solver.solve(model)
print('Total operation costs:', pyo.value(model.obj))

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