简介这是一份面向电气工程及其自动化专业学生的电力系统分析课程设计资料聚焦潮流计算这一稳态分析的核心内容适合正在完成课程设计、需要理解算法原理并动手编程实现的学习者。资源为单个PDF文档压缩包约376KB正文含概述、计算方法简介、题目解析、程序设计、总结与参考文献及附录等章节结构完整。文档先说明潮流计算在电网规划、运行方式编制与故障分析中的作用再引入牛顿-拉夫逊算法并逐项梳理节点类型划分、待求量确定、节点导纳矩阵形成与潮流方程建立等关键环节。第4章给出程序框图与基于MATLAB的算例计算对电压幅值、相角及功率分布结果进行分析读者可据此掌握从建模、编程到结果校验的完整思路。目前已有1510人学习该资料可作为课程设计报告的写作范本与算法复习参考。1. 潮流计算在课程设计里到底算什么从节点方程到可交付报告课程设计题目一出来写着潮流计算很多人的第一反应是打开 MATLAB 从零敲牛顿法结果卡在雅可比矩阵的偏导符号上三天最后靠抄一份同学的程序交差。其实电力系统分析的课程设计考的不是推导能力而是一条完整链路把给定的支路参数和节点负荷整理成节点导纳矩阵选一种迭代方法把非线性功率方程组解到收敛再回头核对功率平衡、网损和电压是否越限最后把过程写成别人照着能复现的报告。这条链路里每一步都有明确的输入输出任何一步对不上后面全会歪。标题里的潮流计算本质是在给定网络拓扑、支路阻抗和节点注入功率的前提下求各节点电压幅值与相角进而算出支路功率和网损。河南工业大学的课程设计通常给出一个几节点到十几节点的算例要求用 MATLAB 编程实现并给出迭代收敛过程。适合两类人看一类是第一次做课程设计、需要一个能跑通的最小闭环的人另一类是程序能跑但结果和参考答案对不上、想知道错在哪的人。下面按数学模型、牛顿法实现、PQ 分解法与排错、报告出图四段推下去。2. 潮流计算的数学模型与 MATLAB 环境准备先把方程写清楚再动手写代码否则后面调试会无从下手。潮流计算要解的是节点功率平衡方程对每个节点 i注入功率等于该节点电压与所有相邻节点导纳作用的结果。把节点电压写成 V_i Vm_i·exp(j·θ_i) 的形式展开后功率方程变成 Vm 和 θ 的非线性三角函数组合节点数稍微一多就没法解析求解只能迭代。2.1 节点导纳矩阵从支路参数表到 Ybus 的组装课程设计给的原始数据一般是一张支路表每行包括首末端节点、电阻、电抗、对地充电电纳如果有变压器还要给变比。组装 Ybus 的规则固定对角线元素是接在该节点上所有支路导纳之和加上对地电纳非对角线元素是两节点间支路导纳取负号。下面是可用的 MATLAB 函数。function Ybus build_ybus(branch, nb) % branch 每行: [i j R X Bc kT]i、j 为 1 起始节点编号 % R、X 为标幺值Bc 为线路总充电电纳kT 为变比(无变压器填 0) Ybus zeros(nb); for k 1:size(branch,1) i branch(k,1); j branch(k,2); R branch(k,3); X branch(k,4); Bc branch(k,5); kT branch(k,6); y 1/(R 1i*X); % 串联导纳 if kT 0 % 普通线路 Ybus(i,i) Ybus(i,i) y 1i*Bc/2; Ybus(j,j) Ybus(j,j) y 1i*Bc/2; Ybus(i,j) Ybus(i,j) - y; Ybus(j,i) Ybus(j,i) - y; else % 变压器支路i 侧接变比 Ybus(i,i) Ybus(i,i) y/(kT^2); Ybus(j,j) Ybus(j,j) y; Ybus(i,j) Ybus(i,j) - y/kT; Ybus(j,i) Ybus(j,i) - y/kT; end end end逻辑说明普通线路按 π 型等值两端各摊一半充电电纳互导纳取负串联导纳。变压器支路把阻抗折算到非变比侧自导纳多出 kT² 因子互导纳除以 kT。参数上最容易错的是 kT 填成实际变比还是标幺变比——课程设计里通常直接给标幺变比如果给的是电压比要先除以基准电压比再填进去。注意Ybus 必须是严格对称的复矩阵。跑完这段代码先执行norm(Ybus - Ybus.)结果应当接近零不为零就说明某条支路的互导纳只写了一个方向。2.2 功率方程与 PQ、PV、平衡节点的约束差别节点类型决定了哪些量已知、哪些量待求这是潮流计算的核心分类也是初学者最容易混的地方。三类节点的已知量、待求量和对应方程数见下表。节点类型已知量待求量提供的方程典型设备PQ 节点P、QVm、θΔP、ΔQ负荷母线PV 节点P、VmQ、θΔP发电机母线平衡节点Vm、θP、Q无主调频厂一个 n 节点系统里设 PQ 节点 p 个、PV 节点 v 个平衡节点 1 个则未知量总数是 2p v正好等于能列出的方程数 2p v。这个计数关系是检查自己程序是否漏方程的最快办法如果雅可比矩阵不是方阵一定是节点类型数统计错了。课程设计里经常把某个发电机节点写成 PV但又把它的 Q 给了固定值这种自相矛盾的输入会让迭代直接发散。2.3 MATLAB 环境、目录结构与标幺值基值约定不需要什么工具箱基础 MATLAB 就够。建议目录分成数据文件、函数文件、脚本文件三层case5_data.m存支路表和节点表build_ybus.m、newton_pf.m存函数run_pf.m做顶层调度。这样改算例只动数据文件不用翻主程序。标幺值方面统一选一个基准功率常见 100 MVA各电压等级选各自的基准电压阻抗标幺值按 Z_base V_base²/S_base 折算。如果老师给的支路参数已经是标幺值直接抄进去如果是欧姆值先查清楚对应电压等级的基准电压再折算这一步不做直接算结果会差几个数量级而且看起来还挺「像个结果」很容易蒙混过关到答辩。3. 牛顿-拉夫逊法潮流计算的 MATLAB 实现牛顿法的思路是把非线性方程在当前点做一阶泰勒展开用线性方程的解作为修正量反复迭代。收敛速度是平方级的正常算例三四次就到 1e-8这是它成为课程设计默认方法的原因。代价是每次迭代都要重新组装雅可比矩阵代码量比 PQ 分解法大一截但换来的是对病态算例更强的适应力。3.1 雅可比矩阵各块偏导数的推导与对照表雅可比矩阵按四个子块组织H 是 ΔP 对 Δθ 的偏导N 是 ΔP 对 ΔVm 的偏导M 是 ΔQ 对 Δθ 的偏导L 是 ΔQ 对 ΔVm 的偏导。记 G、B 为 Ybus 的实部虚部θ_ij θ_i − θ_j各元素公式如下表对角元要用到节点当前注入功率 P_i、Q_i。子块非对角元 (i≠j)对角元HVm_i·Vm_j·(G_ij·sinθ_ij − B_ij·cosθ_ij)−Q_i − B_ii·Vm_i²NVm_i·(G_ij·cosθ_ij B_ij·sinθ_ij)P_i/Vm_i G_ii·Vm_iM−Vm_i·Vm_j·(G_ij·cosθ_ij B_ij·sinθ_ij)P_i − G_ii·Vm_i²LVm_i·(G_ij·sinθ_ij − B_ij·cosθ_ij)Q_i/Vm_i − B_ii·Vm_i对照表容易记混的是 M 块比 H 块多一个负号N 和 L 的非对角公式结构相同但一个用 cos 配 sin、一个用 sin 配 cos。写代码时建议先按公式逐个填别图省事用统一表达式符号错一个整轮迭代都跑偏。3.2 迭代主循环、收敛判据与不平衡量更新把上面的公式落成代码主循环如下。function [Vm, Va, hist] newton_pf(Ybus, Sbus, Vm, Va, pq, pv, tol, maxit) nb numel(Vm); hist []; idxP [pq; pv]; idxQ pq; % 有 P 方程 / Q 方程的节点 for it 1:maxit V Vm .* exp(1i*Va); Sinj V .* conj(Ybus * V); % 当前注入功率 dP real(Sbus - Sinj); dQ imag(Sbus - Sinj); mis [dP(idxP); dQ(idxQ)]; % 不平衡量 if max(abs(mis)) tol, break; end hist(end1) max(abs(mis)); % 记录收敛过程 G real(Ybus); B imag(Ybus); nP numel(idxP); nQ numel(idxQ); H zeros(nP); N zeros(nP,nQ); M zeros(nQ,nP); L zeros(nQ); for a 1:nP i idxP(a); for b 1:nP j idxP(b); th Va(i) - Va(j); if i j H(a,b) -imag(Sinj(i)) - B(i,i)*Vm(i)^2; else H(a,b) Vm(i)*Vm(j)*(G(i,j)*sin(th) - B(i,j)*cos(th)); end end for b 1:nQ j idxQ(b); th Va(i) - Va(j); if i j N(a,b) real(Sinj(i))/Vm(i) G(i,i)*Vm(i); else N(a,b) Vm(i)*(G(i,j)*cos(th) B(i,j)*sin(th)); end end end for a 1:nQ i idxQ(a); for b 1:nP j idxP(b); th Va(i) - Va(j); if i j M(a,b) real(Sinj(i)) - G(i,i)*Vm(i)^2; else M(a,b) -Vm(i)*Vm(j)*(G(i,j)*cos(th) B(i,j)*sin(th)); end end for b 1:nQ j idxQ(b); th Va(i) - Va(j); if i j L(a,b) imag(Sinj(i))/Vm(i) - B(i,i)*Vm(i); else L(a,b) Vm(i)*(G(i,j)*sin(th) - B(i,j)*cos(th)); end end end dx [H N; M L] \ mis; % 修正量 Va(idxP) Va(idxP) dx(1:nP); Vm(idxQ) Vm(idxQ) dx(nP1:end); end end参数说明Sbus是节点注入复功率注意发电机节点要填正、负荷节点填负平衡节点的 Sbus 随便填不影响结果因为它的方程不参与。pq、pv是节点编号列向量测试时先打印length(pq)length(pv)1 nb确认没漏节点。tol一般取 1e-8maxit取 10 到 20 都够。3.2.1 收敛判据怎么设与初值怎么给判据用不平衡量最大绝对值小于 tol 就可以不用改成范数课程设计不需要那么讲究。真正影响收敛的是初值平启动所有 Vm 取 1.0θ 取 0在大多数算例上都能收敛如果发散先把所有负荷节点的 Vm 设成 0.95 再试。另一个坑是 PV 节点的 Vm 在迭代中不能被更新代码里idxQ只包含 PQ 节点PV 节点的 Vm 始终保持设定值这一点写错会让结果偏移到莫名其妙的位置。3.3 五节点算例跑通与结果对照准备一个 5 节点算例1 号平衡节点 Vm1.05、θ02、3 号是 PV 节点4、5 号是 PQ 节点支路参数填成标幺值。顶层脚本长这样。run(case5_data.m); % 载入 nb、branch、bus、Sbus Ybus build_ybus(branch, nb); Vm ones(nb,1); Va zeros(nb,1); Vm(slack) 1.05; Vm(pv) bus(pv,2); % 用给定值覆盖初值 [Vm, Va, hist] newton_pf(Ybus, Sbus, Vm, Va, pq, pv, 1e-8, 20); V Vm .* exp(1i*Va); S_calc V .* conj(Ybus * V); fprintf(迭代次数 %d末次不平衡量 %.3e\n, numel(hist)1, hist(end)); disp(table((1:nb)., Vm, rad2deg(Va), real(S_calc), imag(S_calc)));跑完后第一件事是看迭代次数正常情况下是 3 到 5 次如果超过 10 次还没收敛说明要么数据有错要么某个节点的 Sbus 符号反了。第二件事是把S_calc和Sbus逐节点相减除平衡节点外差值应当在 1e-6 以内。第三件事是算平衡节点的实际注入功率这个值不是输入给的是迭代的结果也是网损核算的入口。4. PQ 分解法、收敛排错与结果校验牛顿法代码写对之后老师往往会追问一句「能不能用 PQ 分解法再算一遍」或者在报告里要求对比两种方法的收敛速度。PQ 分解法是牛顿法在高压电网 R≪X 条件下的简化版本把雅可比矩阵拆成两个常数矩阵迭代中不再重新组装单次迭代快但迭代次数多节点数越多优势越明显。4.1 PQ 分解法的 B′ 与 B″ 矩阵构造PQ 分解法的核心是两个常数矩阵B′ 用于修正相角B″ 用于修正电压幅值。教材上的做法是 B′ 只取支路电抗的倒数并忽略电阻和对地支路B″ 取节点导纳矩阵的虚部。工程实现里直接用 Ybus 的虚部更省事代码如下。idxP [pq; pv]; idxQ pq; Bp -imag(Ybus(idxP, idxP)); % B 矩阵取负号放在方程里 Bpp -imag(Ybus(idxQ, idxQ)); % B 矩阵 Vm ones(nb,1); Va zeros(nb,1); for it 1:30 V Vm .* exp(1i*Va); Sinj V .* conj(Ybus * V); dP real(Sbus - Sinj); dQ imag(Sbus - Sinj); if max(abs([dP(idxP); dQ(idxQ)])) 1e-8, break; end Va(idxP) Va(idxP) Bp \ (dP(idxP) ./ Vm(idxP)); % 先修相角 V Vm .* exp(1i*Va); % 用新相角更新功率 Sinj V .* conj(Ybus * V); dQ imag(Sbus - Sinj); Vm(idxQ) Vm(idxQ) Bpp \ (dQ(idxQ) ./ Vm(idxQ)); % 再修幅值 end逻辑说明每次迭代分两半走先解相角修正量并立刻更新电压相量再用更新后的功率解幅值修正量。参数上 B′ 和 B″ 都是实数对称矩阵可以用\直接解规模小的时候不需要稀疏处理。注意这里两个矩阵都取了负号如果把负号去掉迭代会朝反方向修正表现为不平衡量越迭代越大。4.2 不收敛时的排查顺序与常见误用牛顿法和 PQ 分解法不收敛的原因基本集中在输入侧按下面顺序查能覆盖九成情况。现象常见原因检查方式第一次迭代就不平衡量巨大Sbus 单位是 MW 没转标幺打印max(abs(Sbus))应在个位数迭代几次后发散某条支路 R、X 填反或符号错查abs(imag(1./(R1i*X)))不平衡量震荡不下降平衡节点编号没从 idxP 里排除确认slack不在pq、pv中收敛但电压畸变PV 节点 Vm 被迭代改掉了检查idxQ是否只含 PQ 节点功率永远配不平节点负荷总和与发电机出力差太多算总注入看是否为负的网损量级常见误用还有两种。一是把所有节点都设成 PQ结果平衡节点也被要求满足功率方程方程组无解二是把变比 kT 按实际电压比填成 110/10.5 这种数值标幺变比应该是 1.05 一类的量级填错会让变压器支路的导纳大出几个数量级。4.3 结果校验功率平衡、网损与电压限值收敛不等于正确课程设计报告里必须有校验环节。三个必做的检查节点功率偏差、全网功率平衡、电压是否落在合理区间。V Vm .* exp(1i*Va); S_calc V .* conj(Ybus * V); err Sbus - S_calc; err(slack) 0; % 平衡节点不参与校验 fprintf(最大节点功率偏差: %.3e p.u.\n, max(abs(err))); S_loss sum(S_calc); % 全网注入之和即网损 fprintf(网损: %.3f j%.3f p.u.\n, real(S_loss), imag(S_loss)); fprintf(电压范围: %.4f ~ %.4f p.u.\n, min(Vm), max(Vm));逻辑说明平衡节点的注入功率是待求量校验时把它剔除。全网所有节点注入功率之和在数值上等于网损实部为正说明有功损耗合理虚部反映无功分布。电压范围一般要求在 0.95 到 1.05 之间超出就说明某条线路无功缺额大需要在报告里说明原因而不是硬改数据。这三项都过了再拿程序结果和老师给的参考答案逐节点对比通常只会在第四位小数上有差异。5. 课程设计报告出图出表与答辩常见追问的应对报告里最能拉开差距的不是代码而是把迭代过程可视化和把结果解释清楚。收敛曲线用半对数坐标画因为不平衡量是指数下降的线性坐标上看不出平方收敛的特性。5.1 出图出表把收敛曲线和电压分布画进报告在牛顿法函数里把每轮的max(abs(mis))存进hist画图脚本如下。figure; semilogy(1:numel(hist), hist, -o, LineWidth, 1.2); xlabel(迭代次数); ylabel(最大不平衡量 / p.u.); grid on; title(牛顿-拉夫逊法收敛特性); figure; bar((1:nb)., Vm, 0.5); hold on; plot([1 nb], [0.95 0.95], r--); plot([1 nb], [1.05 1.05], r--); xlabel(节点编号); ylabel(电压幅值 / p.u.);第一张图要标出从第几轮开始不平衡量低于 1e-8第二张图把上下限画成虚线一眼看出哪些节点越限。表格方面节点结果表给节点号、Vm、θ(度)、P、Q 五列就够支路结果表给首末端、有功、无功、损耗四列数值统一保留四位小数。5.2 答辩追问的高频问题与准备方式老师问得最多的几个点为什么平衡节点不能设成 PQ雅可比矩阵某块的物理解释是什么牛顿法和 PQ 分解法各迭代了几次、时间差多少网损的实部为什么是正的。前两个考理解建议把 H 块对角元−Q_i − B_ii·Vm_i²里的 Q 和 B 拆开讲后两个考数据把两次运行的时间用tic/toc记下来放进报告网损为正在于线路电阻消耗有功这一句写清楚就够。一个容易被忽略的加分点把负荷增大到 1.5 倍再跑一次记录收敛迭代次数的变化说明重载下电压下降、雅可比矩阵条件数变差。这种小实验花不了十分钟但报告的技术含量会明显高一档。另一类是改用不同的初值做敏感性对比平启动和热启动各跑一次把不平滑收敛的曲线截下来分析通常比堆十页公式更受欢迎。最后一个实操建议所有代码和数据文件按算例分包命名case5、case14各自独立主脚本只改一行文件名就能切换。答辩现场临时被要求换算例演示时这一行改动就是全部的响应时间。本文还有配套的精品资源点击获取