如何在Python中实现微电网双层优化时间序列问题？

一、先确认KKT条件的正确性 #

双层优化转化为单层的核心是内层问题的KKT条件必须完整，针对你提到的产消者最小化问题，完整的KKT条件应包含4部分：

• 梯度条件：内层目标函数的梯度等于各约束梯度与对应拉格朗日乘子的线性组合（等式约束乘子无符号限制，不等式约束乘子非负）
• 原始可行性：内层所有约束（等式/不等式）均满足
• 对偶可行性：不等式约束对应的拉格朗日乘子≥0
• 互补松弛：不等式约束的松弛量与对应乘子的乘积为0（即约束激活时乘子非零，约束非激活时乘子为0）

你提供的KKT条件需要对应以上四点逐一核对，重点注意：如果内层是最小化问题，梯度条件应为拉格朗日函数对所有内层变量的偏导等于0。

二、Python工具选择（替代YALMIP/Gekko） #

Gekko对双层优化的原生支持有限，但Python有两个更适合的工具：

1. Pyomo #

开源建模框架，支持多种商用/开源求解器（Gurobi、Ipopt、CBC等），适合构建带时间序列的复杂约束模型，手动添加KKT约束的灵活性极高，文档社区完善。

2. CasADi #

专注于最优控制和非线性优化的工具，可自动推导内层问题的KKT条件，无需手动编写梯度表达式，对时间序列动态优化的支持极佳，适合大规模问题求解。

三、编码核心步骤（以Pyomo为例） #

1. 定义变量体系 #

• 上层变量：时间序列维度的p_sell(t)（售电价）、p_dg(t)（DG收购价）、P_grid(t)（电网购电量）
• 内层变量：产消侧的DeltaP_dr(t)（需求响应调整量）、电池充放电量等
• KKT乘子：为每个内层约束定义对应乘子，如等式约束的乘子无边界，不等式约束的乘子设为非负

2. 嵌入KKT约束到上层模型 #

将内层问题的KKT条件作为约束添加到上层模型中，针对时间序列的每个时间步循环处理：

• 原始可行性：添加产消者的功率平衡、电池容量、需求响应上下限等约束
• 梯度条件：手动推导内层目标对各变量的偏导，等于约束偏导与乘子的乘积和
• 对偶可行性：设置不等式约束的乘子变量域为NonNegativeReals
• 互补松弛：对不等式约束g(x) ≤ 0，添加lambda_var[t] * g(x)[t] == 0；若求解器不支持严格互补，可改用大M法或0-1松弛变量

3. 上层目标与求解 #

上层目标为最小化总成本（购电成本+DG收购成本-售电收入），选择支持非线性/混合整数的求解器（如Gurobi处理混合整数，Ipopt处理非线性）求解。

四、处理KKT的lambda变量 #

• lambda是对偶决策变量，直接在Pyomo中定义为Var类型，根据约束类型设置边界：

  # 不等式约束的乘子（非负）
  model.lambda_batt = Var(time_steps, domain=NonNegativeReals)
  # 等式约束的乘子（无边界）
  model.lambda_balance = Var(time_steps)
  

• 互补松弛条件的处理：
  • 若使用支持互补约束的求解器（如Ipopt），可直接写：

    model.complementarity = Constraint(time_steps, rule=lambda m, t: m.lambda_batt[t] * (m.batt_capacity[t] - m.SOC[t]) == 0)
    

  • 若求解器不支持，改用大M法：

    model.slack = Var(time_steps, domain=Binary)
    model.comp1 = Constraint(time_steps, rule=lambda m, t: m.lambda_batt[t] <= 1e4 * m.slack[t])
    model.comp2 = Constraint(time_steps, rule=lambda m, t: (m.batt_capacity[t] - m.SOC[t]) <= 1e4 * (1 - m.slack[t]))
    

    其中1e4为大常数，需根据问题规模调整。

五、简化示例代码（Pyomo） #

from pyomo.environ import *
import pandas as pd

# 读取时间序列数据（光伏出力、负荷、电网电价）
data = pd.read_csv('microgrid_ts_data.csv')
time_steps = data.index.tolist()

# 初始化模型
model = ConcreteModel()

# ---------------------- 上层变量 ----------------------
model.p_sell = Var(time_steps, domain=NonNegativeReals, bounds=(0.5, 1.2))  # 售电价范围
model.p_dg = Var(time_steps, domain=NonNegativeReals, bounds=(0.3, 0.8))    # DG收购价范围
model.P_grid = Var(time_steps, domain=NonNegativeReals)                     # 电网购电量

# ---------------------- 内层变量与KKT乘子 ----------------------
model.DeltaP_dr = Var(time_steps, bounds=(-0.2*data['Load'], 0))  # 需求响应仅削减负荷
model.P_batt_charge = Var(time_steps, domain=NonNegativeReals)
model.P_batt_discharge = Var(time_steps, domain=NonNegativeReals)

# KKT乘子：等式约束（功率平衡）的乘子无边界，不等式约束（电池充放电）的乘子非负
model.lambda_balance = Var(time_steps)
model.lambda_charge = Var(time_steps, domain=NonNegativeReals)
model.lambda_discharge = Var(time_steps, domain=NonNegativeReals)

# ---------------------- KKT约束（内层问题） ----------------------
# 1. 原始可行性：功率平衡
def power_balance_rule(m, t):
    return data['PV'][t] + m.P_batt_discharge[t] - m.P_batt_charge[t] - data['Load'][t] + m.DeltaP_dr[t] == 0
model.power_balance = Constraint(time_steps, rule=power_balance_rule)

# 2. 梯度条件：内层目标对各变量的偏导 = 约束梯度×乘子
# 内层目标：min (Load-DeltaP_dr)*p_sell - PV*p_dg + 0.1*(charge+discharge)
def kkt_gradient_dr(m, t):
    return -m.p_sell[t] + m.lambda_balance[t] == 0
model.kkt_dr = Constraint(time_steps, rule=kkt_gradient_dr)

def kkt_gradient_charge(m, t):
    return 0.1 - m.lambda_charge[t] + m.lambda_balance[t] == 0
model.kkt_charge = Constraint(time_steps, rule=kkt_gradient_charge)

# ---------------------- 上层目标函数 ----------------------
def total_cost_rule(m):
    return sum(
        m.P_grid[t] * data['Grid_Price'][t] + 
        m.p_dg[t] * data['PV'][t] - 
        m.p_sell[t] * (data['Load'][t] - m.DeltaP_dr[t]) 
        for t in time_steps
    )
model.total_cost = Objective(rule=total_cost_rule, sense=minimize)

# ---------------------- 求解模型 ----------------------
solver = SolverFactory('ipopt')  # 若需混合整数，改用gurobi
result = solver.solve(model, tee=True)

# ---------------------- 无优化场景对比 ----------------------
model_no_opt = ConcreteModel()
# 固定电价，不考虑需求响应
model_no_opt.p_sell_fixed = Param(time_steps, initialize=0.8)
model_no_opt.p_dg_fixed = Param(time_steps, initialize=0.5)
# 省略无优化模型的约束与目标构建，求解后对比成本与可再生能源利用率

六、CasADi自动推导KKT的替代方案 #

如果手动推导梯度过于繁琐，CasADi可自动生成内层问题的KKT约束，大幅减少编码工作量：

import casadi as ca
import numpy as np

# 时间序列参数
T = 24
PV = np.random.rand(T)*100
Load = np.random.rand(T)*150
Grid_Price = np.random.rand(T)*0.7

# 上层变量
p_sell = ca.MX.sym('p_sell', T)
p_dg = ca.MX.sym('p_dg', T)
P_grid = ca.MX.sym('P_grid', T)

# 构建内层问题的KKT约束
upper_cons = []
inner_vars_all = []
for t in range(T):
    # 内层变量
    DeltaP_dr = ca.MX.sym(f'DeltaP_dr_{t}')
    charge = ca.MX.sym(f'charge_{t}')
    discharge = ca.MX.sym(f'discharge_{t}')
    inner_vars = ca.vertcat(DeltaP_dr, charge, discharge)
    inner_vars_all.append(inner_vars)
    
    # 内层目标与约束
    obj = (Load[t] - DeltaP_dr)*p_sell[t] - PV[t]*p_dg[t] + 0.1*(charge+discharge)
    cons = ca.vertcat(
        PV[t] + discharge - charge - Load[t] + DeltaP_dr,  # 功率平衡（等式）
        -charge,  # charge ≥0 → -charge ≤0
        -discharge,  # discharge ≥0 → -discharge ≤0
        DeltaP_dr + 0.2*Load[t]  # DeltaP_dr ≥-0.2Load → DeltaP_dr+0.2Load ≥0 → 取负为≤0
    )
    
    # 自动生成KKT条件
    nlp = {'x': inner_vars, 'f': obj, 'g': cons}
    kkt = ca.nlpsol('kkt', 'ipopt', nlp, {'ipopt': {'print_level':0}})
    # 添加KKT可行性与互补松弛约束
    upper_cons.append(kkt['g'] == 0)
    upper_cons.append(kkt['lam_g'][:1] == 0)  # 等式约束乘子无符号，直接满足梯度条件
    upper_cons.append(kkt['lam_g'][1:] >= 0)  # 不等式约束乘子非负
    upper_cons.append(ca.mul(kkt['lam_g'][1:], kkt['g'][1:]) == 0)

# 上层目标
upper_obj = sum(P_grid[t]*Grid_Price[t] + p_dg[t]*PV[t] - p_sell[t]*(Load[t]-inner_vars_all[t][0]) for t in range(T))
upper_x = ca.vertcat(p_sell, p_dg, P_grid, *inner_vars_all)

# 求解上层问题
solver = ca.nlpsol('solver', 'ipopt', {'x': upper_x, 'f': upper_obj, 'g': ca.vertcat(*upper_cons)}, {'ipopt': {'print_level':5}})
sol = solver(lbx=0)  # 上层变量非负

# 提取优化结果
p_sell_opt = sol['x'][0:T]
DeltaP_dr_opt = sol['x'][3*T : 4*T]

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