微电网调度中的NSGA-II多目标优化：成本与碳排放的权衡

简介面向电力系统多目标调度优化场景这份资源提供了基于非支配排序遗传算法NSGA-II求解微电网调度问题的完整Python与Matlab实现。需要说明的是资源内容覆盖微电网中风电、光伏、储能设备与电动汽车的协调调度并同时考虑经济成本、可再生能源利用率和供电可靠性等多个相互冲突的目标函数最终输出帕累托最优解集帮助决策者权衡不同运行策略。压缩包内共6个文件其中4个Python脚本负责非支配排序、拥挤距离、交叉变异操作以及复杂约束下十进制与二进制编码的调度模型实现另2个Matlab脚本提供仿真对照与结果分析方便跨语言验证算法性能。整个资源包仅27KB代码紧凑适合快速上手多目标优化在微电网调度中的应用。目前已有2059人学习下载资源中给出了从模型构建、目标函数定义到算法迭代的完整思路是电力系统优化领域性价比很高的学习案例。1. 微电网调度的冲突从哪来NSGA-II 能给出什么微电网调度要解决的核心问题是未来 24 小时或 48/96 个时段内每台柴油机发多少电、储能何时充放电。给定光伏和负荷预测曲线之后经济性与低碳排放这两个目标并不一致午后光伏充足时给储能充电是省钱的但电池循环损耗会折进成本晚高峰多开柴油机虽然缺电风险小但燃料费用和 CO2 排放同时上升。单目标加权法把两个目标乘上权重后拍成一个标量权重一旦定死优化器只会在这一组权重对应的折衷点上取样你拿不到完整的权衡信息。想看清楚“更低碳的方案到底贵多少”就得反复调权重、反复跑优化。NSGA-II 不转单目标而是直接把“成本最低”和“排放最少”作为并列目标用非支配排序把种群划分成多层 Pareto 前沿再用拥挤度距离保证解在目标空间里散开。一轮优化得到一个覆盖多个折衷点的解集合决策者可以事后根据电价、政策、电池状态从中挑一条。对微电网调度这种变量多、约束多、目标只有两三个的问题它在 Python 里手写不过几百行在 MATLAB 里也能直接用优化工具箱的gamultiobj跑同一套模型。下文先讲清楚 NSGA-II 的机制再把调度模型写成两个版本代码最后给出验证解质量的具体做法。2. NSGA-II 运行机制与微电网调度的对应关系2.1 支配关系决定了“谁更值得留下”两个解 A 和 B各自带着两个目标值成本 f1、排放 f2。如果 A 的 f1 和 f2 都不比 B 差并且至少一个严格更小就说A 支配 B。在微电网调度场景里解 A 可能是“柴油机出力更小、储能深度充放”的方案解 B 是“柴油机出力大、储能少动”的方案柴油机出力直接影响燃料成本和排放所以这两类解之间存在大量支配关系。判断完后不被任何个体支配的解构成第一层前沿去掉第一层后再做同样的判断得到第二层依次类推。快速非支配排序的核心是维护两个数组被哪些个体支配的列表S[i]以及支配它的个体数n[i]。n[i]归零时个体进入当前层。下面是一段可以直接复用的 Python 代码def dominates(ind1, ind2, n_obj): better False for i in range(n_obj): if ind1[i] ind2[i]: return False elif ind1[i] ind2[i]: better True return better def fast_non_dominated_sort(fvals): # fvals: (N, n_obj)返回每一层的个体索引列表 N len(fvals) n_obj fvals.shape[1] S [[] for _ in range(N)] n [0] * N fronts [[]] for i in range(N): for j in range(N): if i j: continue if dominates(fvals[i], fvals[j], n_obj): S[i].append(j) elif dominates(fvals[j], fvals[i], n_obj): n[i] 1 if n[i] 0: fronts[0].append(i) k 0 while fronts[k]: nxt [] for i in fronts[k]: for j in S[i]: n[j] - 1 if n[j] 0: nxt.append(j) k 1 fronts.append(nxt) return fronts[:-1]这段代码把“谁支配谁”完整计算了一遍复杂度是 O(N²·m)N 是种群规模m 是目标数。对微电网调度这种 N 在 100~200、m2 的问题单代排序耗时可以忽略如果以后目标数涨到 5 个以上就要换成基于角度的分解方法这里不展开。2.2 拥挤度距离保证前沿不聚成一团同一层内的个体互不支配不能只用支配关系排序需要第二个指标拥挤度距离。做法是把同一层的个体按每个目标分别升序排列相邻两个个体在该目标上的差值做归一化累加到当前个体上边界个体直接赋无穷大防止被淘汰。这个“边界赋无穷大”的细节很关键。它保证种群始终覆盖目标空间的两个极端——调度问题里就是“极端省钱”和“极端低碳”两个方向。如果没有这一条演化后期会把计算资源全部集中到前沿中段两端曲线缺失决策者就看不到全貌。拥挤度的 Python 实现如下def crowding_distance(fvals, front): # front 是某层个体索引列表 dist {i: 0.0 for i in front} for m in range(fvals.shape[1]): order sorted(front, keylambda i: fvals[i, m]) dist[order[0]] float(inf) dist[order[-1]] float(inf) norm fvals[order[-1], m] - fvals[order[0], m] if norm 1e-12: continue for k in range(1, len(order) - 1): dist[order[k]] (fvals[order[k1], m] - fvals[order[k-1], m]) / norm return dist归一化分母是当层该目标的最大值减最小值。如果某一代所有个体在成本目标上几乎一样norm接近 0说明前沿在该维上已经没有扩展空间此时跳过这一维避免除零。2.3 交叉变异算子与调度曲线编码的适配NSGA-II 本身不限制编码方式但微电网调度里决策变量是连续功率适合实数编码。二进制编码把连续量离散化染色体长度膨胀交叉变异后产生大量违背爬坡约束的向量还需要额外修复收益很小。实数编码配合模拟二进制交叉SBX子代在可行域内连续移动能够保留“相邻时段出力接近”这种曲线形状特征。我一般这样设置参数算子/参数常用设置在调度问题中的作用实数编码决策变量连续取值柴油机出力、储能功率是连续量SBX 交叉交叉概率 0.9分布指数 15~20子代曲线接近父代保持负荷跟随节奏多项式变异变异概率 1/n分布指数 20微调单个时段出力扰动幅度可控锦标赛选择规模 2按“层号小优先层号相同拥挤度大优先”比较交叉概率设到 0.9 是因为调度曲线之间存在大量相似结构子代从两个父代各取一部分“时段片段”组合能保留白天的光伏充电视图和晚高峰的柴油机补电策略。变异概率设为1/nn 是决策变量数48 个变量时大约每个个体变异 1 个时段不会把整条曲线打散。2.4 精英保留与一代一代的演化父代 N 个、子代 N 个合并成 2N先按非支配层逐层填充新种群同一层内按拥挤度从大到小填充。这样上一代的 Pareto 最优解无论如何都不会被直接淘汰演化后期约束难满足但目标优秀的个体也有机会保留。对微电网调度这种大量约束来自爬坡和储能 SOC 的问题这一条比单目标遗传算法的“最优个体直接复制”更关键它保留的不是一个解而是一组解。3. 微电网调度模型的数学化决策变量、双目标与约束3.1 决策变量按时间序列展开以 24 小时为一个调度周期、步长 1 小时。为让模型可复现我做一个明确简化单台柴油发电机加一组储能光伏和风电按“不可调度电源”处理预测曲线给定不放进决策变量。决策变量 x 的长度是 48前 24 维是 t 时刻柴油机出力P_DG(t)后 24 维是储能净功率P_ES(t)正值表示放电、负值表示充电。这种编码方式直接对应功率平衡约束任意时刻负荷等于光伏、风电、柴油机、储能四者之和不多引入机组启停的整数变量。如果以后要加外网购电只需要把变量扩成 72 维多出的 24 维作为购电功率并在目标函数里加上分时电价项即可。3.2 两个目标函数运行成本与碳排放成本目标由三部分组成柴油机燃油成本、设备运行维护成本、储能循环退化成本。燃油成本用二次函数近似是因为柴油机在低负载区油耗率高二次项能体现“低载不经济”的特性f1 Σ [ a * P_DG(t)^2 b * P_DG(t) c k_om * |P_ES(t)| ]其中a0.02, b0.35, c0.8是拟合系数k_om0.01是储能单位充放电量的维护成本。储能退化成本用|P_ES(t)|的线性项近似更精细的做法是把 SOC 和放电深度做成二维查表但对 NSGA-II 来说目标函数多一个非线性查表会显著拖慢适应度评估。排放目标只统计柴油机f2 e_dg * Σ P_DG(t)e_dg0.85 kgCO2/kWh。光伏、风电和储能本身不产生直接排放如果模型里加入外网购电还需要按电网边际排放系数折算成第二项。3.3 约束条件与常见的两种处理思路约束类型表达式处理方式功率平衡P_PV P_WT P_DG P_ES P_L保留为等式约束或用罚函数柴油机出力限值0 ≤ P_DG(t) ≤ 50 kW直接设 lb/ub储能功率限值|P_ES(t)| ≤ 50 kW直接设 lb/ub爬坡约束|P_DG(t1) - P_DG(t)| ≤ 15 kW/h罚函数或修复SOC 动态SOC(t1) SOC(t) - η·P_ES(t)·Δt/E_cap迭代计算周期始末 SOCSOC(24) 与 SOC(0) 偏差在 ±0.02 内容差型等式约束功率平衡是等式约束gamultiobj对等式约束的支持不算友好我一般把它转成一个不等式对不平衡量在 ±1 kW 内视为满足即|P_PV P_WT P_DG P_ES - P_L| ≤ 1。爬坡约束和 SOC 约束同理把严格等式放宽成带容差的不等式遗传算法的可行域不再是一条流形搜索效率会好很多。3.4 MATLAB 目标函数文件怎么组织给gamultiobj用的目标函数必须返回一个行向量每个元素是一个目标。这里的x就是 48 维决策向量function f microgrid_obj(x) P_DG x(1:24); P_ES x(25:48); % 正值放电负值充电 a 0.02; b 0.35; c 0.8; k_om 0.01; fuel sum(a * P_DG.^2 b * P_DG c); om k_om * sum(abs(P_ES)); f(1) fuel om; e_dg 0.85; f(2) e_dg * sum(P_DG); end建模时把P_DG和P_ES分开切片操作比把两套变量混在一起写更清晰。后面加外网购电时只需要扩展索引范围不需要改前面两行的含义。4. Python 手写 NSGA-II算子、主循环与约束处理4.1 最小可用的算子实现第 2 章已经给出非支配排序和拥挤度这里补上交叉变异和主循环。SBX 交叉的关键是分布指数eta_c它控制子代偏离父代的程度值越大子代越接近父代。多项式变异的分布指数同理。代码import numpy as np def sbx_crossover(p1, p2, eta_c15): u np.random.random(p1.shape) beta np.where(u 0.5, (2.0 * u) ** (1.0 / (eta_c 1)), (1.0 / (2.0 * (1.0 - u))) ** (1.0 / (eta_c 1))) c1 0.5 * ((1 beta) * p1 (1 - beta) * p2) c2 0.5 * ((1 - beta) * p1 (1 beta) * p2) return c1, c2 def poly_mutate(child, lb, ub, eta_m20, p_m0.02): mask np.random.random(child.shape) p_m delta np.random.random(child.shape) r np.where(delta 0.5, (2.0 * delta) ** (1.0 / (eta_m 1)) - 1.0, 1.0 - (2.0 * (1.0 - delta)) ** (1.0 / (eta_m 1))) child child mask * r * (ub - lb) return np.clip(child, lb, ub)sbx_crossover的beta在 0 附近时子代会坍缩到父代均值附近beta很大时子代会偏向两个父代的边界。微电网调度里傍晚负荷爬坡段的出力曲线边界往往是决策的敏感区SBX 比均匀交叉更容易探索到这些边界位置。主循环固定五步初始化种群 → 评估目标 → 锦标赛选择父代 → 交叉变异生成子代 → 合并后按非支配排序和拥挤度选出下一代。锦标赛选择里比较规则是“层号小优先层号相同拥挤度大优先”def tournament_selection(pop, fronts, crowd_dist, k2): idx np.random.choice(len(pop), sizek, replaceFalse) best idx[0] for i in idx[1:]: rank_i next(r for r, front in enumerate(fronts) if i in front) rank_best next(r for r, front in enumerate(fronts) if best in front) if rank_i rank_best: best i elif rank_i rank_best and crowd_dist[i] crowd_dist[best]: best i return pop[best]crowd_dist是字典键是全局个体索引。这段选择逻辑每次要从fronts里线性查找在种群 200、代数 300 的量级下不是瓶颈但写成rank {ind: r for r, front in enumerate(fronts) for ind in front}的字典后查找会更快一点。4.2 约束处理不要只挂罚函数罚函数是最快的做法把功率不平衡量、爬坡越限量、SOC 偏差量分别乘上大系数加到 f1 上。但系数给得太大Pareto 前沿会被“推”向目标值很大的区域给得太小大量不可行解混进前沿。我在手写版本里倾向于混合策略功率平衡允许 ±1 kW 偏差超出部分罚1000 * 偏差²柴油机爬坡不做罚函数而是产生子代后直接钳制P_DG(t1) clip(P_DG(t1), P_DG(t)-15, P_DG(t)15)SOC 末端偏差罚500 * max(0, |SOC_end - SOC_start| - 0.02)²。爬坡约束用钳制而非罚函数是因为它只在相邻时段之间起作用前期不满足后期也可能被修复钳制直接生成可行结构几乎不损失搜索能力。SOC 末端偏离用带容差的罚函数因为从数学上看它是一个通过储能运行轨迹累加出来的状态约束没法在单个时段上修复只能靠罚函数引导。这部分约束处理就在主循环里evaluate_individual函数内完成罚函数值直接加在第一目标上。两个目标的量级差异要留意成本目标在几百到几千元范围排放目标是几十到几百 kg罚函数统一加到成本目标上系数选择可以粗一些。4.3 种群规模与演化代数的经验值参数推荐范围说明种群规模100~20048 维决策变量太小覆盖不了完整前沿演化代数200~400看前沿端点是否还在移动交叉概率0.9调度曲线相似结构多交叉是最主要搜索手段变异概率1/48 ≈ 0.02平均每个个体变异一个时段SBX 分布指数15~20越大子代越接近父代多项式变异分布指数20与 SBX 搭配保持局部微调能力如果跑 300 代后前沿端点还在移动说明代数不够或者变异概率太低导致探索能力不足。前端点移动的判断方法是把每隔 50 代的第一层fval最小值打出来看两个目标的最小值是否还在下降。端点不动不代表前沿已经收敛还要盯中段点的分布是否均匀这一步等第 5 章验证部分再说。5. 用 MATLAB gamultiobj 跑同一模型并验证前沿质量5.1 gamultiobj 入口与三个常用选项MATLAB 优化工具箱里的gamultiobj本身就是 NSGA-II 的商用实现。目标函数用第 3 章的microgrid_obj再写一个非线性约束函数入口是这样options optimoptions(gamultiobj, ... PopulationSize, 120, ... MaxGenerations, 300, ... ParetoFraction, 0.5, ... Display, iter); nvars 48; lb [zeros(1, 24), -50 * ones(1, 24)]; ub [50 * ones(1, 24), 50 * ones(1, 24)]; [x, fval] gamultiobj(microgrid_obj, nvars, ... [], [], [], [], lb, ub, microgrid_constr, options);ParetoFraction是最终保留的前沿个体占种群的比例设 0.5 表示 120 个个体里保留约 60 个非支配解。设太高会让中段点位很密但两端缺失设太低则前沿形状粗糙按 0.4~0.6 之间调。gamultiobj的约束函数需要返回[c, ceq]我建议全部用ceq[]把 SOC 末端约束写成带容差的不等式对放进cfunction [c, ceq] microgrid_constr(x) P_ES x(25:48); E_cap 200; eta 0.95; SOC zeros(25, 1); SOC(1) 0.5; for t 1:24 SOC(t1) SOC(t) - P_ES(t) / E_cap * eta; end c [abs(SOC(25) - 0.5) - 0.02]; ceq []; end注意eta放在 SOC 递推里写法不同这里假设充电效率等于放电效率功率为正表示放电SOC下降方向与P_ES不一致时需要先确认符号约定再在目标函数和约束函数里保持同一个符号体系。5.2 前沿可视化与选解的三个实用技巧跑完后立刻画fval散点图横轴是成本纵轴是排放。三个判断点前端点是否和单目标优化的最小值接近中段点是否大体均匀分布是否出现前沿中空的大缺口。前端点可以用scatter(fval(:,1), fval(:,2), filled)直接看。选解时如果拿不准挑哪个我一般看两个基准解成本最小的那个以及排放最小的那个。分别取[~, i_min] min(fval(:,1))和[~, i_min2] min(fval(:,2))再把这两个解的调度曲线画在时域图上。晚高峰时段柴油机出力差多少、储能 SOC 曲线差多少一眼就能看出权衡代价。这个对比步骤不需要额外调参但能快速检验算法有没有保证两个极端方向都覆盖到。本文还有配套的精品资源点击获取