做电力系统题目绕不开潮流计算。前两天整理之前写的两节点系统Gauss-Seidel法MATLAB程序发现这个例子虽然模型简单但它把潮流计算里最关键的节点导纳矩阵、PQ节点电压迭代、复功率符号约定这些东西全串起来了——很适合作为入门的第一块拼图。这篇文章就把这个两节点系统的代码从原理到实现拆开讲一遍包含可直接复现的MATLAB脚本、迭代公式的推导、正负号陷阱和常见坑。适合正在学电力系统分析的学生也适合想快速捡回潮流算法的工程师。花半小时跟着敲一遍效果比单纯看书强得多。1. 两节点系统里到底在算什么项目背景与核心思路1.1 潮流计算的本质从“电流平衡”到“功率平衡”潮流计算也叫潮流分布是电力系统分析里最基础的一块内容。所谓潮流指的是电网在不同运行方式下各节点的电压幅值、相角以及各条支路的有功和无功分布。给定发电机出力和负荷需求求解出全网的电压状态和功率流动这就是潮流计算要干的事。很多教材把它描述成“解非线性方程组”听起来很吓人但其实核心模型并不复杂每个节点都满足基尔霍夫电流定律也就是电流平衡方程再把电流用复功率和电压表示就变成了功率平衡方程。为什么不能像普通电路分析那样直接列欧姆定律一次出结果因为电网里的负荷绝大多数给定的是功率不是电流。而功率和电压之间是乘积关系电压一变电流跟着变功率又要满足给定值这就形成了一个非线性方程组。直白点说你没法把电压“解”出来只能猜一个初始电压然后按照功率偏差反复修正直到稳定。Gauss-Seidel法就是最经典的迭代修正方法之一。它不追求一步到位而是先给所有未知电压设初值然后一个个节点轮流更新每更新一个就立刻用等电压变化足够小的时候就认为收敛了。可以拿供水管网做个类比电网里的发电机像水厂出口压力已知用户像各个家庭用水量已知但每个用户处的实际水压是多少取决于整张管网里所有人一起用水后的综合结果。你没法一眼看出每个水龙头的水压只能先假设一个压力然后根据用水量修正再全管网算一遍反复几次才能逼近真实值。潮流计算里的电压幅值和相角就是那个“水压”。1.2 为什么拿“两节点系统”练手两节点系统在电力系统教材里可能显得太小但潮流计算真正要处理的核心难点它一个都不少。这个系统只有两条母线母线1设成平衡节点电压幅值和相角都给固定母线2设成PQ节点给定有功和无功电压幅值和相角都是未知量。再加上一条线路阻抗、一个负荷就构成了一个完整的最小可求解潮流模型。别看系统小它能验证的东西很多。节点导纳矩阵怎么构建、复功率的符号约定是什么、Gauss-Seidel迭代初值怎么给、收敛判据怎么设、结果怎么用手算校验这些在两节点系统里都能一条龙跑通。等到后续接触IEEE 5节点、14节点甚至更大系统时你会发现整个迭代框架没变只是Y矩阵变大、节点类型变多核心公式还是同一个。所以两节点系统不是玩具而是一台显微镜能把潮流计算的细节全部放大给你看。另一个原因是可手算校验。两节点系统的方程只有两个实部和虚部各一个完全可以用笔算验证MATLAB结果。做数值程序最怕的就是“跑出来一个数但不知道对不对”两节点系统恰好给了你一个手算对照的机会。1.3 节点类型与PQ节点的角色潮流计算中节点通常分成三类理解它们的区别是读懂代码的前提。节点类型已知量待求量典型元件平衡节点电压幅值、相角有功、无功参考节点、主发电机PQ节点有功、无功电压幅值、相角负荷节点、恒功率设备PV节点有功、电压幅值无功、相角可调无功的发电机节点本项目只有两个节点节点1是平衡节点电压固定为1∠0°标幺值下这相当于把系统参考电压和参考相角都定死了。节点2是PQ节点负荷给定为吸收有功1.0pu、无功0.5pu我们的任务就是算出节点2的电压幅值和相角。PQ节点的“P”和“Q”代表节点注入的有功和无功是给定的。这里有一个非常容易混淆的点负荷吸收功率负号怎么写等下在迭代公式里我会专门讲。因为没有PV节点也不需要处理发电机无功修正Gauss-Seidel法的实现可以保持最纯粹的形式对初学者很友好。1.4 系统参数与导纳矩阵构建我用的系统参数如下参数符号数值说明线路阻抗Z120.02 j0.06 pu电阻0.02电抗0.06平衡节点电压V11∠0°标幺值负荷功率S_load1.0 j0.5 pu节点2吸收节点2注入功率S2spec-1.0 - j0.5 pu负荷的负值电压初值V21∠0°迭代起点收敛容差tol1e-8电压变化量线路阻抗用标幺值忽略对地导纳。那么线路导纳是Y12 1/Z12 1/(0.02j0.06)复数值大致为5 - j15。节点导纳矩阵就是Y [ Y12, -Y12 -Y12, Y12 ]这个两节点系统的Y矩阵是2×2。如果线路有对地导纳需要对角元加上对地支路导纳的一半但在最简模型里通常忽略不影响理解算法主线。为什么要用标幺值因为电网实际电压、功率、阻抗数量级差很多直接用有名值容易让数值计算出现病态。统一折算到标幺值后电压都在1附近、功率都在个位数迭代收敛判据也更好设。2. 高斯-赛德尔迭代公式是怎么一步步“挤”出来的2.1 从节点导纳方程到GS迭代公式先写节点电流方程。对任意节点 i有I_i Y_ii * V_i sum_{j≠i} Y_ij * V_j节点注入复功率定义为S_i V_i * conj(I_i)所以I_i conj(S_i / V_i) conj(S_i) / conj(V_i)。把电流方程两边整理一下把Y_ii * V_i单独拎出来Y_ii * V_i conj(S_i) / conj(V_i) - sum_{j≠i} Y_ij * V_j于是得到V_i的显式表达式V_i (conj(S_i) / conj(V_i) - sum_{j≠i} Y_ij * V_j) / Y_ii这就是Gauss-Seidel潮流的迭代种子。右边仍然含有V_i本身所以这不是纯粹的代数求解而是固定点迭代每轮用当前V_i计算出一个新的V_i然后立刻替换旧值参与下一个节点的计算。这种“算一个换一个”的更新顺序正是Gauss-Seidel和Jacobi方法的本质区别。Jacobi要用上一轮所有旧值GS用本轮已经更新的值因此GS一般收敛更快内存占用也更少。对两节点系统节点i2时求和项只有j1公式直接变成V2 (conj(S2spec) / conj(V2) - Y21 * V1) / Y22这个式子就是MATLAB代码中最核心的一行。2.2 正负号陷阱负荷到底是正还是负这一节必须单独拿出来讲因为太多新手卡在这里。潮流计算里的S_i在公式推导中默认是“注入功率”发电机向节点注入功率为正负荷从节点吸收功率等价于负注入。如果负荷实际吸收1j0.5那么参与迭代的S2spec必须写成-1-j0.5而不是1j0.5。很多人第一次写直接把负荷正数代进去结果发现节点2电压越迭代越高甚至超过平衡节点电压跑到1.05以上。这不是迭代发散而是你把负荷写成了发电机系统相当于在一个没有电源的节点上凭空注入了功率电压当然会升。后续再检查Y矩阵、收敛判据都没有用因为根子就在符号。更稳妥的写法是先把负荷定义成S_load 1.0 1i*0.5然后代码里统一写成S2spec -S_load。这样不仅代码可读性好也降低了符号错误的概率。如果你用的是实际有名值还要注意功率基准换算让S_load先除以基准容量变成标幺值再取负号。养成“注入功率为正”的思维习惯对以后写多节点潮流、牛拉法潮流都有好处。2.3 收敛判据与松弛因子迭代是否结束要用判据来判断不能单纯看循环次数。常见的判据有两类一是电压变化的模小于容差也就是相邻两次迭代的|ΔV|足够小二是功率失配小于容差也就是把当前电压代回功率方程算出的注入功率和给定功率之间的误差足够小。电压判据简单直观代码里最常用但严格来说电压变化小不代表功率一定平衡所以工程上更推荐加上功率校验。两节点系统规模小两种判据都可以跑后面我会给出带功率校验的工程版代码。Gauss-Seidel的收敛速度是线性的收敛过程中电压变化会比较均匀地衰减。为了加速可以在更新时使用松弛因子V2_new V2_old alpha * (V2_cal - V2_old)alpha大于1时叫超松弛可以加快收敛alpha小于1叫欠松弛适合系统比较刚、容易振荡的情况。经验上配电网潮流、教学示例里alpha取1.2到1.6比较常见但并不是越大越好。当alpha接近2甚至超过2时原本收敛的迭代也可能直接发散因为固定点迭代的谱半径限制了超松弛的上限。在两节点系统里可以试一下alpha1.7和alpha1.9前者可能很快收敛后者可能电压原地打转这个对比能直观体会收敛域的概念。2.4 GS法的定位和局限Gauss-Seidel法在潮流计算史上有过辉煌时期早期计算机内存极小GS法占用内存少、编程简单是大型电网潮流计算的主力算法。后来牛顿法的雅可比矩阵迭代策略逐渐成熟收敛速度远快于GS于是大规模输电系统潮流基本都转向牛拉法及其改进形式。但GS法并没有完全退出一些辐射型配电网、三相不平衡潮流、以及牛拉法初值生成里GS仍然有用武之地。更重要的是它最适合教学。GS和牛拉法的差异可以简单对比一下算法每次迭代计算量收敛速度对初值的要求适用范围Gauss-Seidel低线性较慢宽松小型系统、配电网、教学Newton-Raphson高需算雅可比二阶很快较敏感要接近解大型输电系统两节点系统用GS计算几十步以内肯定能收敛完全感受不到收敛慢的缺点。但如果你直接用这个代码去跑IEEE 30节点可能要好几百步甚至上千步到时候就能理解为什么实际工程更偏爱牛拉法了。3. MATLAB代码实现每一步都能直接抄的作业3.1 最简可运行版本下面这段代码是完整可运行的MATLAB脚本核心思路和前面公式一一对应。我刻意去掉了花哨的图形界面只保留迭代主体方便你把注意力放在算法本身。%% 两节点系统Gauss-Seidel潮流计算核心版 clear; clc; % 线路阻抗和节点导纳矩阵标幺值忽略对地导纳 Z12 0.02 1i*0.06; Y12 1 / Z12; Y [Y12, -Y12; -Y12, Y12]; % 平衡节点节点1 V1 1.0 1i*0.0; % 1∠0° % PQ节点节点2负荷吸收 1.0j0.5 pu S_load 1.0 1i*0.5; S2spec -S_load; % 注入功率 负负荷 Y22 Y(2,2); Y21 Y(2,1); % 迭代初值 V2 1.0 1i*0.0; tol 1e-8; kmax 100; alpha 1.0; % 松弛因子 fprintf(迭代步 |V2|(pu) 相角(deg) ΔV\n); for k 1:kmax V2old V2; V2cal (conj(S2spec)/conj(V2old) - Y21*V1) / Y22; V2 V2old alpha * (V2cal - V2old); dV abs(V2 - V2old); fprintf(%3d %10.6f %10.5f %.3e\n, ... k, abs(V2), angle(V2)*180/pi, dV); if dV tol fprintf(收敛于第%d步\n, k); break; end if k kmax warning(达到最大迭代次数仍未收敛); end end fprintf(\n最终结果\n); fprintf(V2 %.6f j%.6f pu\n, real(V2), imag(V2)); fprintf(|V2| %.6f pu, 相角 %.6f°\n, abs(V2), angle(V2)*180/pi);代码里最关键的一行是V2cal (conj(S2spec)/conj(V2old) - Y21*V1) / Y22;它跟公式里的(conj(S2spec)/conj(V2) - Y21*V1)/Y22完全一致。这里用conj(V2old)而不是V2old是因为复功率公式里电流等于共轭功率除以共轭电压少了共轭符号结果会完全不对这一点特别容易踩坑。3.2 带日志与功率校验的工程版如果只是为了交作业上面的核心版够了。但实际调试程序时我更建议加上日志记录和功率校验。日志记录能让你看到迭