fsolve求解电力系统潮流计算：MATLAB实现与工程技巧

简介：电力系统分析中，潮流计算是网络规划与运行的基础，其本质是求解一组节点功率平衡的高阶非线性方程。MATLAB优化工具箱中的fsolve作为通用非线性方程组求根器，只需将节点导纳矩阵Ybus与PQ、PV、松弛节点的物理约束映射为F(x)=0的标准形式，即可获得节点电压幅值与相角。该方案不依赖商业潮流软件，适用于教学、科研及中小规模电网离线分析。实现路径涵盖Ybus构造、残差函数封装、fsolve算法选择与参数整定，并辅以三节点系统的完整可复现脚本，同时讨论负荷展布、PV节点无功越限等工程化技巧。

1. 潮流计算到 fsolve 求解器之间，差的是方程而不是电学

电力系统潮流计算的本质，是在已知发电机有功出力、负荷功率和网络拓扑结构的条件下，反求电网中每个节点的电压幅值与相角。因为节点电压的真实值和收敛解之间相互缠绕，整张网络的描述最终落到一组高阶非线性联立方程上。教科书里规定的牛顿—拉夫逊法需要自行推导雅可比矩阵、逐列更新、处理稀疏存储，工程量不小。而 MATLAB 优化工具箱中的 fsolve 是一个通用非线性方程组求根器，它把雅可比矩阵构造、迭代更新和收敛判据全部封装在底层，用户只需要把电网方程写成标准的 F(x)=0 形式，就可以直接获得数值解。

把潮流计算问题交给 fsolve，有两个关键工作要做：一是把网络参数组织成节点导纳矩阵 Ybus，二是把 PQ 节点、PV 节点、松弛节点的物理约束映射成未知量和残差函数。本文面向在 MATLAB 环境下从事电力系统分析的工程师，给出从 Ybus 构造、残差函数封装、fsolve 参数配置到结果校验的完整实现路径。全文不依赖任何商业潮流计算软件包，脚本可复现，参数含义逐一说明。

2. 搭建 Ybus 导纳矩阵与节点功率方程

把一个电气网络交给 fsolve 之前，数学形式必须准确。潮流计算使用标幺值系统，全网电压、功率、阻抗都折算到统一基准之下，这样方程中的系数都在可接受的数值量级内，也是 fsolve 收敛顺利的前提条件之一。后续会解释单位不一致导致的收敛问题，这里先把网络拓扑翻译成 Ybus 矩阵。

2.1 Ybus 的构造方法

一条连接节点 f 和节点 t 的线路，设串联阻抗为 z = r + jx，其支路导纳 y = 1/z。根据电路理论，这条线路对 Ybus 的贡献是：自导纳 Y(f,f) 和 Y(t,t) 各加上 y，互导纳 Y(f,t) 和 Y(t,f) 各减去 y。考虑到线路对地电容和变压器的励磁支路，若存在并联导纳 y_sh，通常在两端节点各加上 y_sh/2。用 MATLAB 函数实现如下：

function Y = buildYbus(branch, n) % branch: Nx4 数组，每行是 [首端节点, 末端节点, 电阻标幺值, 电抗标幺值] % n: 系统总节点数 Z = branch(:,3) + 1j * branch(:,4); y = 1 ./ Z; Y = zeros(n, n); for k = 1:size(branch, 1) f = branch(k, 1); t = branch(k, 2); Y(f,f) = Y(f,f) + y(k); Y(t,t) = Y(t,t) + y(k); Y(f,t) = Y(f,t) - y(k); Y(t,f) = Y(t,f) - y(k); end end

这里把复数运算直接用 MATLAB 的complex类型承载，Ybus 整体是复数矩阵，后续计算有功和无功时需要分别取实部和虚部。函数处理的是一般性线路，没有计入变压器变比和移相器，教学算例和常规小规模电网足够用。若面对抽头变压器，需要对两端自导纳乘上变比平方，这里暂不展开。

2.2 PQ 节点、PV 节点、松弛节点的方程数量分配

在潮流计算中，节点按已知量的不同被分成三类。松弛节点承担全网功率差额，电压幅值和相角都是给定值，通常在系统中只保留一个，作为其他节点相角的参考点。PQ 节点给出有功注入和无功注入，电压幅值和相角是未知数，所以它对应两条方程。PV 节点给出有功注入和电压幅值，无功注入由系统决定，未知量只有相角，对应一条方程。

由上述规则，系统的未知量总数为 2×PQ节点数 + PV节点数，方程数量也必须精确匹配这个数目。常见错误是在 PV 节点上同时写有功和无功方程，导致方程数多于未知量数，fsolve 必然无法收敛。正确的映射是：PQ 节点提供有功误差和无功误差两条残差，PV 节点只提供有功误差一条。为了在程序中方便索引，我习惯用一个结构体或者多列数组统一存放节点信息：

% bus 数组格式: [节点号, 类型, P注入标幺值, Q注入标幺值, V幅值给定] % 类型: 1=松弛, 2=PQ, 3=PV bus = [1, 1, 0, 0, 1.00; 2, 2, -1.2, -0.5, 1.0; 3, 3, 0.6, 0, 1.02]; PQ_idx = find(bus(:,2) == 2); PV_idx = find(bus(:,2) == 3); n_pq = length(PQ_idx); n_pv = length(PV_idx);

注意 P 注入和 Q 注入的符号约定。潮流计算中，节点注入功率定义为发电功率减去负荷功率，也就是说负荷节点的注入功率是负值。这个约定必须写清楚，否则结果符号完全反相，松弛节点的功率也会算错。

2.3 节点功率平衡方程的极坐标展开形式

节点注入功率与电压、导纳之间存在确定关系。设节点 i 的电压相量为 V_i∠θ_i，节点 i 和节点 j 之间的导纳为 Y_ij = G_ij + jB_ij，那么节点 i 的注入有功表达式为：

P_i = sum_j V_i V_j (G_ij cos(theta_i - theta_j) + B_ij sin(theta_i - theta_j))

注入无功表达式为：

Q_i = sum_j V_i V_j (G_ij sin(theta_i - theta_j) - B_ij cos(theta_i - theta_j))

对于 PQ 节点，fsolve 需要找到一组电压幅值和相角，使得上面两式的计算值等于给定的 P_i 和 Q_i。残差函数返回的应该是给定值与计算值的差，不是计算值本身。写成向量后，fsolve 迭代使这个残差向量二范数趋近于 0。

2.4 将节点映射写入残差函数

残差函数必须接收一个向量 x，并从 x 中解出各节点的电压幅值和相角，再计算残差。以下函数是我常用的模板，它把 x 的排列顺序设计为：前 n_pq 个元素是 PQ 节点相角，接下来 n_pq 个元素是 PQ 节点电压幅值，最后 n_pv 个元素是 PV 节点相角。这个排列顺序与 fsolve 的初值 x0 必须完全一致。

function F = power_balance_residual(x, Y, bus, PQ_idx, PV_idx) n = size(Y, 1); theta = zeros(n, 1); V = zeros(n, 1); % 松弛节点 slack_idx = find(bus(:,2) == 1); theta(slack_idx) = 0; V(slack_idx) = bus(slack_idx, 5); % PQ节点：相角从 x 前段取，幅值从 x 中段取 theta(PQ_idx) = x(1:length(PQ_idx)); V(PQ_idx) = x(length(PQ_idx)+1:2*length(PQ_idx)); % PV节点：相角从 x 后段取，幅值使用给定值 theta(PV_idx) = x(2*length(PQ_idx)+1:end); V(PV_idx) = bus(PV_idx, 5); G = real(Y); B = imag(Y); F = zeros(length(x), 1); % 功率计算 for i = 1:n if bus(i,2) == 2 || bus(i,2) == 3 P_calc = 0; Q_calc = 0; for j = 1:n theta_ij = theta(i) - theta(j); P_calc = P_calc + V(i)*V(j)*(G(i,j)*cos(theta_ij) + B(i,j)*sin(theta_ij)); Q_calc = Q_calc + V(i)*V(j)*(G(i,j)*sin(theta_ij) - B(i,j)*cos(theta_ij)); end if bus(i,2) == 2 row_p = find(PQ_idx == i); F(row_p) = bus(i,3) - P_calc; F(length(PQ_idx) + row_p) = bus(i,4) - Q_calc; elseif bus(i,2) == 3 row_pv = find(PV_idx == i); F(2*length(PQ_idx) + row_pv) = bus(i,3) - P_calc; end end end end

代码逻辑上先重建电压向量，再按节点类型计算残差。注意 PQ 节点和 PV 节点在 bus 数组中的编号与它们在 PQ_idx 中的位置不是一回事，find用来建立这两者之间的映射。这个函数的时间复杂度是 O(n^2)，对几百个节点的系统毫无压力，但对上千节点的大型系统，可以考虑用稀疏矩阵配合向量化功率计算函数Ybus * V，这里不再展开。

3. fsolve 的算法选择与关键参数整定

MATLAB fsolve 提供 trust-region-dogleg、trust-region-reflective、levenberg-marquardt 三种算法。潮流计算中残差函数平滑、规模有限，默认的 trust-region-dogleg 通常表现最好，它结合了牛顿法的快速收敛与信赖域方法的稳定性。

3.1 算法特性比较

trust-region-dogleg 适合中小规模的非线性方程组，利用雅可比矩阵的 QR 分解构造迭代方向，整体数值稳定性好，默认选项下大多数潮流算例都能一次收敛。levenberg-marquardt 在残差函数接近线性时有优势，但对初值更敏感，而且求解过程中涉及正规方程的条件数问题，在节点多、导纳矩阵稀疏性强的系统中表现不如 dogleg。trust-region-reflective 主要面向带边界约束的问题，潮流方程本身没有变量边界，除非要额外嵌入电压上下限约束，否则没有必要选用。

3.2 TolFun、TolX、MaxIterations 的工程推荐值

opts = optimoptions('fsolve', ... 'Algorithm', 'trust-region-dogleg', ... 'Display', 'iter', ... 'TolFun', 1e-10, ... 'TolX', 1e-12, ... 'MaxIterations', 200, ... 'MaxFunctionEvaluations', 2000, ... 'StepTolerance', 1e-8);

TolFun控制残差向量二范数的停止阈值，一般取 1e-8 到 1e-10。取值太小会让求解器在达到收敛边界后继续空转，太大则结果不够精确。TolX是状态变量变化量的阈值，当连续两次迭代的解变化小于该值时停止，这里取 1e-12 比较保守，能够避免因为电压幅值更新过小而提前退出。MaxIterations通常设 100 到 300 次，牛顿法正常情况下几十次迭代就能达到机器精度，如果超过 200 次还不收敛，多半是网络参数有问题，此时继续增加迭代上限毫无意义。

3.3 匿名函数句柄传参模式

残差函数除了 x 以外还需要 Ybus 和节点数组，而 fsolve 的调用接口只接受 f(x)，要把额外参数绑定进去。最常见做法是匿名函数捕获工作区变量：

F_handle = @(x) power_balance_residual(x, Y, bus, PQ_idx, PV_idx); [x_sol, fval, exitflag] = fsolve(F_handle, x0, opts);

调用后[x_sol, fval, exitflag]三个返回值的含义分别对应解向量、最终残差和退出标志位。exitflag=1 表示收敛到了解，exitflag=0 表示迭代次数达到上限，exitflag=-1 或 -2 表示算法因某些数值原因终止，这时必须检查残差函数是否定义连续，或者更换初值和算法。

4. 从头验证的一个潮流问题：三节点系统的完整步骤

理论叙述再充分，也不如一个完整算例来得直观。这里以一个三节点系统为例，包含一个松弛节点、一个 PQ 负荷节点、一个 PV 发电节点，足以覆盖所有节点类型和收敛路径。

4.1 节点与支路参数设定

三个节点通过两条支路连接：支路1-2的串联阻抗为 0.02+j0.06，支路1-3的串联阻抗为 0.03+j0.09，支路2-3的串联阻抗为 0.04+j0.12，所有阻抗都是标幺值。节点数据如下表所示：

节点1作为平衡节点吸收全网功率差额，节点2作为负荷节点从电网吸收有功和无功，节点3带有发电机，有功注入为 0.8，电压幅值设定为 1.02。注意所有支路都没有对地导纳，因此并联导纳项为 0。

4.2 从初值到求解的完整脚本

branch = [1, 2, 0.02, 0.06; 1, 3, 0.03, 0.09; 2, 3, 0.04, 0.12]; bus = [1, 1, 0, 0, 1.00; 2, 2, -1.0, -0.5, 1.0; 3, 3, 0.8, 0, 1.02]; n = 3; Y = buildYbus(branch, n); PQ_idx = find(bus(:,2) == 2); PV_idx = find(bus(:,2) == 3); n_pq = length(PQ_idx); n_pv = length(PV_idx); % 平启动初值：相角为0，PQ节点电压幅值为1 x0 = zeros(2*n_pq + n_pv, 1); x0(n_pq+1:2*n_pq) = 1; % PQ节点幅值取1.0 opts = optimoptions('fsolve', ... 'Algorithm', 'trust-region-dogleg', ... 'Display', 'final', ... 'TolFun', 1e-10, ... 'TolX', 1e-12); F_handle = @(x) power_balance_residual(x, Y, bus, PQ_idx, PV_idx); [x_sol, fval, exitflag] = fsolve(F_handle, x0, opts); % 分解结果 theta = zeros(n, 1); V = zeros(n, 1); slack_idx = find(bus(:,2) == 1); theta(slack_idx) = 0; V(slack_idx) = bus(slack_idx, 5); theta(PQ_idx) = x_sol(1:n_pq); V(PQ_idx) = x_sol(n_pq+1:2*n_pq); theta(PV_idx) = x_sol(2*n_pq+1:end); V(PV_idx) = bus(PV_idx, 5); % 输出关键结果 fprintf('节点 2 电压幅值: %.6f pu, 相角: %.6f rad\n', V(2), theta(2)); fprintf('节点 3 电压幅值: %.6f pu, 相角: %.6f rad\n', V(3), theta(3));

这个脚本的关键点在于x0的构造顺序必须与残差函数中对 x 的拆解顺序完全一致。如果 PQ 节点在 bus 数组中的位置和它在 PQ_idx 中的位置不一致，disp出来的结果会串号，但 fsolve 本身不会报错，这种 bug 排起来很闹心。一个稳妥的做法是先打印PQ_idx和PV_idx内容，和 bus 表逐一核对。

4.3 结果合理性判据

拿到 x_sol 之后，先用fval检查残差是否接近零向量。如果fval的最大绝对值超过 1e-6，说明没有收敛到高精度解，此时即使 exitflag=1，结果也只能作为粗略参考。其次，松弛节点的有功逸出了系统的不平衡量，有功全网应当满足发电等于负荷加网损，节点1的有功流出方向和数值可以用支路潮流计算公式验证。

多数情况下，节点2的电压幅值会比设定值 1.0 略低，节点3电压幅值等于 1.02 的给定值。如果节点2电压大于 1.0，且节点3无功输出极端，说明网络参数或符号约定可能出了错。另一点值得注意，PV 节点的无功注入在初始数据里没有给定，收敛后要从残差里反推电力值，这个值必须在发电机无功上下限范围内，否则 PV 节点建模有误，下一节会详细讨论。

5. 工程化适配：初值策略、PV节点无功越限与节点类型转换

基础算例跑通之后，将脚本迁移到自己的网格中会遇到更现实的问题：负荷重载带来的收敛振荡、PV节点无功功率越限导致解不可行，以及大型网络中的计算效率和变量尺度问题。这一节把工程中常见的三个坎逐一拆开。

5.1 负荷展布法解决重载收敛振荡

当系统负荷接近极限传输功率，即潮流方程存在多个解或仅有临界解时，平启动的 fsolve 常陷入迭代发散或两个解之间振荡。最直接的解决办法是负荷展布法：把目标负荷按比例分成 0.2、0.4、0.6、0.8、1.0 五个阶段逐级求解，每阶段以上一阶段的解作为初值。这个方法在数学上等价于沿着连续参数路径跟踪解曲线，能避开鞍结点附近的收敛陷阱。

lambda_list = 0.2:0.2:1.0; x_current = x0; for lam = lambda_list bus_lam = bus; % 按比例缩放PQ节点的P和Q给定值 for i = 1:n if bus(i,2) == 2 bus_lam(i,3) = bus(i,3) * lam; bus_lam(i,4) = bus(i,4) * lam; elseif bus(i,2) == 3 bus_lam(i,3) = bus(i,3) * lam; end end F_handle = @(x) power_balance_residual(x, Y, bus_lam, PQ_idx, PV_idx); [x_current, ~, ef] = fsolve(F_handle, x_current, opts); if ef <= 0 warning('负荷水平 %.1f 处未收敛, 停止计算', lam); break; end end

lambda_list的步长并不是越小越好。步长过小会导致总迭代次数翻好几倍，步长过大则失去了延拓的意义。实用经验是 0.2 的步长足够应对绝大多数重载场景，只有在接近静态电压稳定极限时才把末两段步长加密到 0.05。

5.2 PV 节点无功越限的处理

PV 节点的物理含义是发电机恒压控制。但发电机无功输出存在上下限，若收敛结果中 PV 节点无功越限，这个节点的约束就不再可行，必须转为 PQ 节点且 Q 固定在极限值。实际工程中处理流程是：先按 PV 节点求解一次多个潮流断面，如果某个 PV 节点无功持续越上限，则将其类型改为 PQ 节点并给 Q 赋上限值。这相当于放松了该节点的电压控制能力，更接近真实物理世界中的发电机失磁场景。在执行节点类型转换时，需要注意未知量数量和残差函数长度都变了，x0 要做相应的增删，否则 fsolve 直接报错。这个步骤不应该看作求解失败，而是一种有价值的边界信息：它告诉你这个网络中哪些节点缺乏无功支撑。

5.3 变量尺度与雅可比矩阵状况的改善

当节点电压标幺值在 1.0 左右、相角在 0.1 弧度量级时，变量量级差异不大，fsolve 的默认尺度处理足够稳。如果直接使用有名义值，电压动辄 10kV 级、功率动辄 100MW 级，雅可比矩阵各列元素差距巨大，线性代数求解误差急剧放大。所有工程节点数据在馈入潮流计算前必须完成标幺化。使用 MATLAB 的diag预调节是不可取的，最干净的方法是直接使用基准值换算。在所有数据统一到 1 pu 附近的前提下，fsolve 内建的自动缩放机制已经足够应对大多数问题。

另外一个容易被忽略的细节是导纳矩阵对角元是否含有零。若某个节点没有任何支路连接，它对应的 Ybus 对角线是 0，雅可比矩阵奇异，fsolve 必然报错。这就是为什么隔离子网必须与主网相连，否则即使程序逻辑正确，数学上依然无解。

6. 收敛验证的三大指标与一个周期性求解技巧

判断一次 fsolve 潮流计算是否可靠，不能只看返回标志。工程实践中有三条验证路径，可以筛掉绝大多数隐性错误。

第一是残差范数检查：用返回值 fval 计算norm(fval, inf)，与设定的 TolFun 阈值对比，如果最大残差仍在 1e-6 以上，说明解只是局部稳定点，没有达到期望精度。第二是全网功率平衡检查：把解出的节点电压代回功率方程，重新计算全网 P 和 Q 的总注入，加上网损之后，松弛节点的功率应正好补足差额。第三是支路潮流反向验证：从电压解和 Ybus 出发计算每条支路的潮流，检查它是否与节点注入连续性方程匹配。这三种检查用 MATLAB 向量化操作可以在数行内完成：

V_phasor = V .* exp(1j * theta); I_inj = Y * V_phasor; S_inj = V_phasor .* conj(I_inj); % S_inj的实部与虚部应分别等于各节点的给定注入功率 residual_power = max(abs([real(S_inj) - bus(:,3) '; imag(S_inj) - bus(:,4) ']));

代码中I_inj = Y * V_phasor是节点注入电流的向量化表达式，任意母线节点注入功率等于电压乘以电流共轭。这段验证逻辑的好处是彻底绕开了残差函数本身的实现，用完全独立的计算方法复核解的正确性。

至于周期性求解技巧，在多断面连续潮流计算场景中，逐次调用的初值可以从上一次解直接搬过来用。比如研究日内 24 小时负荷曲线时，把 t 时刻的解作为 t+1 时刻的初始猜测，通常只需要两次迭代就能收敛到新平衡点。而如果每天都从平启动开始，除了浪费大量计算资源，还可能在某些临界负荷水平处无谓发散。合理的初值利用不仅是性能优化，更是让 fsolve 沿功率参数空间平滑追踪解轨迹的关键手段。把这一招和前面的负荷展布法结合，等于同时拥有了速度和稳健性，也是我处理数十个连续断面算例时最常用的一套完整方案。

本文还有配套的精品资源，点击获取