前几天被一位师弟问起连续功率流（Continuation Power Flow，CPF）在MATLAB里怎么实现，他想复现一篇IEEE-14节点系统的电压稳定性分析。我顺手把以前做过的流程整理成了可复用的脚本，从读取IEEE-14标准数据、构建导纳矩阵，到写出预测-校正主循环、绘制P-V曲线，全程不到三百行。今天把这套东西完整拆开讲，重点说清楚每一步为什么这么做，以及实际跑起来会踩到哪些坑。
CPF解决的核心问题是：当系统负荷沿某个方向持续增长时，电压能不能撑住、能撑到哪一步。它输出的是负荷增长因子λ与节点电压幅值的关系曲线，也就是常说的P-V曲线，曲线末端的λ值就是系统的静态电压稳定裕度。IEEE-14节点系统因为规模适中、数据公开、节点类型和变压器配置齐全，经常被拿来做算法验证平台。MATLAB做这个事情很顺手，矩阵运算、迭代调试、数据可视化都在一个环境里完成，省去来回切换工具的时间。
我用的是较新的MATLAB版本，但下面这段代码没有用任何特殊版本特性，R2016a之后的版本应该都能直接跑。基础较弱的朋友也不用担心，我会把牛顿拉夫逊潮流、雅可比矩阵、SVD求切线这些关键点都用大白话讲清楚。
1. 思路拆解：连续功率流到底在算什么
1.1 普通潮流在电压崩溃点为什么会失败
先看普通潮流做了什么。给定发电机出力和负荷功率，求解各节点电压幅值和相角，本质是求解一组非线性方程F(x)=0，其中x由非平衡节点的相角和PQ节点的电压幅值组成。这种问题用牛顿拉夫逊法迭代求解，每一步都要求解潮流雅可比矩阵J。
当负荷持续增加时，运行点会沿着P-V曲线向上移动，系统的雅可比矩阵行列式逐渐趋近于零。在接近电压崩溃点的时候，J变得病态，牛顿法的平方收敛特性严重退化：你计算出的修正量越来越大，迭代却总是不收敛，最后直接发散。也就是说，普通潮流程序在电压崩溃点附近就“罢工”了。
但工程师和研究者恰恰最关心崩溃点在哪、裕度有多大。连续功率流就是为了解决这个问题而提出的：它不把某个负荷水平下的潮流解当作唯一目标，而是把负荷增长因子λ当作一个新变量，把问题变成追踪一整条解曲线。因为多了λ这个维度，即使原潮流雅可比接近奇异，增广系统依然可以继续求解，这就是CPF能“越过”普通潮流发散点的根本原因。
1.2 预测校正框架：从基态解到P-V曲线的完整闭环
CPF的求解框架可以概括成八个字：预测、校正、步长、循环。
已知当前运行点(x_i, λ_i)，先用切线法预测下一个点的近似位置，得到初值(x_i^pred, λ_i^pred)；然后用牛顿法在这个初值附近迭代，把它拉回到精确曲线上，得到(x_{i+1}, λ_{i+1})；接着根据收敛情况调整步长，继续下一轮预测校正。整个循环从λ=0的基态潮流解出发，一直推进到越过电压崩溃点为止。
这里最关键的一步是“补方程”。扩展潮流方程F(x,λ)=0的方程个数是m，但未知量是x加λ共m+1个，比方程多一个自由度。必须额外补一个参数化方程G(x,λ)=0，把系统变成m+1个方程、m+1个未知数。这个参数化方程选得好不好，直接决定整条曲线能不能顺利追踪下来，后面我会展开讲。
预测用切线方向，校正用牛顿法，两者配合起来，每步走的其实是P-V曲线切线上的一小段，再修正回曲线。这有点像用直线段去逼近一条弧线，只要步长控制得当，逼近误差就能控制在可接受范围内。
1.3 为什么选IEEE-14和MATLAB这个组合
IEEE-14节点系统是IEEE标准测试系统家族里很经典的一个算例。它含有14条母线、5台发电机、20条支路，节点类型覆盖了平衡节点、PV节点和PQ节点，还带三绕组变压器和并联支路。规模比IEEE-9大，能体现算法对较复杂系统的适应性；比IEEE-30、IEEE-118小很多，迭代一次只要几毫秒，方便反复调参和验证。
从算法验证的角度看，IEEE-14有大量公开的文献结果可以对比。你在P-V曲线上画的某条曲线是否合理、λ_max落在什么范围，都能找到参考值，这对新手调试程序特别重要。如果一上来就直接跑IEEE-300节点，数据又长、调参又慢，出了问题根本不知道是自己程序写错还是系统太复杂导致的。
MATLAB的选型理由就更直接了：矩阵运算是强项，CPF里每一步都要解线性方程组、做SVD分解；画图工具很成熟，P-V曲线几行代码就能出图；脚本语言调试方便，断点随便打，变量区域直接看矩阵内容。对于这种科研验证型任务，MATLAB确实比C++或Python更省时间。
2. 数据准备：把IEEE-14系统读入MATLAB
2.1 标准数据格式与手动录入矩阵
IEEE-14的数据通常有三种来源：IEEE官网公开的CDF文本文件、MATPOWER库里的case14.m、以及老教材附录里的参数表。无论哪种来源，核心信息都是三张表：母线负荷表、发电机参数表、支路阻抗表。
我在教学场景下习惯先用矩阵手动录入，让初学者彻底搞清楚每一列的含义。以母线负荷表为例，最常用的格式是每行对应一条母线：第1列母线编号，第2列节点类型（1表示PQ，2表示PV，3表示平衡节点），第3列有功负荷Pd，第4列无功负荷Qd，第5列电压幅值初值，第6列电压相角初值。IEEE-14的完整母线数据可以写成下面这样：
% IEEE-14 母线负荷数据（基准容量 100 MVA） % [母线编号 类型 P负荷(MW) Q负荷(Mvar) 电压初值(p.u.) 相角初值(deg)] bus = [ 1 3 0.0 0.0 1.060 0; 2 2 21.7 12.7 1.045 0; 3 2 94.2 19.0 1.010 0; 4 1 47.8 -3.9 1.019 0; 5 1 7.6 1.6 1.020 0; 6 2 11.2 7.5 1.070 0; 7 1 0.0 0.0 1.062 0; 8 2 0.0 0.0 1.090 0; 9 1 29.5 16.6 1.056 0; 10 1 9.0 5.8 1.051 0; 11 1 3.5 1.8 1.057 0; 12 1 6.1 1.6 1.055 0; 13 1 13.5 5.8 1.050 0; 14 1 14.9 5.0 1.036 0; ];注意第4号母线的无功负荷是-3.9 Mvar，这是数据源给的原始值，表示该节点实际上是容性负荷。遇到这种带负号的数不要掉以轻心，也不要当成笔误删掉，IEEE-14原始数据里就是这样的。
发电机参数表里需要关注的列是：所在母线编号、有功出力、无功出力、无功出力上限Qmax、无功出力下限Qmin、端电压设定值。IEEE-14的五台发电机分布在5个节点上，写成矩阵大致如下：
% [母线编号 有功(MW) 无功(Mvar) Qmax(Mvar) Qmin(Mvar) 电压设定(p.u.)] gen = [ 1 232.4 -16.9 10 0 1.060; 2 40.0 42.4 50 -40 1.045; 3 0.0 23.4 40 0 1.010; 6 0.0 12.2 24 -6 1.070; 8 0.0 17.4 24 -6 1.090; ];这里提醒一句：不同渠道拿到的IEEE-14数据，Qmax和Qmin可能有细微差别，因为部分文献为了研究场景做过调整。你以自己下载的数据源为准即可，关键是程序里别写死。
支路表格式是每条支路一行，包含首端母线、末端母线、电阻R、电抗X、对地电纳B、变压器变比tap。IEEE-14有20条支路，包括普通输电线路和变压器支路，非变压器支路的tap填0即可。由于数据行数较多，我建议直接从标准数据文件读入或者复制MATPOWER的case14定义，这里先给出格式和数据片段：
% [首端 末端 R(p.u.) X(p.u.) B(p.u.) tap] branch = [ 1 2 0.01938 0.05917 0.0528 0; 1 5 0.05403 0.22304 0.0492 0; 2 3 0.04699 0.19797 0.0438 0; 2 4 0.05811 0.17632 0.0340 0; 2 5 0.05695 0.17388 0.0346 0; 3 4 0.06701 0.17103 0.0128 0; 4 5 0.01335 0.04211 0.0000 0; 4 7 0.00000 0.20912 0.0000 0.978; 4 9 0.00000 0.55618 0.0000 0.969; 5 6 0.00000 0.25202 0.0000 0.932; 6 11 0.09498 0.19890 0.0000 0; 6 12 0.12291 0.25581 0.0000 0; 6 13 0.06615 0.13027 0.0000 0; 7 8 0.00000 0.17615 0.0000 0; 7 9 0.00000 0.11001 0.0000 0; 9 10 0.03181 0.08450 0.0000 0; 9 14 0.12711 0.27038 0.0000 0; 10 11 0.08205 0.19207 0.0000 0; 12 13 0.22092 0.19988 0.0000 0; 13 14 0.17093 0.34802 0.0000 0; ];支路数据最容易被忽略的是变压器变比。如果忘记处理tap，基态潮流解会直接偏掉，你能从电压幅值明显异常上看出来。IEEE-14里变比不在1.0的变压器支路有好几条，处理时必须把tap映射到导纳矩阵里。
2.2 构建导纳矩阵Ybus的完整代码
有了bus、gen、branch三张表，下一步就是构建节点导纳矩阵Ybus。这个矩阵的物理意义是节点电压和注入电流的关系，潮流计算的所有关键运算都要用到它。构建规则并不复杂：先把每条支路的串联导纳算出来，y = 1 / (R + jX)；线路的对地电纳B平均分到首端和末端；变压器支路按变比折算到两侧。
下面这个函数是我反复精简过的一个版本，对IEEE-14这种20条支路的小系统完全够用：
function Y = build_Ybus(bus, branch) nb = size(bus, 1); Y = zeros(nb, nb); nl = size(branch, 1); for k = 1:nl f = branch(k, 1); t = branch(k, 2); R = branch(k, 3); X = branch(k, 4); B = branch(k, 5); tap = branch(k, 6); if tap == 0 tap = 1; end y = 1 / (R + 1i * X); % 串联导纳分别加到首端和末端 Y(f, f) = Y(f, f) + y + 1i * B / 2; Y(t, t) = Y(t, t) + y / tap^2 + 1i * B / 2; % 互导纳 Y(f, t) = Y(f, t) - y / tap; Y(t, f) = Y(t, f) - y / tap; end end注意变压器支路的R通常是0，X往往比较大，这时y = 1/(0 + jX) = -j/X，表示纯感性支路，这在物理上是合理的。构建完成后可以检查一下Ybus矩阵：对角线元素是各节点关联支路导纳之和，非对角线元素是对应支路互导纳的负值，对称性要满足Y(i,j) = Y(j,i)。
2.3 用MATPOWER加载case14做交叉验证
如果你手头装了MATPOWER工具包，加载IEEE-14系统只需要一行代码：mpc = loadcase('case14')。它返回一个结构体，里面包含bus、gen、branch、baseMVA四个关键字段，正好对应我们手动录入的三张表和基准容量。用MATPOWER的runpf(mpc)可以直接算出基态潮流，和自己手写的牛顿法结果对比电压幅值误差。
我强烈建议读者在跑CPF之前，先做一次这种交叉验证。原因是：基态潮流是一切的基础，如果基态都算不对，后面连续功率流算出来的P-V曲线一定有问题。交叉验证找到的问题通常集中在两类：一类是数据录入错误，比如某条支路的电抗符号写反了；另一类是单位错误，比如忘了把MW除以基准容量换成标幺值。IEEE-14的基准容量是100 MVA，所有功率都要除以100，这一点很容易漏。
3. 连续功率流算法落地：从预测到校正的MATLAB实现
3.1 扩展潮流方程与负荷增长方案定义
连续功率流的第一步，是给潮流方程加一个负荷增长因子λ。最常用的增长方案是比例增长，也就是所有负荷的有功和无功都按基态值的比例增加。用公式表示就是：
Pd_i(λ) = Pd_i_0 + λ × dPd_i
Qd_i(λ) = Qd_i_0 + λ × dQd_i
如果取dPd_i = Pd_i_0、dQd_i = Qd_i_0，那么λ=0对应基态负荷，λ=1对应基态负荷的两倍。这种定义方式直观，文献里也最常用。
确定负荷增长方向后，把潮流方程写成扩展形式F(x, λ)=0。x由非平衡节点的相角θ_pvpq和PQ节点的电压幅值V_pq组成，方程个数m = (n-1) + n_pq。对应的功率残差方程是：
P方程：对除了平衡节点外的所有节点，Pg_i - Pd_i(λ) - Pcalc_i = 0
Q方程：对PQ节点，Qg_i - Qd_i(λ) - Qcalc_i = 0
这里Pcalc和Qcalc是由当前电压算出的节点注入功率，也就是：
Pcalc_i = V_i × Σ(V_j × (G_ij×cosθ_ij + B_ij×sinθ_ij))
Qcalc_i = V_i × Σ(V_j × (G_ij×sinθ_ij - B_ij×cosθ_ij))
在增长过程中，我采用一个简化方案：只有负荷按比例增长，所有发电机的有功出力保持基态值不变，功率不平衡量由平衡节点自动承担。这样dF/dλ的表达式会非常干净：对每个P方程，偏导是-dPd_i；对每个Q方程，偏导是-dQd_i。程序里组装的dFdλ向量就是把这些值按方程顺序排成一列。
3.2 预测步：扩展雅可比的SVD切线
预测步的目标是求当前运行点处的切线方向。把潮流雅可比J和dFdλ并排拼起来，得到扩展雅可比矩阵J_ext = [J, dFdλ]，它是一个m行、m+1列的矩阵。切线方向t满足J_ext × t = 0，也就是t属于J_ext的零空间。由于J_ext是短胖矩阵，零空间至少有一维，这个零空间方向就对应P-V曲线在当前的切线。
求解零空间最稳的办法是奇异值分解SVD。MATLAB里svd函数一步到位，取最小奇异值对应的右奇异向量就是零空间方向：
[~, ~, Vmat] = svd(J_ext); t = Vmat(:, end); % 取最小奇异值对应的右奇异向量 if t(end) < 0 t = -t; % 确保负荷增长方向为正 end代码最后一步的符号判断很重要。SVD回来的向量方向是任意的，如果不强制规定，可能某一步切线方向突然反了，P-V曲线画出来就会来回折叠。判定规则是让λ分量大于0，也就是从基态出发沿负荷增长方向推进。
拿到切线方向t后，预测点就是当前点沿t方向走一个步长ds：
x_pred = x + ds * t(1:end-1); lambda_pred = lambda + ds * t(end);对IEEE-14系统，初始步长ds取0.05到0.1之间一般没问题。后面会讲自适应步长的调整策略，这里先用固定值跑通主循环再说。
3.3 校正步：局部参数化和弧长参数化选型
校正步要做的事情是：从预测点(x_pred, λ_pred)出发，用牛顿法把它拉回精确解曲线。因为扩展方程F(x,λ)=0只有m个方程、m+1个未知数，必须选一个变量固定下来。局部参数化的做法是：挑选预测切线的m+1个分量中变化幅度最大的那个变量，比如第k个状态变量x_k，强制它在整个校正过程中保持不变，等于预测值x_k_pred。
为什么选变化最大的分量？原因很简单：在P-V曲线中，这个变量沿曲线方向变化最快，把它固定下来，得到的校正方程数值条件最好，牛顿迭代不容易发散。如果随便固定一个变化很小的变量，可能出现校正方程病态，导致迭代失败。
校正迭代的方程写成扩展形式：
H = [F(x, λ); x_k - x_k_pred] = 0
它对应的雅可比矩阵是：
JH = [J, dFdλ; 0_{1×m}, 0; 但第k个位置设为1]
用MATLAB组装：
JH = [Jac, dFdlam; zeros(1, m+1)]; JH(end, k) = 1; dx = -JH \ H; x = x + dx(1:end-1); lambda = lambda + dx(end);牛顿迭代一直重复，直到H的无穷范数小于某个容差，比如1e-10。对IEEE-14来说，一般3到6次迭代就能收敛，如果超过10次还没收敛，多半是步长太大或者参数化变量选得不好。
另一种更稳健但稍微复杂一点的选择是弧长参数化，校正方程变成：
(x - x_pred)^T (x - x_pred) + (λ - λ_pred)^2 = ds^2
这个方程的意思是校正后的点必须落在以预测点为圆心、半径为ds的超球面上，对曲线曲率变化的适应能力更强。我个人的建议是：教学代码先用局部参数化，因为简单直观、代码量少；等理解了原理，再升级到弧长参数化做工程应用。两种方式的区别可以简单类比为：局部参数化是沿着某个坐标轴方向修正，弧长参数化是沿着弧线方向修正。
3.4 步长自适应与P-V曲线绘制
固定步长跑CPF最大的问题是效率与鲁棒性难以平衡。步长太大，校正步容易发散；步长太小，整条曲线推进很慢，浪费算力。自适应的基本思路是观察校正步的牛顿迭代次数：迭代次数少说明曲线平滑，可以放大步长；迭代次数多说明曲线弯曲剧烈，需要缩小步长。
我在程序里用的是这个简单策略：
if nit <= 4 ds = min(ds * 1.2, 0.2); elseif nit >= 8 ds = max(ds * 0.5, 0.005); end其中nit是牛顿校正步用掉的迭代次数。上限0.2和下限0.005是经验值，对于IEEE-14系统配合局部参数化，这个区间足够用。如果换用弧长参数化，步长上限可以适当放大到0.5左右，因为弧长校正的鲁棒性更好。
主循环全程把每一轮的λ和关键节点的电压幅值存下来，最后就可以画P-V曲线。电压崩溃通常先在负荷最重的区域暴露出来，IEEE-14系统里母线9、母线14往往是电压下降最明显的节点。绘图代码非常简单：
figure; plot(lambda_trace, V_trace_bus14, 'b-o', 'LineWidth', 1.5); xlabel('负荷增长因子 \lambda'); ylabel('母线14电压幅值 (p.u.)'); grid on; title('IEEE-14系统连续功率流P-V曲线');如果想把多条曲线叠在一张图里，可以同时记录母线9、10、14的电压，用不同颜色区分。曲线尾部电压急剧下降、越走越陡的那一段，就是系统快要失稳的部分。
4. 常见问题与避坑实录
4.1 基态潮流不收敛怎么排查
CPF跑不动，十有八九是基态潮流就出了问题。最常见的有三类：第一类是单位没换算，负荷功率没有除以基准容量100 MVA，导致所有标幺值都大了100倍，牛顿法从初值开始就找不到解；第二类是支路数据录入错误，特别是某条支路的电抗X填成了电阻R的位置，导纳矩阵完全变形；第三类是节点类型设置错误，比如把PQ节点误设为PV节点，导致方程数和未知数不匹配。
最有效的排查方法就是我前面说的交叉验证法。先打印基态潮流求出的节点电压幅值，和IEEE-14的标准解对比。如果所有电压都在0.95到1.1的范围附近，说明数据基本正确；如果某些电压变成1.5甚至负值，回到数据表一张张检查吧。另外，牛顿法迭代初始值也是有讲究的，简单起见所有非平衡节点相角初值设0，所有PQ节点电压幅值初值设1，这对IEEE-14足够用了。
4.2 校正迭代发散的处理
校正步发散是CPF新手遇到最多的问题。表现是：预测点算出来很正常，但牛顿法迭代越走越偏，最后H的范数爆炸。这时候首先检查步长ds，把ds从0.1降到0.05甚至0.02。小步长会让预测点离真实曲线更近，牛顿法从更好的初值出发，收敛概率大幅提高。
如果缩小步长还没用，就要怀疑参数化变量选得不对。按我前面的规则，应该选取切线方向中变化最大的分量，但有时候程序里max(abs(t))会选到λ分量本身，导致校正方程固定的是λ而不是某个状态变量，这在接近分岔点时特别容易出问题。一个简单有效的改法是：只从状态变量即x的分量里选，不把λ纳入候选。
还有一种情况是数值差分的步长选择问题。如果用数值雅可比，差分步长h取1e-6是个比较折中的值。h太大，差分会引入较大截断误差；h太小，浮点舍入误差会主导结果。如果发现雅可比矩阵出现奇异的对角线，可以尝试把h改为1e-7到1e-5之间多试几次。
4.3 PV转PQ、无功越限的坑
IEEE-14的发电机无功出力都有上限和下限，比如母线2的发电机关无功上限50 Mvar、下限-40 Mvar。在CPF过程中，随着负荷增长，部分发电机的无功出力会持续上升，一旦触及Qmax，这台发电机就失去了电压调节能力。此时必须把节点类型从PV转为PQ，在该节点上固定无功出力等于Qmax，把电压幅值从固定值变为待求变量。
这个转换如果不做，程序会出现一个非常明显且迷惑人的现象：雅可比矩阵子块的行列式突然接近零，牛顿迭代发散，但数据看起来没什么问题。转换的实现方式不复杂：维护一个动态节点类型向量，每轮迭代后检查所有PV节点的无功出力，如果越限就更新类型标签和方程结构。注意转换是单向的，在静态电压稳定分析中，一般不允许从PQ转回PV，因为发电机一旦达到无功极限，电压升高也不能恢复调节能力。
4.4 怎么判定已经越过电压崩溃点
CPF循环不能无限跑下去。判定越过崩溃点的标准，是观察预测步求出的切线向量t的最后一个分量，也就是λ方向的切线分量。在P-V曲线的上半支，负荷增长因子随曲线推进而增加，所以t的λ分量始终为正。当曲线越过鞍节点分岔点进入下半支后，λ开始回落，此时t的λ分量会变成负值。
程序里只需要加一条判断：
if t(end) < 0 break; end这样循环在越过崩溃点后立刻停止，并且上一个循环点对应的λ值就是系统的静态电压稳定裕度λ_max。需要强调的是，因为步长是离散的，程序输出的λ_max比真实值会略偏大或略偏小，如果你要做精细的裕度计算，可以在检测到t(end)<0后，把步长缩小，再从上一个点重新推进几次，让λ_max的估计更精确。这个操作有点像二分法逼近极限值，对精度要求高的场景很值得加。
另外，在P-V曲线下半支，电压已经失稳，实际物理系统不可能运行在那里，所以画图的时候一般只画上半支就足够了。如果你看到文献里的P-V曲线画出了完整的“鼻子形”曲线，那是为了展示算法能穿过分岔点，不代表系统可以运行到那个区域。
个人实操体会
最后分享一点我在跑这个算例时最深的感受：调试CPF的顺序极其重要，先跑通基态潮流，再上线CPF，不要试图一步到位。先用固定小步长跑一条粗糙的P-V曲线，确认曲线的形状和大致范围合理，再去调自适应步长、优化参数化方式。我见过不少同学一上来就追求大步长和高级算法，结果连基态数据都没核对清楚，浪费了大把时间在排查低级错误上。
另外，代码写出来之后，一定要做一次结果合理性检查，比如基态的电压幅值、平衡节点的出力是否在合理范围、曲线最大λ值是否落在已有文献的区间里。这个习惯能帮你提前筛掉很多数据问题。后续如果想把项目扩展下去，可以尝试改负荷增长方向看不同场景下的裕度变化，或者接入柔性输电设备研究对电压稳定性的改善效果，主循环逻辑基本不用动。
