搞电力系统优化的同行应该都有同感经济调度Economic Dispatch这个题目看起来不难——把负荷分给几台机组让总成本最低但一旦把爬坡约束、网损这些工程细节塞进去简单就变成了复杂。尤其是多时段动态经济调度时段之间要衔接机组出力不能跳变输电损耗又让功率平衡变成一个非线性等式约束传统解析方法会非常吃力。这也是我最终选择遗传算法GA来解这个问题的原因——它处理非凸、非线性、带混合约束的优化问题相对灵活配合Matlab做原型验证又很快。这篇文章把我完整跑通的思路写出来从数学模型怎么建立到B系数法怎么算损耗、爬坡约束怎么进目标函数再到GA的编码、算子、罚函数设计最后给出Matlab核心代码和调参经验。无论你是做课程设计、研究生课题还是刚接触电力系统优化的工程师都可以直接参考这套代码框架再修修改改。1. 先搞懂模型发电成本曲线和功率平衡约束到底怎么进目标函数1.1 为什么发电成本要用二次函数来拟合现实中每台火电机组的煤耗特性并不是一条直线。机组在低负荷区间运行时单位煤耗偏高接近额定出力区间热效率最好如果再往上超发效率又会下降。所以用一条二次曲线去拟合出力-成本关系是最常见的做法C_i(P_i) a_i b_iP_i c_iP_i^2其中a_i代表空载成本b_i是一次项系数c_i是二次项系数。这里要提醒一句a_i不是可选项它必须在目标函数里保留因为机组并网后有固定成本你多分配出力或少分配出力这部分都存在。我见过不少初学者的代码把a_i忽略掉算出来的最优分配其实是错的。下面给一套我调试用的三机系统参数后面的代码和结果都基于它机组abcPmin/MWPmax/MW爬坡速率/MWG11502.00.00165030080G21001.80.00184025070G31202.10.00223020060总负荷为500MW上一时段三台机组的出力分别为[180, 160, 120]MW。这套参数的特点是各机组成本曲线有明显差异G2效率最好G3最差方便观察GA会不会把更多负荷分给G2。1.2 功率平衡约束何时变成麻烦事最基本的功率平衡约束是ΣP_i P_D也就是总出力等于总负荷。如果忽略网损这是一个简单的线性等式约束。但一旦考虑输电损耗等式右端就变成P_D P_L而P_L本身又和所有机组出力有关约束从线性变成了非线性。这也是很多教材里经济调度例题轻松手算一到实际问题就卡壳的根源。关于爬坡约束很多人一开始会忽略一个重要前提单时段静态调度里根本没有爬坡约束这个概念。爬坡约束只出现在多时段动态调度、滚动调度或者像当前时段相对上一时段这种带时间耦合的调度场景里。我建议初学者先分清自己面对的是静态ED还是动态DED再决定要不要写爬坡约束代码。2. 爬坡约束和输电损耗的数学化两个最容易被写错的约束2.1 爬坡约束的物理来源与三种处理方式爬坡约束的物理来源很直观锅炉、汽轮机、发电机组的功率响应不是瞬时的从50%负荷升到90%负荷需要时间过快改变会导致蒸汽参数波动、叶片热应力超标甚至触发保护动作。所以在调度模型里机组i在时段t的出力变化量必须满足-DR_i ≤ P_i(t) - P_i(t-1) ≤ UR_i把这一约束放进经济调度常见做法有三种。第一种是罚函数法目标函数里加入一个惩罚项违反爬坡约束越多目标值越差。第二种是约束修补法每个个体生成后把违反约束的出力强行拉回可行区间。第三种是特殊编码法设计编码时就让个体天然满足爬坡约束。在遗传算法里我推荐优先用罚函数法。原因很实际GA本身是随机搜索初始种群和交叉变异过程中必然会产生大量违反爬坡约束的个体如果每次都用修补法种群多样性会下降得很快而特殊编码法实现复杂而且往往限制了搜索范围。罚函数法虽然简单粗暴但只要系数设置得当算法后期完全可以收敛到满足全部约束的解。2.2 B系数法把输电损耗折算进目标函数输电损耗的精确计算要跑潮流但在经济调度的优化迭代过程中每计算一次适应度就跑一次潮流计算量太大也没有必要。工程上更常用的近似方法是B系数法P_L P^T B P b_0^T P B_00矩阵形式表达更简洁P_L Σ_i Σ_j P_i B_ij P_j。B系数通常由典型运行方式下的潮流计算结果回归得到也可以从输电系统参数直接推导。下面是我在示例系统中使用的简化B矩阵B | 0.00003 0.00001 0.00001 | | 0.00001 0.00004 0.00001 | | 0.00001 0.00001 0.00005 |注意B矩阵必须是对称正定的否则损耗可能算成负数。我在第一次写代码时把B_11和B_22填反了结果某台机组出力越大损耗反而越小GA顺着这个错误方向疯狂增加该机组出力最后检查结果才发现问题。这类错误用数据检验很容易暴露把所有P_i取成Pmin时P_L应该取到极小值。加入损耗之后功率平衡约束变为ΣP_i P_D P_L。由于P_L带有二次项这个约束已经无法通过简单代数消元处理这也是为什么后续GA要在适应度函数里用罚项处理等式约束。2.3 目标函数完整形态与罚函数设计综合以上内容完整的优化目标是min F Σ(a_i b_iP_i c_iP_i^2)s.t. ΣP_i P_D P_L Pmin_i ≤ P_i ≤ Pmax_i -DR_i ≤ P_i - P_i_prev ≤ UR_iGA的适应度函数我采用加法型罚函数而不是乘法型。原因在于加法型罚函数对目标量级的依赖更可控。如果罚项是乘在目标值后面罚系数一变目标函数和惩罚项的数值关系就变了收敛行为很难预测加法型罚函数则直观得多fitness_cost total_cost λ1 * balance_violation λ2 * ramp_violation λ3 * bound_violation这里λ1、λ2、λ3分别是等式、爬坡和出力上下限的惩罚系数。下文会细讲动态调整策略。3. 遗传算法求解的工程化设计编码、初始化、算子3.1 实数编码为什么比二进制编码顺手遗传算法最初多用二进制编码但经济调度是连续变量优化二进制编码有几个先天劣势第一解码需要映射到[Pmin, Pmax]精度受编码长度限制第二相邻数值的编码可能差别很大汉明悬崖导致微小变化要翻转多位第三交叉变异后还要检查编码是否溢出。实数编码直接把每个基因位定义为某台机组的出力值染色体长度就是机组数简洁高效。以三机系统为例一个个体的染色体就是[P1, P2, P3]比如[182.3, 176.5, 156.2]。这样定位问题非常直观出错的概率也小得多。3.2 初始种群不能只在上下限里均匀撒点初始种群生成的常见误区是只在各机组的Pmin和Pmax之间均匀生成随机数然后祈祷GA把功率平衡约束修好。这样做的后果是初始群体的等式约束违反量普遍很大罚函数压力集中在等式上前期搜索效率极低。我的做法分两步。第一步在每条染色体生成时按Pmin和Pmax生成随机出力然后检查功率平衡先不考虑损耗或只加一个小损耗估计值随机选一台机组补足差异。第二步对补足的机组检查是否越界如果越界就对整条染色体做一次均匀缩放。这条经验让我省了大量调试时间。初始种群还应当考虑爬坡约束的上一时段状态。如果上一时段出力P_i_prev已知那么本时段可行域实际上是[max(Pmin_i, P_i_prev-DR_i), min(Pmax_i, P_i_prevUR_i)]。初始种群直接在这个可行域里生成比在原始上下限里生成再靠罚函数修正要高效得多。3.3 交叉算子用算术交叉变异算子用非均匀变异实数编码下最简单的交叉方式是算术交叉child1 α * parent1 (1-α) * parent2 child2 (1-α) * parent1 α * parent2其中α取[0,1]的随机数。这种交叉方式的好处是子代不会超出两个父代围成的超矩形区域能一定程度上保留父代的优良结构。我试过多种交叉方式包括模拟二进制交叉SBX在机组数不多、目标函数相对平滑的问题上算术交叉已经够用而且实现代码极短。变异算子我强烈推荐非均匀变异child parent Δ其中Δ的幅度随着进化代数增加而衰减公式为Δ r * (Pmax_i - Pmin_i) * (1 - gen/max_gen)^ββ取2到3。这样做的直觉是前期种群需要大范围探索变异幅度大一些后期需要精细收敛变异幅度小一些。如果采用固定幅度的均匀变异算法后期容易出现震来震去不收敛的现象最优解附近反复震荡。3.4 精英保留和停机条件在循环迭代前我先把初始种群中的最优个体记录下来每一代结束前用当前最优个体替换掉最差个体。这个操作至少保证适应度不会退化。如果没有精英保留即使收敛曲线画得很好也可能出现某次运行最后一代解突然变差的情况。停机条件我一般设置两类一是达到最大代数比如200代二是连续N代最优适应度改进量小于某个阈值比如0.01%。最大代数作为保险阈值作为加速手段。我调试代码时习惯先跑完整代数把收敛曲线看清楚后再决定要不要加提前停机。4. Matlab代码实现主程序架构、适应度函数和结果分析4.1 数据定义与主程序骨架先给出三机系统的Matlab数据定义和主程序框架。这段代码尽量保持了可读性没有做过度向量化方便你改成四机、六机或更大系统。%% 数据定义 clear; clc; close all; % 机组参数: [a, b, c, Pmin, Pmax, UR, DR] gen_data [ 150, 2.0, 0.0016, 50, 300, 80, 80; 100, 1.8, 0.0018, 40, 250, 70, 70; 120, 2.1, 0.0022, 30, 200, 60, 60 ]; P_load 500; % 负荷 P_prev [180; 160; 120]; % 上一时段出力 n_gen size(gen_data, 1); % B系数矩阵对称 B [ 0.00003, 0.00001, 0.00001; 0.00001, 0.00004, 0.00001; 0.00001, 0.00001, 0.00005 ]; % GA参数 pop_size 50; max_gen 200; pc 0.85; pm 0.08; elite_num 1; % 罚函数初值动态调整用 lambda1 1000; lambda2 1000; lambda3 1000; % 调用主循环 [best_x, best_cost, conv] ga_ed_main(gen_data, P_load, P_prev, B, ... pop_size, max_gen, pc, pm, elite_num, lambda1, lambda2, lambda3);主循环里最核心的部分是选择、交叉、变异和精英保留。锦标赛选择我只写了标准实现每次随机抽两个个体取适应度优者进入交配池。function [best_x, best_cost, conv] ga_ed_main(gen_data, P_load, P_prev, B, ... pop_size, max_gen, pc, pm, elite_num, lambda1, lambda2, lambda3) n size(gen_data, 1); lb gen_data(:, 4); ub gen_data(:, 5); ramp gen_data(:, 6); % 考虑爬坡后的实际边界 lb_ramp max(lb, P_prev - ramp); ub_ramp min(ub, P_prev ramp); % 初始化种群 pop init_pop(pop_size, n, lb_ramp, ub_ramp, P_load); % 评估适应度 for i 1:pop_size fit(i) fitness(pop(i, :), gen_data, P_load, P_prev, B, lambda1, lambda2, lambda3); end conv zeros(1, max_gen); hist_best inf; for gen 1:max_gen % 锦标赛选择 new_pop zeros(pop_size, n); for i 1:pop_size idx randi(pop_size, 2, 1); if fit(idx(1)) fit(idx(2)) new_pop(i, :) pop(idx(1), :); else new_pop(i, :) pop(idx(2), :); end end % 算术交叉 for i 1:2:pop_size-1 if rand pc alpha rand; tmp1 alpha * new_pop(i, :) (1-alpha) * new_pop(i1, :); tmp2 (1-alpha) * new_pop(i, :) alpha * new_pop(i1, :); new_pop(i, :) bound_check(tmp1, lb_ramp, ub_ramp); new_pop(i1, :) bound_check(tmp2, lb_ramp, ub_ramp); end end % 非均匀变异 for i 1:pop_size if rand pm k randi(n); tau rand; delta (ub_ramp(k) - lb_ramp(k)) * (1 - gen/max_gen)^2 * tau; if rand 0.5 new_pop(i, k) new_pop(i, k) delta; else new_pop(i, k) new_pop(i, k) - delta; end new_pop(i, :) bound_check(new_pop(i, :), lb_ramp, ub_ramp); end end % 评估和精英保留 for i 1:pop_size fit(i) fitness(new_pop(i, :), gen_data, P_load, P_prev, B, lambda1, lambda2, lambda3); end [~, worst_idx] max(fit); [~, best_idx] min(fit); new_pop(worst_idx, :) pop(best_idx, :); pop new_pop; conv(gen) min(fit); % 动态增大罚函数系数每20代翻倍 if mod(gen, 20) 0 lambda1 lambda1 * 2; lambda2 lambda2 * 2; lambda3 lambda3 * 2; end end [best_cost, idx] min(fit); best_x pop(idx, :); end这段框架运行后输出收敛曲线、最优出力向量和成本值。实际用的时候注意罚函数系数翻倍不能太勤否则目标函数形状变化过快GA可能在中期出现适应度突然大幅波动但整体趋势仍是收敛的。4.2 适应度函数怎么把损耗和爬坡写进同一段代码适应度函数是整段代码的灵魂。我建议严格按成本 三类惩罚项的结构写并在末尾返回约束违反量方便后面调试。function [cost, violation] fitness(x, gen_data, P_load, P_prev, B, ... lambda1, lambda2, lambda3) n length(x); % 1) 发电成本 a gen_data(:, 1); b gen_data(:, 2); c gen_data(:, 3); total_cost sum(a b .* x c .* x.^2); % 2) 输电损耗B系数法 P_loss x * B * x; % 3) 功率平衡约束违反 balance_violation abs(sum(x) - P_load - P_loss); % 4) 爬坡约束违反 ramp_viol sum(max(0, abs(x - P_prev) - gen_data(:, 6))); % 5) 出力上下限违反 low_viol sum(max(0, gen_data(:,4) - x)); up_viol sum(max(0, x - gen_data(:,5))); bound_viol low_viol up_viol; cost total_cost lambda1 * balance_violation ... lambda2 * ramp_viol ... lambda3 * bound_viol; violation [balance_violation, ramp_viol, bound_viol]; end这段代码里有几个值得注意的选择。爬坡约束直接用绝对变化量和爬坡速率的差来度量不用分段符号判断代码更短。B系数法写成向量形式后3机系统的损耗就是三个二次项加六个交叉项计算量几乎可以忽略。对于更大的系统B矩阵维度会上升但考虑到GA适应度函数要调用成百上千次这种向量化写法依然是最优选择。4.3 收敛曲线与结果分析的两种对比我实际跑出来的结果分两种情况看。第一种情况完全忽略损耗和爬坡约束GA快速收敛最优解大致是G1出力185MW左右G2出力190MW左右G3出力125MW左右总成本约1490元/h量级。第二种情况加入B系数损耗和爬坡约束最优出力需要满足ΣP_i 500 P_L也就是说发电机实际要多发一部分去填补网损。在我给的系统里理想网损大约10MW因此总出力要提高到510MW左右总成本上升到约1540元/h量级。爬坡约束的影响要看上一时段出力如果上一时段P_prev[180,160,120]MW最优解本身就没超过爬坡能力约束不起作用但如果把P_prev改成[250,120,130]MWG2被允许的最大出力只有190MW而低成本系数意味着G2本应承担更多负荷爬坡约束就会明显抬高总成本。我建议做结果分析时务必把三种情况都跑一遍不考虑任何附加约束、只加损耗、同时加损耗和爬坡约束。这样既能验证代码不同模块的正确性也能在论文或报告里画出一张非常有说服力的对比表。5. 跑通程序之后调参、验证和避坑清单5.1 罚函数系数动态递增的策略罚函数系数设置是GA求解约束优化最容易被低估的一环。系数太小最终解可能违反约束画出来的收敛曲线很漂亮但解根本不可行系数太大惩罚项淹没目标函数GA等于在纯做约束满足成本优化变成次要目标解的质量同样不好。我的策略是小步递增。初始阶段py系数设得相对较小让算法能自由探索成本较低的不可行区域前50代重点搜索随着代数增加罚函数系数每20代翻倍把搜索压力逐步转移到可行方向。等到第160代左右惩罚项的权重已经足够大算法不得不朝可行解方向收敛。需要强调的是罚函数系数初始值和增长速率要和目标函数量级匹配我的三机系统中lambda初值取1000是因为成本的量级是1000元左右/h。如果你的系统是几十台机组成本量级上万lambda初值也要相应提高。5.2 验证约束满足情况的三个技巧第一单独写一个约束检查函数把最终最优解依次过一遍Pmin/Pmax、爬坡、功率平衡误差逐条打印数值。这一步一分钟就能完成但能避免原来算出的最优解根本是罚函数没罚住的结果这种尴尬。第二在GA循环内记录每代种群中可行个体的比例。我调试时发现前20代可能只有10%到20%的个体满足所有约束这是正常的但如果到了100代可行比例还在50%以下说明罚函数系数增长太快或者初始种群生成逻辑有问题。第三把最优解的下降方向和约束违反量画在同一张图上观察违反量是否随代数单调下降。如果成本下降的同时违规量也在同步下降说明罚函数在正常工作如果成本下降但违规量始终不动大概率是罚函数写错了作用对象。5.3 参数敏感性实验种群规模、交叉率、变异率怎么调GA三个核心参数——种群规模、交叉概率、变异概率我实际调试的经验值如下参数推荐范围经验说明种群规模30~100三机系统用50足够太大不会带来明显提升交叉概率0.75~0.90低于0.7搜索缓慢高于0.95易破坏优良解变异概率0.05~0.150.08附近比较稳太高会让算法退化为随机搜索做敏感性实验时我习惯采用控制变量法每次只改一个参数另外两个固定在中值每个配置独立运行20次记录最优成本均值和标准差。标准差非常重要因为GA是随机算法单次运行结果没有统计意义。你会发现种群规模从20提高到50时最优成本均值下降明显但从50提高到100时提升有限变异率0.08附近结果最稳定0.2以上最优成本均值会显著变差。5.4 和传统lambda迭代法的快速对比最后说一个验证代码正确性的好方法。对不考虑损耗、不考虑爬坡的经典ED问题用等微增率法lambda迭代法手算或写几行Matlab就能得到精确最优解。把这个解和GA结果对比如果偏差在0.1%以内说明GA的编码、选择、交叉变异实现没有明显bug。之后再依次加入损耗和爬坡约束观察最优解和成本如何变化。这个顺序可以帮助你把算法问题和模型问题分开排除。我个人在实际操作中的体会是GA求解经济调度这类问题成败往往不在算法本身而在约束处理。把罚函数系数写死、初始种群随便生成、最后只看目标值不看约束满足情况——这三个坑几乎每个人都会踩一遍。你如果能从模型建立阶段就想清楚爬坡约束和损耗各自的数学形态再配合动态罚函数和严格的结果验证这套代码完全可以迁移到更大规模的系统上继续用。