储能参与现货-调频市场的双层交易决策与Matlab实现

储能参与现货电能量-调频辅助服务市场的双层交易决策研究(Matlab代码实现)最近在复盘储能参与电力市场的项目发现不少同学一上来就盯住“现货价差套利”不放却忽略了调频辅助服务这条同样重要的收益路径。现在很多省份已经允许储能同时参与现货电能量市场和调频辅助服务市场两个市场的出清逻辑不同、结算规则不同储能要在两个市场之间分配容量、申报价格本质上是一个典型的双层决策问题上层是储能运营商的收益最大化下层是电力市场的出清优化。用Matlab做双层建模与求解是我个人比较推荐的一条技术路线Yalmip Gurobi组合基本能覆盖大部分研究场景。这篇文章我会把整个思路、模型设计、代码实现框架和踩坑记录都整理出来适合正在做储能市场决策、调频辅助服务、双层优化相关课题的研究生和工程师参考拿到思路就能自己动手改。1. 储能参与现货电能量-调频辅助服务市场的思维起点1.1 为什么储能要把两个市场放在一起决策储能参与单一现货市场时收益模型其实比较简单低电价时段充电高电价时段放电赚的是峰谷价差。但实际项目里储能设备通常具有快速响应的优势调频辅助服务市场对这类资源的奖励往往比单纯峰谷套利更可观。问题也随之而来储能容量有限如果一部分容量用于调频另一部分用于电能量套利两者之间就产生了竞争关系。更麻烦的是两个市场不是独立运行的。储能申报电能量市场的充放电曲线会影响现货出清价格储能申报的调频容量和调频价格又会影响调频辅助服务市场的出清结果。反过来两个市场的出清价格会决定储能的最终收益。这个闭环用传统单层优化描述不了必须用双层模型储能作为领导者Leader先做容量分配与报价决策市场出清作为跟随者Follower在给定报价下优化出清结果。我见过不少论文把这两个市场拆分优化先算现货最优充放电再算调频最优申报最后简单相加。这种做法在理论上不严谨因为调频备用容量的机会成本来自电能量市场而电能量市场的充放电计划又受到调频约束的影响两个问题天然耦合。把一个耦合系统拆成两个独立的子问题计算出来的“最优策略”基本不是全局最优。1.2 现货电能量市场和调频辅助服务市场的运行逻辑现货电能量市场目前普遍以日前市场加实时市场为主储能主要参与日前市场按节点边际电价LMP结算。市场出清的核心任务是安全约束经济调度也就是在满足负荷、网络安全约束的前提下让发电成本最小化。储能作为可充电可放电的灵活资源在出清模型中通过充放电变量参与功率平衡。调频辅助服务市场的机制各省级有差异一般有两种一种按调频容量付费一种按调频里程付费更多的地区是“容量里程”组合。储能申报调频容量后在调度周期内需要实时跟随AGC指令所以储能的可放电功率会被调频容量占用这部分功率不能再参与电能量市场的放电报价。在双层模型里两个市场的耦合约束主要体现为储能总功率约束充放电功率与调频备用容量不能超过额定功率储能能量SOC约束调频动作会消耗能量所以调频容量会影响SOC变化竞价策略耦合储能申报电能量价格和调频价格会影响市场出清价格进而影响收益。1.3 双层交易决策到底在优化什么用一句话概括上层决策储能“怎么用”——每时段申报多少充放电功率、申报什么价格、预留多少调频容量下层市场出清决定储能“被怎么分配”——在给定申报值的情况下储能实际上中标多少、获得什么结算价格。上层目标函数是储能的总收益最大化包括电能量收益和调频收益同时扣除电池循环损耗成本或者折旧成本。下层目标函数是市场出清目标通常是最小化系统总购电成本或者最大化社会福利。上下层之间的关联通过储能申报的报价曲线和市场中其他机组报价共同参与出清来实现。这类模型最常见的形式是MPEC带均衡约束的数学规划而下层线性规划问题可以通过KKT条件或者强对偶理论等价转换到上层约束中再用大M法或SOS1处理互补松弛条件最终转化为一个混合整数线性规划MILPMatlab配合YalmipGurobi就能求解。这也是我接下来要重点展开的实现路径。2. 下层市场出清模型与上层决策目标的建模细节2.1 下层市场出清的数学表达做储能参与市场的研究我建议先把下层市场出清模型写熟练。以一个简化但保留核心机制的现货出清模型为例目标函数是最小化系统总成本目标函数 系统总成本 常规机组发电成本 储能放电成本或报价成本 调频容量购买成本 调频里程成本常规机组发电成本可以用二次函数或者分段线性函数近似。为了把下层模型保持为线性规划通常用分段线性逼近或者直接用线性报价曲线。储能放电成本可以设为一个很小的正数代表电池放电的边际损耗。系统约束至少包括功率平衡约束每个时段的负荷等于常规机组出力加储能放电功率减去储能充电功率机组出力上下限约束储能充放电功率上限约束其中放电上限要扣除调频容量调频容量平衡约束系统调频需求等于各机组和储能中标调频容量之和调频容量上限约束线路潮流约束如果采用直流潮流模型。下层模型用线性规划表示时变量包括各机组出力、储能充放电功率、储能调频容量、调频里程中标量等。所有约束都是线性的KKT条件可以严格写出。2.2 上层储能运营商的收益函数上层的核心是收益表达式我实际建模时用的收益公式如下总收益 电能量市场收益 调频容量收益 调频里程收益 - 电池充放电损耗成本其中电能量市场收益 各时段放电功率 × 现货节点电价 - 各时段充电功率 × 现货节点电价调频容量收益 各时段调频中标容量 × 调频容量出清价格调频里程收益 各时段调频里程中标量 × 调频里程出清价格。注意现货电价、调频容量价格、调频里程价格在单层优化里都是外生参数但在双层模型里它们是下层市场出清的对偶变量是内生变量。这就是双层模型“难”的本质上层决策变量出现在下层约束中下层对偶变量又出现在上层目标函数里。2.3 两个市场之间怎么通过容量约束耦合容量耦合是最容易写错的地方我重点说一下。储能系统额定功率为P_max假设储能参与电能量放电和调频容量申报那么在任意时段t必须满足P_dch(t) R(t) P_max其中P_dch(t)是放电功率R(t)是调频容量。充电功率单独由P_ch(t)表示P_ch(t) P_max同时充电和放电不能同时进行 P_ch(t) P_dch(t) P_max我最初建模时图省事只写了P_dch(t) R(t) P_max和P_ch(t) P_max结果出现充电和放电同时为正的荒谬结果。后来加了防同时充放电的约束才正常。在MILP中这个约束需要引入二进制变量但在很多学术简化模型里也可以用P_ch(t) P_dch(t) P_max来近似虽然严格说这个约束在物理意义上不是完全等价于“不能同时充放电”但在目标函数激励下同时充放电没有经济意义所以可以接受。SOC耦合约束也很关键SOC(t1) SOC(t) (P_ch(t) * eta_ch - P_dch(t) / eta_dch - beta * R(t)) * dt这里的beta表示调频容量对SOC的等效消耗系数。严格说调频动作是双向的上下调节会部分抵消但考虑到调节过程存在能量损耗简化模型中用一个线性系数折算到SOC变化量上是论文里很常见的近似。2.4 为什么必须用双层而不是加权单层我经常被问到一个问题既然下层出清是线性规划能不能直接把两个市场的目标函数加权成一个单层优化理论上如果没有上下层之间的策略性交互单层确实可以做。但储能的市场力体现在它可以调整申报策略来影响出清结果这时候单层优化就无法表达“储能预测市场出清结果并响应”的逻辑。举个例子储能知道夜间申报低价充电会导致自己小额度中标白天申报略低于边际机组价格可以保证放电优先中标这种策略性报价行为天生就是Stackelberg博弈储能是领导者市场是跟随者。单层优化等价于把市场价格当常数根本反映不出储能报价对市场出清价格的反馈影响。这就是双层模型的价值所在。3. 双层模型的KKT转换与求解思路3.1 KKT条件的工程化处理下层是线性规划满足Slater条件所以最优解等价于其KKT条件。把下层KKT条件加入上层模型得到MPEC然后用大M法将互补松弛条件线性化就可以交给Gurobi求解MILP。具体来说下层每个不等式约束都会对应一个对偶变量和一个互补条件。互补条件的标准形式是0 对偶变量 ⊥ (约束右端项 - 约束左侧表达式) 0这个“⊥”符号表示两者至少有一个为零。要转化为线性条件引入二进制变量和大M常数对偶变量 M * z (约束右端项 - 约束左侧表达式) M * (1 - z)M的取值很关键。取太小会把最优解排除在外取太大会导致数值不稳定。我的经验是先求解一次不含互补条件的松弛问题观察对偶变量和约束松弛量的取值范围再据此设置M一般取10^3到10^4之间比较稳妥。3.2 Yalmip建模与Gurobi求解的代码框架我在Matlab里用Yalmip建模根据是否引入了二进制变量分别调用gurobi或求解器处理MILP。下面的代码框架是我调试过的版本逻辑完整可以直接在此基础上改算例。%% 参数定义 % 系统参数 T 24; % 调度时段数 load_demand [40;35;...]; % 各时段负荷24x1 gen_max [80;50;...]; % 常规机组容量上限 gen_cost [25;30;...]; % 常规机组边际成本 % 储能参数 Pmax 20; % 额定功率 MW Emax 100; % 额定容量 MWh eta_ch 0.95; % 充电效率 eta_dch 0.95; % 放电效率 SOC_init 0.2 * Emax; % 初始SOC SOC_min 0.1 * Emax; SOC_max 0.9 * Emax; %% 变量定义 % 上层决策变量 P_ch sdpvar(T,1); % 充电功率 P_dch sdpvar(T,1); % 放电功率 R sdpvar(T,1); % 调频容量申报 % 下层变量 P_gen sdpvar(T,size(gen_cost,1)); % 机组出力 R_freq sdpvar(T,1); % 系统调频中标变量简化为单资源 % 对偶变量由Yalmip自动生成通常在KKT转换时手动写出严格的双层MPEC求解比较繁琐我建议分两步第一步写清楚下层的拉格朗日函数第二步用KKT替换下层约束。Yalmip里有一个功能是可以用kkt(Constraints,Objective)直接得到KKT系统但实际用下来发现变量多的时候容易出错我反而更喜欢手动写出KKT条件逻辑更清楚也方便调试。手动写KKT的一个代码片段如下%% 下层KKT条件 % 对下层变量P_gen求导 Stationarity [...]; % 导数表达式 % 互补条件线性化 M 10000; % 大M常数 z binvar(size(ConstraintMatrix,1),1); % 二进制变量 Complementarity []; for i 1:num_constraints Complementarity [Complementarity, dual_var(i) 0]; Complementarity [Complementarity, constraint_residual(i) 0]; Complementarity [Complementarity, dual_var(i) M * z(i)]; Complementarity [Complementarity, constraint_residual(i) M * (1 - z(i))]; end3.3 强对偶法的替代方案KKT大M法是经典做法但有一个明显缺点引入了大量二进制变量求解速度会随着节点数和时段数增加而迅速恶化。我在做24时段、96时段的算例时明显感到求解时间从几秒上升到几十秒甚至几分钟。一个替代方案是强对偶法。利用强对偶定理可以用下层原始目标值等于对偶目标值来替代互补松弛条件。强对偶等式不引入二进制变量模型规模小很多得到的是一个非线性优化问题需要精心处理双线性项对偶变量乘上层变量。如果上层价格变量只出现在对偶目标中可以采用线性化的技巧处理。实际工程中我还用过一种迭代逼近方法先固定上层申报求解下层出清根据下层出清价格回代求解上层优化反复迭代。这个方法不保证收敛到全局最优但算得快适合大规模系统初步筛选方案。做学术研究建议还是用KKT强对偶的严格方法工程预研可以考虑迭代法。4. 算例设计与Matlab实现过程4.1 算例系统的构建我用的算例是一个简化火电-储能测试系统包含2台常规火电机组、1台储能、1个聚合负荷节点。这种简化足够验证储能参与两个市场后的行为变化又不至于让下层出清的KKT条件过于复杂。系统参数设置如下参数数值说明火电1额定功率80 MW边际成本 25元/MWh火电2额定功率50 MW边际成本 30元/MWh储能额定功率20 MW充放电共用储能额定容量100 MWhSOC范围 10%-90%充/放电效率95%线性简化处理系统调频需求5 MW每个时段固定调频里程系数0.2容量到里程的折算负荷曲线设置为典型的两峰一谷形状这样储能可以在低谷充电、两个高峰放电同时在高峰时段预留部分容量参与调频能明显看出容量分配的变化。4.2 仿真结果的关键观察点调好代码后我最先关注三个输出各时段储能充放电曲线各时段调频容量申报量两个市场分别贡献的收益。实际上储能参与调频后电能量市场的充放电量会明显减少。以我上面参数的系统为例如果不参与调频储能会在两个高峰时段各放电约15MWh参与调频后因为每个时段都要预留5MW的调频容量放电能力被压缩高峰放电量降到了10MWh左右。但总收益反而是增加的因为调频容量收益补偿了电能量套利减少的损失总收益提升约15%到20%。这个结论符合直觉储能响应速度快的天生优势在调频市场上比单靠峰谷价差更有价值。我还对比了不同电价场景。如果峰谷价差超过0.6元/kWh电能量套利收益占比会上升如果价差小于0.4元/kWh储能更倾向多分一些容量给调频市场。这说明双层模型的自适应能力是有价值的固定比例分配容量的经验方法是不可靠的。4.3 代码运行效率优化算例规模不大时双层模型很快但想扩展到IEEE 30节点或者96时段求解效率就成了硬伤。我常用的优化手段一是减少二进制变量。互补条件里的二进制变量数量 下层不等式约束数。如果某些不等式约束在最优解中不可能取等号可以在建模时就剔除避免白白增加变量。二是对大M做灵敏度分析。先跑一个初步结果观察各互补条件实际的对偶变量最大量级再调整大M值。大M从10^6降到10^3后Gurobi的求解速度能提升30%以上。三是对称性消除。如果下层出清模型中存在多个同质机组可以用聚合机组替代减少约束数量。我实际测试过不加这些优化24时段模型求解耗时约45秒做了变量削减和M调整后降到8秒左右。对需要批量跑蒙特卡洛或灵敏度分析的研究任务这个提升非常关键。5. 常见问题与排查技巧实录5.1 求解器报“无法求解非线性模型”这个是我最常遇到的问题。用Yalmip建双层模型时如果直接把下层对偶变量和上层变量相乘模型就变成了非线性Gurobi无法直接处理。报错通常提示“Nonlinear constraints detected”。解决办法是先看目标函数里是否有双线性项。上层收益中的价格对偶变量乘以储能功率上层变量就是双线性项。如果价格是市场出清价格这个乘法必然存在。这时要么用强对偶方法消去双线性项要么在KKT框架下把价格变量对应的互补约束线性化后再用目标函数中的“收益”替换为线性表达式。一个实用技巧如果坚持要用Gurobi求解MILP那么在目标函数里不要出现对偶变量与上层变量的乘积考虑用下层最优目标的对偶等价表达式替代收益这样可以把双线性项消掉。5.2 模型可行但结果不符合物理直觉比如储能莫名其妙地在高电价时段充电或者在低电价时段放电。我排查的顺序是检查SOC约束是否写反了充放电方向检查P_dch和P_ch是否用了同一个功率上限导致同时充放电检查调频容量R是否错误地减了充电功率上限检查对偶变量符号是否与约束方向匹配。其中对偶变量符号问题最隐蔽。下层约束是“”和“”两种方向对应不同符号的对偶变量一旦写反求出“价格”可能为负储能就会反向套利。我的建议是先在无储能报价参与的情况下用Matpower或者通用优化求解器跑一遍下层出清检查节点电价是否为正且随负荷升高而升高确认下层模型正确后再加入上层模型。5.3 Matlab代码报错排查代码层最常见的错误集中在Yalmip变量维度不匹配、sdpvar拼接成大矩阵时索引越界、以及kkt函数对包含二进制变量的约束不适用。我在写代码时习惯每个约束组合都用size()打印维度确认左端和右端维度一致再append到总约束集。另外遇到“Index exceeds array bounds”时第一时间不是看循环变量而是看Yalmip的变量矩阵维度因为很多表达式展开后维度超出了预期。关于Yalmip本身的安装不少同学卡在setup没反应。这个组件下载依赖网络可以手动下载对应版本放入工具箱路径再用addpath(genpath(...))加载。Matlab版本太旧也容易出兼容问题建议至少在R2020a以上跑YalmipGurobi我自己在R2022b和R2023a环境里都测试过没有遇到兼容性问题。5.4 大M取值不当导致的数值问题大M的取值不能拍脑袋。M太小模型切掉了真实的最优解结果总是与原问题不一致M太大Gurobi求解时会出现numerical trouble甚至报“infeasible or unbounded”。我的排查方法先写出下层约束的残差范围比如功率平衡约束的残差可能是机组出力累计值通常不超过500。对应的大M取5000就足够调频容量约束残差不超过50大M取500就行。不要全局统一一个M最好每个互补条件单独设置一个差别不大的M能显著改善数值稳定性。另外使用SOS1约束替代大M也是常见做法Gurobi对SOS1支持良好且不需要手动调M只是建模时要Yalmip中手动指示变量组成SOS1集合。6. 项目扩展与经验总结储能参与两个市场的双层交易决策我做完一遍的最大感受是模型结构不难难在把市场机制转换成数学表达式时不要失真。国内各省电力市场的规则还在快速迭代比如有的省份现货市场采用节点边际电价出清调频市场又有性能指标和里程补偿系数这些都直接影响储能参与策略。做研究时可以先抓主要矛盾用简化模型验证方法再逐步叠加地区规则细节。如果后续想扩展我有几个建议方向引入不确定性的两阶段或分布鲁棒优化考虑新能源出力波动对现货价格的影响多储能主体博弈上层从单领导者变成多领导者模型从MPEC变成EPEC求解复杂度明显上升在约束中加入调频里程的动态响应过程把现在简化的线性折算改成更精细的动态模型结合储能电池寿命模型把循环老化成本内生到上层收益函数中这样储能频繁充放电的行为会被自动抑制结果也更有工程参考价值。我个人的习惯是把Matlab代码模块化数据参数文件、下层模型文件、上层模型文件、求解主文件分开写算例调整只需要改参数文件。遇到模型改版比如从确定性改成鲁棒优化只需要在约束层加场景变量不需要大规模改动主程序。做储能市场决策方向的课题建议一开始就按这个结构来否则反复调试会非常消耗时间。希望这篇拆解能帮你绕开我踩过的坑直接跑通自己的储能双层交易决策模型。