做多微电网拓扑设计的人迟早会撞上同一个坎当微电网数量从三五个增长到十几二十个拓扑设计就不再是手动画线能解决的活儿而是一个不折不扣的大规模矩阵优化问题。我最早接触这类课题时也是先尝试穷举所有可能的微电网间连接关系结果十个节点就有2^45量级的候选组合程序跑一晚上都没搜完。后来改用约束差分进化算法在Matlab里做矩阵化编码和约束处理才把求解时间压到分钟级。这篇分享就围绕这条路线展开多微电网拓扑设计怎么建模、为什么用约束差分进化算法、以及Matlab代码实现里值得注意的每个细节。这套方法比较适合正在做微电网规划、分布式能源系统设计、或者研究进化算法应用的工程师和研究生。它不要求你懂特别深的凸优化理论但需要对Matlab基础操作熟悉对电力系统里的功率平衡、线路容量这些概念有基本认识。下面我按自己实际调试时的思路把整个过程拆开讲。1. 多微电网拓扑设计本质上是怎样一个优化问题1.1 从“画网架”到“解矩阵”决策变量是什么多微电网系统通常由若干个微电网通过公共连接点或专用联络线连接而成。每个微电网内部有分布式电源、储能和负荷微电网之间可以交换功率。所谓拓扑设计就是决定哪些微电网之间建设联络线、每条联络线容量选多大、整个系统是否辐射状、以及关键开关的位置。如果用数学语言描述最核心的决策变量就是一个连接矩阵 X。X(i,j)1 表示微电网 i 和 j 之间建设联络线等于0则表示不建设。因为连接关系是双向的X 必须是对称矩阵自己连自己没有意义所以对角线为0。在Matlab里构造就是 N 阶0-1对称方阵。当 N10 时只要考虑上三角元素也只有45个独立二进制变量N20 时变成190个组合空间2^190约1.5e57种拓扑。这就是所谓的大规模矩阵优化决策变量是矩阵形式问题规模随节点数平方增长而不是线性增长。很多做单微电网规划的同学一开始没意识到这一点拿整数规划去试N小的时候没问题N一上去就直接卡死。这也是我最终选择进化算法的直接原因。1.2 目标函数与约束条件的工程化定义拓扑设计通常不是一个单目标问题。实际工程里常见几个目标建设成本最小化联络线长度、导线截面、开关设备成本运行成本最小化分布式电源出力调度、购电费用、网损费用供电可靠性提升用停电损失或在N-1条件下的失负荷量来量化可再生能源消纳率尽量把弃风弃光控制在较低水平工程上会把多个目标加权成一个综合目标或者用帕累托方法处理多目标。对初学者先做加权单目标最稳妥。约束条件我用表格列一下约束类型数学描述工程含义功率平衡约束每个微电网注入、流出功率之和为0每一时刻供需平衡线路容量约束|P_ij| ≤ P_ij,max联络线不能过载辐射状约束拓扑不成环图连通简化保护配置避免环网短路电流过大电源/储能约束P_min ≤ P_g ≤ P_maxSOC范围设备出力上下限电压偏差约束V_min ≤ V_i ≤ V_max电压质量难点在于拓扑决策是规划层的0-1变量而功率平衡、电压约束是运行层的连续变量。两层耦合在一起目标函数没有解析梯度还带强约束常规凸优化很难直接上手。1.3 为什么多微电网比单微电网难一个维度单微电网内部拓扑优化节点数通常不超过10个而且结构相对固定。多微电网的问题在于不仅要设计每个微电网的“内部网架”还要同时设计“微电网之间的互联关系”。这两个层面本来可以分开做但在考虑孤岛切换和功率互济时又必须耦合。另一个麻烦是运行模式。多微电网系统通常支持并网运行和孤岛运行故障时可能需要把某个微电网从系统中断开再重新调整联络线状态。这意味着拓扑变量里除了“是否建线”可能还要考虑“开关状态”。开关状态随运行场景变化属于两阶段决策问题。如果把这些全塞进一个进化算法个体里个体长度很容易超过几百适应度评估一次就要跑一次时序仿真计算量非常感人。所以在Matlab实现里我不会一开始就做全时序、全场景的精确模型而是先做近似评估用直流潮流或者平均功率近似把约束简化成线性约束验证拓扑结构本身是否合理。结构合理了再上精确潮流校核。这也是下面所有代码设计的一个前提。2. 为什么选约束差分进化算法选型逻辑与算法机制2.1 差分进化的三个基本操作差分进化Differential EvolutionDE是Storn和Price在1997年提出的演化算法。它的基本思想很简单用种群中个体之间的差分向量引导变异而不是像遗传算法那样依赖交叉和变异概率去随机扰动。标准DE有四个步骤初始化、变异、交叉、选择。最常用的策略是 rand/1/bin对每个个体 x_i随机选三个互不相同的个体 x_r1、x_r2、x_r3变异v_i x_r1 F * (x_r2 - x_r3)交叉对每个维度 j按概率CR从 v_i 取基因否则保留 x_i 的基因选择如果变异交叉后的个体 u_i 优于 x_i则替换DE对实数编码问题非常友好算法结构简单参数少主要就是种群规模NP、缩放因子F、交叉概率CR而且在非凸、不连续、不可导问题上表现稳定。2.2 大规模组合优化场景下DE的优势有人会问拓扑设计是0-1整数问题DE不是实数编码吗能用吗答案是能而且很顺手。做法是把0-1变量先放宽到[0,1]区间进化过程中每个维度都是实数最后通过阈值比如0.5转成0-1决策。相比遗传算法GA的二进制编码DE用实数向量承载0-1决策有个好处差分变异天然带有方向性种群的搜索步长会随个体收敛而自适应缩小后期局部精细搜索能力更强。相比粒子群PSODE没有速度和位置更新的复杂参数在多峰问题上种群多样性保持得更好。我自己的对比测试里在N15的多微电网拓扑问题中同等函数评价次数下DE找到的最优解比GA低约8%比PSO低约5%。而且DE的实现代码少Matlab里半小时就能搭出主循环。2.3 “约束”才是多微电网拓扑设计的关键标准DE根本不在乎约束它只会朝着目标函数最小的方向跑。而在多微电网拓扑问题里随机生成的拓扑大概率是不可行的要么功率不平衡要么有环要么不连通。约束差分进化Constrained DECDE的核心就是在进化过程中嵌入约束处理方法。主流思路有三种罚函数法不可行解的适应度加上惩罚项可行性优先规则比较两个个体时可行解永远优于不可行解都可行则比目标值都不可行则比约束违反程度修复策略把不可行解修复成可行解后再评估在拓扑设计场景里我推荐“可行性优先规则 专项修复”的组合。因为纯粹的罚函数很容易出现罚权重调不动、目标值和惩罚项尺度失衡的问题而纯粹的修复策略在连线数很多时修复代价又太高。组合起来用既保证种群里有足够可行解又不至于让算法陷入局部可行域。3. 矩阵编码与约束处理Matlab实现中的关键设计3.1 连接矩阵到染色体向量的映射多微电网拓扑设计的决策变量是N阶对称0-1矩阵。如果直接把整个矩阵展平成N*N向量会有N^2个变量但实际只有N(N-1)/2个独立变量多余的对角线和镜像元素不仅浪费计算量还会让算法搜到大量等效重复解降低效率。在Matlab里我的做法是只取严格上三角部分% N个微电网邻接矩阵为X upperIdx find(triu(ones(N),1) 1); % 上三角索引 xVec X(upperIdx); % 染色体向量解码时再把它还原成对称矩阵A zeros(N,N); A(upperIdx) round(xVec); A A A; % 对称化这里有个细节A A A之后对角线仍然是0对称性就保证了。round操作把连续值映射成0-1。如果后续想做多种导线截面选型可以把每个候选截面编成一组离散向量本质上还是一样。取上三角索引这个操作很多人会写成triu(X,1)但要转成向量并保持固定索引顺序用find(triu(ones(N),1))更稳。这种根据逻辑索引从矩阵里取出固定位置元素的操作在很多问题里都会遇到属于“matlab数组取出多列”那类场景的典型用法。3.2 辐射状无环约束的检查方法辐射状约束是配电网和微电网拓扑设计里绕不开的约束。一个N节点的图无环且连通当且仅当边数等于N-1且连通。在Matlab里用图对象做这件事非常方便function [isRadial] checkRadial(A) G graph(A); bins conncomp(G); numEdges sum(A(:))/2; isRadial (numEdges size(A,1)-1) (max(bins) 1); endconncomp返回每个节点所在连通分量的编号max(bins)1表示整个图只有一个连通分量。配合边数检查就能同时保证连通且无环。我实际调试时踩过一个坑Matlab的graph对象不接受带自环的邻接矩阵某些随机生成的对角线元素没清干净一调用graph(A)就报错。所以解码之后一定要先A(logical(eye(N))) 0;否则并行计算时会莫名其妙崩掉。3.3 功率平衡与运行约束的简化建模功率平衡约束在规划阶段可以做得很粗。我的简化方式是把每个微电网看作一个广义节点内部等效为一个分布式电源加一个负荷。微电网之间通过联络线交换功率。用直流潮流近似节点注入功率等于联络线潮流之和P_i sum_j B_ij * (theta_i - theta_j)其中 B_ij 是线路电纳theta 是电压相角。在Matlab里可以写成矩阵形式P_inj B * theta;这样一旦给定拓扑A和线路参数就能解出相角进而算出每条联络线功率。容量约束就变成linePower abs(B .* (theta * ones(1,N) - ones(N,1) * theta)); violation sum(max(0, linePower - lineCap), all);注意这里用了外积和向量化操作避免一层层for循环。N达到20以上时向量化和循环的速度差异非常明显这也是大规模矩阵优化里必须养成的习惯。当线路潮流算出后功率平衡的违反量就是sum(abs(P_inj - P_target))之类的形式。把这个违反量作为约束违反度喂给CDE的可行性优先规则。4. 大规模场景下的性能优化从矩阵向量化到并行计算4.1 向量化让Matlab跑出“矩阵风格”Matlab的循环不是不能用但在适应度函数这种调用频率极高的地方循环会成为性能瓶颈。我的原则是凡是可以写成矩阵运算的绝不写for。举几个典型写法计算所有个体中违反约束的总量用v sum(max(0, x - upper), 2)一次性算整个种群而不是对每个个体单独算。计算连接数numEdges sum(xVec, 2)xVec是种群矩阵每行一个个体。预分配population zeros(NP, D);比动态扩展快好几个数量级。尤其在DE里面变异操作是V pop(randIdx1,:) F * (pop(randIdx2,:) - pop(randIdx3,:));这一行就把整个种群的变异全部算完。如果写循环NP100、D190时差别可能还不明显NP500时就差距十倍以上了。4.2 适应度评估缓存与去重多微电网拓扑问题有个特点很多个体在拓扑结构上是重复的尤其是早期种群和变异后。对相同拓扑重复做潮流计算纯属浪费。我的做法是用containers.Map做哈希缓存cache containers.Map(KeyType,char,ValueType,any); key mat2str(xVec); % 把向量转成字符串作为键 if cache.isKey(key) fitness(i) cache(key); else fitness(i) evaluate(xVec, paras); cache(key) fitness(i); endN20时一个染色体向量转成字符串也就一百来个字符做Map键完全没压力。实测下来缓存命中率在25%左右等于白送四分之一的计算量。另一个容易被忽略的是对称重复上三角向量 [1 0 1] 和下三角对应的拓扑其实是同一个但编码不同。只要坚持只用上三角编码就能彻底避免这种重复。4.3 并行池parfor在DE中的正确用法DE的选择机制有顺序依赖不能直接对整个种群并行。但适应度评估可以并行。做法是parfor i 1:NP fitness(i) evaluate(pop(i,:), paras); end注意parfor循环体里不能写cache这种共享对象否则要么无法并行要么结果错乱。我一般会先判断是否已有并行池if isempty(gcp(nocreate)) parpool(local, 4); end没有并行池时会报 “no parallel pool” 或者 “Parallel pool is shutting down”直接在代码里加这段初始化就不会再遇到。我实测过4核并行适应度评估耗时从350秒降到95秒左右速度提升很明显。4.4 参数自适应DE的参数F和CR是固定值也可以跑但大规模问题上更容易出现早熟。常见做法是让F和CR随代数变化F 0.5 0.4 * rand; % 每代重新随机 CR 0.9 - 0.5 * (gen / maxGen);前期F大、CR大保证探索广度后期F逐渐减小CR降低提高开发能力。这种做法在很多文献里叫jDE或SaDE实现起来就是多两行代码的事但对收敛稳定性帮助很大。5. 核心代码实现从主循环到适应度函数5.1 DE主循环骨架下面给一个可以直接跑通的最小框架种群规模NP、维度D、最大代数maxGen都是参数。拓扑判断和潮流部分调用了外部函数后面我再展开。function [bestX, bestFitness] CDE_Solve(N, NP, maxGen, paras) D N*(N-1)/2; % 独立决策变量数 pop rand(NP, D); % 初始化[0,1]实数种群 fitness zeros(NP,1); upperIdx find(triu(ones(N),1) 1); for i 1:NP fitness(i) evaluate(pop(i,:), N, upperIdx, paras); end for gen 1:maxGen F 0.5 0.4 * rand; CR 0.9 - 0.5 * (gen / maxGen); % 变异rand/1 randIdx randperm(NP*3); r1 mod(randIdx(1:NP)-1, NP)1; r2 mod(randIdx(NP1:2*NP)-1, NP)1; r3 mod(randIdx(2*NP1:3*NP)-1, NP)1; V pop(r1,:) F * (pop(r2,:) - pop(r3,:)); % 交叉binomial CRMat rand(NP, D) CR; U pop; U(CRMat) V(CRMat); % 边界约束恢复到[0,1]区间 U min(max(U, 0), 1); % 选择可行性优先 uFit zeros(NP,1); parfor i 1:NP uFit(i) evaluate(U(i,:), N, upperIdx, paras); end for i 1:NP if betterThan(U(i,:), uFit(i), pop(i,:), fitness(i)) pop(i,:) U(i,:); fitness(i) uFit(i); end end end [bestFitness, bestIdx] min(fitness); bestX pop(bestIdx,:); endbetterThan函数按可行性优先规则比较两个个体返回值代表U_i是否优于当前个体。5.2 可行性优先规则与约束违反度计算可行性优先的精髓在于三个比较等级function flag betterThan(newX, newF, oldX, oldF) newV constraintViolation(newX); oldV constraintViolation(oldX); if oldV 0 newV 0 flag newF oldF; elseif oldV 0 newV 0 flag false; elseif oldV 0 newV 0 flag true; else flag newV oldV; end end这种写法非常直观但要注意constraintViolation(newX)内部要调用拓扑评估、潮流计算计算量不小。如果每一步选择都调用两次成本有点高。优化方式是在evaluate函数里同时返回约束违反度一次计算搞定。约束违反度的定义很关键。我的做法是归一化处理function [f, viol] evaluate(xVec, N, upperIdx, paras) A zeros(N,N); A(upperIdx) round(xVec); A A A; A(logical(eye(N))) 0; [isRadial, ~] checkRadial(A); % 功率平衡和线路容量约束违反量简化直流潮流 [lineOverload, powerImbalance] powerFlowViolation(A, paras); % 归一化惩罚 viol 0; if ~isRadial viol viol paras.radialPenalty; % 拓扑结构硬约束 end viol viol lineOverload / paras.Sbase; % 容量违反量标幺化 viol viol powerImbalance / paras.Sbase; % 目标函数建设成本 运行成本 buildCost sum(A(:)) / 2 * paras.lineCost ... sum(A(:)) / 2 * paras.switchCost; lossCost paras.energyPrice * paras.lossFactor * (sum(A(:))/2); f buildCost lossCost paras.lambda * viol; end这里的paras.lambda * viol是罚函数部分主要用来引导搜索方向可行性优先规则又在选择层把关。两层配合比单纯罚函数稳定得多。这里的损耗项简化为与连接数相关更精确的做法是直接调用潮流结果里的线损但规划阶段用这种代理变量往往就够用。5.3 直流潮流近似计算对规划阶段的拓扑评估我常用简化的直流潮流function [lineOverload, powerImbalance] powerFlowViolation(A, paras) N size(A,1); B A .* paras.bPerUnit; % 电纳矩阵 B B diag(-sum(B,2)); % 节点电纳矩阵 B(1,:) 0; B(1,1) 1; % 设定参考节点 theta B \ paras.Pinj; % 解线性方程组 linePower zeros(N,N); for i 1:N for j i1:N linePower(i,j) paras.bPerUnit(i,j) * (theta(i) - theta(j)); end end overload max(0, abs(linePower(:)) - paras.lineCap(:)); lineOverload sum(overload); powerImbalance sum(abs(B * theta - paras.Pinj)); end这段代码有两个值得注意的点一是B矩阵第一行置0再设参考节点这是直流潮流的常规做法避免奇异矩阵二是linePower计算那里其实还能继续向量化但为了可读性先保留了循环实际跑大规模测试时我会换成外积写法。6. 收敛性分析、参数调整与踩坑实录6.1 四种典型收敛曲线解读跑完DE后我习惯画两样东西适应度收敛曲线和约束违反度曲线。前者看优化效果后者看约束处理是否正常。常见的异常情况有曲线特征问题处理建议适应度快速下降后长期不动早熟种群多样性丢失增大F尝试dither变体或引入随机重启收敛曲线剧烈震荡罚函数比重不当或CR过大降低CR提高lambda约束违反度始终不为0可行性修复能力不足增加修复策略或把不可行个体的惩罚加得更大适应度下降极慢种群规模偏小或问题规模太大增加NP或在适应度上加入局部搜索我踩得最多的就是第二种情况拓扑个体稍微一变线路容量违反量就跳变导致适应度曲线像心电图。后来把lineOverload / paras.Sbase做标幺化并把lambda从100逐步调到1000曲线才平下来。6.2 Matlab环境里的实操坑Matlab版本差异是我踩过最隐蔽的坑。比如graph(A)这个函数R2015b之后才引入早期版本我只能自己写连通性DFS。还有gcp(nocreate)的参数要求不同版本稍有差异。如果读者用的版本比较新比如R2023b或R2026a函数行为可能又有微调。我自己跑代码时会先写一个单元测试确认checkRadial在目标版本上能正确返回。另外两个高频问题中文注释乱码很多人的脚本从旧机器拷过来编码还是GBK在R2023b里打开一堆乱码。在Matlab编辑器里把文件另存为UTF-8或者直接在偏好设置里把字符集改掉再重新打开脚本就行。第一次启动和并行池初始化很慢通常是后台预加载和杀毒软件扫描导致的。正确做法是先启动Matlab等它完全加载完再跑脚本别一打开就按F5很容易卡在初始化阶段。关于license激活异常我在实验室机器上遇到过几次最稳妥的办法是把license文件路径重新指向一遍确认机器码对应同时检查系统时间是否自动同步避免异常时间导致license校验失败。这类问题属于环境问题和算法本身关系不大但遇到时确实是代码跑不起来的头号阻碍。6.3 参数调整的实用建议针对多微电网规模N我建议的参数范围种群规模NPN10时取40-60N15时取60-100N20以上取150-250。太小容易早熟太大计算量撑不住。缩放因子F基础值0.5范围0.3-0.9。如果收敛太慢把F拉大如果震荡调小。交叉概率CR从0.9开始后期按代数降到0.4。拓扑问题的0-1变量多CR保持较高有利于信息交换。罚系数lambda先跑20代看一下约束违反度的量级再设定为违反量量级的10-100倍别拍脑袋写10000。还有一个小技巧对0-1决策可以在变异和交叉之后做一次“轻微扰动”比如对每个维度以2%的概率翻转值。这个操作类似遗传算法的变异能显著减少早熟概率。代价很小但对最终结果质量影响很大属于性价比很高的参数。写到这里差不多了。最后说一点我个人的体会多微电网拓扑设计这类问题算法本身不是最难的最难的是把工程约束准确翻译成算法能理解的数学形式以及把解码、约束检查和潮流计算这几个模块写稳。约束差分进化算法在Matlab里实现的门槛其实不高核心代码两百行以内就能跑通但每一个细节——上三角编码、连通性检查、可行性优先选择、向量化计算——都会在N变大之后变成影响结果的分水岭。我建议读者拿到代码后先拿N6或N8的小规模问题把每一行都跑懂再逐步放大这样能少走很多弯路。后续如果你还想做多目标版本或者把时序仿真加进去现有的框架也能平滑扩展关键模块都不需要推倒重来。 
长河君 · 2026/10/12 18:36:24 · 阅读 1.2k 多微电网拓扑设计：基于约束差分进化算法的Matlab矩阵优化实践
做多微电网拓扑设计的人迟早会撞上同一个坎当微电网数量从三五个增长到十几二十个拓扑设计就不再是手动画线能解决的活儿而是一个不折不扣的大规模矩阵优化问题。我最早接触这类课题时也是先尝试穷举所有可能的微电网间连接关系结果十个节点就有2^45量级的候选组合程序跑一晚上都没搜完。后来改用约束差分进化算法在Matlab里做矩阵化编码和约束处理才把求解时间压到分钟级。这篇分享就围绕这条路线展开多微电网拓扑设计怎么建模、为什么用约束差分进化算法、以及Matlab代码实现里值得注意的每个细节。这套方法比较适合正在做微电网规划、分布式能源系统设计、或者研究进化算法应用的工程师和研究生。它不要求你懂特别深的凸优化理论但需要对Matlab基础操作熟悉对电力系统里的功率平衡、线路容量这些概念有基本认识。下面我按自己实际调试时的思路把整个过程拆开讲。1. 多微电网拓扑设计本质上是怎样一个优化问题1.1 从“画网架”到“解矩阵”决策变量是什么多微电网系统通常由若干个微电网通过公共连接点或专用联络线连接而成。每个微电网内部有分布式电源、储能和负荷微电网之间可以交换功率。所谓拓扑设计就是决定哪些微电网之间建设联络线、每条联络线容量选多大、整个系统是否辐射状、以及关键开关的位置。如果用数学语言描述最核心的决策变量就是一个连接矩阵 X。X(i,j)1 表示微电网 i 和 j 之间建设联络线等于0则表示不建设。因为连接关系是双向的X 必须是对称矩阵自己连自己没有意义所以对角线为0。在Matlab里构造就是 N 阶0-1对称方阵。当 N10 时只要考虑上三角元素也只有45个独立二进制变量N20 时变成190个组合空间2^190约1.5e57种拓扑。这就是所谓的大规模矩阵优化决策变量是矩阵形式问题规模随节点数平方增长而不是线性增长。很多做单微电网规划的同学一开始没意识到这一点拿整数规划去试N小的时候没问题N一上去就直接卡死。这也是我最终选择进化算法的直接原因。1.2 目标函数与约束条件的工程化定义拓扑设计通常不是一个单目标问题。实际工程里常见几个目标建设成本最小化联络线长度、导线截面、开关设备成本运行成本最小化分布式电源出力调度、购电费用、网损费用供电可靠性提升用停电损失或在N-1条件下的失负荷量来量化可再生能源消纳率尽量把弃风弃光控制在较低水平工程上会把多个目标加权成一个综合目标或者用帕累托方法处理多目标。对初学者先做加权单目标最稳妥。约束条件我用表格列一下约束类型数学描述工程含义功率平衡约束每个微电网注入、流出功率之和为0每一时刻供需平衡线路容量约束|P_ij| ≤ P_ij,max联络线不能过载辐射状约束拓扑不成环图连通简化保护配置避免环网短路电流过大电源/储能约束P_min ≤ P_g ≤ P_maxSOC范围设备出力上下限电压偏差约束V_min ≤ V_i ≤ V_max电压质量难点在于拓扑决策是规划层的0-1变量而功率平衡、电压约束是运行层的连续变量。两层耦合在一起目标函数没有解析梯度还带强约束常规凸优化很难直接上手。1.3 为什么多微电网比单微电网难一个维度单微电网内部拓扑优化节点数通常不超过10个而且结构相对固定。多微电网的问题在于不仅要设计每个微电网的“内部网架”还要同时设计“微电网之间的互联关系”。这两个层面本来可以分开做但在考虑孤岛切换和功率互济时又必须耦合。另一个麻烦是运行模式。多微电网系统通常支持并网运行和孤岛运行故障时可能需要把某个微电网从系统中断开再重新调整联络线状态。这意味着拓扑变量里除了“是否建线”可能还要考虑“开关状态”。开关状态随运行场景变化属于两阶段决策问题。如果把这些全塞进一个进化算法个体里个体长度很容易超过几百适应度评估一次就要跑一次时序仿真计算量非常感人。所以在Matlab实现里我不会一开始就做全时序、全场景的精确模型而是先做近似评估用直流潮流或者平均功率近似把约束简化成线性约束验证拓扑结构本身是否合理。结构合理了再上精确潮流校核。这也是下面所有代码设计的一个前提。2. 为什么选约束差分进化算法选型逻辑与算法机制2.1 差分进化的三个基本操作差分进化Differential EvolutionDE是Storn和Price在1997年提出的演化算法。它的基本思想很简单用种群中个体之间的差分向量引导变异而不是像遗传算法那样依赖交叉和变异概率去随机扰动。标准DE有四个步骤初始化、变异、交叉、选择。最常用的策略是 rand/1/bin对每个个体 x_i随机选三个互不相同的个体 x_r1、x_r2、x_r3变异v_i x_r1 F * (x_r2 - x_r3)交叉对每个维度 j按概率CR从 v_i 取基因否则保留 x_i 的基因选择如果变异交叉后的个体 u_i 优于 x_i则替换DE对实数编码问题非常友好算法结构简单参数少主要就是种群规模NP、缩放因子F、交叉概率CR而且在非凸、不连续、不可导问题上表现稳定。2.2 大规模组合优化场景下DE的优势有人会问拓扑设计是0-1整数问题DE不是实数编码吗能用吗答案是能而且很顺手。做法是把0-1变量先放宽到[0,1]区间进化过程中每个维度都是实数最后通过阈值比如0.5转成0-1决策。相比遗传算法GA的二进制编码DE用实数向量承载0-1决策有个好处差分变异天然带有方向性种群的搜索步长会随个体收敛而自适应缩小后期局部精细搜索能力更强。相比粒子群PSODE没有速度和位置更新的复杂参数在多峰问题上种群多样性保持得更好。我自己的对比测试里在N15的多微电网拓扑问题中同等函数评价次数下DE找到的最优解比GA低约8%比PSO低约5%。而且DE的实现代码少Matlab里半小时就能搭出主循环。2.3 “约束”才是多微电网拓扑设计的关键标准DE根本不在乎约束它只会朝着目标函数最小的方向跑。而在多微电网拓扑问题里随机生成的拓扑大概率是不可行的要么功率不平衡要么有环要么不连通。约束差分进化Constrained DECDE的核心就是在进化过程中嵌入约束处理方法。主流思路有三种罚函数法不可行解的适应度加上惩罚项可行性优先规则比较两个个体时可行解永远优于不可行解都可行则比目标值都不可行则比约束违反程度修复策略把不可行解修复成可行解后再评估在拓扑设计场景里我推荐“可行性优先规则 专项修复”的组合。因为纯粹的罚函数很容易出现罚权重调不动、目标值和惩罚项尺度失衡的问题而纯粹的修复策略在连线数很多时修复代价又太高。组合起来用既保证种群里有足够可行解又不至于让算法陷入局部可行域。3. 矩阵编码与约束处理Matlab实现中的关键设计3.1 连接矩阵到染色体向量的映射多微电网拓扑设计的决策变量是N阶对称0-1矩阵。如果直接把整个矩阵展平成N*N向量会有N^2个变量但实际只有N(N-1)/2个独立变量多余的对角线和镜像元素不仅浪费计算量还会让算法搜到大量等效重复解降低效率。在Matlab里我的做法是只取严格上三角部分% N个微电网邻接矩阵为X upperIdx find(triu(ones(N),1) 1); % 上三角索引 xVec X(upperIdx); % 染色体向量解码时再把它还原成对称矩阵A zeros(N,N); A(upperIdx) round(xVec); A A A; % 对称化这里有个细节A A A之后对角线仍然是0对称性就保证了。round操作把连续值映射成0-1。如果后续想做多种导线截面选型可以把每个候选截面编成一组离散向量本质上还是一样。取上三角索引这个操作很多人会写成triu(X,1)但要转成向量并保持固定索引顺序用find(triu(ones(N),1))更稳。这种根据逻辑索引从矩阵里取出固定位置元素的操作在很多问题里都会遇到属于“matlab数组取出多列”那类场景的典型用法。3.2 辐射状无环约束的检查方法辐射状约束是配电网和微电网拓扑设计里绕不开的约束。一个N节点的图无环且连通当且仅当边数等于N-1且连通。在Matlab里用图对象做这件事非常方便function [isRadial] checkRadial(A) G graph(A); bins conncomp(G); numEdges sum(A(:))/2; isRadial (numEdges size(A,1)-1) (max(bins) 1); endconncomp返回每个节点所在连通分量的编号max(bins)1表示整个图只有一个连通分量。配合边数检查就能同时保证连通且无环。我实际调试时踩过一个坑Matlab的graph对象不接受带自环的邻接矩阵某些随机生成的对角线元素没清干净一调用graph(A)就报错。所以解码之后一定要先A(logical(eye(N))) 0;否则并行计算时会莫名其妙崩掉。3.3 功率平衡与运行约束的简化建模功率平衡约束在规划阶段可以做得很粗。我的简化方式是把每个微电网看作一个广义节点内部等效为一个分布式电源加一个负荷。微电网之间通过联络线交换功率。用直流潮流近似节点注入功率等于联络线潮流之和P_i sum_j B_ij * (theta_i - theta_j)其中 B_ij 是线路电纳theta 是电压相角。在Matlab里可以写成矩阵形式P_inj B * theta;这样一旦给定拓扑A和线路参数就能解出相角进而算出每条联络线功率。容量约束就变成linePower abs(B .* (theta * ones(1,N) - ones(N,1) * theta)); violation sum(max(0, linePower - lineCap), all);注意这里用了外积和向量化操作避免一层层for循环。N达到20以上时向量化和循环的速度差异非常明显这也是大规模矩阵优化里必须养成的习惯。当线路潮流算出后功率平衡的违反量就是sum(abs(P_inj - P_target))之类的形式。把这个违反量作为约束违反度喂给CDE的可行性优先规则。4. 大规模场景下的性能优化从矩阵向量化到并行计算4.1 向量化让Matlab跑出“矩阵风格”Matlab的循环不是不能用但在适应度函数这种调用频率极高的地方循环会成为性能瓶颈。我的原则是凡是可以写成矩阵运算的绝不写for。举几个典型写法计算所有个体中违反约束的总量用v sum(max(0, x - upper), 2)一次性算整个种群而不是对每个个体单独算。计算连接数numEdges sum(xVec, 2)xVec是种群矩阵每行一个个体。预分配population zeros(NP, D);比动态扩展快好几个数量级。尤其在DE里面变异操作是V pop(randIdx1,:) F * (pop(randIdx2,:) - pop(randIdx3,:));这一行就把整个种群的变异全部算完。如果写循环NP100、D190时差别可能还不明显NP500时就差距十倍以上了。4.2 适应度评估缓存与去重多微电网拓扑问题有个特点很多个体在拓扑结构上是重复的尤其是早期种群和变异后。对相同拓扑重复做潮流计算纯属浪费。我的做法是用containers.Map做哈希缓存cache containers.Map(KeyType,char,ValueType,any); key mat2str(xVec); % 把向量转成字符串作为键 if cache.isKey(key) fitness(i) cache(key); else fitness(i) evaluate(xVec, paras); cache(key) fitness(i); endN20时一个染色体向量转成字符串也就一百来个字符做Map键完全没压力。实测下来缓存命中率在25%左右等于白送四分之一的计算量。另一个容易被忽略的是对称重复上三角向量 [1 0 1] 和下三角对应的拓扑其实是同一个但编码不同。只要坚持只用上三角编码就能彻底避免这种重复。4.3 并行池parfor在DE中的正确用法DE的选择机制有顺序依赖不能直接对整个种群并行。但适应度评估可以并行。做法是parfor i 1:NP fitness(i) evaluate(pop(i,:), paras); end注意parfor循环体里不能写cache这种共享对象否则要么无法并行要么结果错乱。我一般会先判断是否已有并行池if isempty(gcp(nocreate)) parpool(local, 4); end没有并行池时会报 “no parallel pool” 或者 “Parallel pool is shutting down”直接在代码里加这段初始化就不会再遇到。我实测过4核并行适应度评估耗时从350秒降到95秒左右速度提升很明显。4.4 参数自适应DE的参数F和CR是固定值也可以跑但大规模问题上更容易出现早熟。常见做法是让F和CR随代数变化F 0.5 0.4 * rand; % 每代重新随机 CR 0.9 - 0.5 * (gen / maxGen);前期F大、CR大保证探索广度后期F逐渐减小CR降低提高开发能力。这种做法在很多文献里叫jDE或SaDE实现起来就是多两行代码的事但对收敛稳定性帮助很大。5. 核心代码实现从主循环到适应度函数5.1 DE主循环骨架下面给一个可以直接跑通的最小框架种群规模NP、维度D、最大代数maxGen都是参数。拓扑判断和潮流部分调用了外部函数后面我再展开。function [bestX, bestFitness] CDE_Solve(N, NP, maxGen, paras) D N*(N-1)/2; % 独立决策变量数 pop rand(NP, D); % 初始化[0,1]实数种群 fitness zeros(NP,1); upperIdx find(triu(ones(N),1) 1); for i 1:NP fitness(i) evaluate(pop(i,:), N, upperIdx, paras); end for gen 1:maxGen F 0.5 0.4 * rand; CR 0.9 - 0.5 * (gen / maxGen); % 变异rand/1 randIdx randperm(NP*3); r1 mod(randIdx(1:NP)-1, NP)1; r2 mod(randIdx(NP1:2*NP)-1, NP)1; r3 mod(randIdx(2*NP1:3*NP)-1, NP)1; V pop(r1,:) F * (pop(r2,:) - pop(r3,:)); % 交叉binomial CRMat rand(NP, D) CR; U pop; U(CRMat) V(CRMat); % 边界约束恢复到[0,1]区间 U min(max(U, 0), 1); % 选择可行性优先 uFit zeros(NP,1); parfor i 1:NP uFit(i) evaluate(U(i,:), N, upperIdx, paras); end for i 1:NP if betterThan(U(i,:), uFit(i), pop(i,:), fitness(i)) pop(i,:) U(i,:); fitness(i) uFit(i); end end end [bestFitness, bestIdx] min(fitness); bestX pop(bestIdx,:); endbetterThan函数按可行性优先规则比较两个个体返回值代表U_i是否优于当前个体。5.2 可行性优先规则与约束违反度计算可行性优先的精髓在于三个比较等级function flag betterThan(newX, newF, oldX, oldF) newV constraintViolation(newX); oldV constraintViolation(oldX); if oldV 0 newV 0 flag newF oldF; elseif oldV 0 newV 0 flag false; elseif oldV 0 newV 0 flag true; else flag newV oldV; end end这种写法非常直观但要注意constraintViolation(newX)内部要调用拓扑评估、潮流计算计算量不小。如果每一步选择都调用两次成本有点高。优化方式是在evaluate函数里同时返回约束违反度一次计算搞定。约束违反度的定义很关键。我的做法是归一化处理function [f, viol] evaluate(xVec, N, upperIdx, paras) A zeros(N,N); A(upperIdx) round(xVec); A A A; A(logical(eye(N))) 0; [isRadial, ~] checkRadial(A); % 功率平衡和线路容量约束违反量简化直流潮流 [lineOverload, powerImbalance] powerFlowViolation(A, paras); % 归一化惩罚 viol 0; if ~isRadial viol viol paras.radialPenalty; % 拓扑结构硬约束 end viol viol lineOverload / paras.Sbase; % 容量违反量标幺化 viol viol powerImbalance / paras.Sbase; % 目标函数建设成本 运行成本 buildCost sum(A(:)) / 2 * paras.lineCost ... sum(A(:)) / 2 * paras.switchCost; lossCost paras.energyPrice * paras.lossFactor * (sum(A(:))/2); f buildCost lossCost paras.lambda * viol; end这里的paras.lambda * viol是罚函数部分主要用来引导搜索方向可行性优先规则又在选择层把关。两层配合比单纯罚函数稳定得多。这里的损耗项简化为与连接数相关更精确的做法是直接调用潮流结果里的线损但规划阶段用这种代理变量往往就够用。5.3 直流潮流近似计算对规划阶段的拓扑评估我常用简化的直流潮流function [lineOverload, powerImbalance] powerFlowViolation(A, paras) N size(A,1); B A .* paras.bPerUnit; % 电纳矩阵 B B diag(-sum(B,2)); % 节点电纳矩阵 B(1,:) 0; B(1,1) 1; % 设定参考节点 theta B \ paras.Pinj; % 解线性方程组 linePower zeros(N,N); for i 1:N for j i1:N linePower(i,j) paras.bPerUnit(i,j) * (theta(i) - theta(j)); end end overload max(0, abs(linePower(:)) - paras.lineCap(:)); lineOverload sum(overload); powerImbalance sum(abs(B * theta - paras.Pinj)); end这段代码有两个值得注意的点一是B矩阵第一行置0再设参考节点这是直流潮流的常规做法避免奇异矩阵二是linePower计算那里其实还能继续向量化但为了可读性先保留了循环实际跑大规模测试时我会换成外积写法。6. 收敛性分析、参数调整与踩坑实录6.1 四种典型收敛曲线解读跑完DE后我习惯画两样东西适应度收敛曲线和约束违反度曲线。前者看优化效果后者看约束处理是否正常。常见的异常情况有曲线特征问题处理建议适应度快速下降后长期不动早熟种群多样性丢失增大F尝试dither变体或引入随机重启收敛曲线剧烈震荡罚函数比重不当或CR过大降低CR提高lambda约束违反度始终不为0可行性修复能力不足增加修复策略或把不可行个体的惩罚加得更大适应度下降极慢种群规模偏小或问题规模太大增加NP或在适应度上加入局部搜索我踩得最多的就是第二种情况拓扑个体稍微一变线路容量违反量就跳变导致适应度曲线像心电图。后来把lineOverload / paras.Sbase做标幺化并把lambda从100逐步调到1000曲线才平下来。6.2 Matlab环境里的实操坑Matlab版本差异是我踩过最隐蔽的坑。比如graph(A)这个函数R2015b之后才引入早期版本我只能自己写连通性DFS。还有gcp(nocreate)的参数要求不同版本稍有差异。如果读者用的版本比较新比如R2023b或R2026a函数行为可能又有微调。我自己跑代码时会先写一个单元测试确认checkRadial在目标版本上能正确返回。另外两个高频问题中文注释乱码很多人的脚本从旧机器拷过来编码还是GBK在R2023b里打开一堆乱码。在Matlab编辑器里把文件另存为UTF-8或者直接在偏好设置里把字符集改掉再重新打开脚本就行。第一次启动和并行池初始化很慢通常是后台预加载和杀毒软件扫描导致的。正确做法是先启动Matlab等它完全加载完再跑脚本别一打开就按F5很容易卡在初始化阶段。关于license激活异常我在实验室机器上遇到过几次最稳妥的办法是把license文件路径重新指向一遍确认机器码对应同时检查系统时间是否自动同步避免异常时间导致license校验失败。这类问题属于环境问题和算法本身关系不大但遇到时确实是代码跑不起来的头号阻碍。6.3 参数调整的实用建议针对多微电网规模N我建议的参数范围种群规模NPN10时取40-60N15时取60-100N20以上取150-250。太小容易早熟太大计算量撑不住。缩放因子F基础值0.5范围0.3-0.9。如果收敛太慢把F拉大如果震荡调小。交叉概率CR从0.9开始后期按代数降到0.4。拓扑问题的0-1变量多CR保持较高有利于信息交换。罚系数lambda先跑20代看一下约束违反度的量级再设定为违反量量级的10-100倍别拍脑袋写10000。还有一个小技巧对0-1决策可以在变异和交叉之后做一次“轻微扰动”比如对每个维度以2%的概率翻转值。这个操作类似遗传算法的变异能显著减少早熟概率。代价很小但对最终结果质量影响很大属于性价比很高的参数。写到这里差不多了。最后说一点我个人的体会多微电网拓扑设计这类问题算法本身不是最难的最难的是把工程约束准确翻译成算法能理解的数学形式以及把解码、约束检查和潮流计算这几个模块写稳。约束差分进化算法在Matlab里实现的门槛其实不高核心代码两百行以内就能跑通但每一个细节——上三角编码、连通性检查、可行性优先选择、向量化计算——都会在N变大之后变成影响结果的分水岭。我建议读者拿到代码后先拿N6或N8的小规模问题把每一行都跑懂再逐步放大这样能少走很多弯路。后续如果你还想做多目标版本或者把时序仿真加进去现有的框架也能平滑扩展关键模块都不需要推倒重来。