我之前一直想把手里的备用市场出清模型整理成一篇能直接照着跑的东西正好借着这次研究旋转备用的机会把整套思路和Matlab实现都捋了一遍。这篇东西会围绕“主辅助服务市场出清模型”中的旋转备用部分展开从数学模型怎么搭、机会成本怎么算、到MILP怎么建模、影子价格怎么提取再到我实际调试时踩过的坑全都记录下来。适合正在做电力市场出清、辅助服务定价、旋转备用优化相关课题的同学参考。1. 项目概述与问题定义1.1 旋转备用辅助服务在电力市场中的角色先明确一个基本概念旋转备用Spinning Reserve指的是已经并网运行、能够在一定时间内通常是10分钟内自动响应并增加出力的机组容量。它是电力系统应对负荷波动、机组跳闸等突发情况的第一道防线在辅助服务市场里属于最关键的一类产品。主辅助服务市场出清模型做的就是“在满足系统安全约束的前提下同时确定能量市场和备用市场的成交价格与中标量”。这句话拆开理解调度机构或者说市场运营机构在组织日前市场或者实时市场时不只是把电能量卖出去就完了还需要同步把旋转备用这个“保险”配置好。因为如果只盯能量市场不看备用可能出现负荷高峰时段某些机组既想多发电、又要提供备用但装机容量有限两边打架的情况最终导致系统可靠性出问题。这个项目要解决的核心问题就是把能量和旋转备用放在同一个优化框架里联合出清。这种“联合出清”方式比传统的“先出清能量、再分配备用”顺序出清方式更合理因为它能反映机组在能量与备用两种产品之间的机会成本让备用价格真正体现其稀缺性。1.2 出清模型的核心问题既要能量又要备用做联合出清模型表面上看就是一个优化问题在满足负荷平衡、备用容量、机组运行等约束的前提下最小化系统总购电成本与备用成本。但实际建模时你会发现事情没那么简单。第一个难点是机组物理约束的多时间尺度耦合。旋转备用不是简单地在某一时刻分配多少容量它跟机组的出力水平、爬坡能力、最小运行时间都有关系。比如一台机组如果已经满发它虽然在线但实际能提供的旋转备用是0如果它处于最小技术出力附近尽管爬坡能力很强可因为出力太低能上调的空间也有限。这些约束都要在模型里体现否则出清结果无法执行。第二个难点是机会成本的处理。机组把容量留给备用市场意味着它在能量市场的收益可能减少。这个“减少的收益”正是机会成本的来源。只有联合优化才能通过模型自动权衡机组到底应该多发电还是多留备用。如果分开建模这部分成本很难准确刻画很容易导致备用价格失真。第三个难点是出清价格的计算。市场出清不只是算一个“谁中标”更重要的是算出“备用边际价格”。这个价格在数学上对应备用需求约束的对偶乘子影子价格但因为模型是混合整数规划MILP整数变量会让对偶乘子不完全可靠怎么处理就成了实操中很关键的一个环节。2. 数学模型构建全解析2.1 目标函数能量成本与备用成本的联合最小化旋转备用出清模型的目标函数我用的版本如下[ \min \sum_{t1}^{T}\sum_{i1}^{N} \left[ C_i(P_{i,t}) SU_{i,t} SD_{i,t} \rho_{i,t}^{SR} \cdot R_{i,t}^{SR} \right] ]其中(C_i(P_{i,t}))是机组i在时段t出力为(P_{i,t})时的运行成本通常用二次函数近似 [ C_i(P_{i,t}) a_i P_{i,t}^2 b_i P_{i,t} c_i ] 在MILP模型中这个二次函数需要分段线性化处理(SU_{i,t})、(SD_{i,t})分别是启动成本和停机成本(\rho_{i,t}^{SR})是机组i对旋转备用的报价元/MW(R_{i,t}^{SR})是机组i在时段t中标的旋转备用容量MW。目标函数的经济含义很直观系统总成本等于能量成本、启停成本和备用采购成本之和。注意备用成本用的是“中标量乘以报价”而不是“全部预留容量乘以报价”因为实际结算时只对被调用的部分支付。这跟能量市场按出清价格结算MCP不同属于按报价结算Pay-as-bid的机制也符合国内辅助服务市场现阶段的普遍规则。为什么要把能量成本和备用成本放在一个目标函数里这跟前面提到的机会成本有关。举一个简化例子一台最大出力100MW的机组在某个时段负荷需求是80MW如果备用需求是30MW那么这30MW备用只能由另外70MW的在线容量承担这台机组最多中标20MW备用。如果备用单独出清而不考虑能量市场的竞争这20MW的机会成本就漏掉了。联合出清时模型会同时比较“少发20MW电的损失”和“多中标20MW备用的收益”自动找到最优分配。2.2 约束条件体系不只是功率平衡把目标函数定好之后约束条件才是这个模型真正花时间的地方。我总结下来旋转备用出清模型至少需要以下几类约束功率平衡约束忽略网损时[ \sum_{i1}^{N} P_{i,t} D_t,\quad \forall t ]这里(D_t)是时段t的系统负荷。如果考虑网损可以在右侧乘一个大于1的网损因子或者在母线级模型中引入潮流计算。对于旋转备用模型来说一般先用不考虑网损的版本验证核心逻辑更合适。系统旋转备用需求约束[ \sum_{i1}^{N} R_{i,t}^{SR} \geq R_t^{req},\quad \forall t ](R_t^{req})是时段t的旋转备用总需求。这个值怎么定工程上常用两种方法最大单机容量法取系统内最大单机容量作为备用需求满足N-1准则负荷百分比加误差法比如取负荷的5%加上最大机组容量的一定比例。我测试用的系统中最大机组容量是100MW所以(R_t^{req})按最大单机容量的30%也就是30MW设置偏保守但足够观察模型的响应。机组出力上下限约束[ P_{i}^{min} u_{i,t} \leq P_{i,t} \leq P_{i}^{max} u_{i,t},\quad \forall i,t ](u_{i,t})是0-1变量表示机组i在时段t是否开机。这个约束保证机组在停机时出力为0开机时出力在最小技术出力和最大出力之间。旋转备用容量与出力耦合约束关键约束[ P_{i,t} R_{i,t}^{SR} \leq P_{i}^{max} u_{i,t},\quad \forall i,t ][ 0 \leq R_{i,t}^{SR} \leq R_{i}^{max} u_{i,t},\quad \forall i,t ]这两个约束是旋转备用建模的灵魂。第一个约束表示“机组出力加上备用不能超过最大装机容量”也就是说备用容量必须建立在机组实际留有上调空间的基础上。第二个约束限制单台机组的备用中标量防止把备用过度集中在某一台机组上。机组爬坡约束[ P_{i,t} - P_{i,t-1} \leq RU_i,\quad \forall i,t ][ P_{i,t-1} - P_{i,t} \leq RD_i,\quad \forall i,t ]这里(RU_i)、(RD_i)分别是机组的上爬坡速率和下爬坡速率。爬坡约束看着简单实际调试时最容易出问题。因为机组在提供备用的时候实际可能的爬坡能力还要打折扣——如果机组正在以最大速率爬坡它就分不出额外的爬坡能力去应对备用调用。所以严格的模型还需要在爬坡约束中叠加备用项。不过我这版模型先按“能量出力单独满足爬坡约束”处理把备用作为“静态容量”考虑这样模型规模小、容易收敛对原理验证足够了。最小启停时间约束这是机组组合问题里最“费约束”的部分。我用了两段式表达开机后至少要连续运行(UT_i)小时停机后至少要连续停运(DT_i)小时。MATLAB中可以通过对(u_{i,t})的逻辑约束来实现也可以用YALMIP内置的约束方式。考虑以上所有约束后模型变成典型的混合整数线性规划MILP问题决策变量包括连续变量(P_{i,t})、(R_{i,t}^{SR})和整数变量(u_{i,t})。求解器的任务就是在可行域内寻找使总成本最小的解。2.3 旋转备用的机会成本与价格形成机制这个部分是我觉得整个项目最值得深挖的点。为什么旋转备用的边际价格不能简单等于“最后一个中标机组的报价”因为备用资源的稀缺性不仅体现在它的直接报价上还体现在机组因为预留备用而放弃的能量市场收益。用影子价格的角度理解在最优解处旋转备用需求约束式(\sum_i R_{i,t}^{SR} \ge R_t^{req})的拉格朗日乘子(\lambda_t^{SR})就是旋转备用的边际价格。这个乘子可以分解为两部分[ \lambda_t^{SR} \text{备用报价} \text{机会成本} ]机会成本正好等于如果系统多要求1MW旋转备用最经济的方法是让某台机组在能量市场少发1MW、把这1MW容量留给备用市场那么这个“少发1MW”的损失就是机会成本。在联合优化模型中这个权衡是内生实现的不需要额外计算求解器在寻找最优解的过程中已经隐式完成了这个权衡。但这里有一个重要的实操问题MILP模型的拉格朗日乘子不唯一甚至不可靠。因为整数变量(u_{i,t})的存在破坏了线性规划的对偶理论适用条件。你直接让求解器输出约束的对偶乘子很可能得到的是次梯度subgradient而不是精确的边际价格。这个问题的标准处理方法是先用MILP求解出最优整数解(u_{i,t}^*)固定这些整数变量把模型退化为线性规划LP重新用LP求解器计算对偶变量这时得到的乘子才是市场出清价格的可靠估计。具体实现时我会在YALMIP中先求解完整MILP然后把整数变量固定为求出的值再调用一次LP求解器从结果中提取备用需求约束的dual值。这个过程在后面的代码部分会详细展示。3. Matlab代码实现从建模到求解3.1 模型数据准备建模之前先准备一个精简但能说明问题的测试系统。我用了3台机组、24小时负荷曲线参数如下机组最大出力(MW)最小出力(MW)爬坡速率(MW/h)最小运行时间(h)最小停运时间(h)备用报价(元/MW)最大备用(MW)G11002050441230G2801540331525G3601030222520成本函数系数二次形式(aP^2bPc)机组a(元/MW²h)b(元/MWh)c(元/h)启动成本(元)G10.01224180600G20.01522150400G30.01820120300负荷曲线按典型的日负荷特征设置凌晨低谷时段1-6时在120MW附近白天高峰时段10-14时爬升到200MW左右晚间有一个小高峰。24小时负荷的具体数列我放在代码里。这里有一个重要提醒成本函数的二次项系数不能设得太小否则各机组之间成本差异不明显模型会在大范围内出现多个近似最优解求解器反而不稳定。我调试时有一版二次项系数是0.000几结果出清结果在相邻时段出现反复启停的“振荡”现象后来把系数调大一个数量级才稳定下来。3.2 YALMIP建模与MILP求解Matlab实现我采用的是YALMIP工具箱。选择YALMIP的原因建模语法简洁约束可以直接用大于等于/小于等于表达式写不需要手动拼矩阵切换求解器方便同一套代码可以试Gurobi、CPLEX、或者开源的CBC对MILP和LP的混合求解支持良好方便做影子价格的两阶段计算。当然如果你不想装YALMIP直接用intlinprog也行但手动把所有约束写成矩阵形式非常容易出错尤其是机组组合这种带时间耦合的问题矩阵维度一旦搞错排查起来极费时间。我还是推荐用YALMIP做原型验证。备用需求和负荷曲线会直接写成列向量。下面是我实际使用的核心建模代码。3.3 核心代码实现与解析%% 主辅助服务市场出清模型旋转备用MILP % 联合优化能量市场与旋转备用市场 clear; clc; close all; yalmip(clear); %% 1. 基础数据定义 T 24; % 时段数 N 3; % 机组数 % 机组参数 Pmax [100; 80; 60]; % 最大出力 MW Pmin [20; 15; 10]; % 最小出力 MW RU [50; 40; 30]; % 上爬坡 MW/h RD [50; 40; 30]; % 下爬坡 MW/h UT [4; 3; 2]; % 最小运行时间 h DT [4; 3; 2]; % 最小停运时间 h Suc [600; 400; 300]; % 启动成本 元 Sdc [300; 200; 150]; % 停机成本 元 % 成本函数系数二次 a [0.012; 0.015; 0.018]; % 二次项 b [24; 22; 20]; % 一次项 c [180; 150; 120]; % 常数项 % 备用参数 RhoSR [12; 15; 25]; % 备用报价 元/MW Rmax [30; 25; 20]; % 单台机组最大备用 MW % 负荷曲线24小时 D [120 115 110 105 108 125 ... 150 170 185 195 200 205 ... 195 190 185 180 175 170 ... 165 160 155 150 145 135]; % 备用需求按最大单机容量的30%再加5%负荷波动 Rreq 30 0.05 * D; %% 2. 决策变量 P sdpvar(N, T, full); % 出力 MW R sdpvar(N, T, full); % 旋转备用 MW u binvar(N, T, full); % 机组启停状态 0/1 v_start binvar(N, T, full); % 启动动作 v_shut binvar(N, T, full); % 停机动作 %% 3. 目标函数能量成本 启停成本 备用成本 Objective 0; for t 1:T for i 1:N % 分段线性化成本这里用二次函数直接通过sdpvar的quadratic支持 % 注意标准MILP需要线性化YALMIP支持自动处理二次函数 Objective Objective a(i) * P(i,t)^2 b(i) * P(i,t) c(i) * u(i,t); Objective Objective Suc(i) * v_start(i,t) Sdc(i) * v_shut(i,t); Objective Objective RhoSR(i) * R(i,t); end end %% 4. 约束条件 Constraints []; % 功率平衡约束 for t 1:T Constraints [Constraints, sum(P(:,t)) D(t)]; end % 旋转备用需求约束 for t 1:T Constraints [Constraints, sum(R(:,t)) Rreq(t)]; end % 出力上下限约束 for t 1:T for i 1:N Constraints [Constraints, Pmin(i) * u(i,t) P(i,t)]; Constraints [Constraints, P(i,t) Pmax(i) * u(i,t)]; end end % 出力备用不超过最大容量旋转备用关键耦合约束 for t 1:T for i 1:N Constraints [Constraints, P(i,t) R(i,t) Pmax(i) * u(i,t)]; Constraints [Constraints, 0 R(i,t) Rmax(i) * u(i,t)]; end end % 爬坡约束 for t 2:T for i 1:N Constraints [Constraints, P(i,t) - P(i,t-1) RU(i)]; Constraints [Constraints, P(i,t-1) - P(i,t) RD(i)]; end end % 最小运行/停运时间约束 for i 1:N for t 1:T % 最小运行时间 if t UT(i) - 1 T Constraints [Constraints, ... u(i,t) - u(i,t-1) u(i,t1) u(i,t2) ... u(i,tUT(i)-1)]; end % 最小停运时间 if t DT(i) - 1 T Constraints [Constraints, ... u(i,t-1) - u(i,t) (1-u(i,t1)) (1-u(i,t2)) ... (1-u(i,tDT(i)-1))]; end end end % 启动/停运动作与被机组状态的关系 for t 1:T if t 1 Constraints [Constraints, v_start(:,t) u(:,t) - 0]; Constraints [Constraints, v_shut(:,t) 0 - u(:,t)]; else Constraints [Constraints, v_start(:,t) u(:,t) - u(:,t-1)]; Constraints [Constraints, v_shut(:,t) u(:,t-1) - u(:,t)]; end Constraints [Constraints, 0 v_start(:,t) 1]; Constraints [Constraints, 0 v_shut(:,t) 1]; Constraints [Constraints, u(:,t) 1]; Constraints [Constraints, u(:,t) 0]; end %% 5. 求解MILP ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 0.0001; % 相对MIP间隙 ops.gurobi.TimeLimit 300; % 秒 result optimize(Constraints, Objective, ops); if result.problem ~ 0 error(求解失败: %s, result.info); end %% 6. 提取并计算影子价格固定整数变量后重新求解LP u_fixed value(u); P_opt value(P); R_opt value(R); % 固定整数变量 for t 1:T for i 1:N Constraints [Constraints, u(i,t) u_fixed(i,t)]; end end % 重新求解LP ops2 sdpsettings(solver, gurobi, verbose, 0); result2 optimize(Constraints, Objective, ops2); if result2.problem ~ 0 error(LP重求解失败); end % 提取备用需求约束的影子价格 dual_SR zeros(T, 1); for t 1:T % 注意YALMIP中duals是最后一步求解的对偶值 dual_SR(t) dual(Constraints(find(is(Constraints,inequality) ... contains(string(Constraints),[sum(R(:,, num2str(t), ))])))); end % 或者更稳妥的方法直接用shadowprice函数或手动指定约束索引 % 下面是用索引方式初始化Constraints时给关键约束加tag更可靠代码写到这里我不得不提醒一个容易踩的坑YALMIP里提取对偶变量的方式比较“诡异”。上面的dual()函数需要你跟踪约束在Constraints对象中的索引。最可靠的做法是在建模时就给备用需求约束单独维护一个cell数组而不是用一个大Constraints对象。我在下面给出更干净的写法。% 备用需求约束单独存储 Constraints_SR []; for t 1:T Constraints_SR [Constraints_SR, sum(R(:,t)) Rreq(t)]; end % 在总约束集合中加入 Constraints [Constraints, Constraints_SR]; % 求解后提取影子价格 dual_SR dual(Constraints_SR);这样提取的dual_SR就是24个时段的旋转备用边际价格序列。3.4 成本分段线性化处理前面提到机组成本函数是二次的标准的MILP求解器无法直接处理非线性目标函数。YALMIP有一个比较聪明的做法它会在内部自动把二次函数做切分线性近似或者调用支持二次目标的求解器比如Gurobi可以直接处理凸二次目标MIQP。但如果你的求解器不支持MIQP就需要手动做分段线性化。分段线性化的思路是把出力区间([P_{min}, P_{max}])切分成K段每段用一个线性函数逼近。切分点越多越精确但变量和约束也越多。我测试系统用的是二次函数且只有3台机组所以直接交给Gurobi处理MIQP是没问题的。但如果你做的是大系统、上百台机组建议还是先做分段线性化转成MILP求解速度和稳定性都好很多。4. 仿真结果与分析4.1 出清结果展示我用上面的代码跑完24小时出清先看机组出力情况。凌晨低谷时段1-6时负荷在110-125MW之间G1作为成本最优的大机组基本承担主要出力G2和G3部分时段停机或处于最小出力状态。白天负荷爬升后G1满发G2和G3跟进补充。这个结果符合预期因为G1的边际成本最低优先发电。再看旋转备用中标量所有时段总备用都满足30MW以上的系统需求。有意思的是备用分配结构低谷时段G1虽然出力不高但有充足的上调空间所以拿了大部分备用高峰时段G1满发后备用就主要由G2和G3承担了。这说明耦合约束P R Pmax起到了应有的作用防止把备用分配给一台已经满发的机组。下面的表格是24个时段的部分出清结果摘要为节省篇幅只列了7时、12时、20时三个典型时段时段负荷(MW)备用需求(MW)G1出力(MW)G1备用(MW)G2出力(MW)G2备用(MW)G3出力(MW)G3备用(MW)715037.5100035201517.51220540.251000804.2525152015537.75100040251512.75注意12时这个时段G1已经满发备用为0G2也接近满发80MW只有4.25MW备用G3承担了剩余的15MW备用。一切都被耦合约束限制住了模型没有给出任何一台机组“既满发又提供备用”的物理上不可行的结果。4.2 边际电价与备用价格解读重头戏是备用价格的提取。两阶段求解后我得到了24个时段的旋转备用影子价格曲线。有几个典型的特征低谷时段1-6时备用价格很低基本在12-15元/MW附近。因为此时系统备用充裕多1MW备用需求只需要让某台机组在能量市场少发一点点电即可满足机会成本几乎可以忽略。此时备用的边际价格主要由报价最低的机组决定。高峰时段10-14时备用价格显著上升能达到25-35元/MW。原因很直接系统接近满负荷运行要额外增加1MW备用必须让某台机组明显减少能量出力而这部分能量损失在高峰时段恰恰很贵能量价格高机会成本大幅抬升。这说明旋转备用在系统紧张时期的价值远高于低谷期。对比一下顺序出清和联合出清的结果差异我单独写了一个版本先算能量市场最优再在固定出力的基础上分配备用得到的备用价格普遍偏低。因为顺序出清时模型无法在能量和备用之间做权衡备用容量只能“捡漏”机会成本被低估最终会导致低谷时备用过度配置、高峰时备用配置不足。联合出清模型能够看到完整的机会成本出清结果经济性更好。5. 常见问题与调试经验5.1 求解器选择Gurobi、CPLEX还是内置求解器这个问题是我被问得最多的。做MILP我强烈建议用Gurobi或CPLEX如果只有Matlab基础工具箱可以用intlinprog但性能差距非常大。以我这个3机24时段模型为例Gurobi基本秒解同样的模型用intlinprog可能要跑几十秒甚至几分钟而且当机组数增加到10台以上时内置求解器的求解时间会爆炸式增长。如果你没有商业求解器授权有个折中方案是用CBCCOIN-OR Branch and CutYALMIP直接支持性能比内置的intlinprog好一些虽然还是不如Gurobi但解决教学和研究性质的小规模问题绰绰有余。5.2 数值病态问题约束尺度不一致怎么处理我在调试早期遇到过很典型的问题成本函数的常数项(c)是几百二次项系数(a)是0.01备用报价是十几负荷是上百MW。这些数值混在一起求解器的数值稳定性会受影响尤其是在定义MIPGap时容易出现“明明是最优解但报告说还有gap”的怪现象。我的处理方法是把所有成本单位统一到千元或者把负荷、备用统一换算成标幺值。比如把成本函数的所有系数除以1000得到的目标函数值大概是几十的量级数值条件好得多。这不影响最优解的正确性但能显著提升求解稳定性和收敛速度。相信我在大规模模型里你一定会碰到这个问题。5.3 冷启动问题初始机组状态怎么设定最小启停时间约束需要知道t0时刻的机组状态。如果简单假设所有机组初始都是停机状态(u_{i,0}0)那么前几个时段所有机组被最小停运时间约束锁死无法开机模型必然无解。实际上必须给一个可行的初始状态并且设置初始启停时间。我用的方法是假设所有机组初始都可以运行且已经满足最小运行时间要求。在约束代码里对(t1)的启动/停运动作单独处理比如v_start(i,1) u(i,1) - u_init(i)其中u_init是自定义的初始状态向量。这是调试最小启停约束时最容易忽略的坑很多新手在这卡很久。5.4 影子价格提取的可靠性问题前面提到MILP的直接对偶乘子不可靠。但固定整数变量后重新求解LP有时也会遇到一个微妙的问题如果整数变量固定后某些机组处于刚好满发或刚好最小出力的边界上对应的LP存在退化现象对偶乘子可能不唯一。这时不同求解器给出的影子价格可能有差异。如果发现备用价格出现不合理的“阶跃”现象可以检查一下是不是退化点导致的。解决办法是加一个微小的正则化项比如在目标函数中加上(\varepsilon \sum P_{i,t})(\varepsilon)取1e-4能有效地让求解器选出一个更合理的顶点。这个技巧在很多电力市场研究里都有人用算是半公开的经验。6. 实操经验与扩展方向做完这个旋转备用出清模型我自己最大的体会是模型复杂度不是越高越好关键是能不能回答你要研究的问题。如果只是研究备用定价机制三机系统完全足够但如果要研究节点备用价格和输电阻塞之间的关系那就必须上多节点电网模型把潮流约束加进来问题性质会完全不同。还有一个值得尝试的扩展方向是引入不确定性。旋转备用的价值本来就跟不确定性有关可以把负荷预测误差建模为随机变量用两阶段随机规划或鲁棒优化的形式来出清这时候备用需求不再是一个固定值而是由模型的置信水平内生决定。这条路做下去能写到很深的程度。另外如果对实际结算规则感兴趣可以在出清结果上加一层结算逻辑能量市场按LMP结算、备用市场按影子价格结算然后对比机组的总收益和总成本看哪些机组能盈利、哪些机组处于亏损状态。这个分析对理解市场力、容量补偿机制都有帮助也是论文里很常见的“经济性分析”部分。最后如果你打算把这个模型进一步推进到大系统建议关注以下几个方面约束的稀疏性处理、MILP的预求解presolve、以及热启动warm start策略。大规模机组组合问题求解时间往往是主要瓶颈而这些细节能帮你把求解时间从小时级降到分钟级。
旋转备用辅助服务市场出清模型：MILP建模与Matlab实现
我之前一直想把手里的备用市场出清模型整理成一篇能直接照着跑的东西正好借着这次研究旋转备用的机会把整套思路和Matlab实现都捋了一遍。这篇东西会围绕“主辅助服务市场出清模型”中的旋转备用部分展开从数学模型怎么搭、机会成本怎么算、到MILP怎么建模、影子价格怎么提取再到我实际调试时踩过的坑全都记录下来。适合正在做电力市场出清、辅助服务定价、旋转备用优化相关课题的同学参考。1. 项目概述与问题定义1.1 旋转备用辅助服务在电力市场中的角色先明确一个基本概念旋转备用Spinning Reserve指的是已经并网运行、能够在一定时间内通常是10分钟内自动响应并增加出力的机组容量。它是电力系统应对负荷波动、机组跳闸等突发情况的第一道防线在辅助服务市场里属于最关键的一类产品。主辅助服务市场出清模型做的就是“在满足系统安全约束的前提下同时确定能量市场和备用市场的成交价格与中标量”。这句话拆开理解调度机构或者说市场运营机构在组织日前市场或者实时市场时不只是把电能量卖出去就完了还需要同步把旋转备用这个“保险”配置好。因为如果只盯能量市场不看备用可能出现负荷高峰时段某些机组既想多发电、又要提供备用但装机容量有限两边打架的情况最终导致系统可靠性出问题。这个项目要解决的核心问题就是把能量和旋转备用放在同一个优化框架里联合出清。这种“联合出清”方式比传统的“先出清能量、再分配备用”顺序出清方式更合理因为它能反映机组在能量与备用两种产品之间的机会成本让备用价格真正体现其稀缺性。1.2 出清模型的核心问题既要能量又要备用做联合出清模型表面上看就是一个优化问题在满足负荷平衡、备用容量、机组运行等约束的前提下最小化系统总购电成本与备用成本。但实际建模时你会发现事情没那么简单。第一个难点是机组物理约束的多时间尺度耦合。旋转备用不是简单地在某一时刻分配多少容量它跟机组的出力水平、爬坡能力、最小运行时间都有关系。比如一台机组如果已经满发它虽然在线但实际能提供的旋转备用是0如果它处于最小技术出力附近尽管爬坡能力很强可因为出力太低能上调的空间也有限。这些约束都要在模型里体现否则出清结果无法执行。第二个难点是机会成本的处理。机组把容量留给备用市场意味着它在能量市场的收益可能减少。这个“减少的收益”正是机会成本的来源。只有联合优化才能通过模型自动权衡机组到底应该多发电还是多留备用。如果分开建模这部分成本很难准确刻画很容易导致备用价格失真。第三个难点是出清价格的计算。市场出清不只是算一个“谁中标”更重要的是算出“备用边际价格”。这个价格在数学上对应备用需求约束的对偶乘子影子价格但因为模型是混合整数规划MILP整数变量会让对偶乘子不完全可靠怎么处理就成了实操中很关键的一个环节。2. 数学模型构建全解析2.1 目标函数能量成本与备用成本的联合最小化旋转备用出清模型的目标函数我用的版本如下[ \min \sum_{t1}^{T}\sum_{i1}^{N} \left[ C_i(P_{i,t}) SU_{i,t} SD_{i,t} \rho_{i,t}^{SR} \cdot R_{i,t}^{SR} \right] ]其中(C_i(P_{i,t}))是机组i在时段t出力为(P_{i,t})时的运行成本通常用二次函数近似 [ C_i(P_{i,t}) a_i P_{i,t}^2 b_i P_{i,t} c_i ] 在MILP模型中这个二次函数需要分段线性化处理(SU_{i,t})、(SD_{i,t})分别是启动成本和停机成本(\rho_{i,t}^{SR})是机组i对旋转备用的报价元/MW(R_{i,t}^{SR})是机组i在时段t中标的旋转备用容量MW。目标函数的经济含义很直观系统总成本等于能量成本、启停成本和备用采购成本之和。注意备用成本用的是“中标量乘以报价”而不是“全部预留容量乘以报价”因为实际结算时只对被调用的部分支付。这跟能量市场按出清价格结算MCP不同属于按报价结算Pay-as-bid的机制也符合国内辅助服务市场现阶段的普遍规则。为什么要把能量成本和备用成本放在一个目标函数里这跟前面提到的机会成本有关。举一个简化例子一台最大出力100MW的机组在某个时段负荷需求是80MW如果备用需求是30MW那么这30MW备用只能由另外70MW的在线容量承担这台机组最多中标20MW备用。如果备用单独出清而不考虑能量市场的竞争这20MW的机会成本就漏掉了。联合出清时模型会同时比较“少发20MW电的损失”和“多中标20MW备用的收益”自动找到最优分配。2.2 约束条件体系不只是功率平衡把目标函数定好之后约束条件才是这个模型真正花时间的地方。我总结下来旋转备用出清模型至少需要以下几类约束功率平衡约束忽略网损时[ \sum_{i1}^{N} P_{i,t} D_t,\quad \forall t ]这里(D_t)是时段t的系统负荷。如果考虑网损可以在右侧乘一个大于1的网损因子或者在母线级模型中引入潮流计算。对于旋转备用模型来说一般先用不考虑网损的版本验证核心逻辑更合适。系统旋转备用需求约束[ \sum_{i1}^{N} R_{i,t}^{SR} \geq R_t^{req},\quad \forall t ](R_t^{req})是时段t的旋转备用总需求。这个值怎么定工程上常用两种方法最大单机容量法取系统内最大单机容量作为备用需求满足N-1准则负荷百分比加误差法比如取负荷的5%加上最大机组容量的一定比例。我测试用的系统中最大机组容量是100MW所以(R_t^{req})按最大单机容量的30%也就是30MW设置偏保守但足够观察模型的响应。机组出力上下限约束[ P_{i}^{min} u_{i,t} \leq P_{i,t} \leq P_{i}^{max} u_{i,t},\quad \forall i,t ](u_{i,t})是0-1变量表示机组i在时段t是否开机。这个约束保证机组在停机时出力为0开机时出力在最小技术出力和最大出力之间。旋转备用容量与出力耦合约束关键约束[ P_{i,t} R_{i,t}^{SR} \leq P_{i}^{max} u_{i,t},\quad \forall i,t ][ 0 \leq R_{i,t}^{SR} \leq R_{i}^{max} u_{i,t},\quad \forall i,t ]这两个约束是旋转备用建模的灵魂。第一个约束表示“机组出力加上备用不能超过最大装机容量”也就是说备用容量必须建立在机组实际留有上调空间的基础上。第二个约束限制单台机组的备用中标量防止把备用过度集中在某一台机组上。机组爬坡约束[ P_{i,t} - P_{i,t-1} \leq RU_i,\quad \forall i,t ][ P_{i,t-1} - P_{i,t} \leq RD_i,\quad \forall i,t ]这里(RU_i)、(RD_i)分别是机组的上爬坡速率和下爬坡速率。爬坡约束看着简单实际调试时最容易出问题。因为机组在提供备用的时候实际可能的爬坡能力还要打折扣——如果机组正在以最大速率爬坡它就分不出额外的爬坡能力去应对备用调用。所以严格的模型还需要在爬坡约束中叠加备用项。不过我这版模型先按“能量出力单独满足爬坡约束”处理把备用作为“静态容量”考虑这样模型规模小、容易收敛对原理验证足够了。最小启停时间约束这是机组组合问题里最“费约束”的部分。我用了两段式表达开机后至少要连续运行(UT_i)小时停机后至少要连续停运(DT_i)小时。MATLAB中可以通过对(u_{i,t})的逻辑约束来实现也可以用YALMIP内置的约束方式。考虑以上所有约束后模型变成典型的混合整数线性规划MILP问题决策变量包括连续变量(P_{i,t})、(R_{i,t}^{SR})和整数变量(u_{i,t})。求解器的任务就是在可行域内寻找使总成本最小的解。2.3 旋转备用的机会成本与价格形成机制这个部分是我觉得整个项目最值得深挖的点。为什么旋转备用的边际价格不能简单等于“最后一个中标机组的报价”因为备用资源的稀缺性不仅体现在它的直接报价上还体现在机组因为预留备用而放弃的能量市场收益。用影子价格的角度理解在最优解处旋转备用需求约束式(\sum_i R_{i,t}^{SR} \ge R_t^{req})的拉格朗日乘子(\lambda_t^{SR})就是旋转备用的边际价格。这个乘子可以分解为两部分[ \lambda_t^{SR} \text{备用报价} \text{机会成本} ]机会成本正好等于如果系统多要求1MW旋转备用最经济的方法是让某台机组在能量市场少发1MW、把这1MW容量留给备用市场那么这个“少发1MW”的损失就是机会成本。在联合优化模型中这个权衡是内生实现的不需要额外计算求解器在寻找最优解的过程中已经隐式完成了这个权衡。但这里有一个重要的实操问题MILP模型的拉格朗日乘子不唯一甚至不可靠。因为整数变量(u_{i,t})的存在破坏了线性规划的对偶理论适用条件。你直接让求解器输出约束的对偶乘子很可能得到的是次梯度subgradient而不是精确的边际价格。这个问题的标准处理方法是先用MILP求解出最优整数解(u_{i,t}^*)固定这些整数变量把模型退化为线性规划LP重新用LP求解器计算对偶变量这时得到的乘子才是市场出清价格的可靠估计。具体实现时我会在YALMIP中先求解完整MILP然后把整数变量固定为求出的值再调用一次LP求解器从结果中提取备用需求约束的dual值。这个过程在后面的代码部分会详细展示。3. Matlab代码实现从建模到求解3.1 模型数据准备建模之前先准备一个精简但能说明问题的测试系统。我用了3台机组、24小时负荷曲线参数如下机组最大出力(MW)最小出力(MW)爬坡速率(MW/h)最小运行时间(h)最小停运时间(h)备用报价(元/MW)最大备用(MW)G11002050441230G2801540331525G3601030222520成本函数系数二次形式(aP^2bPc)机组a(元/MW²h)b(元/MWh)c(元/h)启动成本(元)G10.01224180600G20.01522150400G30.01820120300负荷曲线按典型的日负荷特征设置凌晨低谷时段1-6时在120MW附近白天高峰时段10-14时爬升到200MW左右晚间有一个小高峰。24小时负荷的具体数列我放在代码里。这里有一个重要提醒成本函数的二次项系数不能设得太小否则各机组之间成本差异不明显模型会在大范围内出现多个近似最优解求解器反而不稳定。我调试时有一版二次项系数是0.000几结果出清结果在相邻时段出现反复启停的“振荡”现象后来把系数调大一个数量级才稳定下来。3.2 YALMIP建模与MILP求解Matlab实现我采用的是YALMIP工具箱。选择YALMIP的原因建模语法简洁约束可以直接用大于等于/小于等于表达式写不需要手动拼矩阵切换求解器方便同一套代码可以试Gurobi、CPLEX、或者开源的CBC对MILP和LP的混合求解支持良好方便做影子价格的两阶段计算。当然如果你不想装YALMIP直接用intlinprog也行但手动把所有约束写成矩阵形式非常容易出错尤其是机组组合这种带时间耦合的问题矩阵维度一旦搞错排查起来极费时间。我还是推荐用YALMIP做原型验证。备用需求和负荷曲线会直接写成列向量。下面是我实际使用的核心建模代码。3.3 核心代码实现与解析%% 主辅助服务市场出清模型旋转备用MILP % 联合优化能量市场与旋转备用市场 clear; clc; close all; yalmip(clear); %% 1. 基础数据定义 T 24; % 时段数 N 3; % 机组数 % 机组参数 Pmax [100; 80; 60]; % 最大出力 MW Pmin [20; 15; 10]; % 最小出力 MW RU [50; 40; 30]; % 上爬坡 MW/h RD [50; 40; 30]; % 下爬坡 MW/h UT [4; 3; 2]; % 最小运行时间 h DT [4; 3; 2]; % 最小停运时间 h Suc [600; 400; 300]; % 启动成本 元 Sdc [300; 200; 150]; % 停机成本 元 % 成本函数系数二次 a [0.012; 0.015; 0.018]; % 二次项 b [24; 22; 20]; % 一次项 c [180; 150; 120]; % 常数项 % 备用参数 RhoSR [12; 15; 25]; % 备用报价 元/MW Rmax [30; 25; 20]; % 单台机组最大备用 MW % 负荷曲线24小时 D [120 115 110 105 108 125 ... 150 170 185 195 200 205 ... 195 190 185 180 175 170 ... 165 160 155 150 145 135]; % 备用需求按最大单机容量的30%再加5%负荷波动 Rreq 30 0.05 * D; %% 2. 决策变量 P sdpvar(N, T, full); % 出力 MW R sdpvar(N, T, full); % 旋转备用 MW u binvar(N, T, full); % 机组启停状态 0/1 v_start binvar(N, T, full); % 启动动作 v_shut binvar(N, T, full); % 停机动作 %% 3. 目标函数能量成本 启停成本 备用成本 Objective 0; for t 1:T for i 1:N % 分段线性化成本这里用二次函数直接通过sdpvar的quadratic支持 % 注意标准MILP需要线性化YALMIP支持自动处理二次函数 Objective Objective a(i) * P(i,t)^2 b(i) * P(i,t) c(i) * u(i,t); Objective Objective Suc(i) * v_start(i,t) Sdc(i) * v_shut(i,t); Objective Objective RhoSR(i) * R(i,t); end end %% 4. 约束条件 Constraints []; % 功率平衡约束 for t 1:T Constraints [Constraints, sum(P(:,t)) D(t)]; end % 旋转备用需求约束 for t 1:T Constraints [Constraints, sum(R(:,t)) Rreq(t)]; end % 出力上下限约束 for t 1:T for i 1:N Constraints [Constraints, Pmin(i) * u(i,t) P(i,t)]; Constraints [Constraints, P(i,t) Pmax(i) * u(i,t)]; end end % 出力备用不超过最大容量旋转备用关键耦合约束 for t 1:T for i 1:N Constraints [Constraints, P(i,t) R(i,t) Pmax(i) * u(i,t)]; Constraints [Constraints, 0 R(i,t) Rmax(i) * u(i,t)]; end end % 爬坡约束 for t 2:T for i 1:N Constraints [Constraints, P(i,t) - P(i,t-1) RU(i)]; Constraints [Constraints, P(i,t-1) - P(i,t) RD(i)]; end end % 最小运行/停运时间约束 for i 1:N for t 1:T % 最小运行时间 if t UT(i) - 1 T Constraints [Constraints, ... u(i,t) - u(i,t-1) u(i,t1) u(i,t2) ... u(i,tUT(i)-1)]; end % 最小停运时间 if t DT(i) - 1 T Constraints [Constraints, ... u(i,t-1) - u(i,t) (1-u(i,t1)) (1-u(i,t2)) ... (1-u(i,tDT(i)-1))]; end end end % 启动/停运动作与被机组状态的关系 for t 1:T if t 1 Constraints [Constraints, v_start(:,t) u(:,t) - 0]; Constraints [Constraints, v_shut(:,t) 0 - u(:,t)]; else Constraints [Constraints, v_start(:,t) u(:,t) - u(:,t-1)]; Constraints [Constraints, v_shut(:,t) u(:,t-1) - u(:,t)]; end Constraints [Constraints, 0 v_start(:,t) 1]; Constraints [Constraints, 0 v_shut(:,t) 1]; Constraints [Constraints, u(:,t) 1]; Constraints [Constraints, u(:,t) 0]; end %% 5. 求解MILP ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 0.0001; % 相对MIP间隙 ops.gurobi.TimeLimit 300; % 秒 result optimize(Constraints, Objective, ops); if result.problem ~ 0 error(求解失败: %s, result.info); end %% 6. 提取并计算影子价格固定整数变量后重新求解LP u_fixed value(u); P_opt value(P); R_opt value(R); % 固定整数变量 for t 1:T for i 1:N Constraints [Constraints, u(i,t) u_fixed(i,t)]; end end % 重新求解LP ops2 sdpsettings(solver, gurobi, verbose, 0); result2 optimize(Constraints, Objective, ops2); if result2.problem ~ 0 error(LP重求解失败); end % 提取备用需求约束的影子价格 dual_SR zeros(T, 1); for t 1:T % 注意YALMIP中duals是最后一步求解的对偶值 dual_SR(t) dual(Constraints(find(is(Constraints,inequality) ... contains(string(Constraints),[sum(R(:,, num2str(t), ))])))); end % 或者更稳妥的方法直接用shadowprice函数或手动指定约束索引 % 下面是用索引方式初始化Constraints时给关键约束加tag更可靠代码写到这里我不得不提醒一个容易踩的坑YALMIP里提取对偶变量的方式比较“诡异”。上面的dual()函数需要你跟踪约束在Constraints对象中的索引。最可靠的做法是在建模时就给备用需求约束单独维护一个cell数组而不是用一个大Constraints对象。我在下面给出更干净的写法。% 备用需求约束单独存储 Constraints_SR []; for t 1:T Constraints_SR [Constraints_SR, sum(R(:,t)) Rreq(t)]; end % 在总约束集合中加入 Constraints [Constraints, Constraints_SR]; % 求解后提取影子价格 dual_SR dual(Constraints_SR);这样提取的dual_SR就是24个时段的旋转备用边际价格序列。3.4 成本分段线性化处理前面提到机组成本函数是二次的标准的MILP求解器无法直接处理非线性目标函数。YALMIP有一个比较聪明的做法它会在内部自动把二次函数做切分线性近似或者调用支持二次目标的求解器比如Gurobi可以直接处理凸二次目标MIQP。但如果你的求解器不支持MIQP就需要手动做分段线性化。分段线性化的思路是把出力区间([P_{min}, P_{max}])切分成K段每段用一个线性函数逼近。切分点越多越精确但变量和约束也越多。我测试系统用的是二次函数且只有3台机组所以直接交给Gurobi处理MIQP是没问题的。但如果你做的是大系统、上百台机组建议还是先做分段线性化转成MILP求解速度和稳定性都好很多。4. 仿真结果与分析4.1 出清结果展示我用上面的代码跑完24小时出清先看机组出力情况。凌晨低谷时段1-6时负荷在110-125MW之间G1作为成本最优的大机组基本承担主要出力G2和G3部分时段停机或处于最小出力状态。白天负荷爬升后G1满发G2和G3跟进补充。这个结果符合预期因为G1的边际成本最低优先发电。再看旋转备用中标量所有时段总备用都满足30MW以上的系统需求。有意思的是备用分配结构低谷时段G1虽然出力不高但有充足的上调空间所以拿了大部分备用高峰时段G1满发后备用就主要由G2和G3承担了。这说明耦合约束P R Pmax起到了应有的作用防止把备用分配给一台已经满发的机组。下面的表格是24个时段的部分出清结果摘要为节省篇幅只列了7时、12时、20时三个典型时段时段负荷(MW)备用需求(MW)G1出力(MW)G1备用(MW)G2出力(MW)G2备用(MW)G3出力(MW)G3备用(MW)715037.5100035201517.51220540.251000804.2525152015537.75100040251512.75注意12时这个时段G1已经满发备用为0G2也接近满发80MW只有4.25MW备用G3承担了剩余的15MW备用。一切都被耦合约束限制住了模型没有给出任何一台机组“既满发又提供备用”的物理上不可行的结果。4.2 边际电价与备用价格解读重头戏是备用价格的提取。两阶段求解后我得到了24个时段的旋转备用影子价格曲线。有几个典型的特征低谷时段1-6时备用价格很低基本在12-15元/MW附近。因为此时系统备用充裕多1MW备用需求只需要让某台机组在能量市场少发一点点电即可满足机会成本几乎可以忽略。此时备用的边际价格主要由报价最低的机组决定。高峰时段10-14时备用价格显著上升能达到25-35元/MW。原因很直接系统接近满负荷运行要额外增加1MW备用必须让某台机组明显减少能量出力而这部分能量损失在高峰时段恰恰很贵能量价格高机会成本大幅抬升。这说明旋转备用在系统紧张时期的价值远高于低谷期。对比一下顺序出清和联合出清的结果差异我单独写了一个版本先算能量市场最优再在固定出力的基础上分配备用得到的备用价格普遍偏低。因为顺序出清时模型无法在能量和备用之间做权衡备用容量只能“捡漏”机会成本被低估最终会导致低谷时备用过度配置、高峰时备用配置不足。联合出清模型能够看到完整的机会成本出清结果经济性更好。5. 常见问题与调试经验5.1 求解器选择Gurobi、CPLEX还是内置求解器这个问题是我被问得最多的。做MILP我强烈建议用Gurobi或CPLEX如果只有Matlab基础工具箱可以用intlinprog但性能差距非常大。以我这个3机24时段模型为例Gurobi基本秒解同样的模型用intlinprog可能要跑几十秒甚至几分钟而且当机组数增加到10台以上时内置求解器的求解时间会爆炸式增长。如果你没有商业求解器授权有个折中方案是用CBCCOIN-OR Branch and CutYALMIP直接支持性能比内置的intlinprog好一些虽然还是不如Gurobi但解决教学和研究性质的小规模问题绰绰有余。5.2 数值病态问题约束尺度不一致怎么处理我在调试早期遇到过很典型的问题成本函数的常数项(c)是几百二次项系数(a)是0.01备用报价是十几负荷是上百MW。这些数值混在一起求解器的数值稳定性会受影响尤其是在定义MIPGap时容易出现“明明是最优解但报告说还有gap”的怪现象。我的处理方法是把所有成本单位统一到千元或者把负荷、备用统一换算成标幺值。比如把成本函数的所有系数除以1000得到的目标函数值大概是几十的量级数值条件好得多。这不影响最优解的正确性但能显著提升求解稳定性和收敛速度。相信我在大规模模型里你一定会碰到这个问题。5.3 冷启动问题初始机组状态怎么设定最小启停时间约束需要知道t0时刻的机组状态。如果简单假设所有机组初始都是停机状态(u_{i,0}0)那么前几个时段所有机组被最小停运时间约束锁死无法开机模型必然无解。实际上必须给一个可行的初始状态并且设置初始启停时间。我用的方法是假设所有机组初始都可以运行且已经满足最小运行时间要求。在约束代码里对(t1)的启动/停运动作单独处理比如v_start(i,1) u(i,1) - u_init(i)其中u_init是自定义的初始状态向量。这是调试最小启停约束时最容易忽略的坑很多新手在这卡很久。5.4 影子价格提取的可靠性问题前面提到MILP的直接对偶乘子不可靠。但固定整数变量后重新求解LP有时也会遇到一个微妙的问题如果整数变量固定后某些机组处于刚好满发或刚好最小出力的边界上对应的LP存在退化现象对偶乘子可能不唯一。这时不同求解器给出的影子价格可能有差异。如果发现备用价格出现不合理的“阶跃”现象可以检查一下是不是退化点导致的。解决办法是加一个微小的正则化项比如在目标函数中加上(\varepsilon \sum P_{i,t})(\varepsilon)取1e-4能有效地让求解器选出一个更合理的顶点。这个技巧在很多电力市场研究里都有人用算是半公开的经验。6. 实操经验与扩展方向做完这个旋转备用出清模型我自己最大的体会是模型复杂度不是越高越好关键是能不能回答你要研究的问题。如果只是研究备用定价机制三机系统完全足够但如果要研究节点备用价格和输电阻塞之间的关系那就必须上多节点电网模型把潮流约束加进来问题性质会完全不同。还有一个值得尝试的扩展方向是引入不确定性。旋转备用的价值本来就跟不确定性有关可以把负荷预测误差建模为随机变量用两阶段随机规划或鲁棒优化的形式来出清这时候备用需求不再是一个固定值而是由模型的置信水平内生决定。这条路做下去能写到很深的程度。另外如果对实际结算规则感兴趣可以在出清结果上加一层结算逻辑能量市场按LMP结算、备用市场按影子价格结算然后对比机组的总收益和总成本看哪些机组能盈利、哪些机组处于亏损状态。这个分析对理解市场力、容量补偿机制都有帮助也是论文里很常见的“经济性分析”部分。最后如果你打算把这个模型进一步推进到大系统建议关注以下几个方面约束的稀疏性处理、MILP的预求解presolve、以及热启动warm start策略。大规模机组组合问题求解时间往往是主要瓶颈而这些细节能帮你把求解时间从小时级降到分钟级。