简介:面向微电网多目标优化调度问题,这份资源整合了基于非支配排序蜣螂优化算法(NSDBO)的完整Matlab代码与配套说明,适用于电力系统、自动化、电子信息工程及数学等专业的学生开展课程设计、期末大作业或毕业设计,也适合科研入门者快速复现多目标调度实验。压缩包共32个文件,容量仅5.14MB,其中11个M脚本构成调度主程序、目标函数及结果绘图模块,17张JPG图片展示Pareto前沿与多目标优化结果,另有模型说明PDF、标准数据表格XLSX和说明文档TXT,便于对照学习数据格式与算法结构。已有184人学习。代码采用参数化编程,关键参数易于修改,注释详细清晰;附带可直接运行的测试数据,一键运行即可复现NSDBO求解微电网调度模型的过程。该模型以经济成本、环境排放和系统稳定性为目标,通过非支配排序筛选Pareto前沿解集,为决策者提供多种可行调度方案。读者可在此基础上替换数据、修改约束或扩展目标函数,迁移至不同微电网算例。整体来看,资源兼具教学与科研参考价值。
1. NSDBO不是玄学:给蜣螂算法补上多目标那一课
微电网做多目标优化调度,最头疼的不是算法跑不起来,而是运行成本和碳排放这两个目标天生对着干:多烧柴油机确实便宜,排放大;多用电储能和可再生能源,绿色了但成本高。单目标优化只能给你一个“最省”或“最绿”的调度点,而基于非支配排序的蜣螂优化算法NSDBO要输出的是整条Pareto前沿,让运行人员按自己的偏好去挑方案。它的做法很直接:把蜣螂优化算法DBO的滚球、跳舞、繁殖、觅食位置更新机制保留下来,再嵌进NSGA-II式的快速非支配排序和拥挤度距离,让原本只会朝一个最优解收缩的算法学会在多个目标之间铺开解集。适合正在做微电网日前调度、储能充放电策略或机组组合的研究生和工程师,Matlab里按“建模—算法改造—主循环—调参—验证”的顺序就能落地。
2. 先把调度问题写成优化问题:目标函数、约束与决策变量编码
2.1 两个目标天生互斥:运行成本与碳排放
微电网的典型结构是光伏、风电、柴油机、微型燃气轮机和蓄电池储能,负荷不足或过剩时再跟大电网交互。光伏和风电基本零边际成本、零直接排放,但它们不可控,所以调度真正要决策的是柴油机、燃气轮机、储能充放电以及从大电网买电的功率曲线。运行成本包括柴油和天然气的燃料费、各设备单位出力的运维费、从电网购电的电费;碳排放包含柴油机和燃气轮机的直接排放,以及外购电对应到火电厂的间接排放。成本低的调度往往把柴油机开到高位,把便宜但高碳的网电买进来,排放自然压不住;要减排就得少烧油、多让储能削峰填谷,可储能循环损耗和购电电价又会把成本抬起来。这就是两个目标互斥的根源,也是为什么要用多目标方法而不是简单加权。
目标函数在Matlab里我习惯写成两个独立输出,不合并成一个适应度值。多目标算法要的就是目标向量本身,合并成加权和等于提前把决策者的偏好写死,Pareto前沿就丢了。下面是评估函数的骨架:
function [cost, emission] = evaluate_objectives(x, data) % x是按24时段摊平的决策变量向量,维度为5*24 % 先reshape成5行24列:柴油机、燃气轮机、储能放电、储能充电、电网交互 X = reshape(x, 5, 24); P_de = X(1,:); % 柴油机出力,kW P_mt = X(2,:); % 微型燃气轮机出力,kW P_bd = X(3,:); % 储能放电功率,kW P_bc = X(4,:); % 储能充电功率,kW P_grid = X(5,:); % 与电网交互功率,正买电负卖电,kW % 运行成本:燃料 + 运维 + 储能折算 + 购电 fuel_cost = sum(data.c_de * P_de + data.c_mt * P_mt); om_cost = sum(data.k_de*P_de + data.k_mt*P_mt + ... data.k_bess*(P_bd + P_bc)); grid_cost = sum(data.price_buy .* max(P_grid,0) + ... data.price_sell .* min(P_grid,0)); cost = fuel_cost + om_cost + grid_cost; % 碳排放:本地机组直接排放 + 购电间接排放 emission = sum(data.e_de*P_de + data.e_mt*P_mt + ... data.e_grid * max(P_grid,0)); end要注意的点:储能充放电功率写成两个独立变量而不是一个带符号的变量,是为了避开“同一时段既充又放”的掩盖问题,后面约束里直接用互斥惩罚处理。电网交互功率保留正负号,买电算成本支出和碳排放,卖电回收一点成本但不抵扣排放。price_buy和price_sell可以做成峰谷分时电价,数据里准备好就行。这个函数是所有后续计算的地基,种群里每个个体都要过一遍,所以尽量别在循环里做重复的矩阵运算。
2.2 决策变量编码:24时段出力矩阵与SOC递推
决策变量怎么编码,直接影响算法收敛速度和约束处理的复杂度。最常见的做法是只把可控设备作为决策变量,光伏和风电按预测曲线直接参与功率平衡,不作为变量。这样决策变量就是柴油机、燃气轮机、储能放电、储能充电、电网交互这5个对象在24个时段的出力,合计5×24=120维。120维对种群规模200左右的NSDBO来说不算大,但已经足够暴露“维度灾难”和边界处理的问题了。
每种设备的上下界不同,柴油机有最小技术出力不能低于30%额定功率,储能充放电有额定功率限制,电网交互受联络线容量约束。所以上下界要用矩阵写清楚,初始种群也要在这个范围内生成:
% 决策变量边界:5种设备 x 24时段 lb = repmat([30; 20; 0; 0; -150], 1, 24); % 下界:柴油机最低30kW等 ub = repmat([100; 80; 60; 60; 150], 1, 24); % 上界:各类设备上限 lb = lb(:); % 摊平成120维列向量 ub = ub(:); % 初始化种群:N个个体,每个120维 N = 200; pop = lb' + (ub - lb)' .* rand(N, 120);储能SOC的递推不放在目标函数里,而是放在约束检查之前单独算。因为SOC是跨时段耦合的,前一时刻的充放电决策会改变后一时段的可用容量。SOC递推公式为SOC(t+1) = SOC(t) + η_c*P_bc(t)*Δt - (1/η_d)*P_bd(t)*Δt,初始SOC一般取0.5,调度结束不强制归位但最好落在合理区间,否则第二天没法连续运行。η_c取0.95,η_d取0.9,这是锂离子电池储能常见的效率范围。SOC计算和约束检查放在一起循环24次即可。
2.3 约束检查:功率平衡与机组出力上下限的Matlab实现
微电网调度这步最容易被罚函数带偏。功率不平衡量、SOC越限、机组出力越界这堆约束,如果用硬约束直接淘汰个体,NSDBO的种群会急剧收缩,早熟;全用大罚函数,算出来的“最优解”可能其实还在违反功率平衡,只是罚得不够狠。我的做法是把所有约束的违例量返回成一个标量,然后在目标函数后面按数量级加权。违例量设计成“归一化的越限距离”,这样比直接加一个固定大数更稳。
function vio = get_constraints(x, data) X = reshape(x, 5, 24); P_de = X(1,:); P_mt = X(2,:); P_bd = X(3,:); P_bc = X(4,:); P_grid = X(5,:); % 1. 功率平衡违例:负荷 + 充电 = 柴油机 + 燃气轮机 + 放电 + 光伏 + 风电 + 电网 balance = P_de + P_mt + P_bd + data.P_pv + data.P_wt + P_grid ... - data.P_load - P_bc; vio_balance = sum(abs(balance)); % 2. 机组出力上下界违例,按越界距离累加 vio_de = sum(max(P_de - data.de_max, 0)) + sum(max(data.de_min - P_de, 0)); vio_mt = sum(max(P_mt - data.mt_max, 0)) + sum(max(data.mt_min - P_mt, 0)); % 3. 储能SOC边界与充放电互斥 soc = 0.5; % 初始SOC vio_soc = 0; for t = 1:24 soc = soc + data.eta_c*P_bc(t) - (1/data.eta_d)*P_bd(t); % 假设Δt=1h,容量以kWh计 vio_soc = vio_soc + max(soc - 0.9, 0) + max(0.1 - soc, 0); vio_soc = vio_soc + 1000 * (P_bc(t) * P_bd(t)); % 充放互斥强惩罚 end % 4. 爬坡约束:柴油机和燃气轮机相邻时段变化量限制 vio_ramp = sum(max(abs(diff(P_de)) - 30, 0)) + ... sum(max(abs(diff(P_mt)) - 25, 0)); vio = vio_balance + vio_de + vio_mt + vio_soc + vio_ramp; end爬坡约束用diff一次算完,注意diff结果长度是23,别对齐错位。SOC的初始值、容量、Δt=1h这几个参数在数据里写清楚,不同微电网的量纲不同,抄代码时最容易在这里翻车。罚函数系数我一般放在评估函数外统一给,比如fit = cost + 1e4 * vio,再用这个fit做非支配排序前的单目标排序参考。但NSDBO主循环里非支配排序其实只看[cost, emission]两个原始目标,罚函数只用来淘汰明显违例的个体和辅助选择参考解,这个概念要分清楚。
3. 把非支配排序嵌进DBO:NSDBO的核心改造步骤
3.1 为什么单目标DBO不能直接跑多目标:Pareto支配的定义
蜣螂优化算法DBO原本是单目标全局优化算法,它的核心机制是让所有个体向当前最优位置收缩,配合繁殖和觅食策略做局部精细搜索。单目标问题只有一个最优解,算法往那个方向压没问题。但多目标优化里存在一组相互冲突的Pareto最优解,没有一个点能同时让成本和排放都最小,直接套DBO的话,整个种群会被拽向某个单一目标占优的区域,Pareto前沿根本铺不开。
所以在NSDBO里,先把“谁好谁坏”的比较规则换掉。两个解A和B,只有当A在所有目标上都不差于B,且至少有一个目标严格优于B时,才叫A支配B。互相都不支配的一组解就是同一层的非支配解。这个判定是后面所有排序的基础,代码里最容易犯的错是支配方向写反。
function d = dominates(a, b) % a、b是目标向量,例如[cost, emission] % 返回true表示a支配b d = all(a <= b) && any(a < b); end注意这里的目标都是最小化问题,所以用<=和<。如果是最大化问题,要全部取反。NSDBO每次比较个体优劣,都用这个函数而不是比较加权和,这一步改变了整个算法的选择压力方向。
3.2 快速非支配排序与拥挤度距离:可直接用的两段代码
快速非支配排序的核心思想是给每个解算两个量:被谁支配(n)和支配谁(S)。第一轮先把n=0的解全部放进第1层;然后遍历这些解的S集合,把其中每个解的n减1,减到0就归入第2层,如此递推。这个逻辑比“两两比较后反复删层”快很多,尤其在种群规模200以上时差距明显。
function [rank, front] = nondominated_sort(obj) % obj: N x M 的目标矩阵,N个解,M个目标(最小化) N = size(obj, 1); n = zeros(N, 1); % 被支配次数 S = cell(N, 1); % 被该解支配的解集合 front = cell(N, 1); % 每一层的解索引 % 第一遍:计算n和S for p = 1:N S{p} = []; for q = 1:N if p == q, continue; end if dominates(obj(p,:), obj(q,:)) S{p} = [S{p}, q]; elseif dominates(obj(q,:), obj(p,:)) n(p) = n(p) + 1; end end end % 第一层:所有n=0的解 F1 = find(n == 0); front{1} = F1; k = 1; while ~isempty(front{k}) Q = []; for i = 1:length(front{k}) p = front{k}(i); for j = 1:length(S{p}) q = S{p}(j); n(q) = n(q) - 1; if n(q) == 0 Q = [Q, q]; end end end k = k + 1; front{k} = Q; end front = front(1:k-1); rank = zeros(N, 1); for i = 1:length(front) rank(front{i}) = i; end end这个函数返回两个东西:每个解的层级rank,以及按层分组的front。层级为1的就是当前Pareto前沿。拥挤度距离是同一层内对每个目标排序后,计算相邻解之间的归一化距离之和,首尾解的距离设为无穷大,因为端点的多样性价值最高、必须保留:
function cd = crowding_distance(obj, F) % F是某一层的解索引集合 Nf = length(F); cd = zeros(Nf, 1); M = size(obj, 2); if Nf <= 2 cd(:) = inf; return; end for m = 1:M [~, idx] = sort(obj(F, m)); % 按第m个目标升序 F_sorted = F(idx); cd(idx(1)) = inf; cd(idx(end)) = inf; fmax = obj(F_sorted(end), m); fmin = obj(F_sorted(1), m); if fmax - fmin > 1e-12 for i = 2:Nf-1 cd(idx(i)) = cd(idx(i)) + ... (obj(F_sorted(i+1), m) - obj(F_sorted(i-1), m)) / (fmax - fmin); end end end end3.3 精英保留与排序层级:NSDBO的选择压力来源
排序函数准备好后,关键问题变成“选择压力怎么传给下一代”。DBO原本的策略是每一代用新位置直接替换旧位置,这样最优解容易丢失。NSDBO的做法借鉴了NSGA-II的精英保留:把父代和子代合并成2N个解,先做非支配排序,再按层级从第1层开始逐个收入下一代,直到第L层装不下时,用拥挤度距离排序只保留该层里最分散的那几个。这样Pareto前沿上的优秀解永远不会被新的差解挤掉,这也是NSDBO比“多目标加权+DBO”解法更值得投入的核心原因——它有严格的理论选择机制,而不是靠权重碰运气。
选择参考解时,我通常取当前第1层里拥挤度距离最大的那个作为x_b,而不是随便取一个非支配解。因为拥挤度最大的点往往位于前沿两端或稀疏区域,把它作为DBO偷窃策略的吸引子,可以有效防止种群朝前沿中间一点塌缩。这个细节很微妙,但直接决定了前沿铺得匀不匀。
4. 求解微电网调度的NSDBO主循环:从蜣螂行为到Pareto前沿
4.1 DBO的四类行为策略:滚球、跳舞、繁殖、偷窃的更新公式
DBO的更新策略不是一整块,而是把种群按比例分成四类角色,每类用不同的位置更新公式。理解这四类策略不用把它当黑匣子,只需知道它们是不同粒度的搜索算子:滚球是全局探索、跳舞是遇到障碍时的随机转向、繁殖是局部收缩精搜、偷窃是围绕全局最优的跳跃式开采。NSDBO里这些公式原样保留,只是“全局最优”不再是一个点,而是每次从当前非支配前沿里挑出来的参考解。
滚球策略分两种情形:无障碍物时按x_new = x + α*k*(x - x_prev) + b*(x - x_worst)更新,其中α在-1和1间随机,k=0.1是偏转系数,b=0.3是风速系数,x_worst是当前代最差解。这个公式里“朝最差解拉开距离”的项其实是DBO里重要的探索来源。有障碍物时改为跳舞行为:x_new = x + tan(θ)*|x - x_prev|,θ是0到π之间的随机角,实现角度越大跳得越远。
繁殖策略是在当前最优位置附近收缩出一个产卵区:Lb* = max(lb, x_b - R),Ub* = min(ub, x_b + R),R是变量区间长度的一半,然后在这个区域内随机生成:x_new = x_b + e1*(x - Lb*) + e2*(x - Ub*),e1、e2是随机向量。觅食策略类似,但收缩范围更小,只在最优个体邻域内做精细搜索。偷窃策略的更新公式是x_new = x_b + S*g*(|x - x_b| + |x - x_opt|),S是常数取4,g是正态分布随机数,这个算子的特点是步长较大,用来跳出局部区域。
4.2 NSDBO主循环框架:非支配排序在哪一步介入
NSDBO的主循环不需要把排序函数嵌进每种策略内部,而是把整个“评估—排序—选参考—更新—再排序”作为一个整体循环。我给出一个可以直接跑通框架的版本,种群规模200、迭代800代,决策变量120维。这个框架去掉具体数据文件也能看懂逻辑,替换数据即可用于自己的微电网案例。
% 主循环参数 N = 200; Tmax = 800; num_var = 120; % 初始化种群 pop = lb' + (ub - lb)' .* rand(N, num_var); obj = zeros(N, 2); for t = 1:Tmax % 1. 评估每个个体的原始目标 for i = 1:N [cost_i, emit_i] = evaluate_objectives(pop(i,:), data); vio_i = get_constraints(pop(i,:), data); if vio_i > 0 cost_i = cost_i + 1e4 * vio_i; % 罚函数只影响参考解排序 emit_i = emit_i + 1e4 * vio_i; end obj(i,:) = [cost_i, emit_i]; end % 2. 非支配排序 + 拥挤度距离 [rank, front] = nondominated_sort(obj); cd = zeros(N,1); for fi = 1:length(front) cd(front{fi}) = crowding_distance(obj, front{fi}); end % 3. 选参考解:第1层拥挤度最大的个体 F1 = front{1}; [~, idx_b] = max(cd(F1)); x_b = pop(F1(idx_b), :); % 4. 按比例分四类策略更新 pop_new = zeros(N, num_var); nb_roll = round(0.45*N); % 滚球+跳舞 nb_breed = round(0.25*N); % 繁殖 nb_forage = round(0.15*N); % 觅食 % 剩余为偷窃策略 for i = 1:N if i <= nb_roll % 无障碍滚球与跳舞混合 pop_new(i,:) = roll_update(pop(i,:), pop(max(i-1,1),:), ... x_b, lb, ub, rand < 0.5); elseif i <= nb_roll + nb_breed pop_new(i,:) = breed_update(pop(i,:), x_b, lb, ub, t, Tmax); elseif i <= nb_roll + nb_breed + nb_forage pop_new(i,:) = forage_update(pop(i,:), x_b, lb, ub); else pop_new(i,:) = steal_update(pop(i,:), x_b, lb, ub); end end % 5. 边界修复:越界变量拉回边界内再随机抖动 pop_new = min(max(pop_new, lb'), ub'); pop_new = pop_new + 0.01 * (ub - lb)' .* randn(N, num_var); pop_new = min(max(pop_new, lb'), ub'); % 6. 合并父代和子代,精英选择下一代 pop_all = [pop; pop_new]; obj_all = zeros(2*N, 2); for i = 1:2*N [obj_all(i,1), obj_all(i,2)] = evaluate_objectives(pop_all(i,:), data); end [rank_all, front_all] = nondominated_sort(obj_all); cd_all = zeros(2*N, 1); for fi = 1:length(front_all) cd_all(front_all{fi}) = crowding_distance(obj_all, front_all{fi}); end pop_next = []; fi = 1; while length(pop_next) + length(front_all{fi}) <= N pop_next = [pop_next, front_all{fi}]; fi = fi + 1; end % 最后一层按拥挤度取前几个 rem = N - length(pop_next); if rem > 0 last = front_all{fi}; [~, idx_cd] = sort(cd_all(last), 'descend'); pop_next = [pop_next, last(idx_cd(1:rem))]; end pop = pop_all(pop_next, :); end主循环里最关键的是第6步。合并父代子代后再排序,保证了精英解一定能进入下一代,这是NSDBO和原始DBO最大的区别,原始DBO是就地更新、没有这个保优环节。罚函数我放在评估目标时直接加到两个目标值上,这样如果某个个体严重违反功率平衡,它在非支配排序中的支配关系会被扭曲,自然就会被淘汰。罚函数系数1e4这个量级需要按目标量纲调整,成本和排放通常都在千这个量级,1e4够用。如果数据量纲大,这个系数要等比放大。
4.3 从Pareto解集选折中解:模糊隶属度函数怎么用
NSDBO跑完得到的是一组非支配解,不是单个调度方案。实际运行时要给调度员一个“折中方案”,常见做法是用模糊隶属度函数在每个目标上打分,再取综合隶属度最大的解作为推荐方案。这样选出来的折中解兼顾两个目标,而不是偏向某一端。对不同偏好场景,调整各目标的权重即可。
function [best_idx, best_mu] = select_compromise(obj_front) % obj_front: Pareto前沿解的目标矩阵,每行[cost, emission] M = size(obj_front, 2); nf = size(obj_front, 1); mu = zeros(nf, M); for m = 1:M fmax = max(obj_front(:,m)); fmin = min(obj_front(:,m)); if fmax - fmin > 1e-12 mu(:,m) = (fmax - obj_front(:,m)) / (fmax - fmin); else mu(:,m) = 1; end end total_mu = sum(mu, 2); % 也可以加权重 [best_mu, best_idx] = max(total_mu); end模糊隶属度选择出的解,对应的决策变量向量就是推荐的24时段调度单。把它和所有Pareto解一起保存,后续做灵敏度分析和方案对比都用得上。需要注意的是,如果要硬约束SOC终值归位,折中解未必满足,需要在所选方案上轻微调整储能最后两时段的充放电曲线,这个在下一章避坑部分细说。
5. NSDBO参数设置与实践避坑:常见问题排查清单
5.1 参数表:种群规模、迭代次数、策略比例与边界收缩建议
参数设置没有绝对标准,但按120维决策变量这个规模,我给一组经过大量测试可复用的参考值。种群规模低于100会出现前沿断档,高于300计算量陡增但前沿质量提升有限。迭代次数500代能出一个可用的轮廓,800代基本稳定,1500代对复杂约束微电网能明显改善边界解的分布。策略比例里滚球占比最大是DBO的特点,它承担全局探索;繁殖和觅食负责局部精搜;偷窃占比最小但步长最大,负责跳出局部。这些比例对多目标版本同样适用,不需要因为加排序就大改。
| 参数 | 建议值 | 说明 |
|---|---|---|
| 种群规模 N | 200 | 120维问题建议不低于150 |
| 最大迭代 Tmax | 800 | 看前沿变化率决定是否继续 |
| 滚球+跳舞比例 | 0.45 | 全局探索主力,不宜低于0.3 |
| 繁殖比例 | 0.25 | 收缩搜索,边界收缩速度影响早熟 |
| 觅食比例 | 0.15 | 局部精细搜索 |
| 偷窃比例 | 0.15 | 大跳跃,比例过高会破坏前沿 |
| 偏转系数 k | 0.1 | DBO原始参数 |
| 风速系数 b | 0.3 | 控制滚球步长 |
| 偷窃常数 S | 4 | 步长缩放,可按边界范围放大 |
| 罚函数系数 | 1e4 | 按目标量纲调整 |
环境方面,这段代码用基础Matlab语法写的,不依赖优化工具箱,R2024a跑通后换2026b也兼容。之前有人问我是不是必须装特定工具箱,其实不用,自己实现非支配排序函数反而少了版本差异的麻烦。
5.2 现象一:非支配排序结果永远只有一层
- 现象:
rank数组全是1,Pareto前沿上挤了200个点,后续层级为空。 - 原因:目标函数里不小心把两个目标合并了,或者罚函数加得过大,导致所有个体的目标向量退化成近似成比例的两个数。比如成本和排放都用同一个燃油量推导出来,它们线性相关,非支配排序自然失效。另外也可能是
dominates函数里用了all(a < b)而不是all(a <= b) && any(a < b),严格小于把所有相等目标值的解都判成了互不支配。 - 解决:先检查目标函数的独立性,把成本和排放的散点图画出来,如果点全在一根直线上说明目标耦合了,调整模型;再检查
dominates的判定条件,确保相等目标值时不会误判。这个坑很隐蔽,因为Matlab不会报错,只会在画前沿时发现所有点都在一条曲线上塌成一团。
5.3 现象二:功率平衡被罚函数带偏,SOC出现尖角
- 现象:最优调度单里功率平衡误差不小,SOC曲线出现锯齿状尖角,甚至一天内从0.1跳到0.9。
- 原因:罚函数系数不够或者约束函数里充放互斥的惩罚量级和功率平衡惩罚不在一个数量级。算法会优先满足“罚得狠”的约束,忽视罚得轻的。SOC尖角则是因为SOC约束只在整数边界检查,没有考虑充放电功率连续变化时可能中间越限。
- 解决:把约束函数据所有项归一化到同一量级。我一般把功率平衡误差用负荷峰值归一化,SOC用容量归一化,充放互斥用额定功率平方归一化。这样罚函数系数可以统一调,不会出现某个约束被忽略的情况。调试时直接打印
get_constraints返回的每个分量,看哪个量级明显偏大,逐个修正比瞎调系数快得多。
5.4 现象三:迭代中后期种群塌缩,Pareto前沿只剩一个解
- 现象:前100代前沿铺得挺好,300代后前沿点越来越少,最后只剩三四个解抱团。
- 原因:偷窃策略的步长太大,繁殖策略的收缩区间R随迭代指数级缩小,两者叠加导致整个种群被快速拽向参考解附近。单目标下这是快速收敛的好事,多目标下却是灾难——它会牺牲前沿的均匀性和扩展性。
- 解决:两个办法配合用。第一,对偷窃策略增加一个随机扰动项,让部分个体不完全朝向参考解;第二,精英选择时保留上一层前沿里拥挤度大的端点解进入下一代,即使它们离参考解很远。我在主循环里加了
pop_new更新后的随机抖动就是为了对抗这种塌缩,抖动幅度0.01倍变量范围,既能扰动又不破坏边界约束。
5.5 现象四:拥挤度排序方向搞反,前沿两端被砍
- 现象:Pareto前沿画出来总是缺两头,中间堆得很密,两端只零星几个点。
- 原因:拥挤度排序用了默认的
sort升序排列,把拥挤度最小的点放进了下一代,而拥挤度最小恰恰是前沿中间最拥挤的区域,端点的inf永远排在最后面,等于每次选择都在主动删除端点解。这个方向问题不报错、不提示,只会在结果图里体现。 - 解决:精英选择的最后一层用
sort(cd_all(last), 'descend')取拥挤度最大的若干个。这个细节我在第4章主循环里已经写对了,抄代码时注意别手滑改成升序。验证方法很简单:打印第1层每个解的拥挤度,看端点是否为inf,如果不是就说明排序或索引有问题。
6. 验证NSDBO不是黑匣子:收敛性对比与结果导出
6.1 HV与IGD指标:给NSDBO打分
NSDBO跑完不能只看前沿图花哨,要量化它好不好。IGD需要真实Pareto前沿做参考,但微电网调度问题的真实前沿未知,所以更实用的是超体积指标HV。HV计算Pareto前沿与参考点围成的超体积大小,不需要事先知道真实前沿,只需给定目标空间的一个参考点,通常取比前沿最差点再差10%的位置。HV越大说明前沿的收敛性和分布性越好。Matlab里用delaunayn求凸包体积即可,前沿点少于30个时算出的HV不稳定,所以确保Pareto解集在30个以上再评价。
6.2 和NSGA-II做同参数对比:读图的三种判断
验证NSDBO值不值得用,最直接的方法是让它和NSGA-II在完全相同的问题、种群规模、迭代次数下对抗。对比时看三件事:一是HV随迭代的变化曲线,NSDBO的HV如果前期涨得快、后期稳定值更高,说明收敛速度和多样性都有优势;二是最终Pareto前沿的覆盖范围,哪个算法能在成本低和排放低两个端点上走得更远;三是达到同样HV需要多少代,这决定实际工程中值不值得多花时间跑NSDBO。有一种结果需要警惕:NSDBO的HV和NSGA-II持平甚至略低,但这不影响使用,因为NSDBO的收敛速度快,同样800代下它的前沿已经稳定,而NSGA-II还在缓慢提升。这时把迭代数拉长到1500再比一次才公平。我在微电网算例上常见的读图结论是NSDBO前300代HV接近最终值的90%,后期端点的延展性略好于NSGA-II,这是由精英保留加拥挤度选择带来的结果。
6.3 把Pareto解集导出成表格:调度员真正要看的调度单
算法验证完,最后要把推荐方案落成调度员能直接用的表格。折中解对应的决策变量还原成5×24矩阵,加上光伏、风电预测曲线和SOC曲线,拼成一张以时间为行、设备出力和状态为列的调度表。表的列依次是:时段、负荷、光伏、风电、柴油机、燃气轮机、储能充电、储能放电、电网交互、SOC、运行成本、碳排放。这张表导出为CSV或Excel,再画一张24小时功率堆叠图,曲线上各设备出力的正负关系一目了然。
做这一步时我发现一个让人哭笑不得的细节:从Pareto解集里选出的折中解,有时SOC末值落在0.3,下一轮调度衔接要先手动补一次储能回充。所以实际工程里我在约束检查里加了SOC末值约束,让末值和初值偏差不超过0.1,否则第二天没法连续调度。这个约束会让前沿往里缩一点,但换来的是可执行性,值这个价。
NSDBO这套流程走到这里,从问题建模、非支配排序改造、主循环实现到参数排查和验证对比,每一步都有明确产出和验收方法。我自己的习惯是每次换一个新微电网算例,先把目标函数和约束函数单独跑一遍,确认单目标优化结果合理,再套NSDBO跑多目标,这样可以大概率避免把建模错误归咎到算法头上。希望这些踩过的坑和调参经验帮到你,少走几圈弯路。
本文还有配套的精品资源,点击获取