最近有朋友找我帮忙做电力系统日前调度方案核心诉求是“考虑安全约束和热备用的机组组合优化调度”并且指定要MATLAB实现。国内公开资料里讲机组组合Unit Commitment, UC的代码不少但很多只是简单的经济调度要么忽略网络约束要么把备用简单设成一个固定值不够工程化。这篇文章把我自己调试通过的一套做法整理出来包括问题建模、MATLAB代码框架、算例对比和踩坑记录给正在做电力系统优化调度、毕业设计或者工程排程的同学一个可以直接参考的落地样本。先说清楚这套代码解决什么问题。机组组合本质上是一个大规模混合整数规划MIP在已知未来24小时负荷预测的前提下决定每台机组什么时候开、什么时候停以及每个时段的出力水平让总成本最小。如果只考虑成本那就成了纯经济调度但现实中调度员还必须保证系统安全——比如某条线路断开后其余线路不过载或者某台大机组突然跳闸后备用容量能顶上。这就引入了安全约束和热备用。热备用通常指旋转备用即已经并网、能快速增加出力的机组容量它要覆盖系统最大单一故障损失N-1原则这是电力系统运行里非常关键的一个工程约束。下面直接进入正题从建模思路到MATLAB实现再到结果分析一步步拆开讲。1. 问题建模与方案设计1.1 先理清UC问题的三个层次初学机组组合的人最容易被各种约束搞晕。我习惯把问题拆成三层第一层是时间耦合约束因为机组启停状态要跨时段传递第二层是空间耦合约束因为系统负荷要由所有机组共同平衡线路潮流也和多个节点的注入相关第三层是安全性与可靠性约束这是进阶内容包括备用容量约束、N-1静态安全约束、爬坡约束等等。套进“考虑安全约束及热备用”这个需求时模型会比普通UC多出两类约束备用约束任意时刻所有在线机组能提供的旋转备用之和必须大于系统要求的备用容量。这里有个细节备用容量通常取“当期预测误差最大机组容量”的某个比例工程上常用“最大单机容量”或“负荷的10%”来定。安全约束主要指线路潮流约束尤其是N-1预想故障后的静态安全校验。如果只在正常运行方式下校核潮流调度结果可能在故障时崩溃所以要在优化中显式加入N-1条件下支路潮流不超过限额的约束。1.2 目标函数与三层约束的数学表达我用的目标函数是最小化总成本包含三部分发电成本、开机成本、停机成本。发电成本通常用二次函数拟合为方便求解我会分段线性化。目标函数可以写为[ \min \sum_{t1}^{T} \sum_{i1}^{G} [C_{i}^{gen}(P_{i,t}) C_{i}^{su} \cdot u_{i,t}^{start} C_{i}^{sd} \cdot u_{i,t}^{stop}] ]其中(C_{i}^{gen}(P)) 可以是二次函数 (a_i P^2 b_i P c_i)在MATLAB中我用sdpvar定义变量再用optimizer调用求解器。为了减少非线性求解负担我直接采用分段线性化的方式而不是都交给整数非线性求解器。约束方面我分成几组功率平衡约束每个时段所有机组出力之和等于该时段负荷。这是等式约束。机组出力上下限约束开机能出力区间 ([P_{min,i}, P_{max,i}])停机出力为0。用状态变量 (u_{i,t}) 乘上上下限来约束。最小启停时间约束机组一旦启动至少要连续运行若干小时一旦停运至少要停若干小时。这个约束是在相邻时段状态之间建立的逻辑关系如果不写解出来的开停机看起来可能“毛刺”很多工程上没法执行。爬坡约束机组相邻时段的出力和备用容量提供能力要受爬坡率限制。这里需要非常小心因为爬坡约束不仅约束有功出力还约束旋转备用能力。很多初稿会漏掉“备用响应时段的爬坡能力”导致算出来的备用实际提供不了。旋转备用约束每个时段所有机组的备用容量之和要大于需求 (R_t)。备用容量本身受机组出力上限和爬坡响应能力的双重限制。安全约束部分我用直流潮流模型近似也就是忽略无功和网损只考虑有功潮流的线性关系。对每个节点和每条支路有[ PF_{l,t} \sum_{node} GSDF_{l,node} \cdot (Injection_{node,t}) ]其中 (GSDF) 是发电转移分布因子可以由阻抗矩阵求逆得到。正常运行工况下要求 (PF_{l,t} \le F_{l}^{max})。N-1校核时对每条支路 (k) 开断用断线后的GSDF重新计算其他支路潮流要求不越限。这种模型写起来很麻烦但MATLAB里用矩阵运算可以批量处理尤其是利用kronecker乘法把时段展开。后面代码部分我会展示这个技巧。1.3 热备用需求到底怎么取热备用旋转备用的容量取值直接影响系统安全性和经济性。取少了故障来临时可能躲不过取多了大量机组必须开着低负载运行成本飙升。我采用的工程化方法是[ R_t \max( \lambda_1 \cdot L_t,\ P_{max}^{largest\ unit} ) ]就是一个备用系数 (\lambda_1) 乘以当前负荷和历史最大单机容量取较大值。实际算例中我用10%负荷和最大单机容量600MW对比当最大单机容量超过负荷10%时就取最大单机容量。这个细节看起来小但在算高风电渗透率场景时影响很大。另外要注意热备用不是“想给多少就给多少”它要通过实际在线机组才能兑现。所以我建模时用的是市场化思路系统要求备用需求 (R_t)而机组提供的备用 (r_{i,t}) 是决策变量满足[ \sum_i r_{i,t} \ge R_t ]而且机组组合状态必须保证备用容量可以被提供。也就是说如果一台机组已经满发出力 (P_{i,t} P_{max,i})它的备用容量 (r_{i,t}) 最多只能为0。完整约束是[ P_{i,t} r_{i,t} \le P_{max,i} \cdot u_{i,t} ]这个约束是所有热备用模型的精髓它把“开着的机组”和“能发出的备用”在物理上绑定起来了。1.4 为什么选择线性规划整数变量的形式我见过有人用遗传算法、粒子群算这个模型我强烈不建议。UC是个标准的MIP问题使用Gurobi、CPLEX、或者开源的CBC求解器都能得到带最优性gap的全局解相比智能优化算法效率和结果稳定性完全不是一个量级。在MATLAB里我推荐Yalmip调用Gurobi代码简洁耦合矩阵自动生成省去手写MPS格式的麻烦。关键点是把所有二次成本函数分段线性化把启停逻辑用线性整数约束表达得到的是MILP混合整数线性规划。MILP有成熟的分支定界算法对于几十台机组、96个时段的问题几秒到几分钟就能收敛到1%以内的gap。2. MATLAB实现核心细节2.1 数据准备别把机组参数写死在代码里我习惯把所有数据放在一个结构体mpc里类似Matpower格式但不直接用Matpower因为UC问题比潮流计算复杂得多。机组数据结构建议包含mpc.gen每列分别表示机组编号、所在节点、(P_{min})、(P_{max})、爬坡率MW/h、最小启/停时间h、初始状态、运行成本系数 (a,b,c)、启动成本、停机成本、备用成本可选。mpc.T调度时段数通常取24小时也可以取96点15分钟一个点。mpc.load每个时段的总负荷MW如果是多节点还要给mpc.busload每个节点的负荷占比。mpc.branch支路数据首端节点、末端节点、电抗、容量。mpc.gsf_full预想故障集对应的GSDF矩阵这个要先离线算好。这里特别提醒爬坡率单位是MW/h如果你取15分钟一个时段需要把爬坡率除以4。我早期犯过这个错算出来的非可行解还以为是求解器的问题其实是从分钟到小时单位的换算错了。2.2 决策变量的设计与Yalmip语法我用Yalmip的binvar定义启停状态 (u_{i,t})用sdpvar定义出力 (P_{i,t}) 和备用 (r_{i,t})。T 24; G length(mpc.gen); % 状态变量 u binvar(G, T, full); % 1表示开机 start binvar(G, T, full); % 1表示该时段由停机转启动 stop binvar(G, T, full); % 1表示该时段由开机转停机 P sdpvar(G, T, full); % 出力 r sdpvar(G, T, full); % 热备用启动和停机状态不能和原状态割裂它们本质上是从 (u) 推导出来的但它帮我们把机组启动成本线性化。约束写法是Constraints []; % 初始状态关联 for t 1:T if t 1 % 对应初始状态用mpc.gen(i,9)表示 Constraints [Constraints, start(:,1) - stop(:,1) u(:,1) - init_status]; else % 状态转移u(t)-u(t-1) start(t)-stop(t) Constraints [Constraints, u(:,t) - u(:,t-1) start(:,t) - stop(:,t)]; end % 保证开机/停机不并存 Constraints [Constraints, start(:,t) stop(:,t) 1]; end这段代码的作用是把机组启停状态变化建模成“一个二进制变量交换”。如果不引入 start/stop启动成本就写不出来。很多初学者会用u(t) - u(t-1) 1表示启动但在负状态变化过程中会出现矛盾所以引入两个独立的0-1变量是标准做法。2.3 备用容量约束的正确写法热备用的核心约束除了和出力上限耦合还要考虑爬坡响应能力。具体来说机组在接到调度指令后15分钟内能增加多少出力就决定它能提供多少旋转备用。我用的约束是% 备用与出力上限耦合 Constraints [Constraints, P r repmat(Pmax, 1, T) .* u]; % 备用与爬坡响应能力耦合 (15min备用响应) RampLimitFactor ramp_rate * (1/4); % 15分钟备用的爬坡能力 Constraints [Constraints, r RampLimitFactor * u]; % 备用需求 Rreq max(0.1 * load, MaxGenCap); for t 1:T Constraints [Constraints, sum(r(:,t)) Rreq(t)]; end注意P r Pmax这个约束代替了单独的P Pmax因为它同时限定了机组可以同时提供的出力和备用之和必须低于物理出力上限。而r RampLimitFactor是备用的动态响应约束代表机组必须能快速爬坡否则即使有容量也不能算作热备用。2.4 安全约束直流潮流的批量实现多节点系统里不仅要保证机组总出力等于总负荷还要保证每个节点的功率平衡以及每条支路不过载。我引入节点注入向量Inj它等于该节点的机组出力减去该节点的负荷。然后通过GSDF矩阵计算支路潮流。GSDF的计算基于直流潮流。对于 (N) 节点系统先构建B矩阵电纳矩阵去掉参考节点求逆得到节点阻抗矩阵 (X)然后[ GSDF_{l,k} \frac{X_{i,k} - X_{j,k}}{x_l} ]其中 (l) 连接节点 (i,j)(x_l) 是支路电抗。在MATLAB里可以直接用Matpower的makeB和makeBdc帮助函数但为了确保理解我手写过一次。安全约束的每小时版本% 节点注入矩阵: gen_node_map表示第i台机组所在节点 Inj zeros(N, T); for g 1:G node gen_node_index(g); Inj(node, :) Inj(node, :) P(g, :); end % 减负荷 Inj Inj - NodeLoad; % NodeLoad 是 N x T 的节点负荷 % 支路潮流 GSDF * Inj Pflow GSFDC * Inj; % L x T % 正常运行约束 Constraints [Constraints, -BranchLimits Pflow BranchLimits];N-1校核需要把预想故障集合里的每个开断支路的GSDF都算出来然后在约束里逐条添加。为了避免把模型规模撑爆我采用“故障筛选”思路先运行正常态模型再检查哪些线路在故障态越限只添加越限线路的对应约束迭代求解两三轮即可。这个思路也叫“安全约束生成”工程上非常有效能大幅减小问题规模是实际调度系统里的成熟做法。2.5 求解配置与冷启动设置求解器我用GurobiYalmip调用方式ops sdpsettings(solver, gurobi, verbose, 2, debug, 0); ops.gurobi.MIPGap 0.001; % 设置1%的最优间隙 ops.gurobi.TimeLimit 300; % 最多求解300秒 result optimize(Constraints, Objective, ops);如果之前没装Gurobi也可以换成cplex、gurobi或cbc。开源的CBC求解器也能跑几十台机组的小模型但求解速度和Gurobi差距明显。我建议学术研究用Gurobi自用学习可以用CBC但最好早点习惯商用求解器的参数设置因为毕业后进企业大概率还是要用。求解后我用value(u)、value(P)、value(r)提取结果。这里还有一个小坑如果用Yalmip优化后的变量直接画图有时会带上后缀需要先执行double()或value()转换。3. 算例分析与结果解读3.1 测试系统设计6节点3机再加一条双回线为了把问题讲透我设计了一个小型测试系统3台机组、6个节点、8条支路负荷以3节点为主。参数如下机组所在节点Pmin (MW)Pmax (MW)爬坡 (MW/h)最小启停时间 (h)a ($/MWh^2)b ($/MWh)c ($)启动成本 ($)G1110060030040.00151050001200G2210040020030.0020124000900G335020010020.0030153000600负荷曲线我设成双峰型早高峰09:00-12:00、晚高峰18:00-21:00最高负荷850MW最低400MW。备用需求取10%负荷与最大单机容量600MW的较大值所以夜间低负荷时备用需求反而是600MW这时必须多开机组为的就是保N-1安全。这是一个非常反直觉的结论普通经济调度根本不会这样安排。3.2 有无安全约束的结果对比我跑了三组场景场景A忽略网络约束和N-1只做机组组合与经济调度。场景B加入正常运行方式下的线路潮流约束。场景C在B的基础上再加入N-1安全校核也就是完整的热备用安全约束模型。结果对比如下场景总成本美元G1平均出力占比G2启动次数线路过载越限次数A498,32039%23B525,44045%30正常态C583,76052%40含N-1从A到C总成本上涨了17%。这个涨幅不是“浪费”而是为了安全支付的必要代价。如果某条关键线路因为故障过载导致连锁跳闸可能造成的停电损失远高于这8万美元的日增成本。所以电网公司宁愿多花这个钱也要保证N-1后的系统稳定。这个对比也说明如果只跑忽略网络约束的UC得到的机组组合是没有实际执行价值的。我见过有人拿一个不考虑线路潮流的模型做新能源消纳分析结果算出来联络线功率超过实际限额整个结论都是错的。3.3 备用需求敏感度分析我把备用系数从5%往上调到25%看总成本变化趋势备用系数总成本美元最低时段开机台数5%510,150210%583,760315%618,440325%705,3004我的体会是备用系数从5%到10%成本上涨较快因为系统被迫让一台小机组全天候开机这台机组可能发电成本很高但它提供的容量价值就是热备用价值。从这个意义上看把备用需求放到目标函数里作为可选项让优化器在“开一台贵机组提供备用”和“多预留一台容量”之间权衡比硬性设置一个系数更经济。这部分可以引入备用的失负荷价值VOLL做进一步分析也是利用这个模型做扩展研究的好方向。3.4 机组组合的Gantt图怎么画调度结果最直观的表现就是机组启停时序我习惯用stairs或bar画成梯形图。代码参考figure; cols lines(G); for i 1:G subplot(G,1,i); stairs(0:T, [init_status(i), value(u(i,:))], LineWidth, 2, Color, cols(i,:)); ylim([-0.1, 1.1]); ylabel([G num2str(i)]); xlabel(小时); end这张图能够一眼看出哪些机组全天基荷哪些峰荷启停。我在写论文或报告时通常还会叠加负荷曲线和备用曲线作为第二张子图这样可以直观说明“备用需求高导致某台机组被迫开机”的原因。4. 常见问题与排查技巧4.1 求解时间过长怎么办MILP的求解时间对变量数量和整数变量数量非常敏感。如果变量数超过几万Gurobi也需要不少时间。我的经验是三步削减规模去掉无关约束很多情况下某些节点负荷很小支路容量很大这部分的安全约束可以提前删除。用灵敏度筛选器先判断。并行化时段如果按时间解耦可以引入拉格朗日松弛思路但太复杂。更简单的做法是减少预想故障集只对关键支路做N-1。设置MIP gap工程上不追求绝对最优1% gap已经可以接受。把MIPGap从0.0001放宽到0.01速度往往提升好几倍。4.2 模型提示不可行或求解结果全是NaN这类问题九成出在约束建模上。最常见的原因包括负荷平衡等式与备用约束冲突比如备用需求大于所有机组的最大可调容量导致无论怎么开都不能满足备用约束。这时需要把备用系数调小或增加机组。最小启停时间约束写反了我写过一版把“必须运行至少4小时”误写成“最多运行4小时”结果最优解是全停机。这个错误很难察觉我后来用随机状态测试才定位。初始状态没有考虑第一时段的启停转移约束如果漏了模型会认为所有机组都能瞬开瞬停进而给出不合理结果。排查方法很简单先用极小的负荷比如只有一台机组测试模型如果还不可行就逐步放开约束。用Yalmip的diagnize或check命令也能帮你定位是哪个约束不可行。% 求解后检查 info yalmiperror(result.problem); fprintf(求解状态: %s\n, info); if result.problem ~ 0 [viol, con] check(Constraints); [maxviol, idx] max(viol); if maxviol 1e-4 fprintf(最大违反约束值: %f\n, maxviol); fprintf(可能违反的约束类型: %s\n, class(con(idx))); end end4.3 数值缩放问题电力系统节点电压、功率、电抗大小相差很大如果不做单位归一化求解器数值稳定性会很差。我的习惯是功率用MW阻抗用标幺值但成本系数里可能有小数那么总成本的数值会特别大十万量级这没问题。关键是支路潮流约束里的GSDF元素可能很小而负荷是几百上千如果直接约束Pflow limit可能因为数值差导致约束被误解。稳妥做法是先把所有GSDF乘以1000让系数保持在0.01到100之间。4.4 Yalmip/Gurobi安装问题很多朋友卡在安装上Gurobi需要license但学术版免费。MATLAB里Yalmip的安装其实就是下载后把文件夹加入路径然后yalmiptest能看到已经安装的求解器。我遇到过Gurobi版本和Yalmip不兼容的情况通常换个内部版本号就能解决。如果报No suitable solver先确认是否setup过。这些属于环境问题排起来有点烦但确实是一次性投入。5. 从模型到实际调度的扩展思考5.1 从单时段到动态备用我上面的模型是“静态热备用”即每个时段的备用需求固定不变。但实际系统中负荷预测误差、新能源出力不确定性会让备用需求动态变化。如果你研究的系统里风电占比高建议把备用模型改成场景法比如生成多个风光出力场景计算CCG机会约束或者鲁棒优化版本。这个方向我在后续文章里会专门写MATLAB里用sdpvar也可以扩展但要注意场景数量别太多否则求解时间会指数上涨。5.2 热备用与储能、需求响应的交互热备用本质上是一种系统灵活性资源。随着新能源渗透率提高储能和可中断负荷也可以提供备用。要扩展的话可以在模型中增加储能设备的充电/放电变量和荷电状态变量并把储能提供的备用容量 (r^{ES}_{t}) 也加入备用平衡约束。这样UC模型就变成了一个多类型资源协调优化模型。MATLAB代码的改动主要是增加一组持续状态变量和能量平衡约束其他核心调度逻辑不变。5.3 学习建议与进阶方向如果你从零开始我的建议是先用这个6节点系统把代码跑通理解binvar和sdpvar的变量类型再用Matpower的IEEE 30节点数据替换。千万别一上来就跑几千台机组的系统模型复杂度会把你淹没。跑通后再往里面加安全约束和热备用——你会发现工程问题的难点从来不是算法本身而是把现实约束转换成数学表达式的功夫。我个人在这套代码上踩过最大的坑就是忽略了“备用爬坡响应”约束。最初我只写了P r Pmax结果算出来的备用虽然满足需求但实际机组根本来不及在15分钟内增加这么多出力。后来加入r RampLimitFactor模型才真正符合热备用的物理含义。如果你只关心优化算法而不了解电力系统很容易漏掉这种物理约束导致结果在纸面上好看、现场没法用。还有一个小技巧跑完模型后一定要把P、r、u这组结果回代到原始平衡方程里验算一下用max(abs(sum(P,1) - load))检查功率平衡的误差。如果误差超过0.01MW说明数值求解有问题可能是某条约束的容差设置太宽松。我通常会把Gurobi的IntFeasTol和FeasibilityTol设成1e-6虽然求解会慢一点但可信度高很多。最后把上面这套模型用在你的实际数据前先把“无安全约束”版本跑通作为基线再逐步叠加这样你也更容易定位每个安全约束带来的成本增量到底是多少。这个“由简入繁逐层叠加”的调试习惯是我在做调度系统项目时最重要的心得分享给你。