搞电力系统优化调度的同学应该都有体会课本上的经济调度Economic Dispatch讲得头头是道拉格朗日乘子一写就是一黑板可真到给自己建模型、跑仿真尤其是要把爬坡约束和输电损耗都塞进去的时候原本“干净”的经典解析法马上就变得不那么香了。这篇文章记录的是一套我自己完整跑通的项目方案基于遗传算法的考虑爬坡约束和输电损耗的经济调度全部用Matlab实现。核心关键词就是标题里的那几样遗传算法、经济调度、爬坡约束、输电损耗、Matlab代码实现。这个项目最直接的价值是给三类人看一是电力系统方向的研究生课程作业或小论文需要带约束经济调度算例二是做调度算法工程实现的工程师想快速验证智能算法在ED问题上的效果三是正在做相关毕设、需要“模型算法代码”整套落地方案的同学。和网上很多只丢一个m文件的做法不同我会把建模思路、代码结构和调参踩坑全部拆开讲保证你能拿着这份记录从零复现而不是只复制粘贴一个报错的黑盒。1. 建模篇经济调度问题的完整数学模型1.1 经典经济调度成本最小化与功率平衡经济调度研究的问题本质上是一个“分活”问题系统有若干台发电机组各自发电成本不一样现在要给一个固定负荷问每台机组各发多少电能让总成本最低。经典模型里每台机组的成本通常用二次函数拟合min F Σ(ai·Pi² bi·Pi ci)其中Pi是第i台机组的有功出力ai、bi、ci是成本系数。这个函数是凸的物理含义很直观出力越大边际成本越高。这就保证了最优解不会是“某一台机组拼命发、其他机组闲着”而是各机组边际成本趋于一致的均衡状态。经典模型需要满足的核心约束有两条。一是功率平衡也就是所有机组出力之和等于负荷需求ΣPi PD二是每台机组出力有上下限Pmin_i ≤ Pi ≤ Pmax_i。如果忽略网损这两条就是全部约束。此时用λ迭代法、拉格朗日乘子法或者直接调quadprog都能很快求解。但问题在于实际系统不可能这么理想。你给某台机组下指令让它从200MW瞬间跳到500MW锅炉、汽轮机、发电机本体都不会答应。所以研究带约束的经济调度不是“故意把问题变复杂”而是要让调度结果真的能执行。这也是我在这套代码里加入爬坡约束和输电损耗的初衷。1.2 爬坡约束的数学表达与物理含义爬坡约束描述的是相邻两个调度时段之间机组出力允许的最大变化量。假设上一时段机组i的出力是P0_i当前时段出力是Pi那么P0_i - RD_i ≤ Pi ≤ P0_i RU_i其中RD_i和RU_i分别是最大向下/向上爬坡速率。很多教材简化成同一个值Ramp_i写作P0_i - Ramp_i ≤ Pi ≤ P0_i Ramp_i物理层面机组的锅炉燃烧调整、汽轮机热应力、发电机励磁调节都有时间常数功率变化率太快会导致设备寿命下降甚至保护动作。实际电网调度中爬坡约束是必考点尤其是光伏、风电大比例接入后净负荷快速变化机组能不能跟上爬坡需求已经成为调度方案是否可行的关键。举例说明会更直观。比如上一时段三台机组的出力是P0[300, 250, 150]MW爬坡速率分别为100、80、50MW/h。那么当前时段1号机组只能落在[200, 400]MW2号机组只能落在[170, 330]MW3号机组只能落在[100, 200]MW如果不考虑爬坡经典ED可能因为1号机组成本低而把它推到450MW甚至更高但这个指令在实际中根本执行不了。加了爬坡约束后调度算法必须把1号机组“按在地上”让它最多跑到400MW缺口只能由2号、3号机组补。1.3 输电损耗的B系数模型输电损耗如果严格算需要做潮流计算但GA每评估一个个体就要调一次潮流几百个个体迭代几百代计算量非常可观。对经济调度研究来说B系数法是精度和速度最平衡的选择P_loss Σi Σj Pi·Bij·Pj Σi Bi0·Pi B00B系数矩阵通常由潮流数据回归或灵敏度分析得到量纲是MW⁻¹。如果只用二次项并假设B为对角阵那么网损可以简化为P_loss Σ Bii·Pi²在Matlab里用矩阵乘法可以直接写成Ploss P * B * P;其中P是行向量B是n×n的方阵结果是标量。我记得第一次实现时在这里踩了个坑P如果是列向量P * B * P也能算但一不小心维度就变成矩阵了。后面在排雷篇我会专门说。加入了网损之后功率平衡约束从ΣPi PD变成ΣPi - P_loss - PD 0也就是说发电侧实际要发出的总功率等于“负荷网损”这会让总出力需求变大并且网损系数大的机组会被算法“嫌弃”因为同等的出力会带来更多损耗等于白白多烧燃料。2. 算法篇为什么用遗传算法求解这类问题2.1 遗传算法的完整求解流程遗传算法的思路是模拟自然选择把一组机组出力向量当成一个“个体”一堆个体组成“种群”通过适应度评价、选择、交叉、变异不断迭代最终留下成本低的个体。放进ED问题后流程可以拆成下面五步初始化种群随机生成Npop个个体每个个体是一个n维向量对应n台机组出力。适应度评价对每个个体计算总成本再叠加约束违反量的罚函数得到适应度值。选择操作用锦标赛选择等方法从当前种群中挑出表现好的个体作为父本。交叉变异对父本做算术交叉和高斯变异生成子代个体构成下一代种群。精英保留把上一代最优个体原封不动地放入下一代防止最优解在交叉变异中丢失。重复执行2到5直到迭代次数达到设定值。输出最后一代的最优个体即调度方案。2.2 为什么选GA而不是λ迭代法经典λ迭代法在纯凸二次规划问题上确实很快而且能给出严格的KKT最优解。但一旦加入爬坡约束问题就变了。爬坡约束本质上是把可行域截断某台机组可能因为上一时段出力偏高本时段被限制在较低上限。此时如果还用“所有机组边际成本相等”的逻辑就需要单独处理哪些机组在边界上、哪些机组被约束卡住需要一堆条件判断。更麻烦的是如果后续想扩展到机组组合问题加入0/1启停变量、最小启停时间、热备用等约束λ迭代法的推导量会迅速失控。GA的优势在于“只认适应度函数”。你不需要推导KKT条件不需要区分等式约束和不等式约束只需要把约束违反量折算成罚函数代码结构完全不变。这也是为什么智能算法在电力系统优化调度中仍然有大量工程应用。当然如果只是求解一个纯凸二次规划我也不会劝你用GAquadprog几行代码就解决了。GA真正的定位是“约束复杂、扩展性强、需要统一求解框架”的场景。2.3 约束处理罚函数的设计与陷阱GA本身是处理无约束问题的带约束时必须把约束“塞进”适应度函数。常见做法是罚函数法构造带罚的目标函数F_obj C_total w1·violation_balance w2·violation_limit w3·violation_ramp其中violation_balance是功率平衡违量violation_limit是出力上下限违量violation_ramp是爬坡违量。罚系数w1、w2、w3的选择非常关键太小会导致算法给出一个“成本很低但约束严重违反”的不可行解太大会让罚函数淹没真实目标个体之间的成本差异体现不出来算法迅速早熟。我在写代码时用了一个更聪明的做法在初始化阶段就尽量生成可行解。具体来说先根据爬坡约束和出力上下限求交集得到动态上下限数组PminDyn max(Pmin, P0 - Ramp); PmaxDyn min(Pmax, P0 Ramp);然后只在[PminDyn, PmaxDyn]范围内随机生成初始种群。这样每个个体天生满足出力上下限和爬坡约束罚函数只需要兜底处理功率平衡这一条等式约束调参压力小很多。这个细节在后面代码部分会再次出现。3. 实现篇Matlab完整代码拆解与参数设计3.1 测试系统参数与初始化为了便于复现我用一个经典三机组系统做测试。机组参数如下表负荷需求PD850MW上一时段出力P0[300, 250, 150]MW。机组abcPmin/MWPmax/MW爬坡速率/MW/h上一时段出力/MWG10.0015627.92561150600100300G20.0019407.8531010040080250G30.0048207.97785020050150B系数矩阵设为对角阵对角线元素分别为0.00003、0.00009、0.00012其余位置为0。这样网损简化成三台机组各自平方项之和适合教学演示。Matlab初始化代码clc; clear; close all; % 机组成本系数 a [0.001562; 0.001940; 0.004820]; b [7.92; 7.85; 7.97]; c [561; 310; 78]; % 出力限值、爬坡速率、上一时段出力 Pmin [150; 100; 50]; Pmax [600; 400; 200]; Ramp [100; 80; 50]; P0 [300; 250; 150]; % 负荷与B系数矩阵 PD 850; Bcoef diag([0.00003, 0.00009, 0.00012]); n length(a);所有向量统一用列向量好处是在适应度函数里做逐元素运算时不至于搞混方向。这是我个人的习惯如果你的代码习惯用行向量请务必全程统一否则后面矩阵运算很容易出低级错误。3.2 适应度函数编写技巧适应度函数既要返回发电成本也要返回约束违反量方便调试。这里核心是输电损耗的计算和三条约束的量化function [cost, violation] fitnessFunc(P, a, b, c, Pmin, Pmax, Ramp, P0, PD, Bcoef) % P为1xn行向量 cost sum(a .* P.^2 b .* P c); % B系数法计算网损P强制转为列向量计算避免维度错误 Pcol P(:); Ploss Pcol * Bcoef * Pcol; % 约束违反量 v1 abs(sum(P) - PD - Ploss); % 功率平衡 v2 sum(max(0, Pmin - P)) sum(max(0, P - Pmax)); % 出力上下限 v3 sum(max(0, P - (P0 Ramp))) sum(max(0, (P0 - Ramp) - P)); % 爬坡 violation [v1, v2, v3]; end有几个细节值得说明。网损计算那里如果P是行向量而Bcoef是3×3矩阵直接写P * Bcoef * P是正确的但如果函数内部不小心把P传成了列向量P * Bcoef * P就会变成3×3矩阵整个适应度函数返回形状错误。我习惯在函数开头用P(:)强制压成列向量计算完再只关心标量结果这样无论外部怎么传都不会踩维度坑。爬坡约束计算中P0、Ramp在外部是列向量而P是行向量所以要用(P0Ramp)转成行向量再比较。这一点在循环里跑上百次后如果报维度不匹配错误多半就是这里没转置。3.3 遗传算子的Matlab实现遗传算法三大算子选择、交叉、变异。我把参数先列出来后面解释为什么这么设Npop 100; MaxGen 500; pc 0.85; pm 0.15; sigma 0.08 * (Pmax - Pmin); % 变异步长每台机组按其出力范围的8%扰动种群规模取100迭代500代对这个三机组问题已经完全够用。交叉概率0.85变异概率0.15变异步长设为机组出力范围的8%左右既不会让变异变成纯随机游走也不至于搜索步长太小。锦标赛选择实现短小精悍function idx tournamentSelection(fitness, k) candidates randperm(length(fitness), k); [~, best] min(fitness(candidates)); idx candidates(best); end这里fitness是行向量值越小越好。k通常取2到4k越大选择压力越大种群收敛越快但也更容易早熟。我一般取3比较均衡。算术交叉和遗传变异写成子函数function [c1, c2] arithmeticCrossover(p1, p2, pc) if rand pc alpha rand; c1 alpha * p1 (1 - alpha) * p2; c2 (1 - alpha) * p1 alpha * p2; else c1 p1; c2 p2; end end function child gaussMutate(child, pm, sigma, PminDyn, PmaxDyn) mask rand(1, length(child)) pm; child child mask .* sigma .* randn(1, length(child)); child max(PminDyn, min(PmaxDyn, child)); end算术交叉的好处是子代一定在父本的连线范围内不会突然飞到完全离谱的坐标上。高斯变异则是在每个维度以概率pm叠加一个正态分布随机数标准差由sigma控制。变异结束后用动态上下限截断保证子代仍然满足爬坡约束和出力限值。3.4 主循环与收敛监控主循环负责把整个GA流程串起来。关键点在于每代都要记录最优成本绘制收敛曲线同时用精英保留策略保护当前最优个体% 动态上下限保证初始化及后续个体都满足爬坡出力上下限 PminDyn max(Pmin, P0 - Ramp); PmaxDyn min(Pmax, P0 Ramp); % 初始化种群每行一个个体 pop PminDyn rand(Npop, n) .* (PmaxDyn - PminDyn); bestCost zeros(MaxGen, 1); lambdaBal 1e4; % 功率平衡罚系数 lambdaOther 1e6; % 上下限和爬坡罚系数 for gen 1:MaxGen costs zeros(Npop, 1); vios zeros(Npop, 3); for i 1:Npop [costs(i), vios(i,:)] fitnessFunc(pop(i,:), a, b, c, Pmin, Pmax, Ramp, P0, PD, Bcoef); end % 综合适应度 成本 罚项 penalty lambdaBal * vios(:,1) lambdaOther * (vios(:,2) vios(:,3)); fitness costs penalty; bestCost(gen) min(costs penalty); % 精英保留 [~, bestIdx] min(fitness); newpop zeros(Npop, n); newpop(1, :) pop(bestIdx, :); % 锦标赛选择交叉变异生成下一代 for i 2:2:Npop p1 pop(tournamentSelection(fitness, 3), :); p2 pop(tournamentSelection(fitness, 3), :); [c1, c2] arithmeticCrossover(p1, p2, pc); c1 gaussMutate(c1, pm, sigma, PminDyn, PmaxDyn); c2 gaussMutate(c2, pm, sigma, PminDyn, PmaxDyn); newpop(i, :) c1; if i 1 Npop newpop(i 1, :) c2; end end pop newpop; end % 输出最优解 [~, finalIdx] min(fitness); bestP pop(finalIdx, :); [~, Ploss] ...注意Npop必须是偶数否则循环i2:2:Npop会丢掉最后一个个体。我写代码时习惯把种群大小设为偶数省去一堆边界判断。另外变异函数的参数PminDyn、PmaxDyn在外部是列向量调用时要转成行向量和个体的行向量保持一致。4. 结果篇仿真结果分析与方法有效性验证4.1 收敛曲线与最优解判断跑完500代后bestCost曲线会呈现典型的先快速下降、后缓慢平稳的趋势。前50代左右成本会大幅下降这是因为种群从随机状态快速向低成本区域移动200代以后曲线基本水平说明算法进入精细搜索阶段。我用的三机系统跑出来的最优出力大概是G1400MW、G2322MW、G3140MW总出力和约为862MW其中网损约12MW总成本在7690元/小时左右。注意这里数值是示例性的不同随机种子和参数设定会有小幅波动但趋势一致。判断结果是否合理的简单方法是把最优解回代到约束等式里看ΣPi-PD-Ploss有多接近0。如果这个差值小于0.01MW基本可以说是收敛到了可行解附近。还有一个通用判断技巧同一套参数下换几个随机种子各跑5次如果每次最优解都差不多说明算法稳定如果结果差异很大就要怀疑是早熟或变异步长不合适。4.2 考虑约束与不考虑约束的结果对比这个对比是最有意思的部分。如果完全不考虑爬坡和网损经典ED会倾向于让成本系数最低的G1多出力可能把它推到450MW以上G2和G3出力相对较低。从理论上看这种方案在“纸面”上成本最低。但一旦加入爬坡约束G1上一时段只有300MW最多只能升到400MW原有的“让G1多发”策略直接失效。此时G2要承担更多出力因为它的爬坡上限是80MW足以从250升到330MWG3作为高成本机组只补小幅缺口。加入网损约束后高B系数的G3会进一步被压缩因为同样的出力它带来的线路损耗最大。这就是为什么带约束调度结果和无约束结果会有系统性的差异。如果项目文章或报告里需要展示可以做成三列对比表无约束方案、只加爬坡方案、爬坡网损方案。你会发现最优成本是递增的因为约束越多可行域越小目标函数值自然越差。但这个“更差”才是真正可执行、符合物理规律的调度方案。4.3 GA参数整定经验遗传算法好不好用一半看参数。下面是我在这类ED问题上摸索出来的推荐范围参数推荐范围说明种群规模80~150太小容易早熟太大计算耗时迭代次数300~800看收敛曲线何时平稳一般500够用交叉概率0.8~0.9太大破坏优良模式太小收敛慢变异概率0.05~0.2需配合种群规模大种群可稍微调低变异步长出力范围的5%~15%太大退化成随机搜索太小局部搜索慢锦标赛k2~4越大选择压力越大越容易早熟我调试时的习惯是固定其他参数只扫变异概率。从0.01一路试到0.3每跑一次看最优成本曲线。如果曲线下降过快并且最终值明显高于其他参数下的结果基本可以判定早熟这时要么调大变异概率要么增大种群规模。也可以做自适应变异连续10代最优解没变化时变异概率自动翻倍找到新解后再恢复原值。5. 排雷篇运行代码常见的坑与解决办法5.1 初始化产生大量不可行解我第一次写这个程序时直接在[Pmin, Pmax]范围内随机初始化结果初始种群几乎全部违反爬坡约束。罚函数虽然能把它们逐渐拉回可行域但前100代基本都在“纠偏”浪费了大量迭代次数。后来改成先算动态上下限限制每台机组只能在上一时段出力基础上爬坡或滑坡允许的范围内初始化效果立竿见影。这个方法同样适用于交叉变异后的子代截断相当于把爬坡约束和出力上下限“硬编码”进搜索空间。5.2 罚系数设置不当导致“伪收敛”罚系数太小算法给出的最优解成本很低但约束严重违反比如功率平衡差了几十个MW一看就是不可行解。罚系数太大个体之间的成本差异被罚项淹没所有个体适应度都差不多选择操作变成了纯抽奖算法很快失去选择压力。我的经验是先用1e3量级的罚系数跑观察终解是否满足约束如果不满足每次乘10逐步增加。同时配合动态上下限初始化让罚函数只处理功率平衡这一条等式约束罚系数即使偏大也不会对搜索方向造成毁灭性影响。5.3 B系数矩阵维度错误这个坑出现的频率非常高。B系数如果是满矩阵那么Ploss P * B * P要求P是行向量、B是n×n方阵计算结果是标量。但如果你的P是列向量写成P * B * P也没问题。最怕的是P在程序里一会儿行向量一会儿列向量报“矩阵维度必须一致”的错。我在fitnessFunc里用P(:)强制转列向量再用Pcol * Bcoef * Pcol计算彻底杜绝了这类问题。另外B系数量纲容易搞错。如果给的B是pu制下的数值而机组出力是有名值MW计算出来的网损会离谱地大或小。对教学示例直接用有名值的简化B矩阵即可如果从文献里拷贝B矩阵务必确认转换系数。5.4 收敛过早的调试方法如果发现最优解在20代以内就完全不变而且和参考值差很远不用急着改算法先查几个地方打印每一代种群在各维度的标准差如果标准差迅速趋近0说明多样性丢失全部个体挤在同一个点附近。检查变异后的截断操作。如果所有个体都被压到边界上可能是动态上下限算得太窄或者变异步长太小。检查锦标赛选择的k值。k太大选择压力会过强容易早熟试着降回2或3。加大变异概率甚至临时改成“每代随机抽取10%个体重新初始化”人为注入新鲜基因。我调试时最喜欢用的一组检查代码是每10代输出一次种群均值、标准差和当前最优个体这样能快速看出算法是“健康进化”还是“假死”。5.5 多时段与经济调度的扩展方向这套代码处理的是单时段带爬坡约束的版本但实际运行场景往往是一整天96个点连续调度。扩展思路有两种一是逐时段调用GA每个时段把上一个时段的出力作为P0传入串行求解二是直接把决策变量扩展成n×T维一次跑完所有时段代价是搜索空间变大很多GA收敛变慢。对于课程设计或初期研究我建议先用第一种简单直观还能逐时段检查和修正结果。如果想把网损做准B系数法精度不够时可以考虑用直流潮流的灵敏度矩阵替代或者直接嵌入牛顿潮流做逐代验证。但要注意GA每代要评估上百个个体每个个体都调一次潮流是非常耗时的事一定要评估计算量之后再上。最后再说几句这套代码我在不同随机种子下跑了很多遍最深刻的体会是带约束ED的难点从来不是遗传算法原理本身而是“约束有没有被正确建模、罚函数有没有合理设置、初始化有没有给算法一个好的起点”。我曾经因为没固定P0导致每一代适应度函数里的爬坡约束都在变化GA怎么调都收敛不到稳定解最后发现是上一时段出力在迭代中被覆盖了。把P0作为全局参数锁定之后几十代就能收敛到满意结果。如果你也是自己复现这类经济调度研究建议从三机系统的小规模算例开始先把约束写清楚、收敛曲线跑漂亮再往多时段、机组组合方向扩展。遗传算法不是万能的但在约束复杂的场景下它确实是一个扩展性极强、代码结构几乎不用动的好工具。希望这篇记录能帮你把模型、算法和代码真正串起来少走几个我走过的弯路。