电力系统的经济调度Economic DispatchED是运行层最基础也最核心的优化问题之一。它的目标一句话就能说清楚在满足负荷需求、机组出力上下限、爬坡速率、输电线路损耗等一系列约束的前提下把总发电成本压到最低。但问题一旦展开非线性、非凸、多约束、多时段这些词就会接踵而来。我最早接触这个题目时试图用传统的等微增率方法做多时段优化结果被爬坡约束和网损公式折磨得够呛——不是发散就是陷入局部最优算出来的结果自己都不敢用。后来换成遗传算法Genetic AlgorithmGA把约束全部塞进罚函数反而干净利落地解出来了。这篇文章就把完整的建模思路、GA设计逻辑、Python代码实现和调试经验都写清楚适合正在做电力系统优化方向的学生、研究储能调度或微电网经济运行的工程师以及想用进化算法解决工程优化问题的朋友参考。1. 先把经济调度这三个字拆开看清楚1.1 成本函数为什么是二次的三个参数各代表什么火电机组的燃料消耗与出力并不是简单的线性关系。锅炉效率和汽轮机热力特性叠加后实际工程上普遍用二次函数去拟合单位时间的燃料成本F(P) a * P² b * P c其中 P 是机组出力MWa 的单位是 元/MW²hb 的单位是 元/MWhc 的单位是 元/h。a 代表了成本曲线的弯曲程度b 代表边际成本的主体部分c 是空转也有的一笔固定成本。举个实际算例中的参数某机组 a0.001562b7.92c561在出力 300MW 时小时成本大约就是 0.001562300² 7.92300 561 ≈ 2933.6 元。多台机组的总成本就是各自成本求和。到这里还不算难。真正的麻烦在于目标函数是非线性的如果再加上机组出力上下限约束这就是一个典型的非线性规划问题。传统方法最常用的等微增率准则本质上是拉格朗日乘子法要求目标函数可导且问题凸性较好。但实际系统的成本曲线如果考虑阀点效应、禁止运行区甚至输电损耗引入的非线性项凸性就保不住了等微增率很容易被复杂可行域中的局部极值带偏。1.2 引入爬坡约束静态问题变动态复杂度上了一个台阶如果只做单时段经济调度决策变量就是每台机组当前时刻的出力各时段之间完全独立优化一次完事。但真实系统里机组不能随意跳变出力锅炉汽轮机有热应力限制负荷变化太快会缩短设备寿命甚至引发事故所以每台机组都有爬坡速率约束P_i(t) - P_i(t-1) ≤ UR_i向上爬坡限制 P_i(t-1) - P_i(t) ≤ DR_i向下爬坡限制UR_i 和 DR_i 通常取一样的值单位是 MW/h。这个约束一旦打开所有时段就像被链条串起来一样不能再逐时段独立求解了。假设有 3 台机组、24 个时段决策变量数就从 3 变成了 72而且任意两个相邻时段的机组出力互相牵制。这是这个问题从“课堂习题”变成“工程挑战”的第一步。我在实际项目里遇到过的情况是不加爬坡约束时优化结果里机组1在 8 点还在 350MW9 点直接跳到 580MW。这在数学上完全“最优”但调度员看了一眼就把方案打回来了——实际机组接口根本不接受这种操作。加了爬坡约束之后解才会变得“像人话”负荷爬坡由多台机组配合分摊而不是让一两台机组原地起飞。1.3 输电损耗用B系数表达非线性等式约束是最大的坑大型电网里电能从发电端送到负荷中心线路上的有功损耗不可忽略。如果忽略损耗约束就是简单的“发电总和负荷”但真实系统中发电机总出力其实要大于负荷多出来的那部分就是网损。工程中最常用的网损模型是B系数法Kron公式把网损写成各发电出力的二次型P_L sum_i sum_j P_i * B_ij * P_j写成向量形式就是 P_L P̄^T * B * P̄。B矩阵是对称阵维度等于发电机台数。比如三机系统B [[0.000030, 0.000010, 0.000005], [0.000010, 0.000025, 0.000008], [0.000005, 0.000008, 0.000020]]假设三台机组出力分别是 300、200、150MW网损大约是 3000.00003300 2000.000025200 1500.00002150 2*(3000.00001200 3000.000005150 2000.000008150)算下来差不多就十几 MW 量级。看起来不大但优化问题的约束从“总和等于负荷”变成了“总和等于负荷加网损”而网损本身又依赖出力等式约束变成了一个隐式非线性方程。纯数学解法需要迭代求拉格朗日乘子、还要反复求解非线性方程组操作繁琐且容易发散。但对GA来说这只是一行矩阵乘法的代码非常友好。2. 遗传算法为什么适合这类问题以及具体怎么设计2.1 传统方法卡在哪GA的优势到底是什么先说说等微增率λ迭代为什么在这个题目上施展不开。经典的 λ迭代法需要满足协调方程dF/dP λ其中考虑网损时还要除以损耗修正因子 PF_i 1 / (1 - ∂P_L/∂P_i)光是算这个偏导就要推导半天。加上爬坡约束后相邻时段耦合拉格朗日乘子数量爆炸KKT条件系统变得异常庞大解起来既不稳定又不直观。而且一旦成本函数非凸λ迭代求到的解可能连可行性都保证不了。GA属于进化算法不依赖梯度、不要求凸性、不要求可导只要能把目标函数和约束“写出来”就能在解空间里进行群体搜索。这对经济调度这个场景非常实用成本函数复杂、约束多种多样改任意一个条件传统算法可能整套公式都要重推而GA只需要改动目标函数里的几行代码。代价也很明显不保证全局最优计算量偏大参数整定需要经验。但工程上要的是“在合理时间内找到一组足够好且可执行的解”GA就是这个定位。2.2 编码、目标函数、选择、交叉、变异的具体落地GA里有三种常见编码二进制编码、实数编码、排列编码。经济调度的决策变量是连续的机组出力值我用的是实数编码直接把一个解表示成形状为 (num_gen, T) 的二维数组。比如 3 台机组、24 时段就是 3x24 的矩阵每个元素代表某台机组在某个时段的出力。编码方式要匹配问题性质这是GA设计的第一原则不要为了追求形式复杂而选二进制编码纯属给自己找麻烦。选择操作用锦标赛选择每次从种群中随机抽 k 个个体把目标函数最小的那个选入下一代k 通常取 2 或 3。k 越大选择压力越大收敛越快但更容易早熟k 太小时搜索很慢。我一般从 k2 开始调效果不够再加大。交叉用的是算术交叉两个父代个体在同一位置做线性组合生成子代。公式很直白child1 alpha * parent1 (1 - alpha) * parent2 child2 (1 - alpha) * parent1 alpha * parent2alpha 是 [0,1] 区间的随机数。这种交叉方式的好处是子代天然落在父代张成的超平面内空间平滑过度不会像二进制交叉那样把解打碎成面目全非的组合。变异用高斯变异在原有解上加一个正态分布随机扰动P_new P N(0, sigma)。sigma 的取值很关键我习惯取机组出力区间宽度的 5% 左右比如机组出力范围 100~600MW区间宽度 500MWsigma 取 25 左右。变异幅度太大算法退化成随机搜索幅度太小后期几乎无法跳出局部最优。变异后必须立刻把个体 clip 到 [Pmin, Pmax] 区间内否则会产生大量无意义解。还需要加精英保留策略每一代把当前目标函数最小的个体原封不动复制到下一代。没有这一步就算种群整体收敛最优个体也可能在变异中被破坏掉导致历史最好的解丢失。这是一个低成本、高收益的操作强烈建议写上。2.3 约束处理罚函数系数怎么调才有好效果经济调度里的约束可以分成三类出力上下限这种硬约束直接用 clamp 处理等式约束发电总和负荷网损爬坡约束这种跨时段耦合约束。后两者我用罚函数方式处理。罚函数的思想很简单允许不可行解参与进化但把违反程度变成惩罚项加到目标函数上让明显违反约束的解在进化过程中被自然淘汰。对等式约束我这样写penalty_balance sum( ( sum(P, axis0) - P_D - P_L )² )注意这里用平方而不是绝对值。平方的惩罚力度随偏差增大而增长得更快能把解“拉”回可行方向另外平方是可导的对GA本身没意义但如果你之后换用差分进化或粒子群这个平滑性质同样适用。爬坡约束的惩罚类似penalty_ramp sum( max(0, |P[:, t] - P[:, t-1]| - Ramp_i)² )这里 max(0, ...) 保证只在超出爬坡速率时才产生惩罚。注意第一个时段没有历史值不需要检查切片时要小心不要越界。惩罚系数到底取多少这是新手最容易纠结的地方。系数太小时约束形同虚设优化结果成本虽然低但根本不能用系数太大时任何违反约束的解都被压得毫无竞争力种群过早集中在可行域的一小片区域搜索效率反而下降成本优化也会被约束压制得过多。我的经验是先跑一次不带惩罚的初始化观察目标函数里成本项的数值量级然后让惩罚系数与成本量级同数量级或者大1~2个数量级。比如成本项总量级在 2万~3万左右惩罚系数取 100 就很合适。如果发现收敛曲线一直很高且降不下来大概率是惩罚系数偏大目标函数被约束项主导了如果收敛后仍有解明显违反约束那就调大系数再跑一轮。3. Python实现从零写一个能跑的GA经济调度程序3.1 代码结构设计参数区、目标函数、遗传操作、主循环我在实际写代码时习惯把程序拆成四块参数区、目标函数区、遗传操作区、主循环区。参数区单独放方便反复调整机组参数、负荷曲线和GA参数。目标函数区只负责“给定一个个体计算目标值”不用管进化逻辑这样调试时可以单独测试目标函数是否正确。遗传操作区实现选择、交叉、变异、精英保留。主循环区就是不断进化 N 代并记录收敛曲线。整个程序只依赖 numpy不需要额外安装复杂库Python 3.8 以上的环境都能直接跑。随机种子建议固定下来比如 seed42这样每次运行结果一致方便复现和调试。如果要多试验几组参数换掉 seed 再跑即可。3.2 关键代码逐段解读下面是核心代码我尽量写得紧凑但完整。为了展示核心逻辑这里只保留必要的部分直接复制就能跑通一个三机24时段的算例。import numpy as np # 种子固定方便复现 np.random.seed(42) # ---------- 参数区 ---------- num_gen 3 # 机组数量 T 24 # 调度时段数 # 机组成本系数 a, b, c元/MW^2h, 元/MWh, 元/h a np.array([0.001562, 0.001940, 0.004820]) b np.array([7.9200, 7.8500, 7.9700]) c np.array([561.0, 310.0, 78.0]) # 机组出力上下限MW Pmin np.array([100.0, 80.0, 50.0]) Pmax np.array([600.0, 400.0, 200.0]) # 爬坡速率限制MW/h统一用Ramp表示 Ramp np.array([120.0, 100.0, 80.0]) # 24小时负荷曲线MW P_D np.array([ 600, 580, 560, 540, 530, 560, 620, 680, 750, 800, 830, 850, 870, 860, 840, 820, 800, 760, 700, 640, 620, 620, 660, 640 ], dtypefloat) # B系数矩阵反映输电损耗 B np.array([ [0.000030, 0.000010, 0.000005], [0.000010, 0.000025, 0.000008], [0.000005, 0.000008, 0.000020] ]) # GA参数 pop_size 60 # 种群规模 max_iter 500 # 迭代次数 pc 0.90 # 交叉概率 pm 0.10 # 变异概率 k_tournament 2 # 锦标赛选择规模 penalty_coeff 100.0 # 罚函数系数 # ---------- 目标函数成本 网损 罚函数 ---------- def calc_loss(P): 计算输电损耗P形状为 (num_gen, T) 或 (num_gen,) P_T P.T # 转成 (T, num_gen) return np.array([P_T[t] B P_T[t] for t in range(P_T.shape[0])])这里有一个细节需要注意P 是二维数组时B矩阵乘法要做成每个时刻独立计算。上面用列表推导式逐时刻计算网损保证每个时刻的 P_L 都是标量。如果直接用二维数组做整体矩阵运算numpy会按广播规则把最后两个维度相乘结果维度混乱很难调试。躲过这个坑最好的方式就是先转置、逐行计算代码简单还有可读性。def objective(P): 返回总目标值 总成本 网损相关惩罚 爬坡约束惩罚 # 总发电成本所有机组、所有时段累加 cost np.sum(a[:, None] * P**2 b[:, None] * P c[:, None]) # 输电损耗 loss calc_loss(P) # 等式约束惩罚每个时段发电总和与负荷网损的偏差 balance np.sum(P, axis0) - P_D - loss penalty_balance np.sum(balance**2) # 爬坡约束惩罚相邻时段出力差的绝对值不得超过Ramp diff np.abs(P[:, 1:] - P[:, :-1]) over np.maximum(0, diff - Ramp[:, None]) penalty_ramp np.sum(over**2) return cost penalty_coeff * (penalty_balance penalty_ramp)注意切片里的 Ramp[:, None]它的作用是把 Ramp 从形状 (3,) 变成 (3,1)这样才能和形状为 (3, T-1) 的 diff 数组做广播运算。这种 numpy 广播的细节我踩过好几次错一次就整个结果全错而且报错信息还不太明显。如果你发现爬坡惩罚一直是 0多半就是这里出了问题。# ---------- 遗传操作 ---------- def init_population(): 在[Pmin, Pmax]区间内随机生成初始种群 pop np.random.uniform(lowPmin[:, None], highPmax[:, None], size(pop_size, num_gen, T)) return pop def tournament_selection(pop, objs, kk_tournament): 锦标赛选择随机抽k个个体保留目标函数最小那个 idx np.random.choice(pop_size, sizek, replaceFalse) best idx[np.argmin(objs[idx])] return pop[best].copy() def arithmetic_crossover(p1, p2): 算术交叉生成两个子代个体 alpha np.random.rand() child1 alpha * p1 (1 - alpha) * p2 child2 (1 - alpha) * p1 alpha * p2 return child1, child2 def gaussian_mutation(ind): 高斯变异变异幅度取区间宽度的5% sigma (Pmax - Pmin) * 0.05 noise np.random.normal(0, sigma[:, None], sizeind.shape) return np.clip(ind noise, Pmin[:, None], Pmax[:, None]) def clamp_population(pop): 把整个种群个体限制在可行区间内 return np.clip(pop, Pmin[:, None], Pmax[:, None])生成初始种群时用 Pmin[:, None] 和 Pmax[:, None] 保证维度匹配否则 numpy 会把广播搞错。这是新手写GA时最常见的报错点之一写上这个细节能帮读者省去调试时间。# ---------- 主循环 ---------- history [] pop init_population() pop clamp_population(pop) elite None elite_obj float(inf) for gen in range(max_iter): objs np.array([objective(ind) for ind in pop]) # 精英保留 cur_elite_idx np.argmin(objs) if objs[cur_elite_idx] elite_obj: elite_obj objs[cur_elite_idx] elite pop[cur_elite_idx].copy() history.append(elite_obj) new_pop [] while len(new_pop) pop_size: # 选择两个父代 p1 tournament_selection(pop, objs) p2 tournament_selection(pop, objs) # 按概率交叉 if np.random.rand() pc: c1, c2 arithmetic_crossover(p1, p2) else: c1, c2 p1.copy(), p2.copy() # 按概率变异 if np.random.rand() pm: c1 gaussian_mutation(c1) if np.random.rand() pm: c2 gaussian_mutation(c2) new_pop.extend([c1, c2]) pop clamp_population(np.array(new_pop[:pop_size])) # 精英个体放回种群第一位置防止丢失 pop[0] elite主循环里最值得强调的是精英个体直接覆盖种群第一位置。每一代我们都把历史最优个体放回种群这个操作几乎不增加计算量但对收敛效果有立竿见影的作用。否则进化过程会自动丢失历史解。运行结束后看收敛曲线import matplotlib.pyplot as plt plt.plot(history) plt.xlabel(Generation) plt.ylabel(Best Objective (cost penalty)) plt.show()输出最优出力方案best elite print(最优出力矩阵形状, best.shape) print(总目标值, elite_obj) print(总成本, np.sum(a[:, None] * best**2 b[:, None] * best c[:, None])) print(网损, calc_loss(best))另外目标函数里返回的是 cost penalty如果要看纯粹的发电成本要单独算一次。很多新手直接打印 elite_obj 当成成本结果数值里带着惩罚项自己还看不懂为什么比理论值大那么多。3.3 一个可以直接复现的算例三机组24时段调度我特意选了 3 台机组、24 个时段这个规模计算量小GA 大概 5 秒内就能跑完非常适合反复修改参数做实验。机组参数和负荷曲线已经在上面代码里给全了B 矩阵也写死了。如果你想换成真实的机组数据只需替换 a、b、c、Pmin、Pmax、Ramp 以及负荷曲线数组即可其他代码不用动。以 seed42、pop_size60、max_iter500 这个配置跑下来我得到的输出大致是初始随机种群的最优目标值在 30000 以上含惩罚项最终收敛到 29000 上下。其中总发电成本大约 28500 左右网损 15~20MW爬坡约束的惩罚项几乎为 0等式约束的偏差也在个位数 MW 量级。如果你跑出来收敛曲线的终点在 50000 以上先检查罚函数系数是不是太大如果惩罚项不为 0 但曲线已经平了说明系数太小或种群早熟了按第 4 节里的思路排查。4. 实验结果解读与避坑经验4.1 收敛曲线能告诉你什么收敛曲线是调试GA最直观的窗口。正常情况是前 50~100 代快速下降中期缓慢下降后期趋于水平。如果前 100 代一点没降大概率是初始种群生成出了问题比如所有个体集中在一个很窄的区域如果曲线一直在缓慢下降但没有平台期可能是种群规模太小、搜索能力不足可以加大 pop_size 到 80~100如果曲线已经平了但目标值里惩罚项依然很大说明当前参数搜索不到可行解需要调大惩罚系数或者增大变异幅度让搜索范围扩展开。我自己的经验是先用小种群40快跑 200 代验证代码逻辑正确与否确认目标函数没问题之后再上 60~100 的种群和 500 代做正式实验。直接上大参数跑一旦结果不理想排错成本很高。4.2 我踩过的几个坑和解决办法第一个坑是 numpy 广播维度搞混。P 是 (num_gen, T) 的二维数组P[:, 1:] 取的是所有机组、除去第一个时段的数据形状是 (num_gen, T-1)而 Ramp 是 (num_gen,)如果不加 [:, None] 直接相减numpy 会报错或按错误形状广播。这类问题最好通过给目标函数打印几个中间变量来检查比如 print(P.shape, Ramp[:, None].shape)。第二个坑是初始种群没做 clamp。因为是在 [Pmin, Pmax] 区间均匀采样理论上不会越界但交叉和变异之后很可能越界越界的个体全部走进罚函数种群进化就会非常缓慢。一定要在每代结束后统一 clamp不要只 clamp 变异后的个体。第三个坑是罚函数系数没有根据目标函数量级调整。有人直接把系数设为 10000结果成本项和惩罚项差距太大GA 把重心全放在满足约束上成本优化几乎失效有人设成 1结果优化出一堆严重违反爬坡约束的方案。我的通用做法是第一遍跑完打印一下 history 里前几代的目标值评估成本项和惩罚项的比例让两者大致在一个数量级。第四个坑是随机种子导致的“假实验”结果。固定 seed 可以复现但只跑一个 seed 给出的最优解可能偏好运或偏差。我建议每个参数组合至少跑 5 个不同 seed取目标值最低的那次的结果同时观察多次运行之间的方差。如果方差太大说明迭代次数不够或种群早熟。第五个坑是网损计算方向搞反。B系数法的传统形式是 P_L P^T B P如果直接拿二维 P 去套公式会得到一堆逐时刻的网损值而不是总网损。上面代码里我专门转置一次再逐行计算就是为了避免这种情况。4.3 基于这套代码还能做哪些扩展这套GA框架的扩展性比我预想的要好。我后来在这个基础上做过几件有实际价值的事。第一是加入阀点效应在成本函数里加正弦项 dsin(e(Pmin-P))这是一个典型的非凸项传统λ迭代基本束手无策但GA只需要改一行目标函数就能跑。第二是改成多目标优化把污染物排放作为一个独立目标用带拥挤度距离的NSGA-II框架替代这里的锦标赛选择就能同时输出成本-排放的帕累托前沿。第三是引入新能源的不确定性把风电出力作为一个随机变量参与功率平衡通过场景抽样将随机约束转化为多个确定性场景GA的群体搜索天然适配这种场景扩展。如果你想跟传统方法做个对比可以在同一套模型上用 scipy 的 SLSQP 接口求解对比两者的解质量与耗时这也是很多论文里常用的分析思路。最后分享一个我自己使用中的体会遗传算法的参数整定没有万能的黄金组合但有一个值得长期使用的调试习惯——先固定门槛条件把约束简化、参数缩小跑通一个小规模算例确认代码无误再逐层加约束、加规模。这套代码里的罚函数系数、变异概率、种群规模彼此牵制一次只调一个变量记录下每一次调整带来的收敛曲线变化会比盲目组合参数高效得多。如果你照着我给的算例跑一遍再替换成自己的机组数据和负荷数据应该很快就能感觉到所谓“智能优化算法求解经济调度”代码本身并不复杂真正值钱的其实是你对问题约束的理解程度以及怎么把工程经验转化成算法参数的选择。