简介：该资源面向电气工程、自动化及数学等专业的大学生与科研人员，提供一套基于非支配排序蜣螂优化算法（NSDBO）求解微电网多目标优化调度的Matlab完整代码，可用于课程设计、期末大作业与毕业设计。微电网调度需兼顾经济成本、能量损耗与排放等多重目标，NSDBO凭借较强的全局搜索与收敛能力，通过非支配排序筛选Pareto前沿，为决策者输出多种可选调度方案。压缩包共32个文件，约5.14MB，含11个m脚本文件（主程序、目标函数、适应度与绘图模块）、17张jpg结果图、2份pdf说明文档、1个xlsx数据表及1个txt文本，结构清晰、便于按模块查阅。代码采用参数化编程，注释详尽，可直接运行并灵活调整参数以适配不同场景。目前已有184人学习，适合希望快速上手微电网多目标优化建模与算法验证的读者参考。

1. 微电网多目标调度为什么需要 NSDBO：从一个 24 小时算例说起

微电网多目标优化调度这件事，真正上手做过的人都知道，难的不是写目标函数，而是当运行成本和碳排放两个目标互相拉扯时，怎么在一天 24 个时段里给风机、光伏、储能、柴发和可切负荷分配功率，让 Pareto 前沿既宽又匀。传统做法是把多目标加权成单目标，权重一变结果就变，调参调到怀疑人生。NSDBO（非支配排序蜣螂优化算法）的思路是把蜣螂优化算法（DBO）的滚球、跳舞、觅食、偷窃四种行为，套进 NSGA-II 的非支配排序和拥挤度框架里，一次跑出一整条前沿。这篇笔记就围绕这个标题，把模型怎么建、NSDBO 怎么落地、Matlab 代码怎么改、参数怎么设、坑在哪，一条线讲透，适合正在做微电网调度、综合能源或群智能算法改进的从业者照着复现。

2. 微电网多目标调度模型：目标函数、约束与 24 时段数据结构

2.1 三个目标函数怎么写才不会被审稿人挑

微电网调度的目标通常不止一个，最常见的是运行成本、碳排放、以及平抑联络线波动。我一般会写成三个目标，但标题里说的是多目标，实际落地时建议先做两目标，跑通后再加第三个，否则前沿维度一高，可视化都费劲。

运行成本目标包含四块：柴发燃料成本、风机光伏运维成本、储能充放电折旧成本、以及从大电网购售电成本。写成离散时段求和：

% 运行成本目标函数（24时段） function f1 = cost_obj(x, data) % x: 决策变量矩阵，行=时段，列=[风机,光伏,柴发,储能充放,可切负荷] N = 24; f1 = 0; for t = 1:N P_wt = x(t,1); P_pv = x(t,2); P_de = x(t,3); P_es = x(t,4); P_load = x(t,5); % 柴发燃料成本：二次函数 a*P^2 + b*P + c f_de = data.a*P_de^2 + data.b*P_de + data.c; % 运维成本 f_om = data.k_wt*P_wt + data.k_pv*P_pv; % 储能折旧 f_es = data.k_es*abs(P_es); % 购售电成本，P_grid>0购电，<0售电 P_grid = P_load + P_es - P_wt - P_pv - P_de; if P_grid >= 0 f_grid = data.price_buy(t)*P_grid; else f_grid = data.price_sell(t)*P_grid; end f1 = f1 + f_de + f_om + f_es + f_grid; end end

逻辑说明：决策变量按列组织，每列对应一个设备，行对应时段，这样后面 NSDBO 的种群初始化直接按矩阵展开。参数说明：data.a/b/c是柴发燃料系数，k_wt/k_pv/k_es是运维和折旧系数，price_buy/price_sell是分时电价，通常峰谷差拉到 3 倍以上才能看出储能套利效果。

碳排放目标类似，柴发和购电是主要来源：

function f2 = carbon_obj(x, data) N = 24; f2 = 0; for t = 1:N P_de = x(t,3); P_es = x(t,4); P_load = x(t,5); P_grid = P_load + P_es - x(t,1) - x(t,2) - P_de; % 柴发碳排放系数 + 电网碳排放系数 f2 = f2 + data.emi_de*P_de + data.emi_grid*max(P_grid,0); end end

第三个目标如果是平抑波动，就写联络线功率的一阶差分绝对值之和，这里不展开，两目标跑通后加一行就行。

2.2 约束条件：功率平衡、储能 SOC 与机组出力上下限

约束是微电网调度里最容易翻车的地方。功率平衡必须每个时段成立，储能 SOC 递推不能越界，柴发有最小出力限制，可切负荷不能超过总负荷的一定比例。

function [g, h] = constraints(x, data) N = 24; g = []; h = []; SOC = zeros(1,N+1); SOC(1) = data.SOC0; for t = 1:N P_wt = x(t,1); P_pv = x(t,2); P_de = x(t,3); P_es = x(t,4); P_load = x(t,5); P_grid = P_load + P_es - P_wt - P_pv - P_de; % 等式约束：功率平衡（用h表示） h = [h; P_grid - (P_load + P_es - P_wt - P_pv - P_de)]; % 不等式约束：柴发上下限 g = [g; data.P_de_min - P_de; P_de - data.P_de_max]; % 储能充放电功率限值 g = [g; -data.P_es_max - P_es; P_es - data.P_es_max]; % SOC 递推 SOC(t+1) = SOC(t) - data.eta*P_es*data.dt/data.cap; g = [g; data.SOC_min - SOC(t+1); SOC(t+1) - data.SOC_max]; % 可切负荷比例 g = [g; P_load - data.load(t); data.load_min - P_load]; end end

参数说明：eta是充放电效率，一般取 0.9 到 0.95；dt取 1 小时；cap是储能容量。SOC 递推里符号要注意，充电时 P_es 为负，SOC 上升，这个符号约定必须在目标函数和约束里保持一致，否则会出现 SOC 越充越低的玄学现象。

2.3 24 时段数据结构与初始化

把风机、光伏预测出力和负荷曲线整理成一个结构体，后面算法调用时直接传：

data.dt = 1; data.cap = 500; % 储能容量 kWh data.SOC0 = 0.5; data.SOC_min = 0.2; data.SOC_max = 0.9; data.P_de_min = 20; data.P_de_max = 200; data.P_es_max = 100; data.a = 0.0005; data.b = 0.3; data.c = 10; data.k_wt = 0.05; data.k_pv = 0.04; data.k_es = 0.02; data.emi_de = 0.65; data.emi_grid = 0.58; data.price_buy = [0.4*ones(1,7), 0.9*ones(1,10), 0.5*ones(1,7)]; data.price_sell = 0.35*ones(1,24); data.load = 300 + 150*sin(linspace(0,2*pi,24)); data.load_min = 0.7*data.load;

这些数值不是标准答案，实际项目里按当地电价和设备参数替换。我一般会先把负荷和风光曲线画出来看一眼，确认峰谷时段和电价峰谷对得上，再往下跑算法。

3. NSDBO 算法实现：非支配排序、拥挤度与蜣螂行为映射

3.1 从 DBO 到 NSDBO：为什么不能直接加权

DBO 本身是单目标优化器，靠滚球、跳舞、觅食、偷窃四种行为更新位置。直接拿它跑加权单目标，权重取 0.5 和取 0.6 得到的解可能差很远，而且一次只能出一个解。NSDBO 的做法是保留 DBO 的位置更新公式，但把适应度评价换成非支配排序和拥挤度距离，种群按 Pareto 等级分层，同一层内按拥挤度选，这样一次运行就能得到分布均匀的前沿。

非支配排序的核心是比较两个解在所有目标上的支配关系：

function [fronts, rank] = non_dominated_sort(costs) % costs: 种群目标值矩阵，行=个体，列=目标数 N = size(costs,1); S = cell(N,1); n = zeros(N,1); rank = zeros(N,1); fronts = {}; for i = 1:N for j = 1:N if i == j, continue; end if all(costs(i,:) <= costs(j,:)) && any(costs(i,:) < costs(j,:)) S{i} = [S{i}, j]; % i 支配 j elseif all(costs(j,:) <= costs(i,:)) && any(costs(j,:) < costs(i,:)) n(i) = n(i) + 1; % i 被 j 支配 end end if n(i) == 0 rank(i) = 1; fronts{1} = [fronts{1}, i]; end end k = 1; while ~isempty(fronts{k}) Q = []; for i = fronts{k} for j = S{i} n(j) = n(j) - 1; if n(j) == 0 rank(j) = k + 1; Q = [Q, j]; end end end k = k + 1; fronts{k} = Q; end end

逻辑说明：S{i}存被 i 支配的个体，n(i)存支配 i 的个体数。第一层是 n=0 的个体，然后逐层剥离。参数说明：目标值矩阵要先归一化，否则量纲差异会让支配关系失真，成本是几百、碳排放是几十，直接比较没问题，但如果加一个波动目标在 0.1 量级，就必须归一化。

拥挤度距离计算：

function cd = crowding_distance(front_costs) [n, m] = size(front_costs); cd = zeros(n,1); for k = 1:m [~, idx] = sort(front_costs(:,k)); cd(idx(1)) = inf; cd(idx(end)) = inf; fmin = front_costs(idx(1),k); fmax = front_costs(idx(end),k); if fmax - fmin < 1e-12, continue; end for i = 2:n-1 cd(idx(i)) = cd(idx(i)) + ... (front_costs(idx(i+1),k) - front_costs(idx(i-1),k))/(fmax - fmin); end end end

边界个体拥挤度设为无穷，保证前沿端点不被淘汰。中间个体按相邻解的目标差归一化累加。

3.2 蜣螂四种行为在调度问题里的映射

DBO 的四种行为对应到微电网调度，本质是搜索步长的不同策略。滚球行为是直线搜索，适合在可行域内快速推进；跳舞行为是遇到障碍后的扰动，对应约束边界附近的探索；觅食行为是向当前最优解靠拢，对应局部开发；偷窃行为是从其他个体附近抢解，对应跳出局部最优。

function Xnew = dbo_update(X, Xbest, iter, maxIter, lb, ub) [N, D] = size(X); Xnew = X; for i = 1:N r = rand; if r < 0.25 % 滚球：沿当前方向前进，带随机扰动 delta = 0.1*(1 - iter/maxIter); Xnew(i,:) = X(i,:) + delta*randn(1,D).*(ub-lb); elseif r < 0.5 % 跳舞：在边界附近小步扰动 theta = rand*pi; Xnew(i,:) = X(i,:) + 0.05*cos(theta)*(ub-lb); elseif r < 0.75 % 觅食：向最优解靠拢 Xnew(i,:) = X(i,:) + rand*(Xbest - X(i,:)); else % 偷窃：从随机个体附近取解 j = randi(N); Xnew(i,:) = Xbest + rand*(X(j,:) - X(i,:)); end Xnew(i,:) = min(max(Xnew(i,:), lb), ub); end end

参数说明：delta随迭代衰减，前期大步探索，后期小步收敛；跳舞的0.05系数控制扰动幅度，太大退化成随机搜索，太小跳不出局部最优。我一般会把四种行为的概率设成 0.25 均分，如果发现收敛太慢，把觅食概率提到 0.4。

3.3 NSDBO 主循环：环境选择与外部存档

主循环把 DBO 更新和非支配排序接起来，每代合并父代和子代，排序后按拥挤度选前 N 个：

function [Xbest, Fbest] = nsdbo(fobj, lb, ub, N, maxIter, D) X = repmat(lb,N,1) + rand(N,D).*(repmat(ub,N,1)-repmat(lb,N,1)); F = zeros(N,2); for i = 1:N, F(i,:) = fobj(X(i,:)); end for iter = 1:maxIter % 选当前前沿中拥挤度最大的作为引导解 [fronts, ~] = non_dominated_sort(F); cd = crowding_distance(F(fronts{1},:)); [~, idx] = max(cd); Xbest = X(fronts{1}(idx),:); % DBO 更新 Xc = dbo_update(X, Xbest, iter, maxIter, lb, ub); Fc = zeros(N,2); for i = 1:N, Fc(i,:) = fobj(Xc(i,:)); end % 环境选择：合并后排序取前 N Xall = [X; Xc]; Fall = [F; Fc]; [fronts, ~] = non_dominated_sort(Fall); newX = []; newF = []; k = 1; while size(newX,1) + length(fronts{k}) <= N newX = [newX; Xall(fronts{k},:)]; newF = [newF; Fall(fronts{k},:)]; k = k + 1; end if size(newX,1) < N cd = crowding_distance(Fall(fronts{k},:)); [~, idx] = sort(cd,'descend'); need = N - size(newX,1); newX = [newX; Xall(fronts{k}(idx(1:need)),:)]; newF = [newF; Fall(fronts{k}(idx(1:need)),:)]; end X = newX; F = newF; end Xbest = X; Fbest = F; end

逻辑说明：每代先选前沿中拥挤度最大的个体作为引导解，避免所有个体挤在同一个区域。环境选择时按 Pareto 等级从高到低填，最后一层用拥挤度截断。参数说明：N一般取 100 到 200，maxIter取 200 到 500，决策变量维度 D 是 24 时段乘以设备数，如果 5 个设备就是 120 维，维度偏高时种群要相应加大。

4. 把 NSDBO 接到微电网模型上：Matlab 代码组织与参数设置

4.1 目标函数包装与决策变量边界

NSDBO 的fobj需要接收一个行向量，返回两个目标值。把第 2 章的模型包装一下：

function F = fobj_wrapper(x, data) N = 24; nDev = 5; X = reshape(x, N, nDev); f1 = cost_obj(X, data); f2 = carbon_obj(X, data); F = [f1, f2]; end

决策变量边界按设备逐个设置，展开成 120 维：

lb = zeros(1,120); ub = zeros(1,120); for t = 1:24 idx = (t-1)*5; lb(idx+1) = 0; ub(idx+1) = data.P_wt_max(t); lb(idx+2) = 0; ub(idx+2) = data.P_pv_max(t); lb(idx+3) = data.P_de_min; ub(idx+3) = data.P_de_max; lb(idx+4) = -data.P_es_max; ub(idx+4) = data.P_es_max; lb(idx+5) = data.load_min(t); ub(idx+5) = data.load(t); end

参数说明：风机光伏上限用预测出力，储能上下限对称，可切负荷下限是保底负荷。边界设置错了，算法再强也跑不出可行解，这一步我一般会单独写个脚本检查每个时段的上下限是否合理。

4.2 约束处理：罚函数还是可行性规则

NSDBO 本身不处理约束，常见做法是罚函数。把约束违反量加到目标上：

function F = fobj_penalty(x, data, pen) N = 24; X = reshape(x,N,5); f1 = cost_obj(X, data); f2 = carbon_obj(X, data); [g, h] = constraints(X, data); vio = sum(max(0,g)) + sum(abs(h)); F = [f1 + pen*vio, f2 + pen*vio]; end

参数说明：pen取 1e4 到 1e6，太小约束压不住，太大目标函数被淹没。我一般先用 1e5 跑一遍，看可行解比例，如果低于 30% 就加大罚系数，或者改用可行性规则：可行解永远优于不可行解，两个不可行解比违反量。

4.3 跑一次完整算例：主脚本与结果可视化

data = init_data(); N = 100; maxIter = 300; D = 120; lb = ...; ub = ...; % 按 4.1 设置 fobj = @(x) fobj_penalty(x, data, 1e5); [Xbest, Fbest] = nsdbo(fobj, lb, ub, N, maxIter, D); % 画 Pareto 前沿 figure; plot(Fbest(:,1), Fbest(:,2), 'o'); xlabel('运行成本'); ylabel('碳排放'); title('NSDBO Pareto 前沿'); % 取折中解 [fronts, ~] = non_dominated_sort(Fbest); cd = crowding_distance(Fbest(fronts{1},:)); [~, idx] = max(cd); x_best = Xbest(fronts{1}(idx),:); X_sch = reshape(x_best, 24, 5);

逻辑说明：跑完后先看前沿形状，如果是一条很窄的线，说明目标之间相关性太强或者种群多样性不够。折中解取拥挤度最大的，代表前沿上分布最稀疏的位置，实际决策时也可以按偏好选端点。

5. 避坑与排查：NSDBO 跑微电网调度最常见的 5 个翻车点

5.1 现象：前沿所有解都挤在一起，看不出 trade-off

原因：目标函数量纲差异太大，成本在几千，碳排放只有几十，非支配排序时成本主导了支配关系。解决：对目标做归一化，或者用自适应归一化，每代按当前种群的最大最小值缩放。

5.2 现象：SOC 越跑越低，最后触底

原因：储能充放电符号约定在目标函数和约束里不一致，或者 SOC 递推公式里效率乘的位置不对。解决：统一约定充电为负、放电为正，SOC 递推写成SOC(t+1) = SOC(t) - eta*P_es*dt/cap，跑之前先用一个固定调度方案验证 SOC 曲线是否合理。

5.3 现象：算法收敛到同一个解，多样性丢失

原因：DBO 的觅食行为概率太高，所有个体都向最优解靠拢。解决：把觅食概率降到 0.2，提高滚球和偷窃概率，或者在更新后加一个高斯扰动，扰动幅度随迭代衰减。

5.4 现象：罚函数法跑出的解违反约束但目标值很好看

原因：罚系数太小，约束违反量乘以罚系数后对目标影响不够。解决：先统计初始种群的约束违反量，把罚系数设成目标值量级的 10 到 100 倍，或者改用可行性规则。

5.5 现象：Matlab 跑了几小时还没出结果

原因：决策变量 120 维，种群 200，迭代 500，每次评价都要循环 24 时段，计算量爆炸。解决：把目标函数向量化，去掉时段循环，用矩阵运算；或者先用 50 种群 100 代跑通流程，再加大参数。我一般会在目标函数里加tic/toc，看单次评价耗时，超过 0.1 秒就要优化。

6. 进阶技巧：用外部存档和自适应参数把 NSDBO 再推一步

跑通基础版之后，如果前沿分布还不够匀，可以加两个东西。第一个是外部存档，把每代产生的非支配解单独存起来，最后合并去重，这样前沿端点不容易丢。第二个是自适应参数，把 DBO 的四种行为概率随迭代调整，前期探索为主，后期开发为主。

% 自适应行为概率 p_explore = 0.5*(1 - iter/maxIter) + 0.1; p_forage = 0.5*(iter/maxIter) + 0.1; % 外部存档 archive.X = []; archive.F = []; % 每代结束后 [fronts, ~] = non_dominated_sort(F); archive.X = [archive.X; X(fronts{1},:)]; archive.F = [archive.F; F(fronts{1},:)]; % 去重：目标值距离小于阈值的只留一个 [~, ia] = unique(round(archive.F,4),'rows'); archive.X = archive.X(ia,:); archive.F = archive.F(ia,:);

验证前沿质量可以用两个指标：超体积（HV）和分布均匀性（Spacing）。HV 越大越好，Spacing 越小越匀。我一般会跑 10 次取平均，如果 HV 波动超过 10%，说明算法稳定性不够，需要调种群大小或者加扰动。

还有一个实际技巧：把 24 时段调度结果画成堆叠面积图，一眼就能看出储能什么时候充、柴发什么时候顶上来。如果储能充放电曲线和电价峰谷完全对不上，那大概率是目标函数里购售电符号写反了，这种血泪经验我踩过不止一次。

最后说个习惯，每次改完参数，先跑 10 代看初始种群分布，别一上来就 500 代，不然等结果的时间够你喝三杯咖啡。希望帮到你。

本文还有配套的精品资源，点击获取