我最早认真去写IEEE 33节点潮流计算程序其实不是为了发论文而是为一个35kV实际电网的电压越限问题发愁。翻了几篇文献公式推了一遍又一遍最后发现真正有用的还是先找一个小测试系统把程序跑通再往实际数据上迁移。IEEE 33节点配电网测试系统就是这种情况下绕不开的经典算例而牛顿-拉夫逊法NR法又是绝大多数教科书默认的潮流算法。这篇我就把用Matlab实现IEEE 33节点潮流计算NR法的完整思路、代码细节、坑点以及一个典型的应用场景都整理出来给正在做配电网研究、电力课程设计或者刚入门写潮流程序的朋友做个参考。这套内容适合谁一是电气工程专业学生拿它做毕业设计或者课程实验二是刚进电网相关单位、需要处理配电网运行分析的工程师三是想快速验证自己算法想法、不想从零推公式的研究人员。NR法是理解电力系统分析的核心算法之一能手动实现一遍后面再看BPA、PSASP这些商业软件的潮流模块心里会踏实很多。1. 先搞懂IEEE 33节点这套测试系统到底长什么样1.1 节点支路结构和基础数据IEEE 33节点系统是一条典型的放射状配电网基准电压12.66kV额定总负荷大约5083.5kW无功负荷2547.3kvar。系统一共33个节点、32条常规支路另外还有5条联络开关支路常开状态合上这些联络开关可以形成不同方式的闭环网络所以它既能模拟放射状运行也能模拟环网运行。在这33个节点里节点1是电源点也就是平衡节点一般把它当成变电站出口母线电压恒定在1.0 p.u.相角为0。其余32个节点全部可以看作PQ节点也就是有功负荷和无功负荷给定的节点。当然如果需要研究分布式电源接入也可以在某个PQ节点上增加发电功率让它在一定范围内按PQ节点处理或者改成PV节点这就是后话了。支路数据的形式一般是首端节点、末端节点、线路电阻欧姆、线路电抗欧姆、是否正常运行。IEEE 33节点系统常见参数是32条正常支路加5条联络支路其中联络支路包括8-21、9-15、12-22、18-33、25-29这几条。这些公开数据在大量文献里都能查到我下面会给出算法实现时数据矩阵的组织方式。1.2 为什么这个算例在配电网研究里这么流行IEEE 33节点系统受欢迎的一个关键原因是规模不大但特征丰富。33个节点、37条支路网络既有主干馈线又有分支馈线支路长度和负荷分布差异很大能模拟出馈线末端电压偏低、局部支路重载这些配电网典型问题。相比IEEE 123节点那种动辄几百条支路的大系统33节点跑起来很快调试也方便。更关键的一点是它的计算结果在大量文献里都有现成参照。你写完程序后把电压分布、网损跟经典文献的结果对一下就能立刻知道程序写没写错。对刚入门的人来说这种“有标准答案”的算例是最好的调试工具比拿实际电网数据上来就调要友好得多。2. NR法的数学底子别被雅可比矩阵吓到2.1 潮流计算本质上是解非线性方程组潮流计算要回答的问题很简单已知各节点的负荷和发电求出每个节点的电压幅值和相角然后算出线路里的功率和损耗。这个问题的核心是一组节点功率平衡方程。对于任意节点I注入的有功功率Pi和无功功率Qi必须满足Pi Ui * sum( Uj * (Gijcosθij Bijsinθij) )Qi Ui * sum( Uj * (Gijsinθij - Bijcosθij) )这里面Ui是节点I的电压幅值θij是节点I和节点J之间的相角差Gij和Bij是节点导纳矩阵中对应元素的实部和虚部。需要注意的是Pi和Qi是净注入功率对于纯负荷节点负荷是吸收功率所以净注入功率是负值有发电机或者分布式电源时再在净注入功率里加上发电量。这个方程组的特点是非线性因为电压乘积和三角函数都在里面。无法直接求解只能用迭代方法逼近NR法就是求解这类非线性方程组最经典的手段。2.2 NR法迭代过程的直观理解NR法的思路可以拿“找函数零点”来类比。如果要求f(x)0先猜一个x0在x0处把f展开成泰勒级数保留线性部分得到修正量再更新x反复迭代直到f足够接近0。潮流问题里的“函数f”是功率偏差量也就是按当前电压计算出来的功率和给定功率之间的差。迭代变量则是节点电压的相角修正量和幅值修正量。每次迭代要做四件事根据当前电压计算各节点功率偏差ΔP、ΔQ形成雅可比矩阵J它由有功对相角、有功对电压、无功对相角、无功对电压四块导数组成解线性方程组 J * ΔX -ΔF更新电压相角和幅值进入下一轮迭代。等功率偏差量小于设定阈值比如10^-8就认为潮流收敛。2.3 配电网计算里为什么选NR法而不是别的方法配电网潮流算法其实很多前推回推法在辐射状网络里非常快高斯-赛德尔法历史最悠久但这些方法在通用性和扩展性上各有局限。比如前推回推法处理多电源环网就很麻烦需要额外改造高斯-赛德尔法收敛速度慢处理PV节点也不方便。NR法的主要优势是收敛速度快二次收敛特性意味着收敛阶段误差按平方速度减小通常迭代5到8次就达到10^-8精度。另外它的通用性最好PV节点、环网、多平衡节点、变压器变比这些情况都能在同一个框架下处理改起来比较直接。有人会提到PQ分解法它在输电网中用得非常成功但有一个隐含前提线路电阻远小于电抗也就是低R/X比。配电网恰恰相反线路R/X比较高有些支路电抗和电阻数量级接近PQ分解法的忽略条件不成立收敛性会受影响。所以我个人在做配电网潮流计算时第一选择还是完整NR法稳通用不怕特殊网络结构。3. Matlab代码逐段拆解从数据准备到迭代收敛3.1 支路和负荷数据怎么组织用Matlab做IEEE 33节点算例第一步是把线路参数和负荷参数放进矩阵。常见做法是用行向量表示每条支路每行是[首端节点编号末端节点编号电阻欧姆电抗欧姆开关状态]1表示闭合0表示断开。我习惯先把节点编号原样保留然后在计算导纳矩阵时再转换成矩阵索引因为Matlab数组索引从1开始IEEE 33节点也是从1编号所以这里刚好能对应上。支路数据片段大致如下% 支路数据序号, 首端, 末端, R(ohm), X(ohm), 开关状态 branch [ 1 1 2 0.0922 0.0470 1; 2 2 3 0.4930 0.2511 1; 3 3 4 0.3660 0.1864 1; 4 4 5 0.3811 0.1941 1; 5 5 6 0.8190 0.7070 1; 6 6 7 0.1872 0.6188 1; 7 7 8 0.7114 1.2351 1; 8 8 9 1.0300 0.7400 1; 9 9 10 1.0440 0.7400 1; % 中间支路省略按经典IEEE 33节点数据补全 ];负荷数据则用一个33x3矩阵每行是[节点编号有功功率(kW)无功功率(kvar)]。节点1作为平衡节点没有负荷其他节点按标准数据填入% 负荷数据节点编号, P(kW), Q(kvar) load_data [ 1 0 0; 2 100 60; 3 90 40; 4 120 80; 5 60 30; 6 60 20; 7 200 100; 8 200 100; % 剩余节点同理 ];这里提醒一句网上流传的数据很多不同版本在个别负荷数值上会有细微差别不影响算法验证但如果你想跟某篇文献的结果严格对比一定要用该文献附录给出的同一份数据否则差一点电压就是另一个样子。3.2 基准值选择和标幺化处理IEEE 33节点算例中线路参数给的是欧姆负荷给的是千瓦和千乏而潮流方程里电压、功率、阻抗需要按统一基准做标幺化。我常用的基准是三相基准功率Sb1MVA基准电压Vb12.66kV。这样基准阻抗Zb Vb^2 / Sb 12.66^2 / 1 ≈ 160.28欧姆基准电流约45.6A基准导纳约0.00624S。标幺化时把电阻除以Zb得到电阻标幺值把负荷功率除以Sb得到功率标幺值。电压初值直接设为1.0 p.u.。注意有些文献为了数值计算方便会把基准功率取100kVA甚至10kVA只要全程序统一就行。最忌讳的是阻抗用了欧姆、功率用了兆瓦、电压用了伏特量纲一乱雅可比矩阵里就全是莫名其妙的数字。3.3 节点导纳矩阵的构建导纳矩阵是潮流计算的核心数据结构。对角元素是自导纳等于与该节点相连支路导纳之和非对角元素是互导纳等于两节点之间支路导纳的负值。支路导纳y 1/(R jX)。用Matlab实现时可以先把所有支路导纳算出来然后循环累加nbus 33; Ybus zeros(nbus, nbus); for k 1:size(branch,1) n1 branch(k,2); n2 branch(k,3); R branch(k,4); X branch(k,5); status branch(k,6); if status 0 continue; % 联络开关断开不计入导纳 end y 1 / (R 1i*X); Ybus(n1,n1) Ybus(n1,n1) y; Ybus(n2,n2) Ybus(n2,n2) y; Ybus(n1,n2) Ybus(n1,n2) - y; Ybus(n2,n1) Ybus(n2,n1) - y; end这里R、X如果还没标幺化需要在前面先除以Zb或者直接用阻抗有名值算导纳也是可以的但后面所有功率和电压必须保持一致。我建议一开始就在构建导纳矩阵之前完成标幺化这样后续更清晰。3.4 NR法主迭代段完整实现下面这段是我批量调试后比较稳定的一个版本。它把NR法分成三块计算功率偏差、构建雅可比矩阵、解修正方程并更新电压。% 参数初始化 baseMVA 1; baseKV 12.66; Zbase baseKV^2 / baseMVA; % 节点类型标志1为平衡节点2为PQ节点3为PV节点 % 这里节点1为平衡节点其余为PQ节点 type ones(nbus,1) * 2; type(1) 1; % 电压初值平启动 V ones(nbus,1); theta zeros(nbus,1); % 给定净注入功率标幺值 Pgiven zeros(nbus,1); Qgiven zeros(nbus,1); for k 2:nbus Pgiven(k) -load_data(k,2) / (1000*baseMVA); Qgiven(k) -load_data(k,3) / (1000*baseMVA); end % 分解导纳矩阵 G real(Ybus); B imag(Ybus); % 迭代参数 tol 1e-8; maxiter 20; for iter 1:maxiter % 计算功率偏差 Pcalc zeros(nbus,1); Qcalc zeros(nbus,1); for i 1:nbus for j 1:nbus Pcalc(i) Pcalc(i) V(i)*V(j)*(G(i,j)*cos(theta(i)-theta(j)) B(i,j)*sin(theta(i)-theta(j))); Qcalc(i) Qcalc(i) V(i)*V(j)*(G(i,j)*sin(theta(i)-theta(j)) - B(i,j)*cos(theta(i)-theta(j))); end end dP Pgiven - Pcalc; dQ Qgiven - Qcalc; % 平衡节点的偏差不参与迭代 dP(1) 0; dQ(1) 0; if max(abs([dP; dQ])) tol break; end % 构建雅可比矩阵 p 1; % PQ节点索引 pq_idx find(type 2); % 这里就是所有非平衡节点 n length(pq_idx); J1 zeros(n,n); % dP/dtheta J2 zeros(n,n); % dP/dV J3 zeros(n,n); % dQ/dtheta J4 zeros(n,n); % dQ/dV for ii 1:n i pq_idx(ii); for jj 1:n j pq_idx(jj); if i j J1(ii,jj) -Qcalc(i) - B(i,i)*V(i)^2; J2(ii,jj) Pcalc(i)/V(i) G(i,i)*V(i); J3(ii,jj) Pcalc(i) - G(i,i)*V(i)^2; J4(ii,jj) Qcalc(i)/V(i) - B(i,i)*V(i); else J1(ii,jj) V(i)*V(j)*(G(i,j)*sin(theta(i)-theta(j)) - B(i,j)*cos(theta(i)-theta(j))); J2(ii,jj) V(i)*(G(i,j)*cos(theta(i)-theta(j)) B(i,j)*sin(theta(i)-theta(j))); J3(ii,jj) -V(i)*V(j)*(G(i,j)*cos(theta(i)-theta(j)) B(i,j)*sin(theta(i)-theta(j))); J4(ii,jj) V(i)*(G(i,j)*sin(theta(i)-theta(j)) - B(i,j)*cos(theta(i)-theta(j))); end end end % 组装修正方程若只有PQ节点J2和J4都要保留 J [J1 J2; J3 J4]; F [dP(pq_idx); dQ(pq_idx)]; % 求解修正量 dX J \ F; dtheta zeros(nbus,1); dV zeros(nbus,1); dtheta(pq_idx) dX(1:n); dV(pq_idx) dX(n1:end); % 更新状态量 theta theta dtheta; V V dV; end这一段需要解释几个容易出错的地方。雅可比矩阵中的对角元素和非对角元素公式不一样关键在于对三角函数求导时只有“本节点”求导才会出现含Qcalc、Pcalc的项。非对角元素则全部来自交叉项的求导。很多新手栽跟头就是对角元素里漏掉了自导纳Gii、Bii相关的项。其次修正方程里“对V求导”的子块在代数上往往和实际电压修正量ΔV配合而不是和ΔV/V配合。两个方案在公式上相差一个对角缩放最终结果是一样的但混用的话雅可比矩阵和右端项就不匹配程序表现就是迭代震荡甚至发散。我上面用的是经典的ΔV形式和公式一一对应。另外求解修正量时尽量用J\F而不是inv(J)*F。33节点系统两者差别不大但以后换到几百上千节点系统inv(J)会慢得让人怀疑人生而且数值稳定性也不如直接求解线性方程组。3.5 收敛判据和防发散处理收敛判据一般用功率偏差的无穷范数也就是dP和dQ中绝对值最大的那个值小于阈值就算收敛。我习惯设10^-8这个精度对工程计算绰绰有余。防发散的手段主要有三种限制迭代次数、增加阻尼因子、放松初值。NR法在平启动初值下通常不会发散但如果系统重载或者含PV节点且无功越限可能会出现迭代步长过大导致震荡。这时候可以加一个阻尼因子alpha更新时取dV alpha * dXalpha取值0.8到0.95牺牲一点收敛速度换取稳定。此外可以设一个“迭代未收敛”的提示。我自己在写程序时会在maxiter之外再加一行代码if iter maxiter warning(NR迭代未收敛建议检查导纳矩阵或初值设置); end别小看这个检查。在批量计算里如果某个场景连续多轮不收敛没有这个提示后面所有结果都会白白跑一遍浪费时间还容易得出错误结论。4. 应用实例在33节点系统上模拟分布式光伏接入的电压影响4.1 场景设计IEEE 33节点算例在实际研究中经常被用来做分布式电源接入分析。这里我设计一个很典型的场景在节点17接入一个分布式光伏电站额定出力1MW功率因数0.95也就是大约发出1MW有功和0.33Mvar无功。节点17在原始系统里处于一条较长分支的末端附近线路压降大是电压最敏感的位置之一。接入光伏之后光伏出力会抵消一部分沿线输送功率线路电流减小压降也随之减小所以节点17的电压会明显抬升。问题是抬升多少会不会抬过头这正好可以用NR法算出来。4.2 修改代码模拟新场景在原有代码的负荷数据基础上把节点17的净注入功率改为Pgiven(17) Pgiven(17) 1.0; % 加1MW发电功率 Qgiven(17) Qgiven(17) 1.0*0.3287; % 加0.3287Mvar无功注意这里的符号负荷是负的发电是正的。所以在原来负负荷的基础上加上发电功率即可。如果光伏逆变器运行在单位功率因数无功可以不加但实际并网光伏通常有功率因数要求我这里是按0.95容性设计来模拟。然后重新运行NR主迭代得到新的收敛结果。对比接入前后各节点电压就能清晰看出光伏对电压的支撑作用。4.3 结果分析和工程结论接入前在不含光伏的原始33节点系统里末端节点18或节点17附近的电压通常已经降到0.94 p.u.左右末端负荷很重时甚至更低。接入1MW光伏后由于本地电源提供了大量有功馈线潮流减小电压被抬升到接近1.0 p.u.甚至更高。具体数值跟负荷水平、光伏出力、线路模型都有关系从我跑出来的趋势看末端电压能抬升0.05到0.1 p.u.这是一个很可观的幅值。由此引出的工程结论是分布式光伏对配电网电压有双向影响。光伏出力大而负荷轻的时候本地电压可能越过上限光伏出力小、负荷重的时候电压又可能偏低。这已经不只是一个潮流计算问题而是配电网运行控制里的电压协调问题。用IEEE 33节点系统和NR法你可以快速评估不同接入位置、不同出力水平对电网电压的影响这也是这套算法在实际工程里的重要用法。除了电压幅值NR法算完以后还能得到各支路功率、网损等结果。比如算一次原始系统的总网损大约在0.2p.u.左右以1MVA为基准接入光伏后网损会明显下降。这些指标对于配电网规划、运行优化都很有意义。5. 踩坑实录新手做IEEE 33节点NR法最容易翻车的几个地方5.1 数据单位搞混电压结果全是“糊”的这是我见过最多的问题没有之一。有人把支路电阻用欧姆负荷用kW电压初值用12.66kV结果导纳矩阵和功率不匹配算出来的电压要么全是1要么直接发散。建议先把所有数据标幺化明确写出Sb和Vb再进入下一步。每算一步可以打印几个关键节点的功率偏差如果初始偏差量级就很离谱那肯定是单位标幺化环节出错了。5.2 雅可比矩阵的负号写反NR法修正方程的标准形式是 J * ΔX -ΔF但很多教材在推导时会直接把雅可比矩阵写成“负的偏差量对变量的导数”也就是把负号吸收进J里面。两种写法都能用关键是代码里的J和F必须配套。如果你发现自己从第二、三次迭代开始数值反复跳动优先检查这里。另外雅可比矩阵中dQ/dV的对角元是 Qcalc(i)/V(i) - B(i,i)*V(i)不是加号。这个减号非常容易漏漏掉之后往往迭代会在五到十次时突然震荡现象比较隐蔽。5.3 节点编号和支路首末端方向IEEE 33节点数据里每条支路有明确的首末端节点。导纳矩阵是对称的首末端互换不会影响结果但功率计算时如果后面要提取支路潮流首末端方向就决定了潮流正方向。建议在数据准备阶段就统一好方向别到画支路潮流图时才发现方向全反了。5.4 稀疏矩阵和大系统扩展33节点系统用满矩阵完全没问题但当你想把它扩展到IEEE 123节点、甚至几千节点配电网模型时满矩阵的雅可比矩阵会占大量内存求解速度也明显下降。这个时候把Ybus、J都改成稀疏存储Matlab代码基本不用大改性能却能提升一个数量级。顺手提供一个小技巧Ybus sparse(Ybus); J sparse(J);这两行在节点规模大了以后性价比极高。写程序时养成用sparse的习惯从33节点就开始用后面迁移到大系统会非常轻松。5.5 常见问题速查表为了方便排查我把常见现象、可能原因和处理办法整理成一张表。现象可能原因处理办法第1次迭代就NaN导纳矩阵有零元素或数据单位不一致检查支路数据、开关状态和标幺化迭代到第5次左右震荡雅可比矩阵公式错误可能是负号写反逐项核对J1、J2、J3、J4的公式收敛但电压全部接近1负荷功率没加入或者Pgiven符号错打印Pgiven确认负荷为负、发电为正收敛速度很慢需要几十次平启动初值不合理或数据量纲不对重新检查初值从1.0∠0 开始节点编号是0开头导致索引报错某些数据表把编号从0开始统一编号1或用节点名到索引的映射表更换联络开关后不收敛环路形成后网络参数变化NR法初始值不当尝试运行网闭环打开Java后检查导纳矩阵是否包含联络支路5.6 关于Matlab环境的几句实在话不少朋友问我Matlab环境的事情尤其是刚上手时下载安装、工具箱配置这些。我的建议很简单这套NR潮流计算代码只要Matlab基础环境就能跑不需要额外工具箱连优化工具箱和Simulink都用不上。如果手边还没有顺手的Matlab环境装一个教学版或者学习版就足够了别为了一个潮流计算去折腾一堆用不到的插件。真正花时间的还是算法本身环境配置只是入场券别本末倒置。6. 程序写完之后怎么确认结果是对的跑到这一步程序输出了一组电压、一组相角怎么知道它算对了最直接的办法是对照文献结果。IEEE 33节点原始算例在大量论文里都有标准电压分布数据末端节点电压一般在0.94 p.u.附近系统总网损比例大约5%到8%之间。如果你的结果和这些量级相差太多优先回头检查数据录入和雅可比矩阵。还有一个简单有效的自检方法是功率平衡校验。潮流收敛后把全部节点的Pcalc加起来再和给定的总负荷加总发电比较差值应该小于收敛阈值。无功功率同理。这个校验不依赖外部数据只依赖你的程序内部一致性每次跑完都看一眼是个好习惯。另外我自己的习惯是在每次迭代时把最大功率偏差打印出来。收敛过程看起来应该是前面几步下降很快后面慢慢趋于平稳。如果看到偏差曲线忽大忽小来回跳基本可以断定程序有问题不用等它跑完就能提前止损。如果你想验证更多场景可以把联络开关依次闭合比如把8-21闭合网络就从纯辐射状变成带环网结构。NR法在环网下同样适用不需要改主算法这正好能帮你验证程序的通用性也能理解为什么配电自动化研究里总喜欢拿33节点做重构实验。7. 从33节点到实际电网还需要补哪些功课33节点跑通以后很多朋友会问我能不能直接把程序拿来算我手头几十条馈线的实际配电网可以但有几个改造点需要留意。实际电网的节点数多、支路类型杂不仅有架空线还有电缆、变压器、分布式电源数据往往来自GIS系统或者生产管理系统格式五花八门。你需要写数据转换接口把外部数据清洗成Matlab的branch和load矩阵。这部分工作有时候比潮流算法本身更耗时。实际电网中会有大量PV节点比如有调压能力的变电站母线、无功补偿装置NR法中需要加入PV节点处理逻辑并在迭代过程中检查无功越限越限后要动态转换成PQ节点。33节点模型里全PQ节点这一步可以不做但到了实际系统躲不开。变压器支路也要特殊处理。配电变压器分接头位置对电压影响很大有载调压变压器在线调节还涉及步长和档位。考虑到这些很多人会选择在NR法中扩展变压器变比作为状态变量或者干脆在每次潮流迭代后单独处理分接头这已经属于配电网状态估计和控制优化的范畴了。我不建议一开始就把这些复杂因素全塞进程序里。先把33节点这个纯交流辐射状网络彻底弄明白再一步步加变压器、加PV节点、加环网每一个增量都对比已知结果验证。这样虽然慢一点但踩坑时你能立刻定位到是哪个功能模块出了问题。最后再分享一个我自己调试时的习惯程序里不要在循环体里用一连串变量名区分不同场景而是把每个场景的计算结果封装成结构体比如result.V、result.Ploss、result.theta批量计算时再用cell数组存起来。看起来只是代码风格问题但当你连续跑上百个场景做分布式电源选址分析时这种组织方式会帮你省掉大量整理数据的时间。IEEE 33节点加NR法是一个很小但很完整的起点把这块地基打牢后面无论是做配电网重构、电压优化还是可靠性分析都会顺手很多。