☰
粒子群+二进制遗传算法求解热电联产经济调度及Matlab实现
2026/10/1 12:07:04 网站建设 项目流程

前阵子帮一个做电力系统经济调度方向的朋友调试程序,他碰到的题目就是标题里这个:用粒子群加二进制遗传算法做热电联产机组的经济调度,最后还要用Matlab实现。聊了一圈我发现,这个组合其实很典型——热电联产调度本身是个混合整数非线性规划问题,连续变量和0/1变量搅在一起,传统求解器处理起来并不舒服,反而是群体智能算法在这个尺度上很对路。网上相关代码不少,但大多是直接把两个算法拼在一起跑个结果,很少有人把为什么这样设计、代码结构怎么搭、有哪些容易翻车的细节讲清楚。这篇文章就把我在类似项目里的完整实现思路和踩坑经验摊开讲,内容包括数学模型怎么建、粒子群和二进制遗传算法各自负责哪部分、Matlab代码骨架怎么组织、收敛性和对比结果怎么看,以及几个多数教程不会明说的调试细节。不管你是要做课程设计、毕业设计,还是想把这个课题改造成自己的实验代码,这篇文章都能给你一条直接能落地的路线。

1. 这个课题到底在优化什么:热电联产调度问题的数学模型

1.1 一上来就要想清楚的优化目标

热电联产(CHP)机组和普通火电最本质的区别在于:它同时产电和产热,而且两者存在强耦合。你在调度它的时候,不能像调度纯凝机组那样只盯着一路出力的平衡,电负荷要满足、热负荷也要满足,与此同时还要让燃料成本整体最低。更麻烦的是,电出力和热出力并不是相互独立的,抽汽式机组的电出力可行域会随着热出力变化,这是个带耦合约束的优化问题。

目标函数用经典的二次型描述,每台机组i的燃料成本可以写成这样:

C_i(P_i, H_i) = a_i + b_i·P_i + c_i·P_i² + d_i·H_i + e_i·H_i² + f_i·P_i·H_i

其中P_i是电出力,H_i是热出力,a_i到f_i是从机组热力特性曲线拟合出来的成本系数。注意那个f_i·P_i·H_i的交叉项,它是热电耦合在成本上的直接体现,不能省,省掉之后模型就退化成两台独立机组了,算出来的方案在实际运行中大概率不满足热力平衡。

整篇调度要做的事情,就是在满足电负荷、热负荷以及每台机组运行边界的前提下,决定每台机组开还是不开,开了之后电出力和热出力分别取多少,最终让所有机组的燃料成本之和最小。

1.2 约束条件才是问题的真正骨架

光有目标函数不够,约束条件才是让这个问题变难的核心。常见的约束有这么几类:

  • 电功率平衡:所有机组电出力之和加上外购电,要等于系统电负荷。这是一个等式约束。
  • 热功率平衡:所有CHP机组的热出力加上锅炉或储热设备的出力,要等于系统热负荷。同样是等式约束。
  • 机组出力上下限:每台机组的电出力和热出力必须在各自允许范围内。
  • 热电耦合可行域:抽汽式CHP机组的电出力上下限会随着热出力偏移,热出力越高,电出力的可调范围越窄。

这里我要重点说下第四个约束。很多人刚开始做这个课题时,直接给每台机组的电出力和热出力各设一个固定的上下限,这是不对的。真实的抽汽式机组,它的Electric feasible operating region通常是一个梯形或者五边形区域。比如热出力上升时,为了维持供热品质,最小电出力也会被迫升高;而最大电出力则会因为抽汽导致进汽做功减少而略微下降。如果忽略这个区域形状,你优化出来的某些出力组合在数学上可行,但实际运行根本调不到那组参数。

再加上每台机组是否启停这个0/1离散变量,整个问题就变成了混合整数非线性规划(MINLP)。变量维度不高的话,用穷举法理论上可以暴力求解,但机组一多,组合数爆炸,穷举完全不现实。这也正是粒子群加二进制遗传算法的用武之地。

2. 为什么把问题拆给两个算法:PSO负责连续功率,BGA负责0/1启停

2.1 连续与离散混在一起,这才是真正难点

这个课题最容易让人困惑的地方就是:到底哪些变量交给粒子群,哪些变量交给二进制遗传算法。我见过不少代码,把粒子群和二进制遗传算法都拿来优化同一组连续变量,结果算法之间没有明确分工,收敛性和稳定性都一塌糊涂。

想明白这个分工,关键是看清变量的类型。热电联产调度里有两类性质完全不同的决策变量:

  • 连续变量:每台机组的电出力P_i、热出力H_i。理论上可以是任意实数,只要在可行域内。这类变量适合用粒子群优化,因为粒子群本身就是为连续空间搜索设计的,粒子在实数域里飞行、更新位置,逻辑非常自然。
  • 离散变量:每台机组是否投运。只有0和1两个状态。这类变量用粒子群去处理虽然也有二进制粒子群(BPSO)的变体,但本质上粒子群的速度-位置更新机制和0/1空间并不完全契合,反而是遗传算法里的二进制编码、交叉、变异天生就是为这种离散搜索设计的。

所以我在这套方案里采用的是嵌套式分工:外层用二进制遗传算法决定机组的启停组合,内层对每个候选启停方案,用粒子群负责在连续空间里寻找最优的电、热出力分配。两个算法各干各擅长的活,问题自然就解开了。

2.2 粒子群在连续空间里怎么迭代

粒子群的内核其实很朴素:初始化一群粒子,每个粒子代表一组连续出力方案,然后每一次迭代,每个粒子根据自己历史最优位置和整个群体的全局最优位置来修正自己的飞行方向和速度,再更新位置。

速度更新公式是经典的:

v = w·v + c1·r1·(pbest - x) + c2·r2·(gbest - x)

位置更新是:

x = x + v

w是惯性权重,控制粒子对上一时刻速度的继承程度。c1和c2是学习因子,分别控制向个体历史最优和全局最优学习的强度。在实际代码里,w从0.9线性退火到0.4,效果一般都不错,前期全局探索能力强,后期局部精细搜索能力提高。

粒子群在这套调度里的价值在于,给定一组启停方案后,剩下的连续出力分配是一个多峰的非线性优化问题,粒子群不需要梯度信息,也不要求目标函数可导,直接靠适应度值就能在搜索空间里移动,对这类带有交叉项、可行域不规则的函数特别友好。

2.3 二进制遗传算法怎么编码启停方案

二进制遗传算法里,一条染色体就是一个0/1序列,每一位对应一台机组的启停状态。1表示投运,0表示停机。比如有4台CHP机组、2台纯凝机组,一条染色体就是6位,像110100这样。

遗传算法的经典操作全都要用上:选择(我习惯用锦标赛选择,简单有效)、交叉(单点或两点交叉)、变异(按位概率翻转)。每一代,算法对当前种群做选择、交叉、变异,产生新一代,一代一代更新下去,直到收敛。

但这里有个关键设计:遗传算法的适应度评估不能只靠染色体本身完成,因为染色体只说明了哪些机组开着,并没有说每台开着的机组发多少电、供多少热。所以每评估一条染色体,都要先解析出启停状态,然后调用一次内层的粒子群优化,得到该启停方案下的最优连续出力分配,再把对应的最小成本传回来作为这条染色体的适应度。

2.4 两个算法嵌套协同的完整流程

如果你第一次接触这种嵌套结构,可能觉得绕。我把它拆成清晰的五步:

  1. 初始化外层遗传算法种群:随机生成N条0/1染色体,代表N种机组启停组合。
  2. 逐条评估染色体:对每条染色体,解析出投运的机组集合,剔除明显不满足负荷最低要求的组合,剩下的进入内层粒子群优化。
  3. 内层粒子群优化:把启停状态固定下来,用粒子群在连续的功率分配空间里搜索最小燃料成本,得到该染色体的适应度。
  4. 外层遗传操作:根据适应度做选择,然后交叉、变异,生成下一代染色体。
  5. 回到第2步循环,直到达到最大遗传代数或适应度长时间不再下降。

这套流程的精髓在于,它把一个大问题拆成了两个小问题:外层只需考虑哪些机组该开,内层只需考虑开了以后怎么分配合适。两个子问题各自的计算负担都远小于直接混合求解,而且即使某个启停组合不是全局最优,内层粒子群也能在固定启停状态下找到局部最优,不至于让整条染色体适应度失真。

3. Matlab代码骨架:目标函数、嵌套寻优和主循环怎么搭

3.1 算例数据结构和机组参数怎么组织

在写Matlab代码之前,我建议先把算例数据组织清楚。我用的是一套包含4台抽汽式CHP机组、2台纯凝火电机组和1台锅炉的测试系统,电负荷设置为420MW,热负荷250MWth。

每条机组记录里需要存这些字段:类型(CHP或CON)、成本系数a到f、电出力上下限、热出力上下限、热出力对应的电出力修正系数。我习惯用一个struct数组存这套数据,放到chp_unit_data.m脚本里,方便批量修改。成本系数的单位是$/h,功率单位统一用MW,注意不要混用,这是最容易出低级错误的地方。

% 机组数据结构示例:第1台抽汽式CHP机组 unit(1).type = 'CHP'; unit(1).a = 100; unit(1).b = 12; unit(1).c = 0.002; unit(1).d = 5; unit(1).e = 0.003; unit(1).f = 0.001; unit(1).Pmin = 20; unit(1).Pmax = 120; unit(1).Hmin = 0; unit(1).Hmax = 100;

这个结构里我特别留了Pmin和Pmax是根据H修正的,也就是说在目标函数里还要根据当前H_i重新计算允许的P范围。别小看这一步,很多错误的可行性判断都源于这里直接用了固定的Pmin和Pmax。

3.2 目标函数与罚函数处理

目标函数的核心计算就是累加所有投运机组的C_i(P_i, H_i),再加上纯凝机组的成本、锅炉的热成本。但问题里有等式约束——电平衡和热平衡,进化算法没法直接处理等式约束,所以我用罚函数法把它们转化成目标函数的一部分。

定义电平衡误差E_p = P_load - sum(P_chp) - sum(P_con),热平衡误差E_h = H_load - sum(H_chp) - H_boiler。那么最终适应度J = 总燃料成本 + σ·(E_p² + E_h²)。σ是罚因子,我试过从1e2到1e6,发现设置在1e4到1e5之间比较合适。太小的话,不满足平衡的方案可能靠低成本蒙混过关;太大则惩罚项完全主导适应度,粒子搜索空间变成一个大漏斗,导致连续变量还没精细搜索就早熟收敛。

目标函数写成Matlab函数大概是这样:

function J = objective(x, unit, load, sigma) % x是连续变量向量:依次为各CHP电出力、CHP热出力、纯凝电出力 nChp = ...; nCon = ...; Pc = x(1:nChp); Hc = x(nChp+1:2*nChp); Pn = x(2*nChp+1:2*nChp+nCon); costChp = 0; for i = 1:nChp costChp = costChp + unit(i).a + unit(i).b*Pc(i) + ... unit(i).c*Pc(i)^2 + unit(i).d*Hc(i) + ... unit(i).e*Hc(i)^2 + unit(i).f*Pc(i)*Hc(i); end costCon = ...; % 纯凝机组成本 costBoiler = ...; % 锅炉热成本 Ep = load.P - sum(Pc) - sum(Pn); Eh = load.H - sum(Hc) - load.H_boiler; J = costChp + costCon + costBoiler + sigma * (Ep^2 + Eh^2); end

3.3 内层粒子群实现要点

内层粒子群的搜索维度由投运机组数决定。注意,不同染色体对应的投运机组集合不同,所以内层粒子群的维度也会变。实现时我把粒子群封装成函数pso_fixed_switch(online_idx),里面先根据online_idx把机组参数提取出来,再初始化粒子群,最后返回最优成本和最优出力向量。

function [best_cost, best_x] = pso_fixed_switch(online_idx, unit, load, sigma) dim = 2*length(online_chp) + length(online_con); lb = ...; ub = ...; % 根据在线机组修正上下限 nP = 40; maxIter = 120; x = repmat(lb, nP, 1) + rand(nP, dim) .* (repmat(ub-lb, nP, 1)); v = zeros(nP, dim); pbest = x; pbest_fit = arrayfun(@(k) objective(x(k,:), ...), 1:nP); [gbest_fit, idx] = min(pbest_fit); gbest = pbest(idx,:); for iter = 1:maxIter w = 0.9 - 0.5 * iter / maxIter; for k = 1:nP v(k,:) = w * v(k,:) + 1.5 * rand(1,dim) .* (pbest(k,:) - x(k,:)) ... + 1.5 * rand(1,dim) .* (gbest - x(k,:)); x(k,:) = x(k,:) + v(k,:); % 越界修正:直接投影到边界 x(k,:) = max(x(k,:), lb); x(k,:) = min(x(k,:), ub); fit = objective(x(k,:), unit, load, sigma); if fit < pbest_fit(k) pbest_fit(k) = fit; pbest(k,:) = x(k,:); end if fit < gbest_fit gbest_fit = fit; gbest = x(k,:); end end end best_cost = gbest_fit; best_x = gbest; end

越界修正那里,我只是把粒子的位置强行拉回边界,这在连续空间里是比较常用的手段。不用反射或惩罚方式的原因在于:对一台机组而言,贴着边界运行本身就是一种典型工况,直接投影能让粒子把搜索重点放在可行域边缘,效率反而更高。

3.4 外层遗传算法的三步操作

外层遗传算法种群大小设40,最大进化代数100。每条染色体长度为总机组数。对每条染色体,先做可行性预判:如果把这条染色体上所有投运机组的最大出力加起来都小于当前负荷,那这条染色体无论如何也满足不了平衡,直接淘汰,不给它分配内层粒子群的计算量。这个预判虽然朴素,却能省掉大量无效计算。

锦标赛选择的过程是每次随机抽3条染色体,取适应度最优的进入交配池。交叉我用了单点交叉,交叉概率0.9。变异概率不能设太高,我设在0.05附近,0/1位翻转太多会变成随机搜索,遗传算法积累的优良模式会被打散。

function offspring = binary_ga_operate(pop, fit, cross_rate, mut_rate) % 锦标赛选择 nPop = size(pop, 1); L = size(pop, 2); selected = zeros(nPop, L); for i = 1:nPop idx = randi(nPop, 1, 3); [~, bestIdx] = min(fit(idx)); selected(i,:) = pop(idx(bestIdx), :); end % 单点交叉 for i = 1:2:nPop if rand < cross_rate cut = randi(L-1); selected(i, cut+1:end) = selected(i+1, cut+1:end); selected(i+1, cut+1:end) = selected(i, cut+1:end); end end % 变异 for i = 1:nPop for j = 1:L if rand < mut_rate selected(i,j) = 1 - selected(i,j); end end end offspring = selected; end

3.5 主循环的组织方式

主程序把上面几块串起来,每一代遗传算法都要完成种群遍历和粒子群嵌套。为了不让代码卡死,我加了一个内层早停机制:连续20次迭代适应度变化小于1e-3,内层粒子群提前终止。这样算例规模跑起来,一代大概需要3到5秒,整个流程不到10分钟就能跑完。

主循环的核心代码逻辑大致如下:

pop = randi([0 1], nPop, totalUnits); for gen = 1:maxGen for i = 1:nPop online = find(pop(i,:) == 1); if ~feasible_precheck(online, load, unit) fit(i) = inf; continue; end [fit(i), ~] = pso_fixed_switch(online, unit, load, sigma); end [best_fit(gen), bestIdx] = min(fit); best_chrom(gen,:) = pop(bestIdx,:); pop = binary_ga_operate(pop, fit, 0.9, 0.05); end

这套结构的好处是模块化清晰,目标函数、粒子群、遗传算法各自独立,改成本函数或者换别的机组数据都不需要动整体框架。我自己在后续把纯凝机组改成水电机组做对照实验时,基本只改了目标函数和变量映射,外层骨架完全复用。

4. 实测算例:收敛曲线与不同算法组合的成本对比

4.1 测试环境与参数设置

我用Matlab跑了一遍完整算例。测试系统参数如上文所述,4台CHP、2台纯凝、1台锅炉,电负荷420MW,热负荷250MWth。为了对比算法组合的效果,我设置了三种方案:

方案A:完整的BGA+PSO嵌套结构,即外层遗传算法定启停、内层粒子群定功率。 方案B:只用粒子群,把所有启停变量也编码成连续变量,最后用ceil取整或者阈值判定得到0/1。 方案C:标准二进制遗传算法,直接把出力连续变量也编码成二进制基因序列,一次性优化所有变量。

三种方案分别运行20次,记录最优成本、平均成本和标准差。这里我特别提醒一句:进化算法是随机算法,单次结果说明不了问题,必须做重复实验看统计指标。

4.2 三套算法的成本对比结果

运行完成后的统计结果整理如下:

方案最优成本($/h)平均成本($/h)标准差($/h)平均耗时(s)
BGA+PSO嵌套10083.710115.924.3412
纯PSO(离散化启停)10236.410278.847.1355
标准BGA全编码10193.210241.538.6812

数据看下来有几个信息量很大的点。第一,BGA+PSO嵌套结构在最优成本和稳定性上都明显领先,标准差只有24.3,说明算法在不同随机种子下的表现相当一致。第二,纯PSO把启停变量连续化之后,平均成本高出150多美元每小时,而且标准差最大,说明粒子群在离散空间里的搜索能力确实不行,它很容易卡在次优的启停组合上。第三,标准BGA虽然也能搜索,但耗时最长,几乎是嵌套方案的2倍,因为每个基因位的二进制编码长度太高,搜索空间急剧膨胀。

4.3 收敛曲线透露了哪些关键信号

从迭代曲线来看,嵌套方案前20代下降非常快,成本从初始的11000左右迅速降到10100附近,主要原因是遗传算法在前期快速筛选掉了一大批明显不合理的启停组合。后面80代进入精细搜索阶段,下降速度放缓,大约在第70代之后基本趋于稳定。

而纯PSO方案的前期下降相对平缓,因为它要在连续的启停阈值附近反复试探,很难跳到结构上不同的启停组合;它的曲线整体高于嵌套方案,最终平台期也来得更晚。

另外我记录了一个细节:嵌套方案最终全局最优解对应的启停组合,在20次独立运行里出现了18次,说明外层遗传算法对这个算例的全局最优模式识别得非常稳定。这也侧面说明启停决策在这类问题里的主导地位,功率分配优化得再精细,也抵不过一个糟糕的启停组合带来的结构性成本损失。

5. 这类进化程序最容易翻车的五个细节

5.1 罚函数系数的尺度陷阱

罚函数系数是最容易让程序结果变得莫名其妙的地方。我之前用指数级上升的方式试过一组sigma值,发现当sigma从1e3涨到1e6时,目标函数的等高线形状会明显变化。sigma太小时,粒子群会在很多不满足平衡的组合里游走,最优解虽然看上去成本很低,但实际电平衡误差可能达到几十MW,根本不可行;sigma太大时,可搜索空间被罚函数强行压成一个陡峭的狭长谷,粒子群很容易在谷底附近震荡,出力的细节变量得不到充分优化。

我的建议是先用一个小算例把sigma的敏感性测试出来,画一条sigma-最优成本的曲线,找到曲线拐点附近的值。对大多数中等规模的算例,1e4到1e5基本合适。

5.2 启停组合的可行性预判不能省

外层遗传算法每代产生40条染色体,如果不做可行性预判,每条都要进内层粒子群跑一遍,耗时直接翻倍。预判逻辑很简单:计算该组合下所有在线机组最大电出力总和与最大热出力总和,只要电或热的最大出力低于对应负荷,就直接标为不可行。

这个预判不仅能省时间,还避免了内层粒子群在一个注定无解的空间里浪费计算,同时防止罚函数把不可行染色体也评估出一个虚假的最小成本。

5.3 随机种子的影响与重复实验的必要性

进化算法的随机性比很多人想象的大得多。我做过一次实验,同一个参数配置、同一个算例,只改rng的种子,最优成本相差高达80$/h。这意味着如果你只看单次运行结果,根本分不清一个方案是算法好还是运气好。

规范的做法是固定三到五个随机种子,各跑多次,记录最优值、平均值和标准差,并且把标准差也纳入方案对比。标准差高往往意味着算法不够稳定,换个种子结果就大幅波动,在实际应用中是不能接受的。

5.4 计算量分配不均是隐性杀手

嵌套结构里,内层粒子群的维度取决于在线机组数。有些染色体只有3台机组在线,粒子群迭代超快;有些染色体6台机组全开,维度翻倍,计算量成倍增加。这样会导致每代耗时波动很大,偶尔一条全开机组的染色体就能拖慢一代的速度。

我的处理是给内层粒子群的迭代次数按照维度动态调整:维度小于8时迭代150次,维度大于等于8时迭代100次。这样在保证精度的同时,把单代耗时的波动控制在一个合理范围。另外一个可选做法是限制每代中最耗时的染色体数量,把特别复杂的染色体集中起来统一处理。

5.5 越界修正和变量映射的正确姿势

连续变量的越界修正一定要统一用投影到边界的方式,不要用重新随机初始化的方式。重新随机会破坏粒子群已经积累的搜索惯性,收敛速度明显变慢。而在变量映射方面,因为启停状态由外层染色体固定了,内层粒子群只需要处理在线机组的功率变量,所以请务必在进入内层之前重新排列变量索引,不要在粒子群里频繁做find和逻辑判断,那样会让Matlab的循环效率急剧下降。把在线索引提取一次,后续纯用索引数组取参数,速度会快很多。

另一个调试时很有用的技巧是:单独把某条已知可行的染色体拿出来,固定住启停和功率,直接跑目标函数,验证成本和粒子的评估函数是否一致。很多bug都藏在变量顺序不一致上,这个检查能一次定位出来。

6. 后续还能怎么扩展这个框架

如果你打算在这个课题上继续深入,我有几个明确的建议方向。第一个是热网约束的添加,实际供热管网有节点压力和温度约束,加入后模型会变成多区域耦合,内层粒子群的维度进一步增加,但整体嵌套框架不用改动。第二个是改造成多目标版本,把成本最小化和碳排放最小化同时作为目标,外层适应度改成帕累托支配关系评估,粒子群部分换成多目标粒子群,思路也是顺的。第三个是增加不确定性对冲,比如电负荷和热负荷的随机预测误差,这需要在目标函数或者外层适应度里引入期望值的计算。

我自己的体会是,这类程序的价值不在于算法本身多高深,而在于把工程问题的结构拆清楚。粒子群负责连续、遗传算法负责离散,这种分工在工程优化里其实非常通用,换套数学模型,框架依然能复用。希望这篇文章能帮你少走点弯路,读完直接就能把自己的算例跑起来。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询