储能参与现货电能量-调频市场的双层决策：从论文到可跑代码的落地路径

简介：这份资源面向电力系统与电力市场领域的研究人员、工程师及储能运营商，围绕储能作为独立市场主体参与现货电能量与调频辅助服务市场的交易决策问题，构建了上层收益最大化、下层市场联合出清的双层优化模型，并借助KKT条件将其转化为单层混合整数线性规划求解，为策略性报价与收益最大化提供可复现的技术路径。资源包共1个docx文件，约53KB，以论文复现文档形式呈现，内含完整的数学建模推导、Python代码实现及结果分析，便于读者理解并复现模型构建、求解与算例验证全过程。目前已有70人学习。读者可从中掌握双层模型的KKT转化思路、储能充放电与荷电状态约束的建模方法，以及调频市场收益占比超80%这一关键结论背后的报价策略逻辑，适合用于储能参与市场的方案设计与政策效益评估。

1. 储能参与现货电能量-调频市场的双层决策：从论文到可跑代码的落地路径

储能电站的收益结构正在发生变化。过去靠峰谷价差套利的模式，在现货市场铺开后逐渐被压缩，而调频辅助服务市场的补偿单价高、调用频次密，成了不少独立储能主体真正赚钱的板块。这篇论文复现资源盯的就是这个场景：储能同时参与现货电能量市场和调频辅助服务市场，怎么报价才能把总收益顶上去。核心方法是用双层优化建模——上层储能定报价、下层市场做联合出清，再用 KKT 条件把双层问题压成单层混合整数线性规划（MILP）求解。资源包里给了完整的 Python 代码和逐段解释，适合做电力市场优化、储能调度策略的工程师和研究人员直接上手复现。算例里调频收益占比超过 80%，这个数字本身就值得拆开看看它是怎么来的。

2. 双层模型怎么搭：上层报价、下层出清与 KKT 转化的完整链路

2.1 为什么非得用双层结构，单层 LP 差在哪

储能参与两个市场时，报价策略和市场出清结果是互相牵制的。储能报什么价，影响它能不能中标、中多少；而市场出清的价格又反过来决定储能的实际收益。这种「你决策、我响应、你的收益取决于我的响应」的结构，本质上是 Stackelberg 博弈，用单层线性规划描述不了。

常见做法是把上层写成储能收益最大化，决策变量是它在两个市场的报价bid_energy[t]和bid_freq[t]；下层写成市场出清问题，目标是系统总成本最小或社会福利最大，约束里包含储能的运行约束和市场出清规则。两层嵌套之后，问题变成带均衡约束的数学规划（MPEC），直接求解器搞不定。

KKT 条件的价值就在这里：下层是一个凸优化问题（线性目标 + 线性约束），满足强对偶条件，可以用 KKT 条件等价替换。替换后，下层的「最优性」变成一组约束塞进上层，整个问题就退化成一个单层的 MILP，Gurobi、CPLEX 这类求解器可以直接吃。

注意：KKT 转化成立的前提是下层问题是凸的。如果下层目标函数里出现非凸项（比如含 binary 变量的目标），转化就不严格了，这是很多复现翻车的第一现场。

2.2 参数初始化：储能物理参数与市场价格预测

代码的第一步是initialize_parameters()，把所有参数塞进一个字典。这一步看着简单，但参数取值直接决定结果合不合理。

def initialize_parameters(): params = { 'T': 24, # 时段数，24小时逐时 'P_max': 50, # 储能最大充放电功率(MW) 'E_max': 200, # 储能最大容量(MWh) 'eta_c': 0.95, # 充电效率 'eta_d': 0.95, # 放电效率 'SOC_min': 0.2, # 最小荷电状态 'SOC_max': 0.9, # 最大荷电状态 'SOC_initial': 0.5, # 初始荷电状态 'lambda_energy': [30 + 5*np.sin(2*np.pi*t/24) for t in range(24)], 'lambda_freq': [15 + 3*np.sin(2*np.pi*t/24 + np.pi/2) for t in range(24)], 'C_bid_energy': 0.1, # 电能量市场报价成本系数 'C_bid_freq': 0.15, # 调频市场报价成本系数 'M': 1e5 # 大M法中的大数 } return params

逻辑说明：P_max和E_max决定了储能的功率和能量边界，50MW/200MWh 对应 4 小时储能系统，是当前独立储能项目的典型配置。eta_c和eta_d取 0.95 是磷酸铁锂的常见 round-trip 效率水平。SOC_min=0.2、SOC_max=0.9是防止过充过放的保护区间。

参数说明：lambda_energy和lambda_freq用正弦函数模拟价格曲线，电能量价格均值 30 元/MWh、调频价格均值 15 元/MWh。这里要注意，实际市场里调频价格通常用容量补偿方式结算，单位可能是元/MW·h 而非元/MWh，复现时如果直接套用论文的数值，收益量级可能对不上。C_bid_energy和C_bid_freq是报价成本系数，代表储能因报价产生的间接成本（比如预测偏差惩罚），取值偏小，主要起正则化作用。

2.3 下层出清约束与储能运行约束的代码映射

build_bi_level_model()函数把上下层变量和约束一次性定义出来。变量分两组：上层是bid_energy和bid_freq，下层是p_energy、p_freq、p_charge、p_discharge、soc以及两个 binary 变量u_charge、u_discharge。

# 上层目标: 最大化储能收益 def upper_objective_rule(model): revenue_energy = sum(params['lambda_energy'][t] * model.p_energy[t] for t in model.T) revenue_freq = sum(params['lambda_freq'][t] * model.p_freq[t] for t in model.T) cost_energy = sum(params['C_bid_energy'] * model.bid_energy[t] for t in model.T) cost_freq = sum(params['C_bid_freq'] * model.bid_freq[t] for t in model.T) return revenue_energy + revenue_freq - cost_energy - cost_freq model.upper_obj = Objective(rule=upper_objective_rule, sense=maximize)

目标函数由四块组成：电能量市场收益、调频市场收益、电能量报价成本、调频报价成本。收益项用的是「市场价格 × 中标功率」，注意这里用的是lambda_energy[t]而不是bid_energy[t]，因为市场按统一出清价结算，储能是价格接受者。

储能运行约束里，最关键的是充放电互斥和 SOC 动态更新：

def charge_discharge_rule(model, t): # 充放电互斥约束 return model.u_charge[t] + model.u_discharge[t] <= 1 def soc_dynamic_rule(model, t): if t == 0: soc_prev = params['SOC_initial'] * params['E_max'] else: soc_prev = model.soc[t-1] return model.soc[t] == soc_prev + model.p_charge[t] * params['eta_c'] \ - model.p_discharge[t] / params['eta_d']

互斥约束用两个 binary 变量之和 ≤ 1 实现，保证同一时段不会同时充放电。SOC 动态方程里，充电时乘以效率eta_c，放电时除以eta_d，这个方向不能搞反——充电是「电网给储能」，效率损失后储能实际得到的能量更少；放电是「储能给电网」，要放出p_discharge的量，储能内部需要消耗p_discharge / eta_d。

功率平衡约束p_discharge[t] == p_energy[t] + p_freq[t]把储能放电功率和两个市场的中标功率绑在一起。这里有个隐含假设：储能只通过放电参与市场，充电从电网取电但不计入市场中标。实际现货市场里储能充电也可以作为负荷参与，如果要把充电成本纳入，需要额外加一项购电成本。

2.4 KKT 转化：对偶变量、互补松弛与平稳性条件

convert_to_single_level()是整份代码里最需要小心的部分。它给下层每个约束配一个对偶变量，然后补上 KKT 的四个条件。

对偶变量定义：

model.dual_energy = Var(model.T, within=NonNegativeReals) # 电能量出清约束的对偶 model.dual_freq = Var(model.T, within=NonNegativeReals) # 调频出清约束的对偶 model.dual_balance = Var(model.T, within=Reals) # 功率平衡约束的对偶 model.dual_cd = Var(model.T, within=NonNegativeReals) # 充放电互斥的对偶 model.dual_charge = Var(model.T, within=NonNegativeReals) # 充电限制的对偶 model.dual_discharge = Var(model.T, within=NonNegativeReals)# 放电限制的对偶 model.dual_soc_dynamic = Var(model.T, within=Reals) # SOC动态的对偶 model.dual_soc_min = Var(model.T, within=NonNegativeReals) # SOC下限的对偶 model.dual_soc_max = Var(model.T, within=NonNegativeReals) # SOC上限的对偶

对偶变量的符号有讲究：等式约束（功率平衡、SOC 动态）的对偶变量是自由实数Reals，不等式约束（出清规则、上下限）的对偶变量是非负的NonNegativeReals。这个对应关系搞错，KKT 条件就不成立。

互补松弛条件：

def comp_slack_energy_rule(model, t): return model.dual_energy[t] * (model.bid_energy[t] - params['lambda_energy'][t]) == 0 model.comp_slack_energy = Constraint(model.T, rule=comp_slack_energy_rule)

互补松弛的含义是：如果报价严格低于市场价格（约束松弛），对偶变量必须为零；如果对偶变量大于零，报价必须恰好等于市场价格。这个「乘积为零」的条件是非线性的，Gurobi 处理这种 bilinear 项时通常需要配合 big-M 线性化，或者直接用 Gurobi 的二次约束能力。代码里直接写了乘积等于零，如果求解器报错，需要改成 big-M 形式。

平稳性条件：

def stationary_p_energy_rule(model, t): return (-params['lambda_energy'][t] + model.dual_balance[t] + model.dual_discharge[t] * params['P_max'] * model.u_discharge[t] == 0) model.stationary_p_energy = Constraint(model.T, rule=stationary_p_energy_rule)

平稳性条件是对下层拉格朗日函数求决策变量偏导后令其为零得到的。这里对p_energy求导，得到-lambda_energy + dual_balance + dual_discharge * P_max * u_discharge = 0。注意dual_discharge * u_discharge又是一个非线性项，实际求解时需要线性化处理。

提示：如果用的是 Gurobi 9.0 以上版本，可以直接处理部分二次约束，但dual * binary这种乘积仍然需要引入辅助变量线性化。常见做法是定义w[t] = dual_discharge[t] * u_discharge[t]，然后加约束w[t] <= M * u_discharge[t]、w[t] <= dual_discharge[t]、w[t] >= dual_discharge[t] - M * (1 - u_discharge[t])。

3. 跑通代码：环境配置、求解器选择与结果解读

3.1 环境依赖与 Gurobi 配置

这份代码依赖numpy、pandas、pyomo、matplotlib和gurobi。Pyomo 是建模层，Gurobi 是求解层。安装顺序建议先装 Gurobi 并拿到 license，再装 Pyomo 和其余包。

# 创建虚拟环境 python -m venv venv source venv/bin/activate # Windows 用 venv\Scripts\activate # 安装依赖 pip install numpy pandas pyomo matplotlib # Gurobi 需要单独安装并配置 license pip install gurobipy

Gurobi 的 license 获取方式有两种：学术 license 免费，商业 license 需要付费。如果手头没有 Gurobi，可以换成 CBC 或 HiGHS，但要注意 CBC 对二次约束和 big-M 的处理能力较弱，互补松弛条件可能需要手动线性化得更彻底。

# 替换求解器为 CBC（开源） solver = SolverFactory('cbc') # 或 HiGHS solver = SolverFactory('appsi_highs')

参数说明：SolverFactory的第一个参数是求解器名称，Pyomo 支持gurobi、cbc、glpk、ipopt等。tee=True会把求解器日志打印到终端，调试时建议打开，能看到约束数量、变量数量、求解时间和 gap。

3.2 求解与结果提取：收益拆解和 SOC 轨迹

solve_and_analyze()函数负责求解和结果提取。核心输出是三个数：电能量市场收益、调频市场收益、调频收益占比。

revenue_energy = sum(params['lambda_energy'][t] * p_energy[t] for t in range(params['T'])) revenue_freq = sum(params['lambda_freq'][t] * p_freq[t] for t in range(params['T'])) total_revenue = revenue_energy + revenue_freq freq_ratio = revenue_freq / total_revenue * 100

逻辑说明：收益按「市场价格 × 中标功率」逐时段累加。调频收益占比超过 80% 的结论，来自调频市场价格虽然均值低（15 元/MWh），但储能中标功率在调频市场分配更多，且调频市场的报价成本系数更高，策略性报价的空间更大。

SOC 轨迹图能直观看出储能一天内的充放电循环。正常情况下，SOC 应该在 0.2×200=40MWh 到 0.9×200=180MWh 之间波动，且首尾 SOC 不宜相差太大——如果末端 SOC 明显低于初始值，说明模型在「透支」储能能量来套利，实际运行中不可持续。

注意：如果 SOC 曲线贴着下限跑，说明 SOC_min 约束在起作用，储能被过度放电。这时候要检查P_max是否设得太大，或者市场价格曲线是否过于陡峭，导致模型倾向于把所有能量在高价时段放完。

3.3 结果合理性校验：三个必须检查的点

跑出结果后，别急着信。三个校验点：

第一，检查对偶变量是否满足互补松弛。如果dual_energy[t] > 0但bid_energy[t] < lambda_energy[t]，说明互补松弛条件被违反，KKT 转化有问题。

第二，检查功率平衡是否逐时段成立。p_discharge[t]应该严格等于p_energy[t] + p_freq[t]，如果出现偏差，说明约束没生效。

第三，检查 SOC 动态是否闭合。把每个时段的soc[t]按动态方程手算一遍，和求解结果对比，误差应该在 1e-6 以内。

# SOC 闭合校验 for t in range(params['T']): if t == 0: soc_expected = params['SOC_initial'] * params['E_max'] \ + p_charge[t] * params['eta_c'] - p_discharge[t] / params['eta_d'] else: soc_expected = soc[t-1] + p_charge[t] * params['eta_c'] - p_discharge[t] / params['eta_d'] assert abs(soc[t] - soc_expected) < 1e-6, f"SOC mismatch at t={t}"

4. 避坑与排查：KKT 转化和求解器报错的五条血泪经验

4.1 互补松弛条件导致求解器卡死或不收敛

现象：Gurobi 跑了几分钟还在 root relaxation，gap 不降，日志里反复出现「numerical trouble」。

原因：互补松弛条件dual * (bid - lambda) == 0是 bilinear 约束，Gurobi 在处理时可能陷入数值困难。尤其是当bid和lambda量级差异大时，乘积项的系数矩阵条件数很差。

解决：把互补松弛改写成 big-M 形式，引入 binary 变量z[t]：

model.z = Var(model.T, within=Binary) def comp_slack_bigm_rule(model, t): return model.bid_energy[t] - params['lambda_energy'][t] >= -params['M'] * (1 - model.z[t]) model.comp_slack_bigm = Constraint(model.T, rule=comp_slack_bigm_rule) def dual_zero_rule(model, t): return model.dual_energy[t] <= params['M'] * model.z[t] model.dual_zero = Constraint(model.T, rule=dual_zero_rule)

这样就把非线性乘积转成了线性约束，求解器处理起来稳定得多。

4.2 平稳性条件里 dual 乘 binary 导致模型非凸

现象：模型报「Q matrix is not positive semi-definite」或「non-convex」错误。

原因：平稳性条件里dual_discharge[t] * u_discharge[t]是连续变量乘 binary 变量，属于非凸项。

解决：引入辅助变量w[t]替代乘积，加三组线性约束：

model.w = Var(model.T, within=NonNegativeReals) def w_upper1(model, t): return model.w[t] <= params['M'] * model.u_discharge[t] def w_upper2(model, t): return model.w[t] <= model.dual_discharge[t] def w_lower(model, t): return model.w[t] >= model.dual_discharge[t] - params['M'] * (1 - model.u_discharge[t])

4.3 大 M 取值不当导致数值溢出或约束失效

现象：求解结果里 binary 变量出现 0.9999 或 0.0001 这种接近边界但不精确的值，或者约束明明该生效却没生效。

原因：M=1e5对于报价和功率的量级来说偏大，容易造成数值缩放问题。Gurobi 默认的整数容差是 1e-6，M 太大时，约束的松弛量可能被淹没在容差里。

解决：根据实际变量范围收紧 M。报价上限不超过市场价格的 2 倍，功率上限不超过P_max，所以 M 取2 * max(lambda_energy + lambda_freq)或2 * P_max就够了，通常 100 到 1000 量级。

4.4 下层问题非凸导致 KKT 条件不充分

现象：模型求解成功，但结果明显不合理，比如储能收益为负，或者中标功率全为零。

原因：如果下层问题里包含了 binary 变量（比如机组组合），下层就不是凸的，KKT 条件只是必要条件而非充分条件。这时候用 KKT 转化得到的单层模型，解可能不是原双层问题的均衡解。

解决：检查下层是否含 binary 变量。如果含，要么把下层松弛成 LP（忽略整数约束），要么改用其他方法（如对角化、迭代最佳响应）。这份代码的下层出清规则是线性的，没有 binary，所以 KKT 转化成立。

4.5 市场价格预测曲线过于理想导致结论不可信

现象：调频收益占比稳定在 80% 以上，但换个价格曲线就完全不一样。

原因：代码里用正弦函数模拟价格，电能量和调频价格完全反相（一个 sin，一个 sin 加 π/2），这种理想化曲线在实际市场里不存在。实际调频价格和电能量价格的相关性更复杂，可能同向也可能反向，取决于系统调频需求和新能源出力。

解决：用真实市场历史价格数据替换正弦曲线。如果拿不到数据，至少构造几条不同相关性的价格曲线做敏感性分析，看看调频收益占比在什么范围内波动。

5. 从复现到落地：用真实价格数据替换正弦曲线并做敏感性分析

论文里的正弦价格曲线是为了展示模型机制，但真要用来指导报价，必须换成真实数据。我一般会从两个渠道拿价格数据：一是电力交易中心公开发布的日前和实时出清价格，二是调频辅助服务市场的历史出清结果。拿到数据后，按 24 时段或 96 时段对齐，替换lambda_energy和lambda_freq。

import pandas as pd def load_real_prices(energy_csv, freq_csv): """从 CSV 加载真实市场价格，替换正弦模拟曲线""" df_energy = pd.read_csv(energy_csv) df_freq = pd.read_csv(freq_csv) # 假设 CSV 有 'hour' 和 'price' 两列 lambda_energy = df_energy.sort_values('hour')['price'].tolist() lambda_freq = df_freq.sort_values('hour')['price'].tolist() assert len(lambda_energy) == 24, "电能量价格时段数不对" assert len(lambda_freq) == 24, "调频价格时段数不对" return lambda_energy, lambda_freq

替换后，重点观察三个指标的变化：调频收益占比、SOC 日均循环次数、总收益对报价系数的敏感度。我习惯做一组敏感性分析，把C_bid_energy和C_bid_freq各取 0.05、0.1、0.15、0.2 四档，交叉跑 16 组，看收益排名是否稳定。

这张表是我跑过的一组典型结果。规律是：报价成本系数越高，储能越倾向于把容量往调频市场倾斜，因为调频市场的报价成本相对收益更「划算」。但这个结论依赖价格曲线形态，换成调频价格波动更大的数据，结果可能反过来。

还有一个容易忽略的点：SOC 末端约束。代码里没有强制soc[T-1] == soc[0]，所以模型可能把初始 SOC 的能量「卖光」来套利。实际运行中，储能每天要留够能量应对次日早高峰，我一般会加一条末端 SOC 约束：

def soc_terminal_rule(model): return model.soc[params['T']-1] >= params['SOC_initial'] * params['E_max'] * 0.8 model.soc_terminal = Constraint(rule=soc_terminal_rule)

这条约束加上后，总收益会降一些，但 SOC 轨迹更可持续。从那以后我每次跑储能优化模型，都强制走一遍「末端 SOC 校验 + 互补松弛校验 + 功率平衡校验」这三步，少一步都不敢把结果往报告里写。希望帮到你。

本文还有配套的精品资源，点击获取