基于Matlab的旋转备用与能量联合出清模型:从线性规划到影子价格解析
手里如果有一套能跑通的旋转备用出清代码很多问题会变得特别具体。比如“旋转备用价格到底是怎么算出来的”“为什么某台机组报了低价却没有中标备用”“能量价格和备用价格之间到底有什么关系”这些在你真正把模型写出来、跑一遍、再改几个参数之后答案会非常清晰。这篇文章我把自己做“主辅助服务市场出清模型研究【旋转备用】”时的完整思路整理出来了用一个简化的单节点系统、24小时场景在Matlab里实现能量与旋转备用的联合出清。代码基于线性规划核心求解器是Optimization Toolbox里的linprog模型兼顾了教学可读性和工程扩展性。适合电力市场刚入门的研究生、刚接触出清模型的工程师以及想快速验证一些市场机制想法但不想从零推公式的同行。1. 为什么单独拆出“旋转备用”来建模1.1 旋转备用在辅助服务市场里是什么角色旋转备用Spinning Reserve指的是已经并网、处于同步运行状态的机组能够在很短时间通常是10到15分钟内增加出力用来应对负荷预测偏差、新能源功率波动或者机组非计划跳闸。它和冷备用最大的区别在于“热备用”机组已经在转只是留了一部分容量不发电随时可以顶上去。正因为这个特点旋转备用的成本不只是机组报价单上的那笔容量费用还包含了“本来可以发电赚钱但没发”的机会成本。所以主辅助服务市场不能把能量市场单独算完再慢慢分配备用而要把“买电”和“买备用容量”放到同一个优化问题里联合出清。这就是“主辅助服务市场出清”的核心逻辑也是这个模型和单纯经济调度最大的区别。1.2 联合出清的经济学直觉请一个人干两件事理解联合出清可以拿招人干活来类比。你有一批任务同时还需要有人随时准备处理突发情况。如果有个人干活便宜但不愿意承担额外职责另一个人干活稍贵但能同时兼职救火那你不能只看“干活价格”选前者因为还需要为“救火能力”单独付费。放到电力系统里机组A能量报价低但容量裕度不够报备用容量也少机组B能量报价高一些但容量大、能同时提供较多备用。如果只按能量报价排序调度A可能满发但备用需求会压到B甚至更贵的机组上导致备用容量成本很高。联合出清的优化目标是最小化“能量成本备用容量成本”的总和最终给出的调度方案往往不是单纯最便宜的机组满发而是在能量和备用之间做整体取舍。1.3 模型边界设定先跑通主干再谈复杂化这个版本我刻意做了一些简化避免一上来就被各种细节淹没单节点系统忽略网络潮流和阻塞所有机组面向同一系统负荷和同一备用需求。所有机组已经处于开机状态暂不考虑启停状态变量。只考虑旋转备用不涉及调频、备用分组、爬坡分档等细分。备用容量按容量付费即按“预留了多少MW”计费不模拟被调用后的电量结算。这些假设对于理解市场出清的骨架完全够用。节点价格、机组组合、多类型备用都是在同一套骨架上扩展出来的先把线性规划版本跑透后面加内容会轻松很多。2. 数学模型从物理约束到线性规划2.1 决策变量与目标函数先定义问题里的决策变量。系统有 ng 台机组、T 个时段P(i, t)机组 i 在时段 t 的出力单位 MW。R(i, t)机组 i 在时段 t 预留的旋转备用容量单位 MW。目标函数是让整个调度周期内的总成本最小min sum_t sum_i [ Cp(i) * P(i,t) Cr(i) * R(i,t) ]Cp(i) 是机组 i 的能量报价元/MWhCr(i) 是备用容量报价元/MW。这里 R(i,t) 是“预留容量”而不是“被调用的电量”所以 Cr 的单位是元/MW而不是元/MWh。这是初学时特别容易混淆的一个点备用市场买的不是电量是随时准备发电的权利。2.2 约束条件每一行都要有物理意义模型的约束条件可以分为四类每一条都能对应到实际物理或市场规则功率平衡约束是硬约束任何时刻所有机组出力之和必须等于系统负荷sum_i P(i,t) Load(t), for all t备用需求约束保证系统有足够的旋转备用应对扰动。考虑到负荷的10%左右通常是比较常见的设定但具体比例取决于系统规模和最大单机容量sum_i R(i,t) ReserveNeed(t), for all t机组出力上下限约束0 P(i,t) Pmax(i) 0 R(i,t) Rmax(i)容量耦合约束是这一模型最容易忘也最重要的一条。机组在时刻 t 的出力与备用容量之和不能超过其最大技术出力P(i,t) R(i,t) Pmax(i)这条约束的意思是一台机组不能一边满发一边宣称自己有大量备用容量因为备用是“还能再多发的那部分能力”。我在调试第一个版本时就是因为漏了这条约束导致结果里出现机组出力200MW、备用还有150MW的荒唐结果总备用需求看着满足了实际物理上根本不可能。爬坡约束考虑机组出力调整速度限制。实际机组每时段能够增加的出力有限P(i,t) - P(i,t-1) RampUp(i) P(i,t-1) - P(i,t) RampDown(i)在简化版本里可以设成对称的即 RampUp RampDown。2.3 出清价格从哪来线性规划的对偶变量模型求解后有两个层面的产出决策变量给出调度计划而对偶变量给出价格信号。在 linprog 的返回结果里lambda.eqlin 对应的就是等式约束的影子价格也就是功率平衡约束和备用需求约束的边际成本。功率平衡约束的影子价格可以理解为该时段系统增加1MW负荷时总成本的变化量它最终会体现在市场出清价格上。备用需求约束的影子价格则对应系统增加1MW备用需求时成本的变化量。用大白话说如果某时段能量出清价格是320元/MWh意味着再增加1MW负荷最优解的总成本会增加约320元。这个价格通常由边际机组的能量报价决定。备用价格同理由备用约束的边际机组报价决定。这里有个很有意思的现象备用价格往往不等于某台机组备用报价单上的数字而是能量和备用联合优化的结果。因为当一台机组同时参与能量和备用竞争时它的机会成本会传导到两个市场里。这个问题在结果分析部分再展开。3. Matlab代码实现从零搭一个可运行的出清模型3.1 场景参数与数据准备先构造一个24小时的测试场景。我选了一个典型的夏季度冬型负荷曲线夜间低谷约510MW午高峰和晚高峰分别冲到800MW以上clear; clc; close all; %% 1. 场景数据 T 24; % 时段数每小时一个点 load_p [560 540 520 510 520 540 600 680 760 ... 800 820 810 780 760 760 770 800 820 ... 800 760 720 680 640 600]; % 24小时负荷单位MW % 旋转备用需求取负荷的10%并按10MW取整 reserve_need floor(load_p * 0.1 / 10) * 10; % 机组参数Pmax(MW), 能量报价(元/MWh), 备用容量报价(元/MW) gen [ 200 280 15 180 310 18 160 350 22 120 400 28 100 450 35 ]; ng size(gen, 1); % 机组数量 Pmax gen(:, 1); % 最大出力 Cp gen(:, 2); % 能量报价 Cr gen(:, 3); % 备用容量报价 % 备用容量上限本模型设为机组最大出力的20% Rmax 0.2 * Pmax; % 爬坡能力简化设为最大出力的30%/小时 Ramp 0.3 * Pmax;备用需求取负荷的10%是常见做法更精细的市场会考虑最大单机容量和新能源预测误差工程上还会做备用容量需求曲线。机组报价我故意分成五个梯队前两台大容量机组报价较低相当于大容量火电后面是中等容量机组最后一台小机组报价很高相当于调峰压力大的时候才被叫到的机组。3.2 目标函数与变量上下界装配变量排列顺序是整个代码最容易出错的点。我采用先P后R的排布x [P(1,1), P(2,1), ..., P(ng,1), P(1,2), ..., P(ng,T), ... R(1,1), R(2,1), ..., R(ng,1), R(1,2), ..., R(ng,T)]索引关系是机组 i 在时段 t 的出力索引(t-1) * ng i机组 i 在时段 t 的备用索引ng * T (t-1) * ng i目标函数系数向量可以按这个顺序组装%% 2. 目标函数系数 f n 2 * ng * T; % 总变量数 f zeros(n, 1); for t 1:T for i 1:ng f((t-1)*ng i) Cp(i); % 能量成本系数 f(ng*T (t-1)*ng i) Cr(i); % 备用容量成本系数 end end虽然用 repmat 和 reshape 可以写得非常紧凑但对不熟悉矩阵展开顺序的人来说容易翻车。循环装配虽然看着笨但每一行都能对着公式核调试成本低很多。等代码跑通了再优化成向量化写法也不迟。变量上下界用 lb 和 ub 表示lb zeros(n, 1); ub inf(n, 1); for t 1:T for i 1:ng ub((t-1)*ng i) Pmax(i); % 出力上限 ub(ng*T (t-1)*ng i) Rmax(i); % 备用容量上限 end end这里出力上限直接取Pmax备用上限取Rmax但 P(i,t)R(i,t)Pmax 这条耦合约束会在 A 矩阵里单独加不能只靠 ub 实现。3.3 约束矩阵装配等式与不等式分开搭等式约束包括功率平衡和备用需求约束一共 2T 行%% 3. 等式约束 Aeq zeros(2*T, n); beq zeros(2*T, 1); for t 1:T for i 1:ng % 功率平衡sum P(i,t) load_p(t) Aeq(t, (t-1)*ng i) 1; % 备用需求sum R(i,t) reserve_need(t) Aeq(Tt, ng*T (t-1)*ng i) 1; end beq(t) load_p(t); beq(Tt) reserve_need(t); end这里备用需求我用的是等式约束。实际工程中更稳妥的是不等式约束加失备用惩罚防止某些极端场景下无可解后面第5节会专门讨论。不等式约束包含三部分%% 4. 不等式约束 % 第一部分:PR Pmax共 ng*T 行 % 第二部分:向上爬坡 P(t)-P(t-1) Ramp共 ng*(T-1) 行 % 第三部分:向下爬坡 P(t-1)-P(t) Ramp共 ng*(T-1) 行 n_ineq ng*T 2*ng*(T-1); A zeros(n_ineq, n); b zeros(n_ineq, 1); row 0; for t 1:T for i 1:ng row row 1; A(row, (t-1)*ng i) 1; A(row, ng*T (t-1)*ng i) 1; b(row) Pmax(i); end end for t 2:T for i 1:ng row row 1; A(row, (t-1)*ng i) 1; % P(i,t) A(row, (t-2)*ng i) -1; % -P(i,t-1) b(row) Ramp(i); % 向上爬坡 end end for t 2:T for i 1:ng row row 1; A(row, (t-2)*ng i) 1; % P(i,t-1) A(row, (t-1)*ng i) -1; % -P(i,t) b(row) Ramp(i); % 向下爬坡 end end这一段是整个代码里最容易出错的部分我调试时曾经把行号算错导致某些约束被覆盖。建议每写完一段就打印 row 的数量核对第一部分的行数应该是 ngT120第二部分是 ng(T-1)115第三部分同样是115总行数350。3.4 调用linprog求解并提取结果模型装配完成后求解只需要一行%% 5. 求解线性规划 options optimoptions(linprog, Algorithm, dual-simplex, Display, iter); [x, fval, exitflag, output, lambda] ... linprog(f, A, b, Aeq, beq, lb, ub, options); %% 6. 结果提取 P reshape(x(1:ng*T), ng, T); R reshape(x(ng*T1:end), ng, T); % 影子价格 energy_price lambda.eqlin(1:T); % 能量出清价格 reserve_price lambda.eqlin(T1:2*T); % 旋转备用出清价格 fprintf(最小总成本: %.2f 元\n, fval);reshape 的坑点在于Matlab是按列填充的所以 reshape(x(1:ng*T), ng, T) 得到的结果矩阵里第 t 列正好就是第 t 时段所有机组的出力这个逻辑和变量排列顺序是对应的。linprog 的返回值里 lambda.eqlin 对应等式约束的影子价格lambda.ineqlin 对应不等式约束的影子价格。如果只需要看系统能量出清价格重点看能量平衡约束对应的对偶变量那一列。4. 出清结果解读产出什么、价格怎么来4.1 出力计划怎么读跑完代码后P矩阵的每一列代表一个时段所有机组的出力。观察几个典型时段会发现很有意思的规律夜间低谷时段比如凌晨2点到5点系统负荷只有510MW左右第1台机组200MW、第2台180MW、第3台160MW刚好几乎覆盖全部负荷第4、5台机组的出力可能是0或者刚好留有小出力。此时边际机组大概率是第3台能量出清价格接近350元/MWh。晚高峰时段20点左右负荷825MW左右5台机组总容量760MW明显不够。如果按这个机组配置模型会产生不可行解这时候需要回看3.1节的假设要么降低负荷峰值要么增加机组容量要么允许系统失负荷并引入惩罚项。实际市场不可能容忍“无解”所以工程模型里普遍会加失负荷变量和失负荷价值VOLL惩罚这也是第5节要重点讲的问题。4.2 旋转备用分配的逻辑备用容量分配结果往往是新手最容易困惑的部分。明明第1台机组报价15元/MW最便宜为什么不是所有备用都给它因为机组同时受 PRPmax 的约束。第1台机组可能已经满发承担能量任务了就算备用报价再便宜也腾不出容量来提供备用。这时候备用会分配给那些出力还没顶到上限、且有剩余容量的机组。这就是联合出清的价值所在能量和备用同时竞价系统自动找到“谁发电、谁留备用”的整体最优组合。如果分开优化先调度能量再分配备用很可能出现备用全部压在昂贵机组上的次优结果。4.3 价格信号为什么备用价格可能很低线性规划的对偶变量给出的价格信号取决于约束的松紧程度。如果某时段备用约束是松弛的也就是可用备用容量远大于需求那么备用价格会非常低甚至接近0。这符合市场直觉供过于求时边际备用成本几乎为零。反之如果备用需求很紧张最后1MW备用必须从报价很高的机组那里挤出来备用价格就会跳涨。至于能量价格则通常由边际机组报价决定。晚高峰负荷接近机组总容量时边际机组可能是第5台报价450元/MWh此时能量价格会整体拉高而备用价格反而不一定高——因为高报价机组已经满发根本无法提供备用。价格信号是整个模型最有价值的输出。建议跑完代码后把 energy_price 和 reserve_price 画在同一个图上观察两条曲线的走势和机组报价之间的关系会明显感受到“市场机制”这四个字是怎么从公式中冒出来的。5. 实操中一定会遇到的坑5.1 linprog报 “No feasible solution” 怎么办这个报错几乎是所有第一次跑通这个模型的人都会碰到的。原因可能有两类第一种是数据本身矛盾比如系统总容量小于峰值负荷。我上面给的机组参数总容量是760MW而晚高峰负荷超过800MW如果照抄代码一定会无可解。解决方法也很直接要么把负荷峰值调低要么加一台大容量机组要么在模型里加入失负荷变量和失负荷惩罚让模型在极端场景下可以“甩负荷”但同时付出极高代价。第二种是备用需求设置得太狠比如备用需求设为负荷的20%加上机组Rmax偏小导致某些时段无论如何都无法同时满足负荷和备用。排查思路是先检查所有时段的总容量是否大于“负荷备用需求”的最大值再看有没有单台机组容量特别大、导致其满发时其他机组备用总和不足的情况。实际工程里备用需求一般要大于最大单机容量防止最大机组跳闸时没有足够备用顶上。5.2 备用电价算出来是0或者低得离谱出现这种情况首先要确认这不是bug而是模型在告诉你“备用供给过剩”。但如果实际市场不可能是0通常是因为模型里缺少备用容量的机会成本机制。前面我用的简化模型没有显式计算机组“发不了电”的机会成本只把Cr当成备用容量报价。实际市场中机组参与备用市场时会把能量市场的收益损失打进报价里报价会高于单纯的成本项。如果你想在模型里体现这一点可以在 Cr 之外增加一台“虚拟备用机组”或者把备用需求和能量需求强耦合更精细的做法是引入情景约束scenario-based constraints模拟备用被调用时的能量市场影响。这些属于进阶内容但理解了影子价格机制后扩展起来思路会很清楚。5.3 矩阵维度口算对不上、索引越界调试过程中最常见的低级错误是 A 矩阵行数和 b 向量长度不一致或者变量索引写错导致某些列没有被赋值。我踩过最深的一个坑是把变量顺序从“按时段排列”改成“按机组排列”后忘记同步更新目标函数系数 f 的顺序。结果模型还能跑但结果完全错乱。解决方案是写一个简单的检查函数在调用 linprog 前验证assert(length(f) n, 目标函数系数长度错误); assert(size(A,2) n, A矩阵列数与变量数不一致); assert(size(A,1) length(b), A矩阵行数与b长度不一致); assert(size(Aeq,2) n, Aeq矩阵列数与变量数不一致);这种断言写起来只要半分钟但能省下大把调试时间。另外建议先用小规模数据试跑比如 T 先取4个小时、3台机组手算验证结果合理后再放大到24小时。5.4 版本兼容性和工具箱问题linprog 属于 Optimization Toolbox。如果你的Matlab环境里没有这个工具箱运行时会直接报“未定义函数或变量”。解决办法有两个一是安装 Optimization Toolbox新版Matlab在命令行输入optimoptions如果提示未定义打开附加功能管理器搜索工具箱安装即可二是改用自写单纯形法或内点法但完全没必要工具箱的 linprog 实现成熟且支持大规模稀疏矩阵。另外在不同Matlab版本里linprog 的接口略有差异。R2017b之后推荐用optimoptions设置参数老版本用optimset。如果代码里出现“Invalid OPTIONS”之类的报错优先检查是不是语法版本不匹配。6. 往工程实用方向的三个扩展6.1 从经济调度改为机组组合这个版本假设机组都处于开机状态。实际市场出清需要同时决定每台机组要不要开机这就要引入整数变量。把连续变量 P、R 之外再增加一个启停变量 u(i,t)∈{0,1}表示机组是否开机。同时增加出力上下限与启停状态的耦合u(i,t)*Pmin(i) P(i,t) u(i,t)*Pmax(i)最小启停时间约束这样整个模型就从线性规划变成混合整数线性规划求解器也从 linprog 换成 intlinprog。变量规模不大时intlinprog 完全够用规模大了就需要搭配商用求解器。6.2 从系统价格扩展为节点价格单节点模型只能给出系统统一的能量价格无法反映网络阻塞。要研究输电阻塞对出清结果的影响需要引入直流潮流模型。把每条线路的潮流表示成节点注入功率的线性组合加上线路传输容量约束。这时候不同节点的功率平衡约束影子价格就是节点边际电价LMP。同一时段不同节点的价格差异直接反映阻塞位置和阻塞严重程度。相比单节点模型这个扩展会让结果分析多出非常多的信息量。6.3 给研究者的三条实战建议第一先跑通小规模算例再去碰大规模数据。小算例可以手算验证能快速定位模型错误。我一开始直接上96时段、上百台机组的数据矩阵一跑就崩排查了三天发现是索引错了一位。第二价格信号比调度结果更值得深挖。很多论文和工程报告的核心图不是各机组出力堆叠图而是出清价格曲线。价格发生跳变的时段往往对应约束从松弛变为紧约束的临界点那是理解市场机制最有价值的地方。第三多做敏感分析。把备用需求比例从10%改成15%看总成本和价格的变化把某台机组的备用报价提高20%看中标结果会不会变。这些实验能帮你快速理解模型里每个参数对结果的影响路径比反复读公式更高效。做这个模型最大的收获不是学会了怎么调linprog而是理解了出清模型本质上是把市场规则翻译成数学约束。所有你认为理所当然的市场直觉都可以在模型里找到对应的约束和价格信号来验证。自己动手搭一遍之后再去看实际市场规则文档会发现那些复杂的条款背后就是这一套优化模型不断加细节扩展出来的。