简介:这份课程设计报告书面向电气工程及相关专业学生与电力系统分析初学者,围绕「两机五节点网络」的潮流计算展开,帮助读者掌握稳态分析中的核心求解方法。内容系统梳理牛顿拉夫逊法与PQ分解法的原理、优缺点及适用场景,并借助MATLAB完成程序编写、雅可比矩阵构建与迭代求解,对比两种算法的收敛速度与计算效率。压缩包仅含1个doc文档,约271KB,篇幅完整,涵盖摘要、原理简介、设计资料与参数、程序框图、MATLAB程序编写与运行结果、PQ法程序、总结及参考文献等章节,目录层次清晰,便于按模块检索与复盘。目前已有94人学习,适合作为课程设计参考模板、算法入门读物或潮流计算程序的对照范例,读者可据此理解电力网络建模流程、初值选取思路与结果验证方法,为后续电力系统分析与优化打下实践基础。
1. 两机五节点网络潮流计算到底在算什么
课程设计里给一张两机五节点接线图,两台发电机分别挂在节点1和节点2,另外三个节点接负荷,要求用牛拉法和PQ法各算一遍,交一份电力系统稳态分析报告书。剥掉报告书的外壳,真正要算的只有两样:每个节点的电压幅值和相角,以及由电压推出来的支路潮流和全网网损。反直觉的地方在于,两种方法收敛之后结果几乎重合,误差落在小数点后五六位,但迭代过程的代价差一个量级——牛拉法每轮要组装并求解一个7阶雅可比矩阵,PQ法把它拆成4阶和3阶两组,而且分解后的系数矩阵在整个迭代过程中保持不变。这篇文章适合正在做电力系统稳态分析课程设计、手里已经有电力系统潮流计算MATLAB环境、但不想只抄一份代码交差的人,从节点导纳矩阵怎么填、雅可比四个子块怎么摆,一直讲到收敛不了的时候先看哪一行。
2. 两机五节点网络的数据建模与节点导纳矩阵组装
2.1 节点类型划分:1 个平衡节点、1 个 PV 和 3 个 PQ
潮流计算里每个节点有四个量:P、Q、V、θ,但每个节点只能给两个已知量,剩下两个靠方程解出来,于是节点被分成三类。两机五节点网络里,节点1接主发电机,让它同时承担全网网损和相角参考,所以定为平衡节点,给定 V 和 θ;节点2接另一台机组,有励磁调节能力、能维持机端电压,所以定为 PV 节点,给定 P 和 V;节点3、4、5是纯负荷节点,有功无功都给定,定为 PQ 节点。
这样一划分,未知量就数得出来了:θ2 到 θ5 共 4 个,V3 到 V5 共 3 个,合计 7 个。牛拉法要解的就是 7 阶线性方程组,PQ法解的是 4 阶和 3 阶两个小方程组。下面这张表就是课程设计里最常出现的一组数据,标幺值基准取 Sb = 100 MVA。
| 节点 | 类型 | Pg | Qg | Pl | Ql | V 初值 |
|---|---|---|---|---|---|---|
| 1 | 平衡 | — | — | 0.00 | 0.00 | 1.05 |
| 2 | PV | 0.40 | — | 0.00 | 0.00 | 1.02 |
| 3 | PQ | 0.00 | 0.00 | 0.60 | 0.30 | 1.00 |
| 4 | PQ | 0.00 | 0.00 | 0.40 | 0.20 | 1.00 |
| 5 | PQ | 0.00 | 0.00 | 0.30 | 0.15 | 1.00 |
提示:PV 节点的 Qg 是待求量,不要提前填数;平衡节点的 Pg、Qg 都要等潮流收敛后由全网功率平衡倒推出来,填进去只会让程序逻辑变乱。
2.2 支路参数表怎么填:R、X、B 与标幺值换算
五节点网络一般给 6 到 7 条支路。题目如果直接给标幺值,照抄即可;如果给的是有名值,要先除以阻抗基准 Zb = U²/Sb。线路的对地充电电纳 B 按整条线路给,组装时要对半分到两端,这是很多人第一次写 Ybus 时容易漏掉的一步。
| 首端 | 末端 | R | X | B |
|---|---|---|---|---|
| 1 | 2 | 0.02 | 0.06 | 0.030 |
| 1 | 3 | 0.08 | 0.24 | 0.025 |
| 2 | 3 | 0.06 | 0.18 | 0.020 |
| 2 | 4 | 0.06 | 0.18 | 0.020 |
| 2 | 5 | 0.04 | 0.12 | 0.015 |
| 3 | 4 | 0.01 | 0.03 | 0.010 |
| 4 | 5 | 0.08 | 0.24 | 0.025 |
这套参数里 3-4 支路阻抗最小,是电气距离最近的通道;1-3 和 4-5 阻抗偏大,属于相对薄弱的联络,潮流结果里这两条支路的电压降落会明显一些。
2.3 用 MATLAB 组装节点导纳矩阵
导纳矩阵的规则很简单:对角线元素等于该节点所连全部支路导纳之和,再加对地充电电纳;非对角线元素等于两节点之间支路导纳的负数;没有直接相连的节点位置为 0。
function Ybus = buildYbus(bus, branch) % bus: 节点表, branch: 支路表 [首端 末端 R X B] n = size(bus,1); Ybus = zeros(n); for k = 1:size(branch,1) i = branch(k,1); j = branch(k,2); z = branch(k,3) + 1j*branch(k,4); % 支路阻抗 y = 1/z; % 串联导纳 yc = 1j*branch(k,5)/2; % 充电电纳对半分到两端 Ybus(i,i) = Ybus(i,i) + y + yc; Ybus(j,j) = Ybus(j,j) + y + yc; Ybus(i,j) = Ybus(i,j) - y; Ybus(j,i) = Ybus(j,i) - y; end end这里yc取一半是关键,如果整条线路的 B 全加到一端,Ybus 就不对称了,迭代次数会异常增大甚至不收敛。运行完之后,Ybus 是一个 5×5 的复对称矩阵,对角线上每一项的虚部通常远大于实部,这是高压电网 R≪X 的典型特征,也是后面 PQ 分解法能成立的物理基础。
3. 牛拉法潮流计算:7 阶雅可比矩阵的组装与迭代
3.1 极坐标形式下的功率不平衡量
节点注入功率用导纳矩阵写出来是 S = V·conj(Ybus·V)。展开成极坐标形式,节点 i 的有功无功分别是:
P_i = V_i · Σ_j V_j (G_ij·cosθ_ij + B_ij·sinθ_ij) Q_i = V_i · Σ_j V_j (G_ij·sinθ_ij − B_ij·cosθ_ij)
其中 θ_ij = θ_i − θ_j,G、B 分别是 Ybus 元素的实部和虚部。给定值和计算值相减就是不平衡量 ΔP_i = Ps_i − P_i、ΔQ_i = Qs_i − Q_i,牛拉法要做的就是让这 7 个不平衡量同时趋近于 0。判据一般取 max|ΔP|、max|ΔQ| 小于 1e-6(标幺值),折成有名值相当于 1e-4 MW 量级,对课程设计足够了。
3.2 雅可比矩阵 H、N、M、L 四个子块
雅可比矩阵按未知量分成四块,行对应不平衡量、列对应待求修正量:
| 子块 | 含义 | 本系统中的阶数 |
|---|---|---|
| H | ∂P/∂θ | 4×4 |
| N | ∂P/∂V | 4×3 |
| M | ∂Q/∂θ | 3×4 |
| L | ∂Q/∂V | 3×3 |
非对角元素(i ≠ j)的表达式是固定的:
- H_ij = V_i·V_j·(G_ij·sinθ_ij − B_ij·cosθ_ij)
- N_ij = V_i·(G_ij·cosθ_ij + B_ij·sinθ_ij)
- M_ij = −V_i·V_j·(G_ij·cosθ_ij + B_ij·sinθ_ij)
- L_ij = V_i·(G_ij·sinθ_ij − B_ij·cosθ_ij)
对角元素用更紧凑的写法,也是实际编程里最省事的写法:
- H_ii = −Q_i − B_ii·V_i²
- N_ii = P_i/V_i + G_ii·V_i
- M_ii = P_i − G_ii·V_i²
- L_ii = Q_i/V_i − B_ii·V_i
注意:这四个对角表达式里的 P_i、Q_i 是本次迭代算出来的功率,不是给定的 Ps、Qs。把它们写成给定值,程序会收敛到一个完全错误的结果,而且迭代曲线看起来很“正常”。
3.3 完整可运行的 MATLAB 牛拉法程序
function [Vm, Va, hist] = nr_pf(Ybus, bus, tol, maxit) n = size(bus,1); type = bus(:,1); % 1=平衡 2=PV 3=PQ Vm = bus(:,6); Va = zeros(n,1); Ps = bus(:,2) - bus(:,4); % 给定有功注入 Qs = bus(:,3) - bus(:,5); % 给定无功注入 ang_idx = find(type ~= 1); % 角度未知量下标 pq = find(type == 3); % 电压未知量下标 na = numel(ang_idx); npq = numel(pq); hist = []; for it = 1:maxit V = Vm .* exp(1j*Va); S = V .* conj(Ybus*V); P = real(S); Q = imag(S); mis = [Ps(ang_idx) - P(ang_idx); Qs(pq) - Q(pq)]; hist(end+1) = max(abs(mis)); if hist(end) < tol, break; end G = real(Ybus); B = imag(Ybus); H = zeros(na,na); N = zeros(na,npq); M = zeros(npq,na); L = zeros(npq,npq); for a = 1:na i = ang_idx(a); for b = 1:na j = ang_idx(b); if i == j H(a,b) = -Q(i) - B(i,i)*Vm(i)^2; else t = Va(i) - Va(j); H(a,b) = Vm(i)*Vm(j)*(G(i,j)*sin(t) - B(i,j)*cos(t)); end end for b = 1:npq j = pq(b); if i == j N(a,b) = P(i)/Vm(i) + G(i,i)*Vm(i); else t = Va(i) - Va(j); N(a,b) = Vm(i)*(G(i,j)*cos(t) + B(i,j)*sin(t)); end end end for a = 1:npq i = pq(a); for b = 1:na j = ang_idx(b); if i == j M(a,b) = P(i) - G(i,i)*Vm(i)^2; else t = Va(i) - Va(j); M(a,b) = -Vm(i)*Vm(j)*(G(i,j)*cos(t) + B(i,j)*sin(t)); end end for b = 1:npq j = pq(b); if i == j L(a,b) = Q(i)/Vm(i) - B(i,i)*Vm(i); else t = Va(i) - Va(j); L(a,b) = Vm(i)*(G(i,j)*sin(t) - B(i,j)*cos(t)); end end end dx = [H N; M L] \ mis; Va(ang_idx) = Va(ang_idx) + dx(1:na); % 修正相角 Vm(pq) = Vm(pq) + dx(na+1:end); % 修正 PQ 节点电压 end end逻辑说明:每一轮先由当前电压算出各节点注入功率,跟给定值相减得到不平衡量,判断是否满足精度;不满足就按四块公式组装雅可比,解出修正量 dx,前 na 个分量加到相角上,后 npq 个分量加到 PQ 节点电压上。PV 节点的电压不参与修正,这正是它电压恒定的体现。参数方面,tol取 1e-6、maxit取 20 通常够用,两机五节点这个规模三到四轮就能收敛到 1e-8 以下。
另外,hist里存的是每轮的最大不平衡量,画成半对数曲线是一条接近直线的下降线,这正是牛拉法二阶收敛的直观证据,报告书里的收敛特性图可以直接用它画。
4. PQ 分解法:把 7 阶方程拆成 P-θ 和 Q-V 两组
4.1 快速解耦的三个假设与 B′、B″ 的抽取规则
PQ 分解法能成立,靠的是高压电网的三条近似:支路电阻远小于电抗,所以 G_ij ≈ 0;相邻节点相角差很小,所以 cosθ_ij ≈ 1、sinθ_ij ≈ 0;无功主要影响电压幅值,有功主要影响相角。三条近似压下来,雅可比矩阵的非对角块 N、M 直接消掉,H 退化成只跟电纳有关,L 也一样。
于是迭代方程变成两组独立的线性方程:
ΔP/V = B′·Δθ ΔQ/V = B″·ΔV
B′ 取非平衡节点对应的行列,B″ 只取 PQ 节点对应的行列。工程实现里两者都从 −imag(Ybus) 里抽,这是最省事也最稳的做法。如果按教材上那种“忽略电阻、只用电抗倒数的支路电纳”来构造 B′,收敛会略快一点,但两机五节点这种规模差别看不出来,反而不容易写错。
注意:B′ 是 4 阶、B″ 是 3 阶,它们的阶数分别等于非平衡节点数和 PQ 节点数,不要机械地都取成 5 阶然后一头雾水地发现矩阵奇异。
4.2 PQ 法主循环的 MATLAB 实现
function [Vm, Va, hist] = fdlf_pf(Ybus, bus, tol, maxit) n = size(bus,1); type = bus(:,1); Vm = bus(:,6); Va = zeros(n,1); Ps = bus(:,2) - bus(:,4); Qs = bus(:,3) - bus(:,5); ang_idx = find(type ~= 1); pq = find(type == 3); Bfull = -imag(Ybus); Bp = Bfull(ang_idx, ang_idx); % 4x4,整个迭代过程不变 Bpp = Bfull(pq, pq); % 3x3,整个迭代过程不变 hist = []; for it = 1:maxit V = Vm .* exp(1j*Va); S = V .* conj(Ybus*V); P = real(S); Q = imag(S); dP = (Ps(ang_idx) - P(ang_idx)) ./ Vm(ang_idx); dQ = (Qs(pq) - Q(pq)) ./ Vm(pq); hist(end+1) = max([abs(dP); abs(dQ)]); if hist(end) < tol, break; end Va(ang_idx) = Va(ang_idx) + Bp \ dP; % 先解相角 Vm(pq) = Vm(pq) + Bpp \ dQ; % 再解电压 end end这里 Bp 和 Bpp 在循环外算好,循环里只做两次前代回代,这就是 PQ 法每轮比牛拉法快的原因。严格实现里应该先更新相角、重新计算无功不平衡量、再解电压修正,这样收敛轮数能少一轮;上面这个简化版把两步合在一轮里,对最终结果没有影响,只是迭代次数会多一到两轮。
初值方面,PQ 法的电压初值一定要用平启动,即 PV 和平衡节点取给定值、PQ 节点取 1.0。如果有人手抖把 PQ 节点电压初值填成 1.05,迭代轮数会明显增加,容易误判成方法本身收敛慢。
4.3 两种算法的迭代次数与结果对照
用上面同一套数据跑,结果大致如下表。牛拉法三轮收敛,PQ 法五轮左右收敛,两者终值在 1e-6 精度下基本一致。
| 方法 | 迭代轮数 | 单轮计算量 | 节点3电压 | 节点5电压 |
|---|---|---|---|---|
| 牛拉法 | 3~4 | 组装 7×7 雅可比并求解 | 见运行结果 | 见运行结果 |
| PQ 法 | 5~6 | 两次固定矩阵求解 | 与牛拉法一致到 1e-6 | 与牛拉法一致到 1e-6 |
两者结果差异主要来自收敛精度,把 tol 都压到 1e-8 之后,节点电压幅值差别通常在 1e-7 量级,这个量级完全可以忽略。如果两种方法算出来差得比较多,比如差到小数点后两位,问题一定出在 B′、B″ 的抽取或者雅可比对角元素的写法上,而不是方法本身。
5. 收敛失败与结果校核:几个能省几小时的排查手法
先说最常见的三个报错。
第一个是矩阵奇异,报错信息一般来自反斜杠那一行。绝大多数情况是平衡节点的行列没有剔除干净,或者 PV 节点被误当成了 PQ 节点塞进 B″ 里。检查方式很简单:打印size(Bp)和size(Bpp),两机五节点出来的必须是 4 和 3。
第二个是不收敛、不平衡量在某个值附近来回跳。先看雅可比对角线里的 P_i、Q_i 是不是错用了给定值,再检查充电电纳是不是被整体加到了一端。还有一个隐蔽的坑:Ybus 组装时如果支路表里出现了重复支路,对角元素会多加一次导纳,迭代会稳定地发散。
第三个是收敛但结果离谱,比如节点电压跑到 0.6 以下。这通常是支路阻抗单位没统一,一半用标幺值一半用有名值。统一成标幺值之后重新跑一遍就能确认。
跑出结果之后,校核比调试更重要。我一般按这个顺序来:
- 把最终电压代回 S = V·conj(Ybus·V),看每个 PQ 节点的计算注入跟给定值是否一致,误差应在 1e-6 以内。
- 统计全网功率平衡:所有发电机有功之和减去所有负荷有功之和,应等于全网有功损耗;无功同理。
- 用支路潮流公式 S_ij = V_i·conj((V_i − V_j)·y + V_i·yc) 逐条算,检查是否有支路潮流方向跟直觉相反、或者末端电压反超首端。
第三步里那个支路潮流公式,注意充电电纳只能算本端那一半,写成整条支路的 B 会让网损偏大,而功率平衡校验又刚好能对得上,很容易被忽略。把这三点都过一遍,报告书里的潮流分布表、网损表、收敛曲线基本就能直接用了。
本文还有配套的精品资源,点击获取