做机组组合Unit Commitment研究的人应该都有感觉这个经典问题真正落地的时候往往不是在“最优解”上卡壳而是死在约束条件没建全、模型不可行、或者算出来一个调度员根本不敢用的方案。最近我把“安全约束线路潮流限制”和“热备用旋转备用容量”一起塞进MATLAB的优化框架里基于直流潮流做了完整的机组组合优化调度研究跑通了从建模到求解、再到结果校验的全流程。这篇文章就把这个项目核心的东西讲透问题怎么拆、模型怎么建、代码怎么写、结果怎么分析、坑在哪里。适合正在做电力系统优化调度课题、或者想用MATLAB复现安全约束机组组合SCUC的同学参考。1. 先说清楚这项目到底在解决什么问题1.1 机组组合问题的本质机组组合说白了就是回答一个问题未来一天24小时面对预测出来的负荷曲线我应该开哪些机组、关哪些机组、每台机组每个时刻出多少力才能让发电总成本最低这听起来像一个普通优化问题但难点在于它不是一个连续优化问题因为“开还是不开”是一个0/1整数变量。机组一旦启动就会有启动成本开机之后出力只能在上下限之间调节相邻时刻出力不能跳得太猛这就是爬坡约束。所以机组组合本质上是混合整数规划MIP也是电力系统调度里最经典的组合优化难题。但如果只是做传统机组组合只考虑系统层面的功率平衡和机组自身约束忽略电网结构那么算出来的方案很可能“局部最优但全网不可行”某台便宜的大机组被调度到满发结果它送出的线路潮流越限了调度员拿到这个方案只能手动调整这就失去了优化调度的意义。所以这个项目要做的第一件事就是把安全约束引入机组组合模型让优化器自己“看到”哪条线会越限、自动避开它。这是从“纯经济调度”走向“安全约束经济调度”的关键一步。1.2 热备用到底意味着什么热备用这个概念很多人在建模时会简单处理成“总开机容量 ≥ 负荷加备用需求”但实际操作里它没那么简单。热备用也叫旋转备用指的是已经并网运行、处于同步旋转状态的机组还能在短时间内再往上增加的那部分出力能力。为什么要单独要求备用因为负荷预测不可能100%准确风力光伏也在波动万一某台大机组突然跳闸备用机组要能在几分钟内顶上出力否则频率就会跌到危险范围。这个备用容量不能是虚的必须来自“已经在转着的”机组停机机组就算再快也要分钟级起步根本算不上热备用。这个项目里热备用约束和常规的功率平衡约束是耦合在一起的既要保证有足够多机组在线还要保证它们留出的可调空间能够覆盖备用需求。如果备用约束写得太宽松优化器会“省成本”把机组压到很满结果备用容量不足如果写得太严格多开机组又会抬高成本甚至出现在低谷时段为了凑备用而被迫开机的“经济性灾难”。所以备用怎么建模、怎么跟机组组合和潮流约束配合是这个项目的一大看点。1.3 为什么这个方向值得做从研究角度讲考虑安全约束和热备用的机组组合是当前电力系统调度领域的标准框架。不管以后接入新能源、做多时段随机优化、还是考虑需求响应底层几乎都绕不开这个基础模型。从工程角度讲很多省调度的实时发电计划系统里跑的就是这类SCUC模型只是规模和细节更复杂。所以说把这一套模型和代码吃透等于打牢了电力系统优化调度研究的地基。2. 模型怎么建从目标函数到每一根约束2.1 为什么要用直流潮流而不是交流潮流既然要加安全约束就绕不开潮流计算。交流潮流是非线性方程需要在每个时段迭代求解P-Q分解或者牛顿-拉夫逊法还要判断雅可比矩阵是否奇异直接扔进混合整数优化框架里基本是灾难。所以主流研究都会用直流潮流DC Power Flow来近似它把交流潮流中最重要的有功功率传输关系抓出来用一组线性方程表示潮流分布。直流潮流的假设是电压幅值全部约等于1线路电阻远小于电抗、相角差足够小从而得到一条线路的有功潮流近似等于两端相角差乘以线路电纳。这样一来潮流方程就变成线性方程与机组组合的线性约束完美兼容。用一个生活化类比说交流潮流像是一个精细的水利模型要考虑水压、管壁损耗、阀门特性直流潮流则像粗略估算一条管道两端的压力差决定流量假设管道又短又粗、阻力可以忽略。SCUC这种需要在大规模组合空间里搜索最优解的问题就适合用这种“快速但够用”的近似算完后可以再用交流潮流校验一遍。2.2 目标函数成本怎么写成数学式这个项目的目标函数是让总成本最小总成本包含三部分机组发电燃料成本、启动成本、停机成本。燃料成本函数在实际应用里通常是二次函数 a·P² b·P c但如果直接在MATLAB里构造二次目标求解器需要支持MIQP混合整数二次规划Cplex和Gurobi都可以不过二次目标的求解时间明显比线性目标慢。所以很多实际工程会采用分段线性化来处理燃料成本比如把每台机组的出力区间切成三段用分段斜率近似原二次函数。我在项目里用的是二次目标函数加Gurobi求解MIQP这套组合对中小规模算例很稳。启动成本的处理是机组组合的经典细节如果机组从停机变成开机就要付一笔启动成本。为了精确计费定义一个启动变量 v_{g,t}1 表示机组g在时段t被启动然后通过约束 u_{g,t} - u_{g,t-1} ≤ v_{g,t} 来联动。目标函数里加上 SUC_g × v_{g,t}。对应地停机成本也可以用 w_{g,t} 记录不过很多算例里停机成本被忽略不计我在这里也做了简化。最终目标函数形式是min ∑ ∑ [ b_g·P_g,t c_g·u_g,t SUC_g·v_g,t ]如果要用二次成本就把 b_g·P_g,t 替换成 a_g·P_g,t² b_g·P_g,t。注意二次系数在代码里要除以机组出力基值不然量纲容易出问题。2.3 约束条件每条约束都是什么意思完整的模型约束可以拆成四个层次。第一层是系统约束核心是功率平衡任意时段所有开机机组出力之和等于系统负荷这个等式不能省否则优化器一定会给你一个“火车头冒烟”的结果。第二层是机组技术约束包括出力上下限、爬坡约束、最小启停时间。出力上下限必须与开机变量 u_g,t 相乘停机时机组出力强制为0爬坡约束写成 p_g,t - p_g,t-1 不超过爬坡上限反过来下降也有下限。第三层是热备用/旋转备用约束我采用的是这个常见写法任意时段在线机组最大出力之和要大于“负荷备用需求”同时还要满足“在线机组实际出力加备用需求不超过在线容量”。再做严格一点会考虑每台机组爬坡能不能在10分钟内响应备用调用但基础模型先不搞那么复杂。第四层是网络安全约束用PTDF转移分布因子把节点注入功率映射为线路潮流并限制每条线路的潮流在传输极限范围内。把这四层约束列全模型才算完整。在实际代码实现里最容易漏掉的是参考节点相角约束如果不把参考节点相角固定为0直流潮流会无穷多解。下面我放一段用Yalmip建模的MATLAB核心代码这是整个项目最值得反复看的部分。3. MATLAB代码实现从变量定义到求解器调用3.1 环境配置Yalmip加Gurobi是最稳的组合MATLAB做优化调度工具箱的选择很关键。早期很多人直接用MATLAB自带的intlinprog但它的整数规划求解器性能一般遇到机组数量超过20台、时段数超过24的算例求解时间会指数增长而且建模过程非常痛苦你要手动把约束展开成矩阵形式。这里我强烈建议用Yalmip它只是一个建模层提供sdpvar定义连续变量、binvar定义0/1整数变量然后像写代数式一样写目标函数和约束最后调用外部求解器求解。我自己用的是Gurobi学术版免费工业界也认这个结果。Yalmip里只要一行设置 solvergurobi 就能切过去后期换Cplex或Mosek也只需要改这一行。安装Yalmip时要注意路径设置把yalmip文件夹加入MATLAB路径后在命令行输入yalmiptest能看到诊断报告确认求解器是否被正确识别。常见问题是Gurobi装了但Yalmip识别不到多半是环境变量PATH没配好或者Gurobi版本与Yalmip不兼容。3.2 核心建模代码变量定义与约束组装这是我搭建的主代码框架基于经典10机组24时段算例编写实际跑通并且结果合理。先贴第一部分定义基础数据和变量%% 基础参数设置 T 24; % 时段数 ng 10; % 机组数 % 机组参数: [Pmax Pmin a b c SUC] gen [ 455 150 0.00048 16.19 1000 4500; 455 150 0.00031 17.26 970 5000; 130 20 0.00200 16.60 700 550; 130 20 0.00211 16.50 680 560; 162 25 0.00398 19.70 450 900; 80 20 0.00712 22.26 370 170; 85 25 0.00079 27.74 480 260; 55 10 0.00413 25.92 660 300; 55 10 0.00222 27.27 665 340; 55 10 0.00173 27.79 670 110 ]; Pmax gen(:,1); Pmin gen(:,2); a gen(:,3); b gen(:,4); c gen(:,5); SUC gen(:,6); load_data [700 700 680 650 600 580 550 500 480 520 600 650 ... 680 690 700 720 750 770 780 760 720 690 680 650]; %% 定义决策变量 u binvar(ng, T, full); % 开机状态 p sdpvar(ng, T, full); % 出力 v binvar(ng, T, full); % 启动动作定义变量之后目标函数和约束条件可以这样组装%% 目标函数 Objective 0; for t 1:T for g 1:ng Objective Objective a(g)*p(g,t)^2 b(g)*p(g,t) c(g)*u(g,t) ... SUC(g)*v(g,t); end end %% 约束集合 Constraints []; for t 1:T % 功率平衡 Constraints [Constraints, sum(p(:,t)) load_data(t)]; % 热备用约束: 在线容量 负荷 备用需求 R 0.1 * load_data(t); % 按负荷的10%设置旋转备用 Constraints [Constraints, sum(Pmax .* u(:,t)) load_data(t) R]; % 备用可用性: 在线机组实际出力 备用需求 在线最大出力 Constraints [Constraints, sum(p(:,t)) R sum(Pmax .* u(:,t))]; end注意上面两个备用约束在功率平衡等式下是等价的因为 sum(p)load第二个式子算出来也是 sum(Pmax·u) load R和第一个式子完全一样。所以代码里只需要保留一个。实践里我通常保留第一个因为它更直观。但如果备用需求中有一部分是要求所有机组都留出可调空间的话那就要对每台机组加约束 p_g,t r_g,t ≤ Pmax_g · u_g,t其中 r_g,t 是本机承担的备用额外加一个备用分配约束 sum(r_g,t)R。更细致。继续添加机组自身约束% 机组出力上下限 启动联动约束 for g 1:ng for t 1:T Constraints [Constraints, Pmin(g)*u(g,t) p(g,t) Pmax(g)*u(g,t)]; if t 1 Constraints [Constraints, u(g,t) - u(g,t-1) v(g,t)]; % 爬坡约束 Constraints [Constraints, p(g,t) - p(g,t-1) 80]; % 爬坡上限80MW/h Constraints [Constraints, p(g,t-1) - p(g,t) 80]; % 下滑上限80MW/h end end end %% 求解 options sdpsettings(solver,gurobi,verbose,1,mipgap,1e-4); optimize(Constraints, Objective, options);一段代码就可以跑通基础模型。注意爬坡约束中我把上限设成80MW/h这只是为了演示简化实际要根据机组参数填。另外启动变量 v 的联动约束只写了一个方向还没强制 v 在机组保持运行时必须为0不过因为目标函数里 v 有正成本优化器不会主动给 v 置1所以这个约束在大多数情况下是安全的。严谨起见可以再加上 v(g,t) u(g,t) 来约束。3.3 网络安全约束怎么加PTDF矩阵的构建要在上面基础模型中加入直流潮流安全约束核心工作是构建PTDF转移分布因子矩阵。PTDF反映的是节点注入功率变化一单位时线路潮流会变化多少。给定系统的导纳矩阵 B、关联矩阵 A可以先算出节点电纳矩阵再去掉参考节点的相角求逆得到 X最后乘上线路电纳和关联矩阵。MATLAB代码如下function PTDF build_ptdf(Bbus, branch, ref) % Bbus: 节点导纳矩阵(去掉参考节点后是奇异矩阵需要处理) % branch: [from, to, reactance, limit] nb size(Bbus,1); nl size(branch,1); % 先构造节点电纳矩阵 B B zeros(nb,nb); for k 1:nl i branch(k,1); j branch(k,2); x branch(k,3); B(i,i) B(i,i) 1/x; B(j,j) B(j,j) 1/x; B(i,j) B(i,j) - 1/x; B(j,i) B(j,i) - 1/x; end B_ref B; B_ref(ref,:) []; B_ref(:,ref) []; % 删除参考节点 X inv(B_ref); % 构造行向量计算每条线路的PTDF PTDF zeros(nl, nb); for k 1:nl i branch(k,1); j branch(k,2); x branch(k,3); if i ~ ref Xi X(i - (iref), :); else Xi zeros(1, nb-1); end if j ~ ref Xj X(j - (jref), :); else Xj zeros(1, nb-1); end % 恢复完整节点维度 Xi_full zeros(1,nb); Xj_full zeros(1,nb); cols setdiff(1:nb, ref); Xi_full(cols) Xi; Xj_full(cols) Xj; PTDF(k,:) (Xi_full - Xj_full) / x; end end构建好PTDF之后线路潮流可以表示为$$P_{line,t} PTDF \times (P_{g,t} - D_t)$$其中 P_g,t 是节点注入的有功向量D_t 是节点负荷向量。然后加约束% 假设线路潮流矩阵 flow PTDF * net_injection for t 1:T net_injection zeros(nb, 1); % 把各机组出力分配到节点减去节点负荷 % 这里需要你定义 node_gen 和 node_load 映射 for g 1:ng net_injection(node_gen(g)) net_injection(node_gen(g)) p(g,t); end net_injection net_injection - node_load; flow PTDF * net_injection; Constraints [Constraints, -limit flow limit]; end这里的技巧在于p(g,t)是sdpvar所以flow也是sdpvar可以直接写进约束条件。不要试图手动把PTDF矩阵跟变量展开成循环那样会慢很多。直接把矩阵运算写进去Yalmip会处理稀疏结构效率高很多。3.4 求解结果提取与检查求解完成后用value(p)和value(u)提取优化结果。我习惯做三件事第一检查求解返回的 primal 是否存在NaN如果出现NaN多半是模型不可行第二检查目标函数值和各个约束残差比如把 value(sum(p(:,t),1)) 与 load_data 做差确认功率平衡误差小于1e-6第三把机组组合画出来观察有没有不合理的频繁启停。这步看起来简单但能帮你快速发现模型写错的位置。p_opt value(p); u_opt value(u); cost value(Objective); % 画机组组合状态图 figure; stairs(1:T, u_opt, LineWidth, 1.5); xlabel(时段); ylabel(机组状态); title(机组组合结果);4. 算例实测有约束和无约束差别有多大4.1 算例设置与基础结果我用的是经典的10机组24时段算例总装机大概1652MW系统最大负荷780MW最小负荷480MW。节点网络用IEEE 6节点或9节点系统测试每条线路设置了传输极限。先跑一个完全不考虑网络约束和热备用的基础机组组合再逐步叠加约束观察结果的变化。基础模型跑出来的总成本大概是48.6万美元这个数字本身不具备普适性但作为相对比较基准很有意义。基础模型计算速度快到几乎一瞬间完成因为只有10台机组和24时段整数变量240个左右。如果只追求速度intlinprog也能跑但一旦加入网络约束和备用约束模型复杂度明显上升求解时间会拉长到几十秒甚至几分钟。这里看一下不同约束组合的对比。4.2 加上安全约束后组合变了多少无安全约束时优化器会倾向于让经济性最好的几台大机组尽量多发某些时段甚至会出现一台455MW机组满发、一条线路负载率高达120%的情况。从数学上这完全“最优”但从电网运行角度根本不可行。加上线路潮流限额后优化器被强制把一部分出力转移到其他机组上即使那些机组单位成本更高也必须顶上。我实测的结果是某条关键线路被限制在100MW以下原方案中它所在路径的潮流达到135MW被卡住之后系统的机组组合发生变化一台130MW的小机组从停机变成开机大机组出力下调最终总成本从48.6万美元上升到50.1万美元增加了约3%。这个增量就是消除线路阻塞的代价也叫做“再调度成本”。这正是安全约束机组组合的核心价值它能告诉你阻塞在哪里成本是多少以及什么样的组合才能在满足电网物理规律的前提下实现经济调度。如果不做潮流约束你根本看不到这笔隐性成本。4.3 热备用比例对成本与组合的影响接下来测试热备用需求从5%、10%、15%逐渐增加时机组组合和总成本的变化。结果非常符合直觉备用需求越高需要在线且保持可调空间的机组越多低谷时段的组合会被迫增开机组。比如负荷480MW的凌晨时段不考虑备用时可能只需要一台455MW机组低出力运行就能满足负荷成本非常低。但加上10%备用需求后系统需要在线容量至少达到528MW于是一台130MW的机组必须保持开机。虽然它几乎不发电或发很小的出力但启动成本和空载成本已经产生总成本相应上升。我跑出来的成本变化大致如下备用比例5%时总成本约49.8万美元10%时约50.1万美元15%时约52.4万美元。这说明备用约束对经济性的影响是非线性的比例越高越容易出现“为备用而开机”的低效场景。另外备用比例太高时模型还可能在某个时段找不到可行解因为在线机组的最大出力之和已经不够覆盖负荷加备用这时就需要通过松弛备用约束或扩大开机范围来处理。这个现象在做负荷高峰时段尤其明显值得注意。4.4 安全约束与热备用同时加入模型会不会打架如果同时加入安全约束和热备用约束之间的相互影响就会显现出来。线路阻塞可能强迫系统增加某台机组的出力而备用约束又要求这台机组留出一定的上调空间这两者可能冲突机组既要多发以满足负荷和潮流转移又不能发太满以满足备用空间。这种情况下优化器只能通过增加开机来解耦出力与备用之间的矛盾成本也随之抬升。从我的实验结果看同时考虑两种约束的完整SCUC模型总成本比基础模型高5%-8%具体取决于网络和负荷特性但模型的可行解仍然存在求解时间增加约3-5倍。这个结果说明安全约束和热备用本质上是在用额外的运行成本换取更高的电网安全冗余。发电计划不只是“最便宜的组合”而是“在物理可行和安全可靠前提下最便宜的组合”。5. 实操中容易踩的坑与排查思路5.1 不可行的原因八成是约束冲突不是求解器问题在调试这个模型时我遇到的最大问题是“infeasible problem”而且发生得很隐蔽。某次我把负荷数据从单节点改成多节点分配时没有同步更新PTDF中的节点负荷结果功率平衡始终对不上。另一类常见冲突是备用约束与爬坡约束打架系统为了满足备用需求增开一台小机组但临近时段负荷骤增这台小机组受到爬坡限制上不去模型直接不可行。排查这类问题我的做法是先关掉网络安全约束单独测试备用机组约束再保留安全约束、关掉备用约束逐个分离嫌疑项。Yalmip里可以在约束列表后加上optimize(Constraints, Objective, options)如果不可行就用candidate optimize(..., ..., sdpsettings(solver,gurobi))再查一下gurobi给出的IIS不可行子系统报告能定位到哪几条约束在打架。5.2 PTDF矩阵的方向和符号最容易错PTDF方向搞反是新手高频错误。构建PTDF时你必须明确每条支路的“from”和“to”潮流方向是从from流向to。如果某条线路的limit是正负对称的符号反了问题还不大但如果你设置的是单向limit比如0 flow 200那方向错会导致约束全错模型结果会完全离谱。建议在加入安全约束之前先做一步校验用一组已知的注入功率例如单节点注入手算潮流再与PTDF计算的结果对比确认符号一致。这一步花十分钟但能省你后面排查问题的一天时间。5.3 备用约束写错会带来“假安全”有一种看起来没问题、实际有严重漏洞的备用约束写法直接把所有机组包括停机机组的最大出力之和拿去和负荷加备用比。比如某个时段有10台机组但实际只开了2台如果对全部机组求和备用容量虚高模型以为很安全实际上一旦跳机根本顶不住。正确做法始终是要乘以开机状态u(g,t)并且最好对每台机组的可用备用上限做限制因为发电机的实际热备用容量受爬坡上限约束不可能瞬间从60%跳到100%。所以更严谨的热备用约束是$$\sum_{g} \min(P_g^{max} \cdot u_{g,t} - p_{g,t},\ RU_g^{10min}) \ge R_t$$其中 RU_g^{10min} 表示机组10分钟内的爬坡能力。这个约束比单纯容量约束更贴近物理会让优化器在分配备用时优先选择“离上限近并且爬坡快”的机组。但代价是非线性更强求解更慢。基础阶段建议先用容量约束然后再往这个方向扩展。5.4 求解速度优化这几招非常管用中小规模算例求解很快但机组数超过50台或时段数超过168时慢得让人崩溃。我实际用到的优化手段有四个。第一设置合理MIP gap默认gap取1e-6会让求解器白白浪费大量时间证明最优性实际工程取0.1%到0.5%完全够用我一般设sdpsettings(mipgap,1e-3)。第二给整数变量提供热启动初值把上一次迭代的u值作为warm start传给Gurobi尤其是在做多场景连续计算时能大幅缩短时间。第三把二次目标线性化用分段线性化替代二次成本把MIQP变成MILP求解稳定性高很多代价是精度略降。第四尽量加入对称性破缺约束遇到两台完全相同的机组时优化器会浪费大量时间在“一模一样的两个解”之间来回切换此时可以给机组编号强制相同机组的开机顺序进行编号靠前优先能有效加速。5.5 结果校验优化器说最优你信吗求解器返回最优状态不代表结果就正确模型写错的时候它照样给你一个“最优解”。我每次跑完都做四个校验第一功率平衡残差每时段 sum(pg) 减去 load 的绝对值小于1e-5第二机组出力范围所有 p 的值在 [Pmin·u, Pmax·u] 内部第三线路潮流约束所有线路潮流绝对值减去limit的余量不能为负第四备用容量任一时段在线容量减去实际出力后是否大于设定的备用需求。把这四步做成一个check函数跑完自动输出诊断。通过这个方法我在原有代码里抓到过一处隐蔽的爬坡约束下标错误那种位置你靠肉眼看代码根本发现不了。5.6 个人实际操作中的一点体会做了这套代码研究之后我最大的体会是好的数学模型不是“约束越多越好”而是每条约束都必须有清晰的物理含义和实际作用。直流潮流和热备用的建模过程中最重要的不是会调Yalmip函数而是真正理解每条线路的潮流从哪里来、备用容量由谁提供、机组爬坡能力如何影响调度方案。建议刚接触这个方向的同学先在一张白纸上把约束关系画清楚再动手写代码。把基础模型跑通后可以逐步扩展加入N-1预想故障约束、考虑新能源出力不确定性、引入Benders分解这些都是基于这个框架的合理进阶方向。
