企业官网建设流程全解析

简介：面向电力系统专业学生与从业人员，该MATLAB脚本实现了基于牛顿-拉弗森法（NR法）的IEEE33节点潮流计算，可求解配电网各节点电压幅值与相角、支路功率及网损，适用于教学演示、课程设计和科研算法验证。资源包内仅有1个m文件，压缩包约2KB，代码结构简洁、无额外依赖，便于直接阅读和二次开发。脚本完整涵盖数据输入、初始状态设置、雅可比矩阵形成、迭代修正、收敛判断与结果输出等核心步骤，并生成各节点电压和线路潮流数据，能够直观呈现NR法求解非线性方程组的全过程。目前已有547人学习下载，适合需要掌握电网潮流计算原理、以IEEE33节点为基准算例开展仿真与算法对比的读者。通过该脚本，读者还可自行修改节点参数或网络拓扑，进一步探索不同负荷、线路阻抗条件下的潮流特性。

1. 为什么潮流计算的默认考题是IEEE33节点，又为什么默认解法是NR法

做配电网分析，从IEEE33节点系统开始最省心。这个33节点配电网算例规模小、标准结果齐全，无论是算网损、看电压分布，还是做光伏接入和故障重构，它都是默认的公共试验场。潮流计算的核心任务是在给定负荷下求出网络电压和功率分布，NR法（Newton-Raphson，牛顿-拉夫逊法）用雅可比矩阵逐轮逼近非线性功率方程，在IEEE33这类小型系统上通常4~8次迭代就能收敛。但配电网线路R/X比大，NR法的收敛半径问题不能忽视，工程里经常要配合潮流计算最优因子法才能保证不翻车。下面直接围绕IEEE33 + NR法这条线，把建模、代码、收敛控制与结果校验一次讲透，适合正在做配电网规划、微电网和分布式电源接入的从业者。

2. 把IEEE33节点变成NR法能吃的输入：拓扑、参数与极坐标方程

2.1 节点-支路表是第一步：IEEE33的拓扑与两套编号

IEEE33节点系统由1个平衡根节点、32个PQ节点和32条运行支路组成，基准电压12.66kV，基准功率10MVA。公开资料里通常还列出5条联络开关支路，默认断开，所以潮流计算时网络是辐射状。全部负荷集中在节点上，总负荷3.715MW+j2.3Mvar，这个总量可以用来核对数据抄写是否正确，任何一张标准IEEE33数据集，负荷相加都应该落在这一组数字附近。

做输入前先确定编号规则。有一版常用数据把根节点编号为0，从0到32共33个节点；另一版从1号开始编号。两种编号在网络拓扑表里的支路起止节点完全不一样，拿到数据后要先把根节点、负荷节点标注出来。我一般先画一张单线图，再在程序里把索引写死，避免后面对应错位。如果只拿一串CSV就开始跑，最容易犯的错就是“看结果有点像但细看全错”。

IEEE33的32条运行支路可以按四条路径来记：主线0-1-2-...-17，共17条；从节点1伸出去的左侧分支1-18-19-20-21，共4条；从节点2伸出去的右侧分支2-22-23-24，共3条；从节点4伸出去的末端分支4-25-26-27-28-29-30-31-32，共8条。合起来正好32条，这也是为什么很多支路表看起来像四块拼起来，而不是一条直线。

这个分支结构有两个作用：一是手填支路时可以按路径逐段检查，二是后面做分布式电源接入时，可以快速判断某个节点属于哪条馈线分支。注意联络开关支路不在上表里，通常文献会额外列出5条，编号在不同资料里有差异，默认都断开。

2.2 极坐标功率方程和雅可比矩阵：NR法计算潮流的三个关键量

NR法在电网潮流里常用极坐标形式。对于节点i，注入功率的实部P_i和虚部Q_i可以写成

P_i = V_i * sum_j V_j (G_ij cosθ_ij + B_ij sinθ_ij)
Q_i = V_i * sum_j V_j (G_ij sinθ_ij - B_ij cosθ_ij)

其中θ_ij = θ_i - θ_j，G和B是节点导纳矩阵的实部和虚部。IEEE33系统里平衡节点0的电压幅值固定为1.0p.u.、相角固定为0，其余32个PQ节点的状态量是电压幅值V和相角θ，总共64个未知量。NR迭代要解的是

[ΔP; ΔQ] = J [Δθ; ΔV]

J是2×2分块雅可比矩阵，四块分别是有功对相角、有功对幅值、无功对相角、无功对幅值的偏导。写代码时按节点填充，对角和非对角的公式要严格区分：

H_ii = -Q_i - B_ii V_i^2
H_ij = V_i V_j (G_ij sinθ_ij - B_ij cosθ_ij)
N_ii = P_i/V_i + V_i G_ii
N_ij = V_i (G_ij cosθ_ij + B_ij sinθ_ij)
K_ii = P_i - G_ii V_i^2
K_ij = -V_i V_j (G_ij cosθ_ij + B_ij sinθ_ij)
L_ii = Q_i/V_i - V_i B_ii
L_ij = V_i (G_ij sinθ_ij - B_ij cosθ_ij)

这些公式就是后面代码里雅可比矩阵的填充依据，不要凭印象写，符号错一个，电压结果就会往错误方向跑。另外要注意，IEEE33原始系统全是PQ节点，没有PV节点；如果后续把某个节点改成分布式光伏并网点，就需要在雅可比矩阵里删去该节点的无功失配行，补上电压幅值约束行，这一步NR法比前推回代法灵活得多。

2.3 为什么NR法仍然值得在IEEE33上用

配电网潮流还有一种更常见的前推回代法，利用辐射状拓扑从末端回推电流、从根节点前推电压，一次往返就能更新一轮，编程简单、收敛也快。但它有个硬前提：网络必须是辐射状，很难处理PV节点和多源网络。NR法没有这个限制，任意拓扑都能处理，雅可比矩阵还天然提供灵敏度信息，后面做DG接入、N-1校核或者电压无功优化，都可以直接复用这个矩阵。

代价是NR法的雅可比矩阵每轮都要重算，33节点规模不明显，到几百节点时计算成本会上升。但对IEEE33这个量级，一次矩阵求解在毫秒级，运算时间完全不是瓶颈，这也是大家愿意拿它当教学算例的原因。真正要花心思的不是“能不能算”，而是“怎么让它在重负荷下稳定收敛”，这一点放到第4章专门说。

3. 用Python跑通33节点潮流计算：一个只用numpy的NR实现

3.1 先准备IEEE33数据：支路阻抗和节点负荷

下面是一份可复制的数据准备代码，把支路表整理成[i, j, R, X]四列，负荷表按节点索引存放P(kW)和Q(kvar)。根节点编号用0，这是最常见的一版编号规则。

import numpy as np # IEEE33节点系统，根节点编号0 # 支路: [送端, 受端, R(ohm), X(ohm)] branch = np.array([ [0, 1, 0.0922, 0.0470], [1, 2, 0.4930, 0.2511], [2, 3, 0.3660, 0.1864], [3, 4, 0.3811, 0.1941], [4, 5, 0.8190, 0.7070], [5, 6, 0.1872, 0.6188], [6, 7, 0.7114, 0.2351], [7, 8, 1.0300, 0.7400], [8, 9, 1.0440, 0.7400], [9, 10, 0.1966, 0.0650], [10, 11, 0.3744, 0.1238], [11, 12, 1.4680, 1.1550], [12, 13, 0.5416, 0.7129], [13, 14, 0.5910, 0.5260], [14, 15, 0.7463, 0.5450], [15, 16, 1.2890, 1.7210], [16, 17, 0.7320, 0.5740], [1, 18, 0.1640, 0.1565], [18, 19, 1.5042, 1.3554], [19, 20, 0.4095, 0.4784], [20, 21, 0.7089, 0.9373], [2, 22, 0.4512, 0.3083], [22, 23, 0.8980, 0.7091], [23, 24, 0.8960, 0.7011], [4, 25, 0.2030, 0.1034], [25, 26, 0.2842, 0.1447], [26, 27, 1.0590, 0.9337], [27, 28, 0.8042, 0.7006], [28, 29, 0.5075, 0.2585], [29, 30, 0.9744, 0.9630], [30, 31, 0.3105, 0.3619], [31, 32, 0.3410, 0.5302] ]) # 负荷: [P(kW), Q(kvar)]，节点0为根节点 load = np.array([ [0, 0], [100, 60], [90, 40], [120, 80], [60, 30], [60, 20], [200, 100], [200, 100], [60, 20], [60, 20], [45, 30], [60, 35], [60, 35], [120, 80], [60, 10], [60, 20], [60, 20], [90, 40], [90, 40], [90, 40], [90, 40], [90, 40], [90, 50], [420, 200], [420, 200], [60, 25], [60, 25], [60, 20], [120, 70], [200, 600], [150, 70], [210, 100], [60, 40] ])

这段代码里的支路表覆盖了主线、左侧分支、右侧分支和末端分支，一共32条运行支路；5条联络开关支路没有放进来，默认开环。负荷数组第0行是根节点，负荷为0；后面每一行对应一个PQ节点的有功和无功，单位是kW和kvar，这些值会在下一小节换算成标幺值。如果你拿到的数据是1号起编，记得先把根节点移到下标0，再填负荷数组，否则后面的雅可比矩阵全是错位。

3.2 构建节点导纳矩阵和功率函数

继续写代码，把网络参数变成标幺值，并构建节点导纳矩阵Y，然后定义功率计算函数。

Vbase = 12.66 # kV Sbase = 10.0 # MVA Zbase = Vbase**2 / Sbase # 欧姆 n = 33 # 节点导纳矩阵，单位是标幺值 Y = np.zeros((n, n), dtype=complex) for f, t, R, X in branch: y = 1.0 / ((R + 1j*X) / Zbase) Y[f, f] += y Y[t, t] += y Y[f, t] -= y Y[t, f] -= y G = Y.real B = Y.imag def calc_pq(V, theta): dtheta = theta[:, None] - theta[None, :] P = V * (G * np.cos(dtheta) + B * np.sin(dtheta)) @ V Q = V * (G * np.sin(dtheta) - B * np.cos(dtheta)) @ V return P, Q

逻辑说明：支路阻抗先除以Zbase得到标幺值，再取倒数得到支路导纳；Y矩阵对角线累加自导纳，非对角线累加互导纳，这是标准节点导纳矩阵规则。calc_pq函数用dtheta构造出所有节点相角差的矩阵，然后用矩阵乘法和电压数组完成全网络有功、无功计算，比用双重循环逐节点累加更短，也更容易发现维度错误。

参数说明：Zbase约等于16.03Ω，也就是12.66kV的平方除以10MVA。Sbase用10.0而不是100，因为IEEE33标准算例指定的就是10MVA，如果按一般输电网习惯用100MVA，所有标幺值都会缩小10倍，结果看起来会像“没收敛”。写程序时把Vbase、Sbase放在数据准备区，后续换算都是同一个常量，避免在calc_pq里临时除。

3.3 NR法主循环：失配量、雅可比矩阵和电压更新

下面是最核心的NR迭代部分，这里先给出不带阻尼的原始版本，方便对比后面加最优因子法的效果。

# 给定功率：负荷取负值 P_spec = np.zeros(n) Q_spec = np.zeros(n) P_spec[1:] = -load[1:, 0] / (Sbase * 1000.0) Q_spec[1:] = -load[1:, 1] / (Sbase * 1000.0) V = np.ones(n) # 平启动，单位p.u. theta = np.zeros(n) # 相角，单位rad tol = 1e-8 # 标幺值失配功率阈值 max_iter = 20 m = n - 1 for k in range(max_iter): P, Q = calc_pq(V, theta) dP = P_spec[1:] - P[1:] dQ = Q_spec[1:] - Q[1:] mis = np.concatenate([dP, dQ]) if np.max(np.abs(mis)) < tol: print("converged, iter=", k + 1) break J = np.zeros((2 * m, 2 * m)) dtheta = theta[:, None] - theta[None, :] for i in range(m): ni = i + 1 for j in range(m): nj = j + 1 if ni == nj: J[i, j] = -Q[ni] - B[ni, ni] * V[ni]**2 J[i, m + j] = P[ni] / V[ni] + V[ni] * G[ni, ni] J[m + i, j] = P[ni] - G[ni, ni] * V[ni]**2 J[m + i, m + j] = Q[ni] / V[ni] - V[ni] * B[ni, ni] else: s = np.sin(theta[ni] - theta[nj]) c = np.cos(theta[ni] - theta[nj]) J[i, j] = V[ni] * V[nj] * (G[ni, nj] * s - B[ni, nj] * c) J[i, m + j] = V[ni] * (G[ni, nj] * c + B[ni, nj] * s) J[m + i, j] = -V[ni] * V[nj] * (G[ni, nj] * c + B[ni, nj] * s) J[m + i, m + j] = V[ni] * (G[ni, nj] * s - B[ni, nj] * c) dx = np.linalg.solve(J, mis) theta[1:] += dx[:m] V[1:] += dx[m:] print("Vmin=", V.min(), " Vmax=", V.max()) print("P_loss(MW)=", (P[0] + P_spec.sum()) * Sbase)

这段代码的收敛判据是max(|ΔP|, |ΔQ|)小于1e-8，这是比较严格的要求；实际工程用1e-6已经足够，因为标幺值1e-6对应的有名功率只有约10W。雅可比矩阵直接用解析偏导填充，64×64规模很小，np.linalg.solve的耗时可以忽略。

如果代码正确，IEEE33原始负荷下落点电压最低值通常在0.90p.u.附近，不会低于0.85；总网损在0.2MW上下。如果你跑出来的最低电压只有0.6p.u.，大概率是负荷没除以Sbase，而不是NR法本身的问题。原始NR法在平启动下通常能收敛，但一旦修改支路参数或加重负荷，就可能出现失配量反复震荡，这正好引出下一章的最优因子法。

4. 收敛不行就上潮流计算最优因子法：原理与代码改造

4.1 为什么配电网会让NR法翻车：R/X比与初值敏感

NR法的收敛性对初值非常敏感。输电网电抗远大于电阻，雅可比矩阵条件数好，平启动就能稳稳收敛；配电网线路R/X比高，IEEE33里很多支路电阻大于电抗，矩阵对角占优性变差。这时候如果末端负荷很重，电压真实解已经偏离1.0p.u.很远，全步长修正很容易越过真实解，表现为失配量先变小再变大，电压出现负值或超过1.1p.u.。

我实际跑过不少33节点改造模型，把某条支路R放大或把末端负荷加到2倍，原始NR法就开始发散。这种发散不是代码bug，而是牛顿法收敛半径问题。遇到这种情况，第一反应不要是去改收敛阈值，那只会得到假收敛；应该给修正量加一个自适应步长，这就是潮流计算最优因子法要解决的场景。

4.2 潮流计算最优因子法是什么：给修正量乘一个标量μ

最优因子法的思想很直观。NR法的修正量Δx来自线性化方程，全步长x+Δx只对线性模型最优，非线性模型下这个步长可能太大。最优因子法把更新改成

x^{k+1} = x^k + μ Δx

其中μ取某个标量，让失配量函数

F(μ) = || f(x^k + μΔx) ||^2

尽量小。当μ=1时就是原始NR法；μ<1时相当于阻尼步长，μ>1时可以加速收敛。精确求μ需要解一个高阶多项式，工程上很少这么做，更常见的做法是用回溯线搜索：先试μ=1，如果更新后的失配量范数比更新前更大，就把μ折半，直到失配量下降或达到折半上限。

4.3 代码改造：加入最优因子后的NR主循环

把上一章的更新部分替换成带线搜索的版本，其余代码不变。

# 在每次求得dx后执行，替换 theta[1:] += dx[:m] 那两行 max_mu_try = 12 mu = 1.0 old_norm = np.linalg.norm(mis) for t in range(max_mu_try): V_can = V.copy() theta_can = theta.copy() V_can[1:] += mu * dx[m:] theta_can[1:] += mu * dx[:m] P_can, Q_can = calc_pq(V_can, theta_can) new_mis = np.concatenate([P_spec[1:] - P_can[1:], Q_spec[1:] - Q_can[1:]]) if np.linalg.norm(new_mis) < old_norm: V, theta = V_can, theta_can break mu *= 0.5 else: print("line search failed at iter", k + 1) break

说明：这段代码先尝试全步长，只有失配量范数确实下降才接受；否则把步长μ不断减半。max_mu_try=12意味着最小步长约为2^-12，足以应对绝大多数重负荷场景。如果连续多次线搜索失败或μ长期小于0.2，说明问题大概率在数据本身，而不是算法，继续压步长只会让收敛变慢。

和原始NR法相比，这个改造增加的计算量在最坏情况下是每次迭代多算12次功率方程，但对IEEE33来说，每次calc_pq都是64维矩阵运算，开销完全可以接受。收敛轮数可能会从4轮增加到6轮，但换来的是不再随机发散，这比那一点额外计算重要得多。

4.4 最优因子法的三个使用注意

第一，不要无脑把所有场景都开成强阻尼。正常负荷下原始NR法4轮就收敛，加了线搜索后如果μ一路保持1.0，不影响结果；但如果初始电压给得特别差，μ连续多轮折半，收敛速度会明显下降。这时可以改用上一轮潮流解作为初值，而不是每次都从1.0p.u.平启动。

第二，最优因子法不能修复雅可比矩阵奇异。如果节点编号断链、支路阻抗填成0，或者某个节点孤立，Ybus对应行全为0，np.linalg.solve依然会报singular matrix。线搜索只处理“方向对但步长太大”，处理不了“方向本身错误”。

第三，μ折半上限要跟收敛阈值匹配。max_mu_try=12对应最小步长约0.000244，如果这个步长下失配量仍然不降，那就不是步长问题，应该停止迭代并检查数据。盲目把上限提高到50，只会多算几十次calc_pq，对结果没有帮助。

5. IEEE33潮流计算避坑清单：编号、基准值、收敛条件与拓扑

5.1 节点编号从0还是1起编，最容易引起结果对不上

现象：跑出来的电压分布趋势对，但具体节点电压和公开参考值差一格，越到末端越乱。原因：同一份IEEE33单线图，有人把根节点标0，有人标1；从1起编的数据里“支路1-2”对应0起编的“0-1”，错位一个节点。解决：写代码前先把源数据的路径画出来，明确根节点索引；我一般直接在支路表注释里写“0是变压器低压侧出口”，再让程序打印每个节点电压，肉眼核对首末端。注意负荷数组也要同步平移，曾经有人只改了支路编号，负荷还留在1号起编的位置，结果总负荷分布完全错了。

5.2 基准值选错：12.66kV/10MVA背后的单位换算

现象：所有电压算出来只有1e-5量级，或者网损几乎为0。原因：把负荷kW直接当成标幺值，或者支路阻抗没有除以Zbase。IEEE33的基准电压是12.66kV，基准功率常用10MVA，Zbase约等于16.03Ω。负荷换算要先除以10MVA，支路阻抗要先除以Zbase。解决：在代码里保留Vbase、Sbase常量，换算写在数组前面；不要在calc_pq函数里临时处理单位。另一个常见误用是把Sbase设成100MVA，结果电压剖面看起来“偏低但不离谱”，这种错误最难排查，建议每次打印总负荷标幺值，确认它等于0.3715+j0.2300。

5.3 收敛条件：只看ΔP不看ΔQ，会得到假收敛

现象：迭代3轮就达到收敛，但Q失配量还停在1e-2量级，无功和网损明显不对。原因：收敛判据只取了有功失配，而PQ节点是P和Q两个方程同时要满足。解决：用max(|ΔP|, |ΔQ|)合并判断，且对标幺值取1e-6以下；如果还关心电压稳定性，再要求最大电压修正量小于1e-5。这一条在重负荷场景下特别重要，因为无功失配往往比有功失配收敛得慢，只看有功会把未收敛的解当成结果。

5.4 5条联络开关不加区分，33节点就变成闭环网络

现象：直接用文献里的完整33节点数据，把5条联络开关也放进Ybus，潮流结果与标准辐射状算例相差很大。原因：IEEE33节点系统原始数据包含5条常开联络支路，默认开环运行；只有做网络重构时才逐个合上。解决：支路表先只放32条运行支路，把5条联络支路单独放一个数组，需要时再追加。如果代码里Ybus是循环遍历branch构建的，那就很容易把联络开关误算进去，建议构建前打印branch.shape，确认是32行而不是37行。

5.5 雅可比矩阵奇异或迭代爆炸，先查数据再怀疑算法

现象：np.linalg.solve报singular matrix，或电压更新后出现负值。原因：常见是节点编号断链、支路R或X填成0、某个节点负荷大到让该节点电压无解，而不是NR法本身的锅。解决：先打印Ybus的行非零个数，确认每个节点都挂在网络上；再用一个把各节点负荷乘0.1的轻载算例跑通，确认算法稳定后再恢复全负荷。这样可以把数据问题和算法问题快速切开，不用靠猜。轻载算例如果也发散，那就是编号或雅可比矩阵代码问题；轻载收敛但满载发散，再去考虑最优因子法和初值。

6. 结果可信度验证：功率不平衡检查与前推回代交叉验证

6.1 三个自检指标

收敛后不要急着把V数组写进报告，先用三个量做自检。第一个是根节点功率与全系统负荷的差值，这个值就是总网损，IEEE33在原始负荷下网损大约在0.2MW量级，如果算出来是几MW，要么是收敛判据太松，要么是负荷单位错。第二个是节点电压范围，正常应在0.90~1.0p.u.之间，最低压节点通常在长分支末端，任意节点超过1.05就说明潮流解有问题。第三个是最后一次迭代的失配量最大值，打印max(|dP|, |dQ|)，应该低于1e-6。这三项可以直接写成一个check函数，每次跑完自动打印，避免人工翻日志。

6.2 用前推回代法做交叉验证

NR法实现了不代表雅可比矩阵没有填错，最稳妥的验证是用另一个算法交叉对比。IEEE33是辐射状网络，前推回代法实现起来很简单：从末端节点开始，用节点功率和当前电压回代支路电流，再从根节点向前推算各节点电压，反复迭代到电压差小于阈值。把两种算法算出的末端电压放在一张表里对比，偏差小于1e-4 p.u.就可以基本确认实现正确。我第一次调通NR法时，就是靠前推回代交叉验证发现了支路表里一条R/X数据填反，比对着参考值排查快得多。

6.3 把最优因子法推广到DG接入场景

如果后面要做分布式电源接入，把某个负荷节点改成PV节点，NR法的雅可比矩阵会多一行电压幅值约束，这时最优因子法的价值更明显：DG出力从0逐步增加到额定值，每一步用上一轮的解做初值，配合μ阻尼，一般不会发散。我现在的习惯是跑任何配电网潮流都默认带最优因子开关，先试一次μ=1，再允许线搜索介入，这比反复调初值省事，也少了很多玄学调参。希望帮到你。

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