做配电网优化研究的人大概率都经历过这种尴尬牛拉法潮流算得挺准可一旦想把潮流方程塞进优化模型里非线性方程一进去求解器直接翻脸初始点稍微不好就开始发散。后来大家都在用Distflow解决这件事。Distflow本质上是支路潮流模型在辐射状网络下的标准表达既能做潮流计算又能写成适合优化求解的约束形式配合线性化处理后配电网重构、DG接入、储能调度这类问题才真正有了工程化的解法。我最近在Matlab里把这套东西从非线性迭代到线性化版本完整过了一遍并且针对不同拓扑的网络做了适配。这篇文章就把整个过程中的推导、代码、误差分析和踩坑记录整理出来适合正在做配电网方向的研究生以及刚接触Distflow、想快速上手Matlab实现的工程师参考。1. Distflow到底在算什么支路潮流模型不是玄学1.1 两条基础定律推出来的四个方程Distflow的源头并不复杂就是基尔霍夫电压定律和电流定律在辐射型配电网上的重新排列。传统牛拉潮流法的变量是节点电压幅值和相角方程写出来之后雅可比矩阵几乎是稠密的节点一多求解复杂度就上去了。Distflow换了一套变量每一条支路的有功P_ij、无功Q_ij、电流幅值平方l_ij再加上节点电压幅值平方V_i^2。对于一条从节点i指向节点j的支路假设节点i是父节点节点j是子节点Baran和Wu在1989年给出的经典形式是P_ij sum(P_jk) p_j r_ij * l_ij Q_ij sum(Q_jk) q_j x_ij * l_ij V_j^2 V_i^2 - 2*(r_ij*P_ij x_ij*Q_ij) (r_ij^2 x_ij^2)*l_ij l_ij (P_ij^2 Q_ij^2) / V_i^2这里P_ij、Q_ij是支路首端功率也就是从节点i流向j的功率p_j、q_j是节点j上的负荷功率r_ij、x_ij是支路阻抗l_ij是电流幅值的平方。前两个方程是节点功率平衡第三个是电压降落方程第四个是电流定义。四个方程里前三个在给定l_ij之后全是线性方程唯一真正“难搞”的是第四个式子它把电流平方和支路功率、节点电压耦合在了一起这也是整个Distflow模型非线性来源的核心。1.2 辐射网的树结构让一切变简单Distflow之所以在辐射状网络里特别好用是因为辐射网的拓扑是一棵树每个节点有且只有一个父节点功率只能沿着唯一的路径从根节点流向末端。这意味着每一个节点上的功率平衡方程里流入项只有一个来自父支路流出项是若干子支路结构非常干净。如果网络带环同一个注入功率就有多条传播路径方程之间的耦合关系立即变复杂。配电网虽然理论上设计成辐射状但实际中联络开关闭合会产生弱环这也是后文要单独讨论适配方案的原因。辐射网的树结构还带来一个计算上的优势可以直接用前推回代法求解。这个方法本质上就是在做Distflow方程的定点迭代不需要计算雅可比矩阵收敛速度快占用内存小特别适合配电网这种“长而瘦”的网络。1.3 非线性项藏在哪心里要有数我刚开始学Distflow时总觉得公式多记不住。后来发现只要抓一个关键点就全通了整个模型其实只有两处非线性。第一处是l_ij的定义也就是电流平方等于功率平方除以电压平方。在优化问题里这个式子代表的是二次等式约束非凸很难直接处理。第二处是电压降落方程里的(r_ij^2 x_ij^2)*l_ij这一项它把l_ij乘上了阻抗本身是线性项没引入新的非线性但因为l_ij依赖功率和电压所以整个方程依然是耦合的。理解了这两处非线性后面的线性化思路就顺理成章了要么把l_ij这个变量消掉要么把含有l_ij的项直接忽略掉代价是损失一部分网损和电压降的精度换来的是整个方程组变成线性。工程上为什么要这么做下一章慢慢说。2. 基础版Distflow的Matlab实现先用前推回代把潮流跑准2.1 数据结构parent-children比邻接矩阵好用在Matlab里实现Distflow第一步不是写方程而是设计网络的数据结构。很多新手喜欢用邻接矩阵存网络后面写方程时又要一层层查索引代码又慢又容易错。我建议直接用两个结构parent数组和children元胞数组。parent(j)表示节点j的父节点编号根节点的parent设为0children{i}是一个数组列出节点i的所有子节点。配合支路阻抗数组Z其中Z(j)表示支路(parent(j), j)的阻抗根节点的Z(1)设为0。这个结构的优势在于前推和回代两个过程都只需要按照节点编号顺序从头到尾扫一遍不需要反复查询邻接关系。对辐射网来说节点从根到末端依次编号之后子节点的编号必然大于父节点这为后面向量化也埋下了伏笔。2.2 前推回代核心代码实现下面这段代码是我实测可以稳定运行的基础版前推回代函数。输入是网络结构和节点注入功率输出是节点电压和支路功率。function [V, S_branch, I_branch, iter] distflow_fbs(parent, children, Z, S, V0, tol, maxIter) % parent(j): 节点j的父节点编号parent(1)0 % children{i}: 节点i的子节点编号数组叶子节点为空数组 % Z(j): 支路(parent(j), j)的阻抗标幺值Z(1)0 % S: 节点复功率注入负荷为正DG出力为负 % V0: 根节点电压通常取1.0 % tol: 收敛阈值 % maxIter: 最大迭代次数 n length(parent); V ones(n, 1) * V0; % 复数电压初始值 I_branch zeros(n, 1); % 支路电流I_branch(j)对应支路(parent(j), j) S_branch zeros(n, 1); for iter 1:maxIter V_old V; % 回代从末端向根计算支路电流 for j n:-1:2 I_load conj(S(j) / V(j)); % 节点注入电流 所有下游支路电流 child_sum 0; for k children{j} child_sum child_sum I_branch(k); end I_branch(j) I_load child_sum; end % 前推从根向末端更新电压 V(1) V0; for j 2:n i parent(j); V(j) V(i) - Z(j) * I_branch(j); end if max(abs(V - V_old)) tol S_branch V(parent(2:n)) .* conj(I_branch(2:n)); return; end end error(前推回代不收敛请检查网络参数或初始值); end注意这段代码用的是电流形式好处是不需要先算网损回代公式里已经天然包含了支路电流对功率平衡的影响。电压更新时直接用V(i)减去阻抗压降Z(j)*I_branch(j)逻辑非常直观。2.3 初值、收敛判据与加速经验前推回代对初值不敏感因为迭代格式本身就是围绕辐射网物理特性的不像牛顿法那样需要好初值保证收敛。把所有电压设成根节点电压通常三到五次迭代就能收敛到10^-8精度速度非常快。但有两个细节必须注意。第一收敛判据用max(abs(V - V_old))虽然简单但遇到某些特殊网络会提前误判。比如个别节点电压幅值特别小的时候绝对误差可能已经很小但相对误差还很大。我在代码里用绝对值是因为配电网电压通常在0.9到1.1之间绝对值误差可以接受但如果你的网络重载比较严重建议改成max(abs(V - V_old) ./ max(abs(V), 1e-6))。第二回代过程里出现了内层for循环遍历children{k}这一步在节点数多时是性能瓶颈。优化方法是把children关系按网络层次分好组然后向量化逐层更新。对于几百节点的配电网for循环版本其实已经足够快没必要过度优化。3. 线性化Distflow把约束矩阵装好就能交给优化求解器3.1 为什么优化模型里不能用完整Distflow做配电网优化时比如最小网损重构、DG出力优化、储能充放电调度需要把潮流方程作为约束放进优化模型里。非线性Distflow公式中有l_ij与功率、电压的二次耦合这类约束会让可行域变成非凸集合。交给fmincon这类局部求解器结果严重依赖初始点交给全局优化器节点规模一上去就完全跑不动。这就是线性化Distflow存在的根本原因。它不是要在精度上取代完整潮流而是提供一种从工程角度“足够精确”的潮流近似让优化问题变成线性约束或二阶锥约束从而交给成熟的商业求解器稳定求解。3.2 LinDistFlow的三个关键近似线性化的方式有很多种最经典的是LinDistFlow核心假设是电压接近额定值线路电流不太大网损相对传输功率可以忽略。具体在公式上做了三处处理第一l_ij近似为0于是有功和无功平衡方程变成P_ij sum(P_jk) p_j Q_ij sum(Q_jk) q_j第二电压降落方程中的(r_ij^2 x_ij^2)*l_ij项直接忽略得到V_j^2 V_i^2 - 2*(r_ij*P_ij x_ij*Q_ij)第三如果需要电压本身而不是电压平方可以在V1附近做一次泰勒展开近似写成V_j V_i - (r_ij*P_ij x_ij*Q_ij)我实际使用时更推荐保留V^2变量的版本。因为在电压偏离1较远的末端节点直接线性化V本身会引入额外误差而保留平方项在后续电压上限约束里也更自然。除非你有特殊需求否则别为了省一个变量去用线性电压版本。3.3 矩阵化组装代码可直接抄线性化Distflow进入优化模型后最关键的一步是把方程组装成稀疏矩阵形式。以变量顺序x [P; Q; u]为例其中P和Q是支路有功和无功u是节点电压幅值平方下面的函数可以直接构建约束矩阵。function [Aeq, beq] build_lindistflow_matrix(n_parent, children, R, X) % n_parent: parent数组n_parent(1)0 % children: children元胞数组 % R, X: 支路电阻和电抗标幺值数组R(1)0, X(1)0 % 返回 Aeq, beq使得 Aeq * [P;Q;u] beq n length(n_parent); nl n - 1; % 支路编号支路b对应节点b1方向为 parent(b1) - b1 % 功率平衡矩阵 B (n x nl) B sparse(n, nl); for j 2:n i n_parent(j); b j - 1; B(j, b) 1; % 流入支路b B(i, b) B(i, b) - 1; % 从支路b流出 end % 电压方程矩阵(nl x nl) 2R, 2X 和 (nl x n) u系数 A_p 2 * diag(R(2:n)); A_q 2 * diag(X(2:n)); A_u sparse(nl, n); for j 2:n i n_parent(j); b j - 1; A_u(b, j) 1; A_u(b, i) -1; end % 组装完整Aeq % 变量块P(nl), Q(nl), u(n) Aeq [B, sparse(n, nl), sparse(n, n); A_p, A_q, A_u]; beq zeros(n nl, 1); end注意功率平衡方程的行向量对应节点j的注入功率约束beq中前n行对应负荷向量p公式为B*P p。实际调用时把节点负荷和分布式电源出力合并到p里即可。之所以用稀疏矩阵是因为配电网节点数动辄几百上千稠密矩阵会让约束的内存占用爆炸。代码里虽然用了for循环构建索引但最终生成的是sparse矩阵交给求解器实际求解时的内存表现很好。3.4 接入linprog、quadprog、YALMIP/CVX的写法矩阵组装好之后接入不同求解器只是换一层调用壳而已。我举个最简单的例子假设目标函数是近似网损最小用quadprog求解% 变量 [P; Q; u] P sdpvar(nl, 1); Q sdpvar(nl, 1); u sdpvar(n, 1); % 近似网损目标sum(R .* P.^2) sum(R .* Q.^2)用标幺值近似忽略电压分母 cost P * diag(R(2:n)) * P Q * diag(R(2:n)) * Q; Constraints [B*P p_load - p_pv]; % 加上电压方程约束根节点电压固定节点电压上下限 Constraints [Constraints; A_p*P A_q*Q A_u*u 0]; Constraints [Constraints; u(1) 1.0]; Constraints [Constraints; 0.95^2 u 1.05^2]; optimize(Constraints, cost);如果你偏好CVX写法几乎一样把sdpvar换成cvx变量optimize换成cvx_begin/cvx_end即可。实际工作中我更倾向直接用YALMIP因为后面如果要加二进制变量做网络重构YALMIP的接口切换非常方便底层求解器可以随时从linprog换成gurobi或cplex。4. 网络拓扑一变就崩多馈线、弱环网、DG节点的适配方案4.1 多馈线/多变电站虚拟根节点技巧实际配电网很少有单一电源。两条甚至多条馈线从不同变电站引出各自带一片负荷有时还有联络开关。遇到这种情况很多人直接在程序里写死一个根节点结果功率平衡方程怎么都对不上。解决办法是加一个虚拟根节点。把原本的多个根节点比如节点1和节点2都连接到虚拟根节点0上用阻抗为0的支路表示。这样整个网络仍然是一棵树物理上多电源并列的问题被转化成虚拟根到各馈线根节点的纯拓扑连接Distflow方程照常成立。要注意的是虚拟根节点的电压设置。当多个实际电源点电压不完全一致时虚拟根相当于一个等值电源点直接设为1.0 p.u.实际变电站出口母线电压再通过各自的节点电压约束去限制。4.2 弱环网联络线的开断与Big-M约束配电网正常运行时是辐射状的但检修或故障转供时联络开关闭合会形成弱环。如果直接在图上做Distflow环网会破坏树结构前推回代不再适用。处理方式取决于你的目标是纯粹潮流计算还是优化。如果是优化问题通常会把联络线作为可选闭合支路用二进制变量z表示开断状态并施加Big-M约束-P_max*(1-z) P_b P_max*(1-z) -Q_max*(1-z) Q_b Q_max*(1-z)这里P_max取一个大数至少要大于该支路可能流过的最大功率。再加上辐射性约束比如闭合支路数等于节点数减一并且网络保持连通就能把“选哪些联络线闭合”变成优化问题的一部分。如果只是给定的一种弱环网状态不想做优化那么用牛拉法或其他通用潮流求解器更省事。我个人的习惯是优化场景用线性Distflow加开断约束精确潮流场景直接上Matpower两种工具互补。4.3 DG节点与储能接入功率平衡右边改一个量含DG节点的网络适配逻辑上是最简单的只需要把功率平衡方程的右侧从纯负荷向量改成负荷减DG出力。对于恒功率模型DG的P和Q是固定的直接改右侧常数。对于DG出力可调的优化问题DG的出力要变成优化变量这意味着功率平衡约束要从Aeqx beq变成Aeqx 其他变量 beq本质上是在变量向量里增加DG有功无功变量并在目标函数里加入对应成本或消纳项。储能节点特殊在它有时间耦合。单时段潮流方程里储能就像一个可正可负的负荷但多时段优化时储能电量与前后时段的关系要靠额外的状态约束描述潮流方程本身不需要改动。很多新手把储能约束写进潮流方程里导致变量数和约束规模爆炸其实没必要。4.4 拓扑自动识别用graph和BFS自动生成树结构前面所有方法都依赖parent和children结构手动维护很麻烦。尤其是多次调用不同网络时我建议直接利用Matlab的graph对象自动生成树结构。下面这个函数可以从边列表生成parent和children。function [parent, children] extractTree(edgeList, n, root) % edgeList: 两列的边表一行表示一条支路两个端点 % n: 节点总数 % root: 根节点编号 G graph(edgeList(:,1), edgeList(:,2), n); T minspantree(G, Root, root); % 辐射网本身就是生成树 parent zeros(n, 1); children cell(n, 1); visited false(n, 1); queue root; visited(root) true; while ~isempty(queue) i queue(1); queue(1) []; nbrs neighbors(T, i); for j nbrs if ~visited(j) visited(j) true; parent(j) i; children{i} [children{i}, j]; queue(end1) j; end end end end这个函数有个好处哪怕你的边表初始是乱序的甚至网络带环minspantree都会生成一个以root为根的生成树保证Distflow可解。当然带环网络做生成树意味着物理上断开了某些支路对应到优化问题里就是联络线开断需要自行决定哪些支路可以断开。5. 线性化误差实测什么情况下结果开始失真5.1 误差来源重载、长线路、高R/X线性化Distflow最核心的近似就是忽略网损和电流平方项。误差不是均匀分布的它和三个因素高度相关负载率、线路长度、线路电阻电抗比。负载率越重的网络流过支路的功率越大电流平方增长网损占总功率的比例上升忽略l_ij带来的误差自然变大。长线路的阻抗大电压降方程里忽略的(r^2x^2)*l项变大。高R/X比的线路则更容易让电压误差方向发生偏移因为无功对电压的影响权重和真实模型不同。另外末端节点的误差普遍比靠近根节点的节点大。这是因为误差在沿着树从根向末端逐级累积每经过一条支路就叠加上一次近似误差。对长链式网络来说末端电压的误差可能比想象中大不少。5.2 一个基准算例的实测结果我以IEEE 33节点标准网络做了一组对比前推回代求解完整Distflow作为基准线性化Distflow用纯Matlab线性方程直接解。在基准负荷下各节点电压幅值误差最大约0.5%到1%这个精度对大多数优化决策来说完全够用。把负荷放大到1.5倍后误差明显上升。靠近末端节点电压幅值误差达到2%左右而网损误差更夸张因为线性化Distflow直接忽略了网损所以用它的功率解去估算系统总网损时误差普遍在5%到10%之间。这也说明一个关键点如果你优化目标本身就是网损最小化线性化模型给出的“最优”网损数值会系统性偏小但最优解对应的开关组合或DG出力策略通常仍然是合理的。5.3 误差不可接受时怎么办如果案例里的电压偏差超过5%线性化模型就不太可靠了。我有两个备选方案。第一个方案是分段线性化或者迭代修正。先用线性化Distflow求一个初始解把解代入完整Distflow计算真实网损再把网损作为等效负荷加回功率平衡方程重新求解。迭代两三轮后精度明显提升代价是计算时间增加但对中等规模网络完全能接受。第二个方案是SOCP松弛。把l_ij保留把等式l_ij (P^2Q^2)/V^2放宽成不等式l_ij (P^2Q^2)/V^2写成二阶锥形式后用Mosek或Gurobi求全局最优解。这个方法在辐射网上被证明松得很紧基本能拿到原始模型的精确解代价是建模复杂度提高变量数目翻倍。对于科学研究或者对精度有硬性要求的场景我会直接选SOCP不再纠结线性化精度。6. 我在Matlab里写Distflow踩过的坑与绕坑方法6.1 节点编号与父子关系错位矩阵全是乱的这是最隐蔽也最致命的坑。根节点不一定是1号节点支路编号和节点编号之间也没有天然对应关系。我早期写build_lindistflow_matrix函数时默认支路b对应节点b1结果换了一个边表顺序完全乱掉的网络后功率平衡矩阵B的第一行都错了后面的电压方程更是完全没有意义。后来我强制规定所有内部函数接收的第一参数统一是parent数组而不是边表。边表只在预处理阶段用一次转成parent和children后后续所有计算都只依赖parent和children杜绝了编号歧义。如果你是从第三方数据读入网络第一步一定是先做数据校验确保每个节点有且只有一个父节点。6.2 电压平方与电压本身从头到尾只用一个线性化Distflow里u V^2和V本身是混用的重灾区。我在一个项目里把电压下限0.95直接写到u上约束写成了u 0.95结果解出来的电压“合格”了实际上真实电压只有0.975左右完全错误。正确写法是u 0.95^2。这不是简单的错误是定义侧的不统一。建模时必须决定整套代码用u还是V作为变量。如果选了u那么所有电压约束、目标函数里的电压项、初始值都要用u表示如果选了V则电压降落方程不能直接用V_j V_i - (rPxQ)因为那是从u形式近似推导的。用u作为变量最后输出结果时再统一开根号反而是最不容易出错的方案。6.3 收敛判据别用绝对误差前推回代里如果只用max(abs(V - V_old)) tol在电压本身就比较大的量级下可能没问题但遇到带储能或DG的模型无功功率很大导致电压变化缓慢时绝对误差判据可能提前停止迭代解还没真正收敛。建议使用max(abs(V - V_old) ./ max(abs(V), 1e-6))这样的混合判据。另外要注意有些网络的迭代过程会出现轻微振荡绝对误差先降到阈值以下随后又反弹。判断收敛要同时满足V的变化量和支路功率的变化量都小于阈值或者直接观察连续两次迭代的误差序列不要因为单次满足条件就停止。6.4 向量化与稀疏矩阵从for循环地狱走出来写Distflow和线性化模型的时候for循环虽然直观但节点数超过1000之后性能会变得很难看。优化模型里的约束矩阵我强烈建议用sparse函数直接构建尽量避免用zeros填稠密矩阵再接sparse。矩阵构建如果循环不可避免可以先收集行索引、列索引、数值的三列向量最后一次性sparse(idx_i, idx_j, val, m, n)。这个技巧在构建几百条支路的大网络时提速效果明显而且内存占用少一个量级。前推回代里的内层for循环也类似。如果网络是多层级结构可以用层次分组的方式先找所有叶子节点统一计算一层再向根推进。大部分配电网分层迭代的加速效果很可观。6.5 与Matpower、YALMIP衔接时要注意的单位问题Matpower的默认数据格式是mpc结构节点导纳矩阵、支路参数用的都是标幺值。如果你从Matpower读网络要注意负荷和发电机容量可能用的是MW/MVar而优化模型里功率变量常用p.u.转换公式是p.u. 实际值 / Sbase。根节点电压、支路阻抗这类标幺值则不需要额外处理。YALMIP和CVX内部默认所有变量都是普通实数约束的单位和建模时保持一致即可。但求解器的数值稳定性对变量尺度敏感功率若用MW单位变量可能是几百甚至上千对求解效率有影响。我统一把功率和电压全部转到标幺制下建模到输出结果时再乘回Sbase这样既方便调试也方便从一套数据里同时得到潮流结果和优化结果。最后再分享一个小技巧无论你用的是完整Distflow还是线性化版本一定要写一个自检函数把解出的支路功率代回功率平衡方程看残差量级。我每搭一套新网络第一步就是跑一次这个自检这比任何调试器都管用。Distflow这套工具本身不复杂复杂的是网络数据的干净程度和你对自己模型定义的一致性。把这两点控制好后面无论网络怎么换代码都能稳稳跑起来。