说到电力市场出清很多人第一反应是电能量市场里那个价高者得、按边际价格统一出清的过程。但实际电力系统运行还要同时保证频率、电压和事故备用这些功能不是发电报价就能自动覆盖的需要单独建立辅助服务市场。旋转备用就是辅助服务里最基础、也最绕不开的一个品种——它要求并网机组在10分钟内能同步增加出力应对负荷波动、新能源出力突变或者机组跳闸。把这部分需求也纳入市场出清就构成了主辅助服务市场出清模型今天这篇就围绕旋转备用这个品种把建模思路、数学约束和Matlab实现一起讲透。如果你正在做电力市场方向的研究、课程设计或者刚开始接触辅助服务出清这篇文章应该能帮你省不少事。我会先拆清楚模型在算什么再给出可以直接复现的Matlab代码YALMIP和原生linprog两个版本都有最后聊聊我在实际调试中踩过的坑。尽量让没有太多优化背景的同学也能跟着跑通。1. 先搞清楚旋转备用出清到底在算什么1.1 从辅助服务到旋转备用电力系统的电能生产和消费必须实时平衡但负荷会波动、风机光伏会受天气影响、机组也可能突然故障。为了保证供电可靠性调度机构需要预留一部分富余发电能力这部分能力就是备用。备用又分很多种按响应速度分有瞬时响应的一次调频、10分钟内的旋转备用、30分钟内的非旋转备用按是否并网分有热备用和冷备用。其中旋转备用Spinning Reserve指的是已经并网、处于旋转状态、能在规定时间内通常是10分钟增加出力的机组容量。它不需要经历启动过程响应最快也是最贵的辅助服务之一。因为在实时运行中它必须时刻“在线”这部分容量即使没有被实际调用也要按照中标容量付费相当于给系统买了一份“随时能增出力的保险”。1.2 出清模型解决什么问题出清模型要回答的核心问题有两个一是哪些机组提供旋转备用各提供多少二是旋转备用的结算价是多少。如果只单独优化备用不考虑电能成本很可能出现最便宜的备用机组其实已经满发、根本没有余量这种荒唐结果。所以现在的主流做法是把电能量市场和旋转备用市场联合出清让同一个机组在电能和备用之间做权衡。举个例子一台机组已经发了90%的额定出力它报的备用价格再低实际也提供不了多少备用容量因为它没有往上调的空间。反过来一台机组为了多卖备用可能故意压低电能量报价把出力压在下限附近留出充足的上调空间。联合出清模型要考虑的就是这种电能备用耦合关系在总购电费用最小的目标下同时决定每台机组的电能中标量和备用中标量。1.3 为什么选Matlab做这套模型Matlab在电力系统领域的使用率一直很高原因很直接矩阵运算方便、画图交互顺手、YALMIP这类建模工具箱让优化问题的写法接近数学公式不用手动拼一大堆A、b矩阵。更重要的是Matlab自带的Optimization Toolbox提供的linprog、intlinprog就可以解决线性规划和混合整数线性规划这套出清模型恰好属于线性规划范畴不需要重型商业软件也能跑。这套模型的数学本质是一个线性规划目标函数是购电费用最小约束条件全部是线性的等式和不等式决策变量是各机组电能出力、旋转备用中标量。规模不大时用linprog完全够用规模大了换Gurobi、CPLEX也只需要在YALMIP里改一行求解器设置。所以我的建议是先用Matlab原生工具把模型逻辑调通再根据数据规模决定要不要上商业求解器。2. 数学模型目标函数与约束条件怎么搭2.1 目标函数购电费用最小化与报价驱动的市场逻辑旋转备用市场出清的经典目标函数是系统总购电费用最小包括电能量费用和备用费用两部分min Σ (c_i × p_i q_i × r_i)其中p_i是机组i的电能量中标出力MWr_i是机组i的旋转备用中标容量MWc_i是机组i的电能量报价元/MWhq_i是机组i的旋转备用报价元/MW。这里要强调一个市场设计逻辑报价来自发电企业在市场申报时提交的成本曲线优化目标是最小化按报价计算的购电费用而不是真实成本。也就是说出清模型是一个价格驱动的问题报价越低机组在中标排序上越有优势但最终所有中标机组都按统一的边际价格结算这就是电力市场里按报价排序、按边际价格出清的基本玩法。在目标函数里把电能和备用加起来正是联合出清的关键。如果只优化电能备用需求就会被忽略如果分两步先出清电能再出清备用就会出现前面说的备用机组没有上调空间的问题。联合优化让模型自己在电能收入和备用收入之间做取舍结果更符合实际运行约束。2.2 四类核心约束的物理意义约束条件是出清模型的灵魂。我把基础的旋转备用出清模型通常包含的约束分成四类每一类背后都有明确的物理意义。第一类是系统功率平衡约束Σ p_i P_wind D这个约束表示全系统所有电源出力之和必须等于负荷。D是系统负荷P_wind是风电等新能源出力。在简化的研究模型里新能源出力通常作为固定值处理不考虑弃风弃光如果考虑新能源参与调峰P_wind也可以变成决策变量加入目标函数。第二类是旋转备用需求约束Σ r_i ≥ R_needR_need是系统所需的旋转备用容量一般按最大单机容量的一定比例或负荷的一定比例确定。工程上常见做法是取系统最大一台机组容量和预测负荷的5%—10%的较大者。这个约束保证了系统在任何一台机组跳闸时剩余机组有能力在10分钟内补上缺口。第三类是机组技术约束Pmin_i ≤ p_i ≤ Pmax_i0 ≤ r_i ≤ Rmax_ip_i r_i ≤ Pmax_i前两个不等式是电能出力和备用中标量的上下限Rmax_i是机组申报的最大备用容量。第三个不等式是电能与备用的耦合约束表示机组同时卖出电能和备用后总占用容量不能超过额定出力上限。这个约束决定了满发机组无法再卖备用的物理事实。第四类是可选的爬坡约束p_i - P0_i ≤ RU_iP0_i - p_i ≤ RD_ir_i (p_i - P0_i) ≤ RU_i其中P0_i是机组当前出力RU_i和RD_i是10分钟内的上调和下调爬坡能力。第二条是旋转备用的核心——机组在预留备用容量后必须仍有能力在备用响应时间内完成上调。这个约束在实时市场出清里非常重要但在简化的日前出清模型中经常被省略。我在代码里会用它但默认打开与否取决于你的研究深度。2.3 边际价格出清价格从哪来出清价格的数学含义是约束对应的对偶变量拉格朗日乘子。在最小化问题中功率平衡约束的对偶变量就是电能量的节点边际价格LMP备用需求约束的对偶变量就是旋转备用的边际价格。理解对偶价格可以从经济学角度想边际价格表示系统多增加1个单位需求时总购电费用的边际变化。如果备用需求从60MW增加到61MW系统需要多买1MW最便宜的可用备用这个新增成本就是备用边际价格。在出清结果里它恰好等于最后一个中标备用机组的报价也就是市场出清价在统一价格机制下。在代码实现里YALMIP的dual函数可以直接提取约束的对偶值linprog则通过返回的lambda结构体获取。这里有个很容易踩的坑linprog对不等式约束统一按A*x ≤ b的形式处理如果你把sum(r) ≥ R_need改写成-sum(r) ≤ -R_need得到的对偶乘子是负值需要取相反数才是真实的备用边际价格。3. Matlab代码实现从YALMIP到原生linprog3.1 算例数据说明我构造一个四台火电机组加一个风电场的简化系统方便展示模型逻辑。系统负荷680MW风电出力80MW旋转备用需求取净负荷600MW的10%即60MW。四台机组的参数如下表。机组Pmin (MW)Pmax (MW)P0 (MW)电能量报价 (元/MWh)备用报价 (元/MW)Rmax (MW)G110040025028020100G2803502003201680G3502001504003060G430150805004540这些参数的设计思路是G1便宜但备用价格不算最低G2电能量价格中等但备用价格最低G3、G4电能量价格高正常情况下只承担最低技术出力。这样设置便于观察联合出清下电能和备用的权衡结果。3.2 用YALMIP建模求解YALMIP是一个Matlab建模工具箱语法和数学表达式非常接近。安装好YALMIP后github上直接下载并添加到路径即可不管底层用哪个求解器代码都不用改。下面是我常用的建模代码%% 主辅助服务市场出清模型旋转备用联合出清YALMIP版 clear; clc; close all; %% 1. 机组参数 % 列含义Pmin, Pmax, P0, 电能报价, 备用报价, Rmax G [ 100, 400, 250, 280, 20, 100; 80, 350, 200, 320, 16, 80; 50, 200, 150, 400, 30, 60; 30, 150, 80, 500, 45, 40 ]; Pmin G(:,1); Pmax G(:,2); P0 G(:,3); c G(:,4); q G(:,5); Rmax G(:,6); % 系统参数 D 680; % 负荷 MW Wind 80; % 风电出力 MW Rneed 60; % 旋转备用需求 MW n length(c); % 机组数 %% 2. 决策变量 p sdpvar(n, 1); % 电能中标出力 MW r sdpvar(n, 1); % 旋转备用中标容量 MW %% 3. 目标函数 Objective c * p q * r; %% 4. 约束条件 Cons []; % 功率平衡火电总出力 风电 负荷 Cons [Cons, sum(p) Wind D]; % 旋转备用需求 Cons [Cons, sum(r) Rneed]; % 机组出力上下限 Cons [Cons, Pmin p Pmax]; % 备用容量上下限 Cons [Cons, 0 r Rmax]; % 电能与备用耦合p r 不能超过机组上限 Cons [Cons, p r Pmax]; % 爬坡约束若研究实时出清可取消注释 % Cons [Cons, p - P0 50]; % 10min上调能力 % Cons [Cons, P0 - p 60]; % 10min下调能力 % Cons [Cons, r (p - P0) 50]; % 上调能力与备用耦合 %% 5. 求解 ops sdpsettings(solver, gurobi, verbose, 1); % 没有Gurobi时可以用 solver, linprog 或直接不指定 sol optimize(Cons, Objective, ops); if sol.problem ~ 0 warning(求解失败%s, sol.info); else P_opt value(p); R_opt value(r); fprintf( 旋转备用联合出清结果 \n); for i 1:n fprintf(机组G%d: 电能出力 %6.2f MW, 备用中标 %6.2f MW\n, ... i, P_opt(i), R_opt(i)); end fprintf(火电总出力: %.2f MW, 风电: %.2f MW, 负荷: %.2f MW\n, ... sum(P_opt), Wind, D); fprintf(备用总计: %.2f MW, 备用需求: %.2f MW\n, ... sum(R_opt), Rneed); fprintf(电能总费用: %.2f 元, 备用总费用: %.2f 元\n, ... c * P_opt, q * R_opt); % 边际价格 LMP dual(Cons(1)); % 功率平衡约束的对偶 - 电能边际价格 RP dual(Cons(2)); % 备用需求约束的对偶 - 备用边际价格 fprintf(电能边际价格(LMP): %.2f 元/MWh\n, LMP); fprintf(旋转备用边际价格: %.2f 元/MW\n, RP); end这段代码有个细节值得说明我建议把数据集成为矩阵G而不是写四个独立的变量这样机组数量改变时不用大改代码。你在做研究时如果要对10台、100台机组做测试只需要把G矩阵扩展就行。3.3 不依赖第三方工具箱的linprog版本YALMIP虽然好用但并不是Matlab自带的。我遇到过不少同学在装YALMIP时被依赖关系、路径设置搞得焦头烂额。所以这里再给一个纯原生版本只用Matlab自带的linprog适合快速验证模型逻辑也能让你更清楚底层到底在解什么样的矩阵。%% 主辅助服务市场出清模型旋转备用联合出清linprog版 clear; clc; close all; %% 1. 机组参数同上 G [ 100, 400, 250, 280, 20, 100; 80, 350, 200, 320, 16, 80; 50, 200, 150, 400, 30, 60; 30, 150, 80, 500, 45, 40 ]; Pmin G(:,1); Pmax G(:,2); c G(:,4); q G(:,5); Rmax G(:,6); n length(c); D 680; Wind 80; Rneed 60; %% 2. 决策变量 x [p1; p2; p3; p4; r1; r2; r3; r4] f [c; q]; % 目标函数系数 %% 3. 等式约束sum(p) Wind D Aeq [ones(1, n), zeros(1, n)]; beq D - Wind; %% 4. 不等式约束 % sum(r) Rneed - -sum(r) -Rneed A1 [zeros(1, n), -ones(1, n)]; b1 -Rneed; % p r Pmax A2 [eye(n), eye(n)]; b2 Pmax; A [A1; A2]; b [b1; b2]; %% 5. 变量上下限 lb [Pmin; zeros(n, 1)]; ub [Pmax; Rmax]; %% 6. 调用linprog求解 options optimoptions(linprog, Display, iter, Algorithm, dual-simplex); [x, fval, exitflag, output, lambda] linprog(f, A, b, Aeq, beq, lb, ub, options); %% 7. 结果整理 if exitflag 0 P_opt x(1:n); R_opt x(n1:2*n); fprintf( 旋转备用联合出清结果linprog\n); for i 1:n fprintf(机组G%d: 电能出力 %6.2f MW, 备用中标 %6.2f MW\n, ... i, P_opt(i), R_opt(i)); end LMP lambda.eqlin; % 电能边际价格 RP -lambda.ineqlin(1); % 备用边际价格注意取负 fprintf(电能边际价格(LMP): %.2f 元/MWh\n, LMP); fprintf(旋转备用边际价格: %.2f 元/MW\n, RP); fprintf(电能总费用: %.2f 元, 备用总费用: %.2f 元, 总费用: %.2f 元\n, ... c * P_opt, q * R_opt, fval); end这里要特别解释lambda.ineqlin那个负号。linprog返回的lambda.ineqlin对应的是A*x ≤ b这类不等式非负量。但原始模型里的备用需求是sum(r) ≥ Rneed是≥约束我把它改写成-sum(r) ≤ -Rneed后对偶乘子会是负值真实的备用边际价格就是它的相反数。这一步看似小但我见过不少人直接把lambda.ineqlin(1)当备用价格输出结果价格是负数怎么看怎么不对。3.4 结果输出与边际价格提取跑完上面任意一个版本都会输出下面这样的结果各机组电能出力、备用中标量、系统总费用以及最重要的两个边际价格——电能LMP和备用边际价格。在YALMIP版本中直接用dual(Cons(1))和dual(Cons(2))就能分别拿到功率平衡约束和备用需求约束的对偶值。要注意的是YALMIP对dual的处理要求约束本身是原始形式且求解成功如果求解器返回的是对偶单纯形法结果dual通常都能正确提取。如果你想让结果更工程化可以把fprintf换成xlswrite或者writematrix把出清结果输出到Excel方便后续做灵敏度分析和图表展示。4. 算例结果解读出清量、出清价和里面的门道4.1 电能和备用联合出清结果运行上面的代码得到的最优结果如下G1电能出力400MW满发备用中标0MWG2电能出力120MW备用中标60MWG3电能出力50MW最低技术出力备用中标0MWG4电能出力30MW最低技术出力备用中标0MW火电总出力600MW加风电80MW正好等于负荷680MW。备用总容量60MW刚好满足备用需求。电能总费用为400×280 120×320 50×400 30×500 185400元。备用总费用为60×16 960元。系统总购电费用186360元。这个结果非常直观地体现了联合出清的逻辑G1虽然是最便宜的电能机组但已经满发剩余备用容量为0所以即使它报了备用价格也没有任何备用中标。G3、G4电能量报价太高只承担最低技术出力同样没有备用中标。所有旋转备用都落在了G2身上——它的备用报价最低而且电能出力只有120MW距Pmax上限还有230MW的余量完全有能力承接这60MW备用。4.2 为什么边际价格是这些值模型输出的电能边际价格LMP320元/MWh因为系统中最后满足负荷的边际机组是G2它的电能量报价就是320元/MWh。G1、G3、G4分别被上下限约束“卡住”不参与边际定价。旋转备用边际价格16元/MW对应G2的备用报价。因为系统只需要60MW备用而G2的备用量还有80MW的上限空间最后一个中标备用兆瓦依然来自G2所以边际价格就是G2的报价。如果备用需求增加到超过G2的备用能力比如需求变成100MW那么G2最多只能中标80MW剩下20MW需要从G1报价20元/MW那里买备用边际价格就会抬升到20元/MW。这里有一个在出清结果分析时非常关键的点边际价格不是所有中标机组的平均报价而是最后那1MW的边际成本。这就是统一出清价的核心含义——所有中标备用都按同一个边际价格结算而不是按各自的报价结算。所以G2虽然只报了16元/MW但它提供的全部60MW备用都按16元/MW结算没有吃亏也正因为这个机制机组才有动力报出自己的真实边际成本。4.3 联合出清与顺序出清的差异为了让你更直观地理解联合出清的价值我对比一下先出清电能、再出清备用的顺序出清结果。顺序出清的第一步只优化电能费用不考虑备用约束。由于G1最便宜G1发到400MWG2发200MWG3、G4分别在50MW和30MW最低点火电总出力恰好600MW。但是此时G1满发G2距Pmax还有150MW余量G3、G4余量分别为150MW和120MW。第二步优化备用时G2的备用报价最低且余量充足G2仍会中标60MW备用。这个结果和联合出清一模一样因为系统里有充足的廉价备用机组兜底。但你把算例改成负荷1200MW风电200MW火电净负荷1000MW情况就完全不同了。联合出清下模型会主动压低G1的出力量避免它满发腾出空间给备用顺序出清下G1满发后用最便宜的G2顶上等第二步买备用时可能发现最便宜的备用机组全都满发或接近满发只能高价买G3、G4的备用总费用明显更高。这就是为什么实际市场大多转向联合出清也解决了两步走的效率损失问题。5. 实操中踩过的坑与排查思路5.1 模型不可行从约束冲突找原因用这套模型最容易遇到的报错就是Problem infeasible。我总结了一下绝大多数不可行问题都出在以下几种情况可以直接对照排查。第一参数自相矛盾。比如Pmin之和大于净负荷也就是说所有机组最低出力加起来都超过负荷需求那功率平衡约束注定无法满足。我见过有人做风电出力很大的场景时忘记把火电出力下限调低结果系统不得不强迫停机但模型里又没有启停变量自然就无解了。第二备用需求设置过大。如果R_need大于全系统可用备用上限模型无解。全系统可用备用上限大概是Σ min(Rmax_i, Pmax_i - Pmin_i)所以如果你把R_need设成200MW而这个系统最大只能提供180MW备用那不管怎么优化都完成不了备用需求。第三爬坡约束和电能出力冲突。设了p - P0 ≤ RU但Pmin和Pmax之间可能根本没有一条从P0出发、在爬坡限制内能到达的路径。做实时出清时初始出力P0往往来自上一时段的出清结果如果上一时段和当前时段的负荷差异太大很容易出现爬坡约束卡死。我的排查顺序一般是先关掉备用需求约束看电能模型是否可行再关掉爬坡约束看是否能找到备用最后逐个加回约束用二分法定位到底是哪条约束在捣乱。YALMIP还可以配合ops sdpsettings(debug, 1)输出详细诊断信息。5.2 求解器选择与结果稳定性Matlab自带的linprog默认算法在高版本里通常是dual-simplex对小规模LP问题完全够用。我跑四机组算例时基本是毫秒级出结果。但如果你把机组规模放大到几百台、约束上万条还是建议用商业求解器。Gurobi和CPLEX是电力市场领域最常用的两个求解器性能稳定对LP和MILP的支持都很好。用YALMIP的话只需要改一行ops sdpsettings(solver, gurobi);如果系统里没装GurobiYALMIP会自动回退到linprog但这可能掩盖一些问题比如在MILP模型里linprog根本解不了整数变量YALMIP会报错。所以正式研究时我建议把求解器当成模型的一部分明确写在论文或报告里不要用默认求解器这种模糊表述。关于数值稳定性我印象最深的是单位不统一导致的病态约束。比如电能报价用元/MWh、备用报价用万元/MW、出力用kW目标函数和约束矩阵的数值量级差了上万倍求解器在计算对偶变量时很容易出现精度问题甚至报错。策略很简单全系统统一量纲我习惯全部用MW和元单位换算在数据导入阶段就完成不在模型里做系数变换。5.3 单位、数值标度和边界工况下面这个表基本涵盖了我这里跑模型时遇到过的典型问题直接贴出来供你对照。现象可能原因排查方式求解器报不可行Pmin之和超过净负荷检查负荷与风电功率的假设是否合理备用边际价格为负把“≥”约束按“≤”处理导致对偶符号错误linprog中对偶取负YALMIP直接用dual出清结果小数泛滥机组参数量纲不统一统一使用MW和元避免混用kW/万kW关掉爬坡约束就能解爬坡能力与Pmin/Pmax区间冲突检查P0是否在当前时段的可行区间内Gurobi求解时间炸裂模型里隐含有整数变量检查是否有逻辑变量必要时调整MIP Gap还有一种边界工况容易被忽略如果R_need恰好等于零模型会退化成一个纯电能量出清问题此时备用需求约束的对偶乘子为0备用边际价格就是0。这在数学上没问题但市场机制里不太可能出现不需要备用的场景所以如果你发现备用价格异常先回去检查R_need是不是被写成了0或太小。6. 往真实市场延伸模型还可以怎么扩展6.1 从单时段到多时段滚动出清前面给的是单时段模型假设各时段相互独立。实际电力市场有日前市场和实时市场日前市场通常以1小时为一个时段优化整个24小时的电能和备用计划。这时候模型要从单时段LP变成多时段LP机组爬坡约束会“跨时段”耦合t时段出力和t1时段出力之间的变化量不能超过爬坡速率乘以时段长度。扩展方式是给决策变量加一个时间下标p变成p(i,t)r变成r(i,t)然后增加时序耦合约束。比如p(i,t) - p(i,t-1) ≤ RU_i × Δt p(i,t-1) - p(i,t) ≤ RD_i × Δt目标函数也变成对所有时段求和。请注意单时段模型里的备用需求是系统总需求多时段模型里每个时段都有一个R_need(t)可以按各时段负荷水平动态设置。6.2 从确定性模型到随机/鲁棒优化基础模型假设负荷和风电出力都是确定的实际中风电预测误差很大一个更贴近工程实际的模型会把风电出力设成随机变量用场景法Scenario-based描述多种可能的风电出力构建两阶段随机优化模型。第一阶段决定火电机组的中标电量和备用容量第二阶段针对每个风电场景计算再调度成本和失负荷惩罚。目标函数变成确定性成本加各场景下期望成本的加权和。这样旋转备用容量就不仅是响应一个固定需求而是真正用来对冲风电不确定性。如果你觉得场景法计算量太大可以做鲁棒优化把风电出力看成在一个不确定集合内波动要求系统在最坏情况下也有足够备用。这类模型一般还是LP或者MILP形式求解器层面并不需要换只是约束构造要花更多心思。6.3 从LP到MILP机组组合联合出清基础模型默认所有机组都处于开机状态只优化出力分配。实际市场里机组是否开机本身就是决策一台机组如果电能量报价太高、备用报价也高最优选择可能是直接停机而不是被迫保持最小技术出力。这就需要引入二进制变量u_i表示机组启停状态模型变成MILP目标函数里加入启动成本出力下限约束变成Pmin_i × u_i ≤ p_i ≤ Pmax_i × u_i备用约束变成r_i ≤ Rmax_i × u_i。求解器从linprog换成intlinprog或者用YALMIP配Gurobi/CPLEX。加入启停变量后模型的复杂度明显上升几十台机器的MILP可能就需要数秒到数分钟。这里有一个工程心得如果你只是研究旋转备用出清机制建议先用LP版本把基本逻辑跑通再加启停变量。LP结果对参数变化非常敏感调参方便能帮你快速建立对市场的直觉直接上MILP一旦模型不可行定位问题会困难很多。我在实际做这个模型时最大的体会是旋转备用出清看起来只是一个LP问题但它的价值全在约束条件的物理含义里。功率平衡约束决定了电能市场的基本盘备用需求约束决定了系统要花多少钱买“保险”而pr≤Pmax这条耦合约束才是联合出清的精髓。把这几条约束想明白了后面无论做多时段、随机优化还是机组组合都是在同一套骨架上的扩展。最后分享一个调试小技巧拿到一段出清代码先别急着跑大数据把机组数改成3台、负荷设成刚好卡在边界上的值手算一遍目标函数和边际价格再去对代码的输出。代码值和手算值对上了说明模型写对了对不上多半是约束符号或对偶方向的问题。这种“小算例对拍”的方法我用了很多年比任何调试工具都好使。