做多微电网规划的人最头疼的往往不是电源选型而是拓扑设计。多个微电网之间通过联络线互联每条线建不建、容量怎么定、怎么保证整个系统既不出现孤岛又不出现环网这些决策的组合数量是爆炸级的。我最近把一个模拟项目X重新用Matlab完整实现了一遍核心算法是约束差分进化算法Constrain Differential EvolutionCDE优化变量直接写成大规模矩阵形式用来求解多微电网的最优拓扑结构。这篇文章就把建模、算法、代码实现和调试过程中的坑一次讲透适合刚接触微电网规划的研究生也适合想在工程上快速搭一版拓扑优化原型的工程师。先说一个总体结论这套方案不是我拍脑袋选的而是实测过几种算法之后觉得“省心”的组合。差分进化不需要求梯度天然能处理混合整数问题矩阵编码让连通性检查、辐射状约束和Matlab向量化计算变得非常顺手再加上基于可行性准则的约束处理机制整个代码的结构很清晰遇到问题也好排查。下面按我的实现顺序展开从问题建模到代码细节再到踩坑实录一步步来。1. 多微电网拓扑设计到底在优化什么1.1 拓扑设计不是一个画图问题而是一个NP难组合优化问题先把问题说清楚。假设系统里有N个微电网节点每个节点内部有分布式电源、储能和负荷。节点与节点之间可以新建联络线从而在某个节点发电富余时把功率送到缺电的节点减少弃电、提高可靠性。拓扑设计的本质就是回答两个问题哪些节点之间建线建线之后系统整体效益最好别小看这两个问题。N个节点之间最多有C(N,2)N(N-1)/2条候选线路每条线路只有“建”或“不建”两种状态整个解空间大小是2^(N(N-1)/2)。N10的时候候选线路最多45条组合数量大约是3.5万亿N20的时候这个数字已经大到无法枚举。更麻烦的是任何一个合法的拓扑还必须满足连通性、辐射状结构、线路容量和功率平衡等约束。综合来看这是一个典型的大规模组合优化问题用穷举法完全不现实这也是差分进化这类元启发式算法能发挥作用的前提。1.2 为什么用矩阵描述拓扑结构图论里多微电网系统天然就是一个图微电网节点是图的顶点联络线是图的边。为了在代码里方便表示和计算我选择用邻接矩阵X来编码拓扑结构。如果系统有N个节点X就是N×N的0-1对称矩阵。X(i,j)1表示节点i和节点j之间建设联络线X(i,j)0表示不建设。这种编码有几个非常实际的好处连通性检查很方便。Matlab内置的graph对象可以直接由邻接矩阵构造然后调用conncomp获取连通分量数量一行代码就能判断整个系统是否处于孤岛状态。辐射状约束可以用“边数N-1且连通”来判定这比在边列表上逐个检查要直观得多。差分进化算法里所有个体的长度都是固定的矩阵可以直接拉成向量参与变异和交叉维度对齐非常自然。我对比过边列表和邻接矩阵两种表示方式直观感受如下表对比维度边列表邻接矩阵存储效率高只存有效边随N的平方增长但可用sparse缓解图算法便利度一般需要额外维护顶点集合高直接用graph、conncomp等内置函数DE维度对齐需要手动映射拉直成向量非常直接可读性一般直观一眼看清节点间关系工程实现时我不会真的用一个N×N稠密矩阵去参与每代计算而是先维护候选线路表用一个M维0-1向量表示哪些候选线被选中在目标函数里再把这个向量组装成稀疏邻接矩阵。这样既能享受矩阵表示带来的便利又不会让内存爆炸。1.3 目标函数和约束条件的完整定义为了让算法能跑起来必须把“方案好坏”变成可计算的数学表达式。模拟项目X里我在目标函数里聚合了三部分成本得到一个年度综合费用年综合费用 线路投资等年值 网损费用 运行费用线路投资等年值可以写成C_inv Σ c_line(i,j) * X(i,j) * CRF其中c_line(i,j)是线路i-j的建设成本CRF是把总建设成本折算到每年的资本回收系数。网损费用采用直流潮流简化模型假设节点注入已知通过节点导纳矩阵求出线路功率再计算损耗。运行费用则粗略按各微电网内非清洁出力的燃料成本来估计。约束条件我拆成了四组拓扑连通性约束整个图必须连通不存在孤立节点或者孤立子微电网。辐射状结构约束边数必须等于N-1且连通等价于要求最终拓扑是一棵生成树。为什么强调这一点微电网群在运行时大多希望采用“闭环设计、开环运行”的模式环网会导致继电保护配合复杂、短路电流增大工程上通常不会让联络线随意成环。线路功率极限约束实际流过线路的功率不得超过该线路的最大传输容量如果线路没有被选中实际功率必须为0。功率平衡约束每个节点的注入功率内部出力减去负荷再减去流出功率之和必须等于0这是潮流计算的基本要求。可能有人会问直接把这些约束写进目标函数里当罚函数不就行了理论可以但工程上纯罚函数法在约束多且强的时候非常容易失效后文我会详细说我是怎么用约束DE的可行性准则来处理这些约束的。2. 约束差分进化算法为什么偏偏选它2.1 差分进化的核心三步变异、交叉、选择差分进化算法最早是为连续优化问题设计的但它的机制决定了它处理混合整数问题也很自然。标准DE的流程非常简单每个个体是一个D维实数向量种群规模设为NP迭代过程围绕三个算子展开。变异操作是DE的灵魂。当前种群中随机挑三个不同的个体 x_r1、x_r2、x_r3生成变异向量v_i x_r1 F * (x_r2 - x_r3)其中F是缩放因子一般在0.5到0.9之间。这个式子的含义就是让每个个体朝向种群中其他两个个体的差分方向偏移相当于用群体自身的多样性引导搜索方向。F越大搜索范围越广但收敛会变慢F太小种群容易早熟。交叉操作将变异向量和目标个体的分量按概率组合。DE常用二项式交叉对每个维度j生成一个[0,1]均匀随机数如果小于交叉概率CR就取变异向量的第j个分量否则取原个体的第j个分量。CR越大种群多样性越高对0-1问题来说CR通常可以设得比较大。选择操作采用一对一竞争。新个体u_i生成后计算它的目标函数值与旧个体x_i比较。如果是无约束问题直接保留目标函数值更优的个体即可。这个流程反复迭代直到达到最大代数或者收敛判据。2.2 约束处理机制Deb可行性准则有约束和无约束的DE最大差异就在选择这一步。如果我直接用“目标函数越低越好”的规则那些严重违反约束但目标函数极小的个体就会大量霸占种群最终搜出一个拓扑上完全不成立、跑不通潮流的方案。我采用的是Deb提出的三条可行性准则也叫约束支配准则。比较两个个体a和b时规则如下如果a和b都满足全部约束谁的目标函数值小谁胜出。如果a可行、b不可行无论b的目标函数多漂亮a都胜出。如果a和b都不可行谁的总约束违反量小谁胜出。这个规则特别直观翻译成大白话就是先要方案合法再谈方案省钱。它不会完全抛弃不可行个体而是给约束违反量小的个体留了生存空间让它们从不可行区域逐步修正到可行区域。在我的模拟算例里这一步对整个算法的稳定性提升非常明显后文会专门展开讲。2.3 与遗传算法、粒子群算法的横向对比既然问题已经明确是组合优化为什么不直接选遗传算法GA或者粒子群PSO我在同一套模拟算例下做过对比简单说说我的感受。算法主要优势主要劣势在拓扑设计中的表现遗传算法天然支持离散编码交叉变异算子设计空间大收敛速度偏慢能用但到后期多样性下降严重粒子群收敛快实现简单容易早熟需要处理速度到位移的离散映射小规模可行大规模容易陷进局部最优差分进化参数少变异机制保留多样性连续问题出身0-1离散需要额外处理与修复策略配合稳定收敛曲线平滑我没说DE绝对碾压其他算法但它在“变量维度高、离散连续混合、约束复杂”这个问题组合下代码实现成本是最低的。尤其它不需要算梯度目标函数里嵌两轮直流潮流计算也无所谓。这对工程原型验证来说很重要调起参来压力小很多。3. 大规模矩阵编码与初始化实操3.1 决策变量的矩阵编码方案先说清楚代码层面的数据结构。假设候选线路有M条网络数据里我会维护三个等长的数组s、t、lineCost分别记录每条候选线的起点、终点和建设成本。决策变量是一个M维的0-1向量xx(m)1表示第m条候选线被选中。返回拓扑矩阵时可以用一行Matlab代码把x组装成N×N的稀疏邻接矩阵X sparse(s(x1), t(x1), 1, N, N); X X X;注意这里用sparse而不是普通矩阵。N30的时候稠密邻接矩阵也就900个元素看着不大但DE每一代要对NP个个体做潮流计算如果代码里不小心把稠密矩阵传给graph或者参与矩阵运算运行时间会成倍增加。用sparse是个好习惯。为了让DE的变异交叉能在这个0-1离散空间上运作我采用“隐变量连续化”策略。真正参与变异交叉的变量v是M维实数向量每次评估个体时通过sigmoid函数映射到(0,1)区间再用0.5的阈值离散成0-1决策向量x double(sigmoid(v) 0.5);这么做的好处是DE那一套差分变异公式可以直接用在实数空间里不需要为离散问题额外设计二进制变异算子。代价是离散后的个体可能不满足拓扑约束所以下一步必须引入修复机制。3.2 初始种群的生成怎么保证每个解合法初始化这一步最容易翻车。如果初始种群完全随机生成0-1向量大量个体既不连通、又成环、边数还乱七八糟算法在可行域边缘反复挣扎收敛效率极低。我在模拟项目X中用的是“随机生成树初始化”策略。思路其实很简单对候选图随机赋一组权值然后用最小生成树算法生成一棵随机生成树把树上所有边对应的决策变量置为1其他边置为0。由于随机生成的边权不同每次调用得到的是不同结构的生成树正好满足种群初始多样性的需求。function x init_from_random_tree(s, t, N) M length(s); w rand(M, 1); % 随机边权 G graph(s, t, w, N); T minspantree(G); % 随机权值下的最小生成树 x zeros(M, 1); % minspantree 返回的T.Edges包含了生成树上的边 for e 1:height(T.Edges) idx find(s T.Edges.EndNodes(e,1) t T.Edges.EndNodes(e,2), 1); if isempty(idx) idx find(s T.Edges.EndNodes(e,2) t T.Edges.EndNodes(e,1), 1); end x(idx) 1; end end这个初始化方法能保证每个初始个体都是连通且无环的生成树后续所有约束的检查都建立在“基础拓扑合法”的前提下省了很多事。3.3 变异和交叉后的拓扑修复策略光有合法初始种群还不够。DE的变异和交叉会在实数空间里乱跳离散成0-1向量之后个体大概率会变成不连通、成环、边数不对的坏结构。所以在每次评估目标函数之前都必须先做拓扑修复把个体拉回到“近似生成树”的状态。修复流程我总结成三步用graph(X)和conncomp检查当前拓扑的连通分量个数。如果连通分量个数为1且边数大于N-1说明有环需要随机断环。Matlab里可以用minspantree或者遍历环上的边随机删除一条直到边数等于N-1。如果连通分量个数大于1说明存在孤岛。随机选择两个不同分量内的节点在两节点对应候选边存在的前提下把这条边加入拓扑直到整个图连通。function x repair_topology(x, s, t, N) M length(x); % 组装当前拓扑 X sparse(s(x1), t(x1), 1, N, N); G graph(X X); comps conncomp(G); % 先保证连通 while numel(unique(comps)) 1 % 选两个不同分量的节点 c1 find(comps 1, 1); c2 find(comps 2, 1); cand find((sc1 tc2) | (sc2 tc1)); if ~isempty(cand) x(cand(1)) 1; else break; end X sparse(s(x1), t(x1), 1, N, N); G graph(X X); comps conncomp(G); end % 再断环至边数N-1 while numel(unique(comps)) 1 nnz(x) N-1 X sparse(s(x1), t(x1), 1, N, N); G graph(X X); % 任意找一条环上边并删除可用findcycle简化 edges table2array(G.Edges); % 这里简化处理随机选一个当前选中的边尝试删除 sel find(x 1); rm sel(randi(numel(sel))); tmp_x x; tmp_x(rm) 0; Xt sparse(s(tmp_x1), t(tmp_x1), 1, N, N); Gt graph(Xt Xt); if numel(unique(conncomp(Gt))) 1 nnz(tmp_x) N-1 x(rm) 0; else % 换一条边再试 sel(sel rm) []; if isempty(sel), break; end end X sparse(s(x1), t(x1), 1, N, N); G graph(X X); comps conncomp(G); end end这段代码我不保证是所有情况下的最优写法但它在我跑的算例里表现稳定。需要提醒的是修复并不能解决线路容量和功率约束那部分要交给DE的选择机制处理。修复只保证拓扑结构上的“形状合法”。4. Matlab代码实现核心模块拆解4.1 主程序框架DE循环与收敛判断整个程序我按模块化思路拆了几个函数数据准备、种群初始化、目标函数评估、拓扑修复、变异交叉、约束支配比较。主程序框架如下控制DE主循环% 数据载入与参数设置 load(network_data_15node.mat); % 模拟算例数据 NP 60; F 0.6; CR 0.9; MaxGen 200; M length(network.s); % 初始化 pop zeros(NP, M); fit inf(NP, 1); cv inf(NP, 1); for i 1:NP pop(i, :) init_from_random_tree(network.s, network.t, network.N); [fit(i), cv(i)] evaluate_individual(pop(i, :), network); end % 迭代 best_fit zeros(MaxGen, 1); waitGen 0; for gen 1:MaxGen for i 1:NP % 变异 r randperm(NP, 3); v pop(r(1), :) F * (pop(r(2), :) - pop(r(3), :)); % 交叉 u pop(i, :); jrand randi(M); for j 1:M if rand() CR || j jrand u(j) v(j); end end % 离散化 修复 u_bin double(sigmoid(u) 0.5); u_bin repair_topology(u_bin, network.s, network.t, network.N); % 评估 [u_fit, u_cv] evaluate_individual(u_bin, network); % Deb准则选择 if is_dominated(u_fit, u_cv, fit(i), cv(i)) pop(i, :) u_bin; fit(i) u_fit; cv(i) u_cv; end end best_fit(gen) min(fit); % 连续停滞判断 if gen 10 abs(best_fit(gen) - best_fit(gen-10)) 1e-4 waitGen waitGen 1; else waitGen 0; end if waitGen 5 break; end end收敛判断我用了最简单的“连续几代最优值变化小于阈值”的规则。实际跑下来15节点算例在80代左右基本稳定提前终止能省不少时间。4.2 目标函数与约束罚函数的向量化写法目标函数是性能瓶颈。忽略向量化对每个个体里每条线路都写一层for循环NP60、M32、MaxGen200的情况下一轮仿真可能要跑好几分钟。我把潮流和目标函数尽量向量化之后同样配置能压到30秒内。直流潮流求解的核心代码如下function [P_line, P_loss] dc_power_flow(X, network) % X: 邻接矩阵(稀疏), 包含当前拓扑 % network.node_power: 各节点净注入(不含线路) N network.N; B network.B_sus; % 节点导纳矩阵虚部(稀疏) theta B \ network.node_power; % 直流潮流解相角 % 线路功率按相角差计算 [s, t] find(triu(X)); % 取出有效线路 P_line (theta(s) - theta(t)) ./ network.xline(sub2ind(...)); % 网损按 I^2*R 近似 P_loss P_line.^2 .* network.Rline / network.Vbase^2; end这里我需要说明直流潮流模型忽略了无功和电压幅值变化但对拓扑寻优的初筛阶段来说足够可靠。真要做详细运行方案评估可以后续再升级成交流潮流架构不变。约束违反度的计算也尽量一条表达式完成。线路容量约束的违反量定义为viol_cap sum(max(0, abs(P_line) - network.Pmax(x1)))功率平衡约束的违反量则直接取直流潮流方程残差因为节点功率不平衡主要体现在相角解算后的注入不匹配。把所有违反量求和就得到总约束违反量cv供Deb准则使用。4.3 仿真算例与结果分析为了验证算法我在模拟项目X中搭建了一个15节点的多微电网算例。每个节点内部包含光伏、风机、储能和常规负荷候选线路32条线路参数和建设成本都随机生成但符合实际数量级。DE参数设置为NP60、F0.6、CR0.9初始种群全部通过随机生成树生成。运行结果很有代表性。初始随机生成树的年度综合费用大约在2350万元左右算法迭代到第80代时降到约1950万元下降幅度约17%。收敛曲线前30代下降非常快之后趋于平缓整个搜索过程没有出现目标值大幅反弹的情况说明可行性准则加修复策略的组合确实让种群稳定在可行域附近。我还运行了同样条件下的标准遗传算法作为对照。GA在90代后目标值仍然比DE最优解高约6%而且最优个体对应的拓扑经常出现孤岛需要额外修复。这个结果当然不绝对但它支持了我在2.3节里的判断DE在这个问题上至少不输GA而且调参负担更小。5. 常见问题与工程调试经验5.1 约束处理太硬种群全是不可行解怎么办这是第一次跑通代码时最容易遇到的情况。具体表现是迭代很多代之后种群里所有个体仍然处于约束违反状态目标函数数值看起来很低但拓扑要么断成好几个小岛要么线路功率全线越限。我踩过这个坑之后的经验是先从三个方向排查。第一个方向是初始种群是否合法。如果初始解本身就不满足拓扑连通性后面怎么修都白搭。用随机生成树初始化可以从根源上解决。第二个方向是修复策略是否真的被执行。我在早期版本里把修复函数放在离散化之后、评估之前但断环那段代码存在死循环风险某些个体根本修不完。后来加了最大迭代保护修复失败就保留修复前状态交给Deb准则自然淘汰。第三个方向是约束违反度的量纲是否失衡。线路功率约束违反量可能是几十兆瓦级别而拓扑约束违反量是0或正整数两者加起来时拓扑约束被完全淹没算法只顾着降功率越限忽略了拓扑合法性。解决方法是先把各约束违反量做归一化再求和。5.2 DE参数敏感F和CR到底怎么调DE参数看起来只有三个但调起来很讲究。网络上很多文章推荐F0.5、CR0.9之类但那是针对标准连续函数测试集的经验放在多微电网拓扑问题上需要修正。我的实际经验是F不要设太小。0.5以下时种群在几十代内就会快速聚集多样性迅速下降结果容易停在某个局部最优拓扑上。F设在0.6到0.8之间搜索过程更稳妥。CR可以适当大一些0.8到0.95都可以因为交叉率越高变异向量的信息保留得越多离散化后个体结构变化越大有利于跳出局部最优。NP的设置更微妙。理论上NP越大越好但在目标函数里要做直流潮流计算NP翻倍直接意味着时间翻倍。我建议NP先取50到80之间跑一轮看收敛曲线。如果曲线在近100代还在持续下降说明搜索能力不足再加大NP或者提高F。我常用的一个偷懒策略是F自适应。前30代用F0.8后面降到0.5初期保持多样性后期加强精细搜索。实测比固定参数稳定不少。5.3 大规模矩阵内存和速度优化当N超过20、候选线路超过100条时优化规模会明显变大。这个规模在Matlab里还不至于内存爆炸但要是不注意代码写法速度会让人崩溃。几个我实测有效的优化习惯所有邻接矩阵都用sparse存储。Matlab对sparse矩阵的图算法支持很完善完全不用为了可读性而逃回稠密矩阵。候选线路预筛。地理距离超过安全送电距离的边直接删掉既缩减变量维度M又避免算法在明显不合理的候选线上浪费时间。目标函数的一次性批量计算。评估NP个个体时不要一个一个在循环里调用潮流求解器而是提前把所有个体的邻接矩阵组装成元胞数组或三维数组在向量化层面统一处理。实在不行就用parfor做并行评估注意每次并行迭代时重置随机种子即可否则所有worker生成的随机数可能完全相同DE的差分变异就失去了意义。缓存不变的网络参数。节点导纳矩阵、线路长度、线路成本这些在迭代中不变的数据全部提前放进struct函数调用时直接传递不要在目标函数里重复计算。最后再分享一点实操体会整套代码跑通之后我最大的感悟是拓扑优化的瓶颈往往不在算法本身而在约束建模是否贴合工程实际。我第一版代码忽略了辐射状约束算法搜出来的拓扑经常带着环网虽然目标函数好看但拿来直接做继电保护整定完全没法用。后来把“生成树初始化修复Deb准则”这三层结构加进去问题立刻少了一大半。这个框架后续扩展空间很大。比如把目标函数里的直流潮流升级成交流潮流或者改成双层优化上层搜拓扑、下层嵌套一个运行调度模型把储能充放电策略也考虑进去。现有代码的模块化设计基本不用大改只要替换目标函数评估模块就行。我接下来就准备往这个方向做一版到时候有结果了再继续分享。