简介一份围绕考虑负荷满意度的微电网运行多目标优化问题的MATLAB代码资源面向电力系统调度与智能优化算法初学者及研究者重点演示NSGA-II二代非支配排序遗传算法在含约束多目标问题中的编码、寻优与结果分析流程。资源包为RAR压缩格式大小约2.32MB页面未给出文件总数内容以.m源码文件为主可在MATLAB中直接运行参考。基于实际算例对参考文献[1]的方法进行复现同时结合Deb的经典NSGA-II算法框架帮助读者理解负荷满意度指标、经济运行目标如何转化为目标函数与约束条件。由于原始模型存在若干可讨论之处资源更适合作为学习非支配遗传算法原理与调试技巧的练习素材微电网调度建模细节仅供参考。目前已有528人学习下载适合需要快速上手多目标优化编程的读者对照验证。1. 负荷满意度与微电网多目标优化复现前先分清解集和目标这个标题的核心不在“多目标”三个字而在“负荷满意度”。很多初读相关论文的人第一版代码就把经济成本、碳排放和满意度加权成一个标量跑完拿一个最优解就交差结果根本解释不了“满意度只降了 5%为什么运行成本能省四分之一”这种问题。因为满意度不是约束它和成本是一对天然冲突的目标压在一起就把这层决策信息抹掉了。把满意度作为第二目标之后微电网调度问题从“找一个解”变成“找一簇解”工作量从建模转移到分析 Pareto 前沿到底把满意度下限放在 0.85 还是 0.95储能循环次数和用户舒适度的交换比是多少靠事后从解集里选。这篇文章按“量化满意度 → 建立双目标模型 → 用 NSGA-II 思路在 MATLAB 里搜前沿 → 处理满意度约束 → 滚动时域扩展”的顺序把这个标题涉及的代码重新立起来。适合正在复现日前优化调度、被多目标算法选型和约束处理卡住的人。2. 负荷满意度量化与双目标优化建模的 MATLAB 索引2.1 负荷满意度的连续量化公式从可削减比例到权重表达切负荷是微电网调度里最常见的满意度牺牲手段所以满意度先要从“切了多少”算起。最直接的量化方式是统计各时段可中断负荷的最大可削减量与实际削减量。设时段 $t$ 的最大可削减量为 $P_{cut,max}(t)$实际削减量为 $P_{cut}(t)$则无权重满意度为$$ S 1 - \frac{\sum_{t1}^{24} P_{cut}(t)}{\sum_{t1}^{24} P_{cut,max}(t)} $$这个式子里的满意度是连续值从 0全切到 1完全不切。它比直接用供需比 $\frac{实际负荷}{需求负荷}$ 更合理因为刚性负荷医院、核心产线根本不在可削减范围内用供需比会把不可控的刚性缺口也归罪于调度策略。如果要区分负荷优先级给每个时段或每类负荷加权重 $w_t$公式写成$$ S 1 - \frac{\sum_{t1}^{24} w_t \cdot P_{cut}(t)}{\sum_{t1}^{24} w_t \cdot P_{cut,max}(t)} $$权重怎么定是复现里的一个关键自由度。常见做法是按负荷价值排序居民负荷权重 0.6商业负荷 0.9工业保底负荷 1.0再按当前时段处于峰平谷来微调。也有论文用模糊隶属度函数把“满意度很好”定义成 $S \geq 0.95$ 时隶属度为 1低于 0.8 后快速下降。这类定义不改变优化框架只改变目标函数曲线形状代码里只需替换这一句% 满意度无权重版本S 的范围 [0, 1] Pcut_max data.Pcut_max; % 24x1各时段最大可削减量 Pcut x(25:48); % 24x1决策变量中的实际削减量 S 1 - sum(Pcut) / max(sum(Pcut_max), eps); % eps 防除零这里Pcut是从决策变量里取出来的data.Pcut_max是算例预设的每时段削减上限来源是历史负荷曲线的一定比例比如取各时段尖峰负荷的 20%。加权重时把sum换成w * Pcut即可注意w要和时段一一对应。2.2 双目标函数设计运行成本与满意度缺失量怎么共存这个标题里“多目标优化”最稳的落地方式是两个目标运行经济性最优、满意度缺失最少。第一个目标度量成本第二个目标度量舒适度损失两者量纲完全不同不能合并只能放在同一坐标系里看 Pareto 折中。运行成本的组成在只含光伏、储能、与大电网购售电的微网里是三项购电费用、售电收益负项、储能充放电循环折算成本。如果算例里还有柴油发电机或微燃机再加燃料费和启停成本。公式表达为$$ C \sum_{t1}^{24} \left( \max(P_{grid,t},0) \cdot \pi_{buy,t} \min(P_{grid,t},0) \cdot \pi_{sell,t} \right) \sum_{t1}^{24} k_{bat} \cdot |P_{bat,t}| $$注意这里不要给切负荷加惩罚项。很多复现代码犯的错误是满意度已经作为第二个目标存在又在成本函数里对切负荷罚钱等于同一个动作被惩罚两次Pareto 前沿会被压得不成形甚至出现“成本差几千元、满意度几乎不动”的假折中。第二个目标选最小化满意度缺失量即 $f_2 1 - S$。前面公式里 $S$ 接近 1 时$1-S$ 接近 0这样两个目标默认都是“越小越好”可以直接喂给 MATLAB 的gamultiobj不需要做取反。两个目标的冲突在分时电价下非常直观峰时电价 1.35 元时削减 100 kW 负荷能省下 135 元代价是满意度下降而谷时电价 0.38 元切负荷几乎省不到钱满意度损失就显得很亏。优化算法会在所有这样的时段组合里搜索最后给出从“高成本高满意度”到“低成本低满意度”的一整条前沿。2.3 微电网运行约束清单与 MATLAB 索引对应约束条件的处理方式直接决定代码能不能收敛。把这组约束分成三类功率平衡等式、储能动态不等式、联络线上下限。下表给出一套典型约束和它在 MATLAB 里的实现位置约束类型数学表达代码处理位置功率平衡$P_{grid} P_{pv} P_{bat} P_{load} - P_{cut}$罚函数进目标 $f_1$SOC 上下限$SOC_{min} \le SOC(t) \le SOC_{max}$非线性不等式c充放电功率$-P_{ch,max} \le P_{bat} \le P_{dis,max}$决策变量LB/UB联络线功率$0 \le P_{grid} \le P_{grid,max}$决策变量LB/UB满意度下限$S \ge S_{min}$非线性不等式c末时 SOC$SOC(24) \ge SOC_{end,min}$非线性不等式c这里最需要注意功率平衡。它在物理上是硬等式但遗传算法处理等式约束收敛极慢尤其 72 维变量、200 个种群时等式容差设置不对会直接弹 “No feasible solution”。常见做法是把功率平衡残差的平方乘以一个大系数加进目标 $f_1$变成软约束balance Pgrid Ppv(:) Pbat - (Pload(:) - Pcut); penalty 1e6 * sum(balance.^2); f1 cost_grid cost_bat penalty;系数 1e6 不是随便拍的。它要远大于正常成本量级万元级保证任何不可行解的成本都被抬高到不可能被选中的程度同时又不至于让目标函数值超出 double 精度范围。算例里成本大概在 1e45e4 元罚项从 1e6 起步能把平衡误差压到 0.1 kW 以内。储能 SOC 的递推是这段模型里最容易写错的地方。充放电效率不是同一个值放电时 SOC 下降要除以放电效率充电时 SOC 上升要乘以充电效率SOC(1) 0.5; % 初始 SOC eta_ch 0.95; eta_dis 0.95; for t 1:24 if Pbat(t) 0 % 放电 SOC(t1) SOC(t) - Pbat(t) * dt / (eta_dis * Cap); else % 充电Pbat 为负 SOC(t1) SOC(t) - Pbat(t) * eta_ch * dt / Cap; end enddt在日前调度里是 1 小时Cap是储能容量kWh。功率单位用 kW、容量用 kWh 时SOC 变化量无量纲正好落在 [0,1]。这段递推不能写成SOC SOC - Pbat*dt/Cap那种统一效率写法虽然能跑但充放电损耗的物理方向不对复现论文数据会对不上。3. 用 MATLAB 优化工具箱的 gamultiobj 构建 NSGA-II 求解框架3.1 多目标冲突下为什么选 NSGA-II 思路而不是线性加权线性加权法把两个目标乘权重再加起来本质上变成单目标优化跑一次只有一组解。权重系数 $\lambda$ 取 0.8 还是 0.3对结果的影响比算法本身还大而且无法解释“为什么这个权重就是对的”。NSGA-II 的思路是一次性维护一整批解靠非支配排序和拥挤距离同时保证收敛性和分布性解之间互相比较最后给出完整 Pareto 前沿。MATLAB 里直接能用的是优化工具箱的gamultiobj它实现的是带约束的 NSGA-II 变体。如果不需要高度定制交叉变异算子工具箱版本比手动写 NSGA-II 稳定得多尤其在约束容差和并行计算的处理上。手动版适合要改编码方式、要输出每代非支配层面的场景但作为一个复现项目先用gamultiobj跑通再在它基础上加自定义算子是最省时间的路径。用gamultiobj之前必须想清楚一个事决策变量怎么排。3 类变量、24 个时段顺序不一致会导致目标函数里索引对不上这是最常见的“代码跑不出合理结果”的原因。下一节给出固定约定。3.2 72 维决策变量编码储能出力、削减量与购电功率决策变量排成一个列向量 $x \in \mathbb{R}^{72}$前 24 维是储能出力中间 24 维是各时段削减量最后 24 维是购电功率。约定变量块MATLAB 索引单位含义上下界Pbatx(1:24)kW储放出力正放电、负充电[-120, 100]Pcutx(25:48)kW可中断负荷削减量[0, Pcut_max]Pgridx(49:72)kW与大电网交互功率正为购电[0, 250]上下界直接传给gamultiobj的LB和UB作用是替代一部分线性约束不用写进约束函数。注意Pbat的下界是 -120、上界是 100因为储能充电功率和放电功率不必对称实际锂电池的 1C 充放对应 600 kWh 容量就是 600 kW这里取 120 是人为限制防止储能过度参与功率调节导致循环寿命成本暴涨。Pcut的上界不是定值而是每个时段各不同所以UB(25:48)要用data.Pcut_max(:)逐维赋值不能直接写标量。Pgrid只允许正方向购电是因为典型算例里微网以消纳光伏和储能为主反向馈电虽可获售电收益但会复杂化功率平衡和电压约束先拿掉后续要加只需把Pgrid的LB改成负值、在成本里加售电项。3.3 gamultiobj 最小可跑主脚本选项设置与调用方式搭建一个最小框架只需要三个文件参数结构体、目标函数、主脚本。下面这个主脚本是完整可跑的骨架所有自定义函数名和参数都对得上% 主脚本 demo_mg_multiopt.m clear; clc; rng(42); data load_case_data(); % 导入算例参数 nvars 72; LB [-120*ones(24,1); zeros(24,1); zeros(24,1)]; UB [ 100*ones(24,1); data.Pcut_max(:); 250*ones(24,1)]; options optimoptions(gamultiobj, ... PopulationSize, 200, ... MaxGenerations, 300, ... ParetoFraction, 0.35, ... ConstraintTolerance, 1e-4, ... FunctionTolerance, 1e-6, ... Display, iter, ... UseParallel, true); [x_pareto, fval_pareto] gamultiobj( ... (x) mg_obj_fun(x, data), nvars, ... [], [], [], [], LB, UB, ... (x) mg_con_fun(x, data), options); scatter(fval_pareto(:,2), fval_pareto(:,1), 8, filled); xlabel(满意度缺失 1-S); ylabel(运行成本 / 元);这段代码里rng(42)固定随机种子保证复现时每次跑出来的前沿不漂移。PopulationSize取 200低于 100 时 72 维解空间覆盖不住前沿会出现明显空洞MaxGenerations取 300在 200 种群下通常 150 代后前沿就稳定了。ParetoFraction是保留在非支配前沿上的种群比例0.35 表示 200 个个体里留 70 个最终解太少会漏掉折中段太多会把非支配层里的劣解也混进来。UseParallel开并行之前要parpool或者至少确认并行工具箱可用否则只在单核上跑。目标函数里尽量把矩阵运算写成向量化并行才能吃到甜头。Display设为iter能实时看到每代的最佳前沿数量和平均代价收敛满不满一目了然。3.4 等式约束在遗传算法里的罚函数替代gamultiobj官方支持非线性约束函数[c, ceq] constraints(x)但对等式约束的默认容差非常苛刻通常设到 1e-6 甚至 1e-8。功率平衡这种带充电效率、负荷波动的等式遗传算法交叉变异后很难精确卡到零结果就是大量个体被判不可行种群退化到只在可行域边缘探索。所以常见的复现做法是把功率平衡从ceq里拿掉改成罚函数项写进目标函数就像 2.3 节里penalty 1e6 * sum(balance.^2)那样。约束函数里只保留不等式SOC 上下限、满意度下限、末时 SOC 下限。这样ceq永远返回[]算法绕开了遗传算法处理等式的天生短板。function [c, ceq] mg_con_fun(x, data) % 只需要 SOC 与满意度不等式功率平衡走罚函数 T 24; Pbat x(1:T); Pcut x(T1:2*T); S 1 - sum(Pcut) / max(sum(data.Pcut_max), eps); soc zeros(T1, 1); soc(1) 0.5; for t 1:T if Pbat(t) 0 soc(t1) soc(t) - Pbat(t) / (data.eta_dis * data.Ccap); else soc(t1) soc(t) - Pbat(t) * data.eta_ch / data.Ccap; end end c [data.soc_min - soc(2:end); soc(2:end) - data.soc_max]; c [c; data.S_min - S]; % 满意度下限约束 c [c; data.soc_end_min - soc(end)]; % 末时 SOC 下限 ceq []; end注意S在这里同时被目标函数和约束函数计算重复计算开销很小但保持两处公式完全一致很重要。曾经见过一处用加权满意度、一处用无权重导致前沿在 $S0.9$ 附近出现一道不自然断裂找了两天才定位到公式不一致。4. 典型算例参数设置与含满意度约束的完整代码实现4.1 典型微电网算例参数与默认曲线复现类代码最怕参数设置和原文章不一致。这里给一套在光伏储能微电网文章里最常见的典型数据单位统一用 kW、kWh、元。储能容量 600 kWh 对应一个中等规模园区微网联络线 250 kW 意味着与大电网交互受限只能部分依赖外购电。参数取值说明光伏额定功率250 kW曲线峰在 12:00~14:00储能容量600 kWh锂电池SOC 范围 [0.1, 0.9]最大充/放电功率120 / 100 kW非对称防循环过深初始 SOC0.5末时要求 ≥ 0.2峰时电价1.35 元/kWh10:00-15:00、18:00-22:00平时电价0.83 元/kWh07:00-10:00、15:00-18:00谷时电价0.38 元/kWh22:00-次日 07:00售电价0.50 元/kWh本算例不用后续扩展可削减上限各时段负荷 20%Pcut_max按此生成储能循环损耗0.02 元/kWh折算到每次充放电data结构体里还需要存dt1、eta_cheta_dis0.95。储能的循环损耗系数在代码里乘的是abs(Pbat)即无论充或放每流过 1 kWh 能量计 0.02 元这是把电池寿命折算成运行成本的最简方式。4.2 分时电价、光伏与负荷数据的 MATLAB 预处理数据预处理直接写进load_case_data()避免在主脚本里放一长串硬编码数组。分时电价按 24 点手动填成一个列向量光伏出力用一个钟形近似生成负荷用一个双峰曲线做基准function data load_case_data() data.dt 1; % 小时 data.eta_ch 0.95; data.eta_dis 0.95; data.Ccap 600; data.soc_min 0.1; data.soc_max 0.9; data.soc_end_min 0.2; data.cost_bat_kwh 0.02; % 分时电价买入谷0.38 / 平0.83 / 峰1.35 price 0.83 * ones(24,1); price(10:15) 1.35; price(18:22) 1.35; price(23:24) 0.38; price(1:7) 0.38; data.price_buy price; data.price_sell 0.50 * ones(24,1); % 光伏出力白天钟形额定250kW t 0:23; data.Ppv round(250 * exp(-((t - 12).^2) / 18)); data.Ppv(data.Ppv 0) 0; % 负荷早晚双峰峰约 280kW Pbase [120 110 105 100 100 110 140 190 240 260 250 ... 230 215 200 210 225 235 245 280 260 220 180 150 130]; data.Pload Pbase; data.Pcut_max round(0.2 * Pbase); % 可削减上限为负荷20% data.S_min 0.85; % 满意度下限约束 end注意光伏出力数组生成后要把负值截断round保证整型数组在后续计算里不会产生小数点位数不一致的问题。负荷曲线和光伏曲线都是示意性数据实际复现时应该替换成目标文章的日曲线替换点只在data.Ppv和data.Pload两行不影响其他代码。4.3 含满意度约束的目标函数完整 MATLAB 实现目标函数是整个代码的核心。前面已经拆过成本组成这里给出完整可运行版本。变量索引、罚函数、满意度公式、SOC 递推全部揉在一起但分段清晰function [f, c, ceq] mg_obj_fun(x, data) T 24; Pbat x(1:T); % 储能出力正放负充 Pcut x(T1:2*T); % 削减负荷量 Pgrid x(2*T1:3*T); % 购电功率非负 % 购电成本只买不卖 cost_grid sum(Pgrid .* data.price_buy(:)); % 储能循环损耗成本 cost_bat data.cost_bat_kwh * sum(abs(Pbat)); % 功率平衡误差 - 罚函数 balance Pgrid data.Ppv(:) Pbat - (data.Pload(:) - Pcut); penalty 1e6 * sum(balance.^2); % 目标1总运行成本含罚项 f1 cost_grid cost_bat penalty; % 目标2满意度缺失量 S 1 - sum(Pcut) / max(sum(data.Pcut_max), eps); f2 1 - S; f [f1; f2]; % 不等式约束SOC 与满意度下限 soc zeros(T1, 1); soc(1) 0.5; for t 1:T if Pbat(t) 0 soc(t1) soc(t) - Pbat(t) / (data.eta_dis * data.Ccap); else soc(t1) soc(t) - Pbat(t) * data.eta_ch / data.Ccap; end end c [data.soc_min - soc(2:end); soc(2:end) - data.soc_max]; c [c; data.S_min - S]; c [c; data.soc_end_min - soc(end)]; ceq []; endf返回一个 2×1 向量gamultiobj会自动识别为两目标。penalty的 1e6 系数让任何不平衡解的成本直接冲到百万级和正常成本的万元级差开两个数量级以上。SOC 递推使用 for 循环而不是向量化因为充放电效率的分支判断让向量化反倒更难读24 步循环的开销在每代 200 个个体下可以忽略。满意度约束data.S_min - S在S0.84时是正值算法判定不可行SOC约束同理。ceq保持空数组所有硬约束已经通过罚函数在其他路径消化掉这里刻意不做二次兜底防止罚函数和ceq对同一个量双重判定导致可行域变形。4.4 gamultiobj 运行选项与报错排查运行主脚本时最常见的报错集中在两类约束导致无可行解以及目标函数返回维度错误。下表是排错对照现象可能原因处理办法输出提示No feasible solutionS_min设太高或 SOC 初值不合理先把S_min降到 0.8 试跑Pareto 前沿只有一个点ParetoFraction太小或种群退化回调到 0.3~0.4加MaxGenerations前沿明显不连续、成两团目标函数里满意度上下限公式不一致核对两个函数的S表达式逐字符一致成本异常大百万级罚项残留检查balance是否真为零减小ConstraintTolerance并行运行变慢目标函数里有带副作用的全局变量移除persistent全部数据传data结构体ConstraintTolerance默认 1e-3对只有不等式约束的问题够用但如果翻代码觉得可行域被侵蚀试着收紧到 1e-4代价是每代收敛速度慢 5% 左右。FunctionTolerance1e-6 是稳定的默认值不用频繁调整。5. Pareto 前沿读取、满意度下限分析与滚动时域改造5.1 从 Pareto 前沿选取调度方案的拐点距离法fval_pareto的每一行是一个非支配解的两个目标值。直接用肉眼从散点图选点不严谨常用做法是“到理想点最短距离法”把两个目标各自的理想最优值作为原点归一化后计算每个解到原点的欧氏距离距离最小的点就是折中最好的方案。fp fval_pareto; fp_norm (fp - min(fp)) ./ (max(fp) - min(fp) eps); dist sqrt(sum(fp_norm.^2, 2)); [~, idx] min(dist); best_x x_pareto(idx, :); Pbat_opt best_x(1:24); Pcut_opt best_x(25:48); Pgrid_opt best_x(49:72);归一化前一定要先看两个目标的量级。成本在 1 万4 万元满意度缺失在 00.2不归一化的话距离完全被成本统治选出来的解还是“最便宜”的。归一化后两个目标的权重天然各占一半这是无偏好选择如果决策者倾向保满意度可以改成0.3 * 归一化成本 0.7 * 归一化满意度缺失公式只改一行。5.2 满意度下限 S_min 对最优前沿形态的影响data.S_min从 0.85 调到 0.95前沿形态会发生一个可预见的移动高满意度区间的成本整体抬升前沿向左压缩最后可能干脆丢失低满意度解。这是约束和目标同时作用的结果。反过来把S_min降到 0.7前沿右端会出现一批极度激进的切负荷方案成本可能比 0.85 时低一半但在实际场景里几乎没有调度员敢用。调试这个参数时建议每次只改data.S_min跑完记录三个指标前沿右端点成本、左端点满意度缺失、折中解对应的总切负荷量。三者放在一张表里就能看出满意度约束是在约束边界起作用还是已经松到不影响优化了。如果 0.85 和 0.75 的前沿完全重合说明在最优点附近算法本来就不会把满意度压到 0.75 以下S_min形同虚设这时候就该把它收得更紧让约束真正参与搜索。5.3 滚动时域改造用热启动把日前优化变成日内优化日前 24 小时优化是静态的但实际运行中光伏预测和负荷预测每 15 分钟都会刷新固定的一组调度指令跑一天根本不现实。滚动时域改造是这套代码最自然的进阶方向把 24 小时切成滚动窗口每步只优化未来 16 小时执行当前时段然后推进窗口。MATLAB 里直接复用已经写好的mg_obj_fun和mg_con_fun只需要在主循环里更新data.Ppv、data.Pload的预测段。关键是让上一轮的 Pareto 解参与下一轮初始种群用热启动避免每步都从头搜索for k 1:48 % 每半小时滚动一次 data_k update_forecast(data, k); % 更新预测窗口 if k 1 options optimoptions(options, ... InitialPopulationMatrix, ... [x_pareto(1:50, :); rand(150, nvars)]); end [x_pareto, fval_pareto] gamultiobj( ... (x) mg_obj_fun(x, data_k), nvars, ... [], [], [], [], LB, UB, ... (x) mg_con_fun(x, data_k), options); apply_control(x_pareto(1, :), k); endInitialPopulationMatrix要求行数不小于种群的 10%这里取前 50 个非支配解加上 150 个随机个体既保留了上一轮优势区域的记忆又给预测误差留了探索空间。实际测试中热启动版本的收敛代数比清空重跑少一半以上尤其当电价时段切换导致最优区域突变时上一轮的优质解仍能提供有效梯度。滚动窗口的时长按S_min的敏感性调满意度约束越紧窗口越短越安全。本文还有配套的精品资源点击获取