做综合能源系统调度这几年,有个感受特别明显:只要题目里带P2X,工作量就翻一倍。原因很简单,P2X把原本“各管各”的电、气、热系统牢牢绑定在了一起,风电光伏出力不确定,电解槽、甲烷化设备、储氢罐、热泵一掺和进来,建模维度直接多出几个数量级。更麻烦的是,实际调度往往不只看经济成本,还得兼顾碳排放和可再生能源消纳,这就成了一个典型的多目标优化问题。我用Matlab实现了一套基于多目标退火算法的求解框架,把含P2X综合能源系统的日前调度问题完整跑通,包括模型构建、算法设计、代码实现和调试经验。这篇文章把整个过程复盘一遍,给同样在做综合能源调度、毕业设计或者实际项目的人一个可以直接参考的模板。
我不打算把文章写成纯理论推导。调度问题最终是要落地的,模型再漂亮,算不出来、算得慢、跑出来的结果不可复现,都是白搭。所以我主要围绕“怎么把这个问题用Matlab实现出来”这条主线讲,中间穿插算法设计的理由和调试时踩过的坑,尽量让你看完之后能直接动手改代码。
1. 调度问题建模:P2X让系统变成了“多能耦合”问题
很多初学者拿到题目第一件事就是找算法、抄代码,这是本末倒置。P2X综合能源系统调度难,难的根本不是在“退火”两个字上,而是难在问题本身的建模和约束处理上。模型建不好,再牛的多目标算法都是在一团乱麻里随机蹦跶。
1.1 P2X设备在调度模型里的抽象方式
P2X的本质就是“电转X”,X可以代表氢气、天然气、热能甚至化工原料。在综合能源系统里,最常见的是电转氢(P2H2)、电转气(P2G)、电转热(P2H)这三条路径。
我在建模时不会把每个设备都详细展开到内部电化学过程,那对调度问题没有意义。调度关心的是:输入多少电,经过多长时间,能产出多少对应形态的能源,以及这个转换过程受到什么限制。
举一个典型例子。电解槽(electrolyzer)的输入是电功率,输出是氢气功率,转换关系可以写成一个带效率系数的表达式:
% 电转氢设备模型 H_P2H2 = eta_P2H2 * P_P2H2; % 其中 P_P2H2 为电解槽消耗功率,H_P2H2 为产氢功率 % eta_P2H2 一般为 0.6~0.75这里有一个关键细节:很多P2X设备的效率并不是常数,它随负载率变化。电解槽在低负载时效率明显偏低,启动过程还有最小稳定运行功率约束。但如果你一开始就把效率曲线做成非线性,后面退火算法的搜索空间会变得非常难处理。我的做法是:第一步用恒定效率做基线版本,把整个算法框架跑通,再逐步把效率曲线替换成一族离散工况点,配合查表或插值函数来处理。
设备层面要抽象的信息大概有这几类:输入输出耦合关系、出力上下限、爬坡速率、最小启停时间、运行维护成本系数、启停成本。建议用结构体统一管理,别散落在一堆全局变量里。我习惯在Matlab里这么组织:
% 设备参数结构体示例 dev.P2G.cap = 500; % 额定功率 kW dev.P2G.eta = 0.62; % 转换效率 dev.P2G.ramp = 100; % 爬坡速率 kW/h dev.P2G.cost = 0.015; % 运维成本 元/kWh dev.P2G.pmin = 30; % 最小运行功率 kW dev.P2G.pmax = 500; % 最大运行功率 kW这套结构体接下来会贯穿整个代码,目标函数、约束检查、邻域生成都要用它。
1.2 目标函数与经济环保指标的数值化
调度问题要优化的目标,工程上最常用的就是两个:运行成本最小、碳排放最小。有时候还会加上“可再生能源弃用率最小”,但弃风弃光惩罚通常可以并入成本目标里,用惩罚项体现。
运行成本一般包含四块:向上级电网购电费用、购买天然气费用、各设备的运行维护费用、启停成本。可以写成:
% 目标1:综合成本 f_cost = f_grid_buy + f_fuel + f_om + f_start_stop; % 目标2:碳排放量 f_co2 = co2_grid + co2_gas;网上购电费用要注意电价是分时的,峰谷电价差异很大。P2X设备往往是大功率用户,把电解槽安排在低电价时段运行,是降低运行成本最直接的手段。这个现象会影响最终调度方案里P2X设备的出力曲线。
碳排放计算也要注意口径。电网购电的排放系数不是固定值,不同时段的电网碳排放强度可能不同,因为电网负荷结构在变化。做研究时可以简化为统一系数,但做工程项目时最好拿到当地调度机构发布的分时排放因子数据。
目标函数的量纲和数量级差异也需要注意。成本是很低的值(几千到几万),但碳排放可能是几十吨到几百吨。这两个目标在算法迭代中的贡献度天然不均衡。我的处理方式是在算法内部对两个目标分别做归一化,用各自的基准值去除,得到一个0到1左右的标量,然后再参与Pareto支配判断和接受概率计算。
常见的目标函数形式:
| 目标 | 单位 | 典型数量级 | 说明 |
|---|---|---|---|
| 综合运行成本 | 元 | 10^4 ~ 10^5 | 含购电、燃气、运维、启停 |
| 碳排放总量 | kg | 10^3 ~ 10^5 | 含电网间接排放和燃气直接排放 |
| 弃风弃光率 | % | 0 ~ 20 | 通常折算成惩罚项并入成本 |
1.3 约束条件里最容易被忽略的几个坑
约束条件我专门拎出来讲,是因为调度问题里80%的bug都出在约束处理上。最容易翻车的有这么几个:
第一,能量平衡约束不是简单的“发电=负荷”等式。加入P2X之后,电功率平衡左边要减去P2X设备的消耗:
% 电功率平衡 P_grid + P_CHP + P_wind + P_pv + P_bat_discharge - P_bat_charge ... == P_load + P_P2H2 + P_EB;热功率平衡类似,P2H设备(电锅炉、热泵)产生的热要并入热力母线:
% 热功率平衡 H_CHP + H_GB + H_EB + H_P2H == H_load + H_charge - H_discharge;氢/气平衡则取决于P2X的具体路径。如果是电转氢,还要考虑储氢罐的充放氢平衡:
% 氢气平衡 H_P2H2 + H_storage_discharge == H_load_h2 + H_storage_charge;第二,储氢罐和储电设备的SOC约束是个“时间耦合约束”。SOC是逐时递推的,初始时刻的SOC和末端SOC必须满足设定值,而中间的SOC变化量约束了前后时段的充放功率。这会导致任意一个时段的决策变化都会沿着时间轴传导,用退火算法做邻域扰动时,如果只改动某个时段的设备出力,很容易在SOC约束上违限。
第三,CHP机组的热电耦合可行域。抽凝式汽轮机或者背压式机组的电出力和热出力之间存在联动关系,不能各自独立设定。常见的做法是用一个“可行域多边形”来描述,而多边形约束是非凸的,这在退火算法的邻域搜索里会让问题变得很棘手。
第四,旋转备用约束。系统需要预留足够的可调节容量来应对负荷预测误差和新能源波动。这个约束通常写成不等式形式,但在P2X参与之后,电解槽削减自身功率也可以视作一种“向下可调容量”,模型里可以纳入。
我的建议是:把所有约束写进一个独立的检查函数里,返回总的违限量。这样既能配合罚函数在算法里使用,又方便调试时单独看是哪个约束在报错。
function viol = checkConstraints(x, data) viol = 0; % 依次检查能量平衡、设备上下限、爬坡、SOC、备用约束 viol = viol + elecBalanceViolation(x, data); viol = viol + deviceLimitViolation(x, data); viol = viol + rampViolation(x, data); ... end2. 多目标退火算法设计:为什么是退火,怎么做Pareto求解
算法选型永远不是越新越好,而是看它和问题的匹配程度。P2X综合能源系统调度这个问题的典型特征是:强非线性、非凸可行域、混合整数变量、多目标耦合。在这个组合下,退火算法有它独特的生存空间。
2.1 为什么把退火用于IES调度
综合能源系统调度的主流解法其实分为两大阵营:一类是数学规划方法,以MILP、MINLP为代表,配合CPLEX、Gurobi这样的商业求解器;另一类是元启发式方法,粒子群、遗传算法、差分进化、退火算法都在里面。
数学规划方法的优点是有全局最优的证明保证,但前提是模型必须转化成可求解的标准形式。P2X设备一旦引入非线性的效率曲线、CHP可行域多边形、储氢的复杂状态方程,模型就变成MINLP,求解器的收敛性和耗时都非常不稳定。做学术研究时可以通过分段线性化把它转成MILP,但这是一项非常繁琐的工作,每改一次模型就要重新推导一次线性化方法。
退火算法的优势在于:它几乎不关心目标函数和约束的数学性质。你给它一个黑箱子,它能给你一个还不错的最优解。工程上,尤其是项目前期做方案比选时,这种鲁棒性是非常有价值的。
此外,退火算法的参数比较少,调参空间比遗传算法和粒子群小很多。GA涉及种群规模、交叉概率、变异概率、选择策略、精英保留个数;粒子群涉及惯性权重、个体学习因子、社会学习因子、速度上限。退火算法核心就是三个参数:初始温度、降温系数、迭代次数。这对工业落地来说极其友好。
2.2 从单目标到多目标:Pareto支配与外部存档
单目标模拟退火的逻辑大家都很熟:以一定概率接受差解,随着温度降低,接受差解的概率逐渐收敛到0。但多目标退火要回答一个核心问题:两个目标相互冲突时,怎么判断“谁更优”?
答案是Pareto支配。成本更小且碳排放更小,那这个解支配另一个解;成本小了但碳排放大了,这两个解互不支配。一组互不支配的解构成Pareto前沿,调度人员可以根据实际偏好从前沿上选择一个折中解。
在我的代码里,是把整个退火过程拆成两条线:搜索线用加权目标来引导,存档线用Pareto支配来维护多样性。具体来说,每次迭代给两个目标随机分配一组权重,把多目标函数转换成单目标值进行退火搜索;但在判断解的质量时,仍然把解放进外部存档中,按Pareto支配关系更新存档集合。
用随机权重的好处是,退火搜索的方向会不断变化,从而有机会探索到不同偏好的解区域,不容易陷入某个单一方向上的局部最优。这也是多目标退火(MOSA)里最常用的做法之一。
需要特别说明的是,存档更新里有个“支配删除”的过程。每加入一个新解,要把它和存档里的所有解做比较。如果新解被某个存档解支配,则丢弃;如果新解支配了某些存档解,则删除那些被支配的存档解;如果新解和存档解互不支配,则加入存档。这个过程写起来逻辑清晰,但要注意存档规模膨胀的问题,后面第4节我会详细展开。
2.3 冷却调度、邻域扰动与接受规则
冷却调度是退火算法的灵魂。温度从高温降到低温的过程,决定了算法前期有多少“探索”能力、后期有多少“开发”能力。
我用的降温公式是经典的等比降温:
T = T0 * (alpha ^ iter);T0取100,alpha取0.95,迭代2000次左右,温度大概降到50。这个温度幅度对综合能源调度问题的目标数量级来说是够用的。但注意,如果目标函数做了归一化,初始温度和终止温度的选择范围会随之变化,一般需要根据目标扰动幅度的统计值来定,我建议做一个一两百行的“温度校准”测试:在初始解附近随机扰动,计算目标值变化量的平均值,T0大概取这个平均变化量的5到10倍。
邻域扰动方面,这个问题的决策变量分两类:连续型(设备出力、SOC)和整数型(机组启停状态)。连续型我用高斯扰动,整数型用随机翻转:
% 连续型变量扰动 x_new(i) = x(i) + sigma * randn(); x_new(i) = max(dev.Pmin, min(dev.Pmax, x_new(i))); % 整数型变量扰动:随机选中一个机组,翻转启停状态 x_new(d) = 1 - x(d);注意两点。第一,sigma不能固定,我按温度动态调整,温度高时扰动幅度大,温度低时扰动幅度缩小,这能明显提升收敛效果。第二,相邻时段的设备出力不要独立扰动,否则很容易破坏爬坡约束。我习惯在邻域生成里引入“跨时段联动”,例如同时扰动连续两三个小时的出力,并给它们加上一个共同的偏移量。
接受规则沿用Metropolis准则:
delta = f_new - f_old; if delta < 0 x = x_new; elseif rand() < exp(-delta / T) x = x_new; else % 保持原解 end其中f_new来自随机加权后的综合目标值。多目标场景下,这个接受判定只对“搜索游走”起作用,存档更新永远是独立于接受判定之外的,这样能保证即使搜索过程被卡住,存档里的Pareto解也不会丢失。
3. Matlab代码实现的完整流程
讲完设计和原理,现在说代码怎么落地。我会按构建顺序从输入数据组织讲到主循环输出,每一步都会给出对应的Matlab实现思路,以及我自己实践中踩过的雷。
3.1 输入数据组织与参数初始化
写调度代码的第一步不是写算法,而是搭数据接口。我习惯把案例数据拆成两个文件:一个存系统参数的load_case.m,一个存负荷和新能源预测数据的load_profile.m。这样后期换算例时只需要改这两个文件,算法代码完全不用动。
% load_case.m 返回结构体 data data.dev.CHP = struct('cap', 800, 'eta_elec', 0.42, 'eta_heat', 0.45, ... 'ramp', 200, 'cost', 0.02, 'pmin', 100, 'pmax', 800, ... 'hmin', 80, 'hmax', 500); data.dev.GasBoiler = struct('cap', 400, 'eta', 0.9, 'ramp', 120, 'cost', 0.01); data.dev.P2H2 = struct('cap', 300, 'eta', 0.65, 'ramp', 80, 'cost', 0.02); data.dev.Battery = struct('cap', 200, 'soc_max', 0.9, 'soc_min', 0.2, ... 'soc_init', 0.5, 'eta_charge', 0.95, 'eta_discharge', 0.95); data.price.grid_buy = repmat([0.35 0.7 1.1], 8, 1); % 峰谷平电价示例这里有个小技巧:把SOC相关的初始化直接放进data结构体,避免后面在多个函数里反复传参。储能初始SOC和期望末端SOC通常设为相等值,比如0.5,这是为了保证一天调度周期的循环可持续性。如果你在做的是多日连续调度,那这个初始值的设定逻辑要另行处理。
负荷和新能源预测数据我用矩阵存,行是时间断面(24行),列是不同母线或电源类型:
% 第1列电负荷,第2列热负荷,第3列氢负荷,第4列风电预测,第5列光伏预测 load_db = csvread('profile_24h.csv');3.2 决策变量编码方式的确定
怎么把调度方案变成一维向量,很大程度上决定了邻域搜索的效率和约束检查的复杂度。我采用的是“按设备、按时间”堆叠,而不是“按时间、按设备”。
% 决策变量 x 的结构: % x(1:24) :燃气轮机输出功率 P_CHP % x(25:48) :燃气轮机热功率 H_CHP % x(49:72) :燃煤电厂购电 P_grid % x(73:96) :P2H2 电功率 % x(97:120) :电锅炉电功率 P_EB % x(121:144) :蓄电池充放电功率(正放负充) % x(145:168) :储氢罐充放氢功率(正放负充) % x(169:192) :蓄电池SOC % x(193:216) :储氢罐SOC按设备堆叠的好处是,在邻域扰动时可以直接对一个设备的整条24小时曲线做局部修改。比如只随机挑某台机组的连续3小时出力做扰动,向量索引范围很清晰。但代价是跨设备做能量平衡检查时,需要先reshape成[时段 × 设备]的矩阵再计算,这样会稍麻烦一些。我不太建议初学者把所有变量混在一个一维向量里不做区分,那样后期debug会非常痛苦。
编码里还有一点:CHP机组的电功率和热功率并不是独立变量。对背压式机组,电热比接近固定;对抽凝式机组,热电可调,但必须在可行域范围内。我建议在邻域生成阶段就把这两者关联起来,先扰动电功率,再根据可行域计算热功率范围,在范围内随机取热功率,这样生成的新解自然满足CHP耦合约束。
3.3 目标函数与约束罚函数的Matlab写法
目标函数我用一个大函数evaluate_solution.m统一处理。输入是决策向量x和data结构体,输出是两个目标值和约束违限量。所有中间变量都放在函数内部,不要把计算结果散落在全局工作区里。
function [f_cost, f_co2, viol] = evaluate_solution(x, data) % 1. 展开决策变量 P_CHP = x(1:24); H_CHP = x(25:48); P_grid = x(49:72); P_P2H2 = x(73:96); P_EB = x(97:120); ... % 2. 计算目标1:成本 cost_grid = sum(P_grid .* data.price.grid_buy) * dt; cost_gas = sum((P_CHP ./ data.dev.CHP.eta_elec + H_GB ./ data.dev.GasBoiler.eta) .* data.price.gas) * dt; cost_om = sum(P_P2H2 * data.dev.P2H2.cost) * dt; % 仅示例 f_cost = cost_grid + cost_gas + cost_om; % 3. 计算目标2:碳排放 co2_grid = sum(P_grid .* data.co2.grid_factor) * dt; co2_gas = sum((P_CHP ./ data.dev.CHP.eta_elec + H_GB ./ data.dev.GasBoiler.eta) .* data.co2.gas_factor) * dt; f_co2 = co2_grid + co2_gas; % 4. 约束检查 viol = check_constraints_all(x, data); end这里要注意一个数值细节:如果用的是逐时步长,乘dt = 1表示时间单位是小时,功率单位是kW,那么能量单位就是kWh。如果时间步长是15分钟,dt = 0.25,所有乘了dt的项才能正确换算能量。很多算法跑出来结果离谱,不一定是算法问题,而是量纲换算错了。
约束罚函数我不会把所有约束的违限量直接相加了事,而是给不同优先级设置不同权重:
viol = viol_balance * 1000 + viol_ramp * 100 + viol_soc * 10 + viol_limit * 10000;能量平衡和设备上下限是最硬性的约束,权重拉满;爬坡约束和SOC约束次之。这样做的目的是让退火算法的搜索过程优先避开那些惩罚量大的区域,而不是在几个约束之间进退两难。
3.4 算法主循环与Pareto存档更新实现
主体算法的循环结构和参数我之前已经介绍过原理,这里直接看关键的代码骨架:
% 参数初始化 T0 = 100; alpha = 0.95; maxIter = 3000; T = T0; x = initSolution(data); archive = []; % 主循环 for iter = 1:maxIter % 随机生成一组权重,用于多目标合并搜索 w1 = rand(); w2 = 1 - w1; x_new = generateNeighbor(x, T / T0, data); [f1_old, f2_old, viol_old] = evaluate_solution(x, data); [f1_new, f2_new, viol_new] = evaluate_solution(x_new, data); % 搜索目标:加权法和罚函数 search_old = w1 * f1_old + w2 * f2_old + lambda * viol_old; search_new = w1 * f1_new + w2 * f2_new + lambda * viol_new; delta = search_new - search_old; if delta < 0 || rand() < exp(-delta / T) x = x_new; end % 更新Pareto存档:只保留可行非支配解 if viol_new <= 1e-6 archive = updateArchive(archive, [f1_new f2_new], x_new); end T = alpha * T; end这里的updateArchive是核心函数。我的实现是:把每个解的二维目标值存成N×2矩阵,同时用一个很小的自定义结构体把对应解向量也存下来。这样后续画Pareto前沿图时,坐标轴非常好办;如果需要进一步分析某个特定折中方案的调度策略,也可以直接从存档里取回完整决策变量。
画Pareto前沿是检查多目标算法最直观的方式:
figure; plot(archive_obj(:,1), archive_obj(:,2), 'o'); xlabel('运行成本(元)'); ylabel('碳排放量(kg)'); title('MOSA求解得到的Pareto前沿'); grid on;4. 算例结果与算法性能分析
代码能跑通只是第一步,结果合不合理、算法性能怎么样,都需要在标准算例上做测试和分析。我会用一个小规模但完整的算例做演示。
4.1 测试场景与参数设置
我用一个典型的小型区域综合能源系统,包含以下组成部分:上级电网(可购电)、CHP机组、燃气锅炉、电锅炉、P2H2(电转氢)、蓄电池、储氢罐,以及24小时电热氢负荷曲线。风电光伏没有直接加进主算例,避免变量太多干扰问题分析。
具体参数如下:
| 设备 | 额定功率 | 效率 | 爬坡速率 | 备注 |
|---|---|---|---|---|
| CHP机组 | 800 kW | 电0.42/热0.45 | 200 kW/h | 热电耦合 |
| 燃气锅炉 | 400 kW | 0.9 | 120 kW/h | — |
| 电锅炉 | 300 kW | 0.95 | 100 kW/h | P2H |
| 电转氢(P2H2) | 300 kW | 0.65 | 80 kW/h | 产氢供氢负荷 |
| 蓄电池 | 200 kWh | 0.95/0.95 | — | SOC 0.2~0.9 |
| 储氢罐 | 50 kg | — | — | SOC 0.1~0.9 |
分时电价按“峰1.1元/kWh、平0.7元/kWh、谷0.35元/kWh”,天然气价设为2.5元/m³,天然气热值取9.7 kWh/m³。碳排放系数:电网购电0.82 kg/kWh,天然气燃烧0.2 kg/kWh。
仿真时段一共24小时,时间步长为1小时,决策变量包括各设备逐时出力、储能SOC、储氢罐SOC,总维度约200维。这个规模在退火算法里属于小到中型,3000次迭代大概一两分钟能跑完。
4.2 收敛性能与Pareto前沿
跑完3000次迭代后,我画了目标函数的收敛曲线和最终的Pareto前沿。退火算法本身的“温度下降”决定了它不会像遗传算法那样有明显的“代际降幅”,而是呈锯齿状波动。因为退火期间会接受部分差解,所以成本曲线不是单调下降的。这一点新手要习惯,不要看到曲线反弹就觉得算法崩了。
判断收敛要看存档里面的非支配解数量是否稳定,以及平均目标值是否主要朝前沿方向推进。我一般会在每100次迭代记录一次“当前解的目标值”,然后在高温度段和低温度段分别统计一步搜索的平均改善幅度。如果低温段的改善幅度已经远小于初始阶段,就说明可以终止迭代。
最终获得的Pareto前沿大致呈一条向左下方凸出的曲线。最左边的点对应最低成本方案,代价是碳排放较高;最右边的点对应最低排放方案,成本明显上升。为什么会有这种趋势?因为在低排放目标的压力下,系统会减少从高碳排放系数的电网购电,转用碳排放相对低的天然气机组或者多利用P2H2与储氢环节进行能源时移,而后者的运行成本往往更高。
4.3 与加权和法与NSGA-II的对比结论
为了验证多目标退火算法的实际效果,我在同一套算例上跑了两个对标算法:一个是传统的加权和方法(单目标退火,固定权重),另一个是经典的NSGA-II。对比结果可以归纳为一张表:
| 算法 | 耗时(秒) | Pareto解个数 | 前沿分布均匀度 | 收敛性 |
|---|---|---|---|---|
| MOSA(本文) | 68 | 45 | 较好 | 优秀 |
| 加权和SA(固定权重) | 42 | 1 | 无多样性 | 优良,但单次只能得到一个点 |
| NSGA-II | 120 | 52 | 中等 | 良好,但种群分布受交叉变异参数影响大 |
加权和法的最大缺陷已经暴露得很明显:固定一套权重只能求出一个解。如果你想生成完整前沿,就得跑很多组权重,耗时换算下来并不比MOSA少,而且权重取值的均匀分布并不能保证前沿上的点均匀分布。
NSGA-II在解集多样性上并不差,但它的参数敏感性更高。我试了不同的交叉概率和变异概率,结果波动比较大,而且耗时比MOSA长。对P2X这种带强耦合约束的问题,NSGA-II的交叉操作经常把CHP的电热耦合关系打散,导致生成大量不可行解,需要额外做修复,进一步拖慢速度。
所以我的结论是:在不追求数学规划全局最优解、更看重快速出可行方案和灵活改模型的场景下,多目标退火算法是一个性价比很高的选择。它不完美,但它稳。
5. 调试与实战中踩过的坑
最后这部分,我把自己在Matlab代码实现过程中踩过的坑梳理出来,希望你能少走一些弯路。
5.1 邻域扰动幅度引起的“假收敛”
我最开始写MOSA的时候,给所有连续型变量固定用sigma = 0.05倍额定功率做高斯扰动。结果算法跑得飞快,前端1000次迭代目标值快速下降,之后完全卡住不动。最后发现问题是sigma太小,后期搜索根本跳出不了局部区域。
但把sigma调大之后又遇到另一个极端:目标值长期不下降,看上去就像随机游走。原因是扰动太剧烈,绝大部分新解都严重违反约束,罚函数把搜索目标带上天,Metropolis准则又不断接受这些差解,最终变成了随机爬坡。
解决方法是“温度自适应步长”。我把sigma与当前温度绑定:高温阶段sigma设为额定功率的0.3倍,随着温度降低线性收缩到0.03倍。这样算法前期大步探索、后期小步精细调整,收敛效果立刻改善。这个调整对调度问题的改善非常明显,我强烈建议你在自己的代码上试一试。
另外还可以加一个限制:记录连续100次迭代中最优目标值的改善幅度,如果低于某个阈值,就对sigma做一个小幅放大,模拟“回火”操作,帮助算法跳出停滞状态。
5.2 罚函数系数和量纲归一化问题
罚函数系数的设置我单独踩过一个大坑。成本数量级是10^4,碳排放数量级是10^3,能量平衡违限量的数量级可能是10^2(kW)。如果直接把viol乘以同一个系数加到搜索目标里,你会发现最后的搜索结果完全被能量平衡约束支配,成本和碳排放根本没有被有效优化。
正确的思路是先归一化。我在代码里把成本和碳排放除以其各自的基准值(比如最大允许成本、历史平均排放量),然后和归一化后的罚函数相加,保证这几个项都在同一个量纲量级下。之后在调试时,我会在目标函数里临时打印出搜索目标的分项构成,观察到底是目标项主导还是罚项主导。
如果你发现罚项一直是0,可行域覆盖率过高,那算法其实没有充分搜索边界区域;如果罚项始终非常大,说明初始解生成得不好,建议先修复初始解生成逻辑。
5.3 Pareto存档膨胀与前沿分布拥挤
存档集合如果不加控制,会随着迭代次数增加不断膨胀。我第一次跑3000次迭代,Pareto存档点达到上千个。画出来看似壮观,但里面大量点实际上非常拥挤,很多解之间的差异可以忽略不计,这对最终决策没有帮助,也拖慢了存档更新时两两支配判断的时间。
我的处理方式是用“网格拥挤度筛选”来修剪存档。把目标空间划分成网格,每个网格内只保留一个代表性解,或者每次修剪时优先移除与相邻解距离最近的解。这样存档数量能控制在30到60个左右,前沿分布也更均匀。
注意:存档修剪不要只在算法结束后做一次,要在迭代过程中定期做,否则迭代后期存档被填满,新解的核心信息容易被原有解淹没。
5.4 Matlab代码性能优化的小技巧
Matlab本身解释执行,循环效率不像C或Python那样有优势,退火算法又是串行迭代,很难直接并行。我把代码性能优化的经验总结成几条:
首先,目标函数和约束检查尽量向量化。我最初用for循环逐时段做电量平衡检查,整个3000次迭代跑下来要接近十分钟。改成矩阵运算后,相同算例只需要一分钟出头。对于24时段的小模型,向量化的收益远大于复杂化带来的调试成本。
其次,避免在迭代循环内反复计算与当前解无关的常量。比如分时电价、排放系数、设备额定值这些,都是只在加载案例时计算一次的,不要放到evaluate_solution里每次重新计算。
再次,如果你在跑更大规模系统(比如多区域IES),Matlab里有几个可用的加速方向:用parfor并行计算多个工况的初始调度方案,只共享存档列表;把目标函数函数体改成编译后的MEX文件,或者用matlabFunction、coder生成C代码;在内存允许的情况下用preallocate预定义所有中间矩阵,不要使用动态增长的数组。
最后,如果只需要结果展示和画图,可以不输出每一代目标值,只记录每次存档更新时的关键指标,能显著降低I/O开销。
我个人在实际项目里最深刻的体会是:多目标退火算法不是一个“开箱即用”的黑盒工具,它出结果容易,出好结果难。真正决定最终效果的,往往不是算法框架本身,而是你对问题的理解程度、约束建模的精细程度,以及那两三个看似不起眼的参数设定。这个项目我最初用两个星期搭出能跑的版本,又花了一个月反复调参、修bug,才让结果稳定到可以写论文和报告的水平。如果你也是第一次接触这块,建议先跑通一个简化案例,确认每一段逻辑都符合物理直觉,再逐步增加P2X设备的数量和复杂度,这样排查问题会轻松很多。