搞过电力系统经济调度的人应该都有体会：如果只是做纯火电机组的负荷分配，那是个经典凸优化问题，很多现成算法都能解；但一旦把热电联产机组加进来，问题立刻变味了。热电联产机组的电出力和热出力之间存在强耦合，可行运行区域通常是一个多边形，再加上机组启停状态这种离散变量，整个问题就变成一个带复杂约束的混合整数非线性优化问题。这也是为什么我在做这个课题时，最终选择了粒子群算法加二进制遗传算法的混合框架，用Matlab从零把整套代码写出来。这篇文章我就把自己在这个项目里的建模思路、代码结构和踩过的坑完整梳理一遍，给正在做类似研究或者准备课程项目的朋友一个能直接参考的流程。

1. 热电联产经济调度的建模：目标函数和机组可行域

经济调度问题的本质，是在满足电力负荷和热负荷需求的前提下，给每台机组分配电功率和热功率，让总燃料成本最低。先谈建模，因为后续所有算法设计、代码结构甚至调参方向，都取决于数学模型怎么定。

1.1 三类机组的目标函数

CHP系统通常包含三类机组。第一类是纯凝汽式发电机组，只发电、不供热，成本函数写作C_e(P) = a_e*P^2 + b_e*P + c_e。这个二次函数里的系数来自机组热力试验数据，实际做研究很多时候直接采用论文公开数据，如果你是自己做工程项目，建议用最小二乘拟合实测的煤耗曲线。

第二类是热电联产机组，也叫抽凝式机组，它可以在一段连续范围内调配电功率和热功率。因为电功率和热功率共享同一个蒸汽循环，成本函数不能简单拆成“发电成本+供热成本”，而必须写成带交叉项的二次函数：C_chp(P,H) = a1*P^2 + b1*P + c1 + a2*H^2 + b2*H + c2 + f*P*H。交叉项f通常为负，反映“抽汽供热后发电效率变化”的耦合关系。研究项目里常看到有人为了省事把f忽略，结果优化解全部落在可行域极端处，导致结果荒谬。这一点务必注意，建模阶段不能偷懒。

第三类是纯热力机组，只供热，成本函数与电功率无关，写作C_h(H) = a_h*H^2 + b_h*H + c_h。

系统目标函数是全部在线机组的成本之和。我在这里特别强调“在线”两个字：如果一台机组被启停决策判为停机，那么它对应的功率变量必须置零，固定成本项c也应该从目标函数中剔除。否则外层遗传算法会看到“多开一台机组，多负担一个固定成本”的信号，反而抑制正常搜索。

另外，启停成本是否计入是一个需要提前想清楚的问题。经典经济调度只看稳态成本，不看启动和停机费用；但一旦引入了二进制启停变量，算法很容易通过频繁启停低效机组来“钻空子”。我的做法是把启动成本做成一个可选项，默认计入，这样外层BGA才会真正比较“多开一台低效机组”和“切换启停状态”的代价。

1.2 抽凝式机组的非线性可行域

抽凝式机组和纯凝汽机组最大的区别在于，它的电出力和热出力不能独立选择，必须落在二维平面上的一个多边形可行域内。这个多边形通常由最小凝汽流量、最大凝汽流量、低压缸最小进汽比例、抽汽压力上限等热工边界决定，顶点坐标可以通过厂家热平衡图或文献获得。

举例说明，一台CHP机组的可行域可能由以下顶点围成：

• 最大电出力、最小热出力点
• 最大电出力、较大热出力点，对应凝汽流量受限的边界
• 中等电出力、最大热出力点，对应抽汽量上限
• 最小电出力、最大热出力点
• 最小电出力、最小热出力点

严谨的做法是把这些顶点定义成一个N×2的矩阵，按逆时针顺序排列。为什么要强调顺序？因为后面做多边形越界判断时，顶点顺序错了，内外部判定就会反过来。我自己就在这个细节上吃过亏，后面第4章会详细说。

编写代码时对这个区域必须做边界判断，否则算法会给出一组位于可行域之外的电热组合，看起来满足功率平衡，实际该机组根本无法运行。

1.3 功率平衡与机组启停约束

除了每种机组的容量上下限，系统层面还有两个强约束：

• 电功率平衡：所有在线机组电功率之和等于系统电负荷
• 热功率平衡：所有在线CHP机组热出力与锅炉/热力机组热出力之和等于系统热负荷

机组启停变量u_i取0或1，同时P_i和H_i必须满足u_i=0时P_i=0、H_i=0。这个约束在优化里等价于把连续变量的可行范围改成[0, u_i*Pmax]，而不是简单地在目标函数里加一个“停机成本惩罚”。

在这些约束中，最难处理的就是CHP多边形。因为它不是简单上下界，而是一个二维凸包约束。喜欢用罚函数的同学要注意，单纯用inpolygon判断点是否在多边形内部，只能给出“是/否”，无法给出越界程度，罚函数梯度对算法毫无帮助。这里面的细节我在第4章单独展开。

2. 为什么选择粒子群+二进制遗传算法的混合框架

2.1 单用PSO处理离散变量的问题

粒子群算法在连续优化上的表现大家有目共睹，代码简单、收敛快、对初值不敏感。但这个案例里存在机组启停变量，这是一个离散组合问题。直接套用经典PSO会非常别扭：位置更新公式里连续速度项对二进制变量没有意义，阈值化的方式又会让粒子频繁在0和1之间震荡，很难稳定收敛。更麻烦的是，启停状态的组合数是2的n次方，纯PSO的那套“向个体最优和全局最优学习”机制，在离散搜索空间里很容易早熟，稍不留神就全部聚集在一个局部组合上，再也跳不出来。

2.2 混合算法的分工：BGA管启停，PSO管功率分配

所以我最终采用了嵌套的混合结构：外层是二进制遗传算法（BGA），负责优化启停组合；内层是粒子群算法，负责在给定启停方案下求连续功率分配的最优解。

BGA的每个个体是一条长度等于机组数量的二进制染色体，每一位代表一台机组的状态：1表示开机，0表示停机。对每一条染色体，先把开机机组筛选出来，然后在内层运行PSO，把电功率、热功率作为连续变量寻优。内层PSO返回的最小燃料成本，就作为外层染色体的适应度值。

选择这种分工的核心原因是：遗传算法有交叉、变异两种组合搜索算子，在离散组合空间里的效果比PSO阈值化做法稳定得多；而连续功率分配又是一个典型的多维、带约束非线性优化问题，正好是PSO的主场。两者各干各擅长的活，不需要把问题硬塞给单一算法。

2.3 整体计算流程

混合算法的完整流程是：

1. 随机生成N条二进制染色体，构成BGA初始种群。
2. 对当前种群中的每条染色体，执行内层PSO：
   • 根据染色体解码出开机机组集合；
   • 在开机机组的功率范围内随机初始化M个粒子；
   • 计算目标函数和罚函数；
   • 更新个体最优、全局最优、速度和位置；
   • 达到PSO最大代数或收敛阈值后，返回最优成本和最优功率分配。
3. 把PSO返回的最优成本作为该染色体的适应度。
4. BGA执行锦标赛选择、两点交叉、位变异，生成新一代种群。
5. 重复步骤2到4，直到达到最大代数或适应度连续多代不再下降。

这种嵌套结构写起来不复杂，但计算量比单一算法大不少。我在第6章会给出几个工程化的加速技巧，比如缓存同一启停组合的内层结果。

3. Matlab代码实现：从数据定义到迭代更新

这一部分直接给可参考的代码骨架。我用Matlab实现时，把所有逻辑拆成了四个部分：系统数据定义、目标函数、内层PSO、外层BGA，这样每一块都比较好单独测试。

3.1 系统数据与算法参数定义

用一个简化测试系统来演示。系统包含2台纯电机组、2台CHP机组、1台锅炉，电负荷300MW，热负荷150MW。

%% 系统数据 systemData.pD = 300; % 电负荷 systemData.hD = 150; % 热负荷 % 纯电机组参数: [a b c Pmin Pmax] electricUnits = [ 0.21 25.8 0.0 20 150; 0.18 24.5 0.0 25 200 ]; % CHP机组多项式可行域，顶点按逆时针排列 chpPolys{1} = [ 30 0; 60 0; 75 30; 45 70; 20 40 ]; chpPolys{2} = [ 25 0; 70 0; 80 35; 50 65; 15 30 ];

当然，代码里还需要记录每个CHP机组的成本系数，包括P平方项、H平方项、交叉项系数，以及锅炉/热力机组的成本参数。为了避免代码过长占用版面，这里不把40个系数全部列出来，重点是要理解数据结构：用一个总的systemData结构体统一保存所有数据，因为外层BGA和内层PSO函数都要访问它。

算法参数单独放一份，方便以后做参数敏感性实验时只改一个地方：

params.popSize = 40; % BGA种群规模 params.maxGen = 50; % BGA最大代数 params.pc = 0.9; % 交叉概率 params.pm = 0.05; % 变异概率 params.nParticles = 30; % 内层PSO粒子数 params.maxIterPSO = 100; % 内层PSO最大迭代次数 params.wMax = 0.9; params.wMin = 0.4; params.c1 = 2.0; params.c2 = 2.0;

3.2 内层PSO的连续变量寻优

内层PSO接收一个启停状态unitStatus，然后只在开启机组的维度上搜索。一个常见错误是把停机机组的维度也放进去，结果粒子不断被罚函数拽回0，浪费大量迭代次数。正确的做法是动态生成变量维度索引：

function [bestCost, bestX] = innerPSO(systemData, unitStatus, params) nUnits = length(unitStatus); % 统计开机机组 activeIdx = find(unitStatus == 1); nActive = length(activeIdx); if nActive == 0 bestCost = 1e10; bestX = []; return; end % 根据开机机组构建上下界向量 lower, upper % 这一步从 systemData 中读取Pmin/Pmax与CHP多边形最小外接框 ... % 初始化粒子 x = zeros(params.nParticles, nActive); v = zeros(params.nParticles, nActive); for i = 1:params.nParticles x(i,:) = lower + rand(1,nActive) .* (upper - lower); end pbestPos = x; pbestVal = inf(params.nParticles, 1); gbestPos = zeros(1,nActive); gbestVal = inf; for iter = 1:params.maxIterPSO for i = 1:params.nParticles cost = chpCostFunc(x(i,:), systemData, unitStatus, activeIdx); if cost < pbestVal(i) pbestVal(i) = cost; pbestPos(i,:) = x(i,:); end end [gbestVal, idx] = min(pbestVal); gbestPos = pbestPos(idx,:); w = params.wMax - (params.wMax - params.wMin) * iter / params.maxIterPSO; for i = 1:params.nParticles v(i,:) = w*v(i,:) + params.c1*rand(1,nActive).*(pbestPos(i,:)-x(i,:)) ... + params.c2*rand(1,nActive).*(gbestPos - x(i,:)); x(i,:) = x(i,:) + v(i,:); x(i,:) = min(max(x(i,:), lower), upper); end end bestCost = gbestVal; bestX = gbestPos; end

这里的chpCostFunc不是只算燃料成本，而是“燃料成本+约束惩罚”，第4章会细说。

3.3 外层BGA的编码、交叉与变异

BGA部分相对简单。种群初始化时我加了保底开机逻辑，后面会解释原因。染色体解码、适应度评估、进化更新的代码如下：

function [bestChrom, bestCost] = outerBGA(systemData, params) nUnits = 5; % 示例 pop = initPopulation(params.popSize, nUnits, systemData); for gen = 1:params.maxGen fitness = zeros(params.popSize, 1); for i = 1:params.popSize unitStatus = pop(i,:); [cost, ~] = innerPSO(systemData, unitStatus, params); fitness(i) = cost; end % 精英保留 [bestCost, idxBest] = min(fitness); bestChrom = pop(idxBest,:); % 锦标赛选择 nextPop = zeros(size(pop)); nextPop(1,:) = bestChrom; for i = 2:params.popSize idx = tournamentSelect(fitness, 2); parent1 = pop(idx(1),:); parent2 = pop(idx(2),:); child = twoPointCrossover(parent1, parent2, params.pc); child = bitFlipMutation(child, params.pm); nextPop(i,:) = child; end pop = nextPop; end end

两点交叉和位变异的实现很直接：

function child = twoPointCrossover(p1, p2, pc) n = length(p1); child = p1; if rand < pc pt1 = randi([1 n-1]); pt2 = randi([pt1+1 n]); child = [p1(1:pt1), p2(pt1+1:pt2), p1(pt2+1:end)]; end end function child = bitFlipMutation(child, pm) mask = rand(size(child)) < pm; child(mask) = 1 - child(mask); end

为什么用两点交叉而不是单点交叉？因为机组启停问题中，整段基因的连续块往往对应一组功能近似的机组，两点交叉能更好地保留和重组这些模块，比单点交叉更不容易把基因片段打碎。

4. 约束处理与调试：罚函数、多边形判定和初始化技巧

4.1 罚函数系数怎么取

约束处理我会优先用罚函数法，而不是失效个体淘汰法。原因很简单：经济调度问题的可行域在约束边界附近占整体搜索空间的比例并不大，如果只保留可行解，粒子群在初始阶段就很难找到足够多的样本点；罚函数则给不可行解一个“靠近可行域”的梯度和压力，引导搜索逐步进入可行区域。

我的罚函数设计如下：

function val = chpCostFunc(x, systemData, unitStatus, activeIdx) % 根据 activeIdx 还原完整机组的功率向量 P = zeros(nUnits,1); H = zeros(nUnits,1); % ... 填入连续变量 % 1) 计算燃料成本 fuelCost = evalFuelCost(P, H, systemData); % 2) 电/热功率平衡惩罚 penBalance = 500 * (abs(sum(P) - systemData.pD) + abs(sum(H) - systemData.hD)); % 3) CHP多边形越界惩罚 penPoly = 1500 * sum(getPolyViolation(P, H, systemData)); val = fuelCost + penBalance + penPoly; end

罚函数系数怎么定是一个经验活。我的经验是：先不加罚函数跑一次，观察正常燃料成本的数量级，再把惩罚系数设在成本量级的100倍以上。如果系数太低，最后给出的解可能“电热不守恒”，结果没法看；如果系数太高，粒子几乎只在可行域边界做微小移动，优化能力反而被压制。理想的情况是，罚函数在不可行区域形成的梯度方向，能把粒子推回可行域，而不是直接淹没成本信号。

4.2 多边形可行域的判定细节

CHP可行域是多边形，Matlab内置的inpolygon判断点在多边形内外很好用，但只能给0/1结果，不能给越界距离。为了让罚函数连续，我实现了一个点到多边形的最短距离函数。

其中核心是点到线段的距离计算：

function d = pointToSegmentDist(pt, A, B) v = B - A; w = pt - A; t = max(0, min(1, dot(w, v) / dot(v, v))); proj = A + t * v; d = norm(pt - proj, 2); end

然后遍历多边形所有边，取最短距离作为越界惩罚量。这里需要特别注意多边形的顶点顺序，必须是顺时针或逆时针统一方向。我一开始以为顶点顺序不影响“并集”，直接用一个乱序表格，结果getPolyViolation返回的越界量忽大忽小，内层PSO一直震荡。后来我画了个图检查顶点连线才意识到，最短距离函数对顶点顺序极其敏感。建议在代码里写一个简单的自检：把所有顶点画出来，肉眼确认是不是一个封闭有序的多边形。

4.3 初始化与“半可行”解的修正

外层GA初始化时，如果完全随机生成0/1染色体，很容易出现“开机总容量不够负荷”的个体。例如电负荷300MW，随机染色体只开了两台小机组，最大出力才180MW，这种染色体无论如何都不可能满足可行性，内层PSO跑得再多也只会返回一个巨大惩罚值。

我在初始化里加了一个保底逻辑：按额定容量从大到小把机组依次打开，直到累计最大电出力不小于电负荷、累计最大热出力不小于热负荷；然后在这个基础上以一定概率翻转部分位，生成多样化个体。这样做能让初始种群基本都具备可行的硬件基础，内层PSO才有意义。

内层PSO的粒子初始化同样不能太随意。我会在生成随机功率后，先按照某个比例缩放到接近电热平衡附近，让粒子起步点在可行域边缘附近，而不是天女散花。这样做能显著加快收敛。

5. 仿真算例：用一个小型CHP系统验证算法效果

5.1 测试系统与参数配置

下面给一个具体的算例，方便你复制复现。系统包含2台纯电机组、2台CHP机组、1台锅炉。负荷为电300MW、热150MW。

各机组的成本数据和可行域我没有全部贴出来，主要是因为这类系数来自公开文献，各家的数据差异挺大；关键是算法框架和参数设置能不能复现。算法参数我按以下配置跑：

5.2 结果对比：纯PSO与BGA-PSO

为了体现混合算法的价值，我做了两组实验：第一组固定所有机组全部开机，只用内层PSO做连续优化；第二组用外层BGA决定启停，再用内层PSO做功率分配。

结果如下表：

从结果看，把低效纯电机组停掉之后，总成本大约下降了12%。这种差距在机组数量更多的系统里会更加明显，因为机组越多，启停组合的决策空间越大，“硬开着一台高煤耗机组”的代价也越高。

这种结果是否可信，要看最后一列成本是不是满足所有约束。我在调试时会额外写一个约束检查函数，把返回的最优解重新算一遍功率平衡误差、多边形越界量，全部打印出来。只要这些误差在允许范围内，才认为这个解是有效的。

5.3 从收敛曲线上能看出什么

跑完实验后，每次迭代都记录当前的适应度。得到的收敛曲线大致有这样的特征：

• 外层BGA前10代成本下降非常快。原因是初始种群经过保底开机，大部分个体都有可行基础，交叉和变异能快速筛选出较优的启停方案。
• 到第20到30代后，适应度下降明显变慢，这时候主要是在做“微调型”搜索，找到一个新的更优启停组合往往需要多次变异。
• 内层PSO在单个启停组合下一般需要30到50次迭代就能稳定到该组合的最优功率分配。

如果看到内层PSO的收敛曲线后期出现“锯齿形跳动”，优先检查速度上限vmax。速度上限设太大，粒子会在可行域边界来回震荡；设太小，又容易早熟。我一般把vmax设为变量范围宽度的10%到20%，效果比较平衡。

6. 我在调参和跑算例过程中积累的几点经验

6.1 嵌套结构的耗时问题与缓存优化

嵌套BGA和PSO最大的现实问题就是计算量。前面算了一下：外层50代、每代40条染色体，每条染色体内层PSO跑100次迭代、每次30个粒子，目标函数评估次数是50×40×100×30等于600万次。再加上多边形投影距离这种计算量偏大的部分，Matlab直接跑会非常慢。

我的加速手段有三个。第一是缓存去重：外层GA的交叉变异会产生很多与父代重复或彼此重复的染色体，用containers.Map把每条二进制染色体的内层优化结果缓存起来，键就是mat2str(unitStatus)，命中直接返回，实测能减少将近一半的重复计算。第二是内层提前终止：PSO连续20次迭代目标值变化小于1e-4就提前跳出，不必跑满100次。第三是先把CHP多边形顶点换算成边向量和法向量，避免在目标函数里重复计算。

缓存代码大致长这样：

cache = containers.Map('KeyType','char','ValueType','any'); key = mat2str(unitStatus); if isKey(cache, key) cost = cache(key); else [cost, x] = innerPSO(systemData, unitStatus, params); cache(key) = cost; end

6.2 早熟和停滞的应对

混合算法跑到后期，最常见的症状是整个BGA种群几乎只剩一两条染色体的变体，多样性严重不足，适应度长期不变化。问题根源在于遗传算法选择的“马太效应”：好个体越选越多，基因越来越趋于一致，交叉算子逐渐失效。

我的做法是加入移民算子：每10代，随机抽取20%个体，按照保底开机逻辑重新生成。这样做既破坏了群体同质化，又保证了新个体具备可行性基础。如果项目对收敛速度要求高，还可以用自适应变异：当种群适应度方差低于某个阈值时，把变异概率临时从0.05上调到0.15。

实测下来，移民算子对这类启停组合问题特别有效，因为它能把一些前期被淘汰的机组重新拉回候选集合，避免算法过早锁定一个并不全局最优的“铁组合”。

6.3 可以继续扩展的方向

这套“BGA管启停、PSO管分配”的框架扩展空间很大。

如果要加入机组爬坡约束和最小启停时间约束，就不能再用一维二进制染色体了，需要在每个基因位额外保存“已连续开机或停机时长”的记忆信息，交叉和变异时做合法性检查，这个改动虽然麻烦，但框架不变。

如果要考虑风电场或光伏并网带来的不确定性，电功率平衡约束变成随机约束，可以结合场景法生成多个风电出力场景，内层PSO对每个场景都要做一次经济分配，外层BGA再用期望成本做适应度。这个思路我在另一个项目里试过，计算量更大，但外层框架完全复用。

如果要做多目标版本，比如同时优化燃料成本和污染物排放，可以把内层PSO换成多目标PSO，外层BGA的适应度改为非支配排序和拥挤度距离。这个扩展在代码结构上比单目标还顺，因为多目标PSO的返回值本身就是一组Pareto解集。

最后分享一点个人体会。做这类研究型代码，最忌讳一上来直接跑最终算例。先把问题缩到最小规模，比如1台纯电机组加1台CHP机组加1台锅炉，用手工枚举穷举验证代码的正确性，再逐步放大。启发式算法框架本身不难，难的是把约束处理得干净、把边界条件都想明白。我自己在CHP多边形投影距离那里就踩过坑，因为顶点顺序写反，整个可行域判定都错了，调了两天才发现。如果你也打算做这个方向，强烈建议先把这些细节调稳了，再去追求更花哨的算法。后面写论文、做对比实验，顺畅程度会是天壤之别。