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

简介这份资源面向电力系统与电力市场领域的研究人员、工程师及储能运营商围绕储能作为独立市场主体参与现货电能量与调频辅助服务市场的交易决策问题构建了上层收益最大化、下层市场联合出清的双层优化模型并借助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 条件等价替换。替换后下层的「最优性」变成一组约束塞进上层整个问题就退化成一个单层的 MILPGurobi、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_min0.2、SOC_max0.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(ruleupper_objective_rule, sensemaximize)目标函数由四块组成电能量市场收益、调频市场收益、电能量报价成本、调频报价成本。收益项用的是「市场价格 × 中标功率」注意这里用的是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, withinNonNegativeReals) # 电能量出清约束的对偶 model.dual_freq Var(model.T, withinNonNegativeReals) # 调频出清约束的对偶 model.dual_balance Var(model.T, withinReals) # 功率平衡约束的对偶 model.dual_cd Var(model.T, withinNonNegativeReals) # 充放电互斥的对偶 model.dual_charge Var(model.T, withinNonNegativeReals) # 充电限制的对偶 model.dual_discharge Var(model.T, withinNonNegativeReals)# 放电限制的对偶 model.dual_soc_dynamic Var(model.T, withinReals) # SOC动态的对偶 model.dual_soc_min Var(model.T, withinNonNegativeReals) # SOC下限的对偶 model.dual_soc_max Var(model.T, withinNonNegativeReals) # 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, rulecomp_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, rulestationary_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 gurobipyGurobi 的 license 获取方式有两种学术 license 免费商业 license 需要付费。如果手头没有 Gurobi可以换成 CBC 或 HiGHS但要注意 CBC 对二次约束和 big-M 的处理能力较弱互补松弛条件可能需要手动线性化得更彻底。# 替换求解器为 CBC开源 solver SolverFactory(cbc) # 或 HiGHS solver SolverFactory(appsi_highs)参数说明SolverFactory的第一个参数是求解器名称Pyomo 支持gurobi、cbc、glpk、ipopt等。teeTrue会把求解器日志打印到终端调试时建议打开能看到约束数量、变量数量、求解时间和 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×20040MWh 到 0.9×200180MWh 之间波动且首尾 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, fSOC mismatch at t{t}4. 避坑与排查KKT 转化和求解器报错的五条血泪经验4.1 互补松弛条件导致求解器卡死或不收敛现象Gurobi 跑了几分钟还在 root relaxationgap 不降日志里反复出现「numerical trouble」。原因互补松弛条件dual * (bid - lambda) 0是 bilinear 约束Gurobi 在处理时可能陷入数值困难。尤其是当bid和lambda量级差异大时乘积项的系数矩阵条件数很差。解决把互补松弛改写成 big-M 形式引入 binary 变量z[t]model.z Var(model.T, withinBinary) 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, rulecomp_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, ruledual_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, withinNonNegativeReals) 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 这种接近边界但不精确的值或者约束明明该生效却没生效。原因M1e5对于报价和功率的量级来说偏大容易造成数值缩放问题。Gurobi 默认的整数容差是 1e-6M 太大时约束的松弛量可能被淹没在容差里。解决根据实际变量范围收紧 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 循环次数C_e0.05, C_f0.05最高约 75%1.8C_e0.10, C_f0.15中等约 82%1.5C_e0.20, C_f0.20最低约 88%1.2这张表是我跑过的一组典型结果。规律是报价成本系数越高储能越倾向于把容量往调频市场倾斜因为调频市场的报价成本相对收益更「划算」。但这个结论依赖价格曲线形态换成调频价格波动更大的数据结果可能反过来。还有一个容易忽略的点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(rulesoc_terminal_rule)这条约束加上后总收益会降一些但 SOC 轨迹更可持续。从那以后我每次跑储能优化模型都强制走一遍「末端 SOC 校验 互补松弛校验 功率平衡校验」这三步少一步都不敢把结果往报告里写。希望帮到你。本文还有配套的精品资源点击获取