潮流计算的计算机算法资料

129 第四章 潮流计算的计算机算法 第一节 概述 潮流计算是电力系统最基本、最常用的计算。根据系统给定的运行条件、网络接线及元件参数，通过潮流计算可以确定各母线的电压（幅值及相角），各元件中流过的功率、整个系统的功率损耗等。潮流计算是实现电力系统安全经济发供电的必要手段和重要工作环节。因此潮流计算在电力系统的规划设计、生产运行、调度管理及科学研究中都有着广泛的应用。 电力系统潮流计算分为离线潮流计算和在线潮流计算。前者主要用于系统规划设计和安排系统的运行方式，后者则用于正在运行系统的经常监视及实时控制。本章主要讨论离线潮流计算问题，它的基本算法同样适用于在线潮流计算。 潮流计算在数学上是多元非线性方程组的求解问题，求解的方法有很多种。自从五十年代计算机应用于电力系统以来，当时求解潮流的方法是以节点导纳矩阵为基础的逐次代入法（导纳法），后来为解决导纳法的收敛性较差的问题，出现了以阻抗矩阵为基础的逐次代入法（阻抗法）。到六十年代，针对阻抗法占用计算机内存大的问题又出现了分块阻抗法及牛顿－拉夫逊（Newton－Raphson）法。Newton—Raphson法是数学上解非线形方程式的有效方法，有较好的收敛性。将N－R法用于潮流计算是以导纳矩阵为基础的，由于利用了导纳矩阵的对称性、稀疏性及节点编号顺序优化等技巧，使N－R法在收敛性、占用内存、计算速度方面的优点都超过了阻抗法，成为六十年代末期以后普遍采用的方法。同时国内外广泛研究了诸如非线形规划法、直流法、交流法等各种不同的潮流计算方法。七十年代以来，又涌现出了更新的潮流计算方法。其中有1974年由B、Stott、O、Alsac提出的快速分解法以及1978年由岩本伸一等提出的保留非线性的高
130 速潮流计算法。其中快速分解法（Fast decoupled load flow）从1975年开始已在国内使用，并习惯称之为PQ分解法。由于PQ分解法在计算速度上大大超过N－R法，不但能应用于离线潮流计算，而且也能应用于在线潮流计算。 本章主要介绍最常用的N—R法和PQ分解两种潮流计算的计算机算法的原理框图及程序。 第二节 潮流计算的基本方程 一、 节点的分类 用一般的电路理论求解网络方程，目的是给出电压源（或电流源）研究网络内的电流（或电压）分布：作为基础的方程式，一般用线性代数方程式表示。然而在电力系统中，给出发电机或负荷连接母线上电压或电流（都是向量）的情况是很少的，一般是给出发电机母线上发电机的有功功率（P）和母线电压的幅值（V），给出负荷母线上负荷消耗的有功功

率（P）和无功功率（Q）。我们的目的是由这此已知量去求电力系统内的各种电气量。所以，根据电力系统中各节点性质的不同，很自然地把节点分成三种类型。 1、PQ节点 这一类节点，我们事先给定的是节点功率（P、Q），待求的未知量是节点电压向量（V、θ）。所以叫“PQ节点”。通常变电所母线都是PQ节点。当某些发电机的出力P、Q给定时，也作为PQ节点。PQ节点上的发电机称之为PQ机（或PQ给定型发电机。在潮流计算中，系统大部分节点属于PQ节点）。 2、PV节点 这类节点给出的参数是该节点的有功功率P及电压幅值V，待求量为该节点的无功功率Q及电压向量的相角θ。 这种节点在运行中往往要有一定可调节的无功电源，用以维持给定的电压值。通常选择有一定无功功率贮备的发电机母线或者变电所有无功补偿设备的母
131 线作PV节点处理。PV节点上的发电机称之为PV机（或PV给定型发电机）。 3、平衡节点 在潮流计算中，这类节点一般只设一个。对该节点，给定其电压值，并在计算中取该节点电压向量的方向作为参考轴，相当于给定该点电压向量的角度为零。也就是说，对平衡节点给定的运行参数是V和θ，因此又称为Vθ节点，而待求量是该节点的P、Q，整个系统的功率平衡由这一节点承担。 关于平衡节点的选择，一般选择系统中担任调频调压的某一发电厂（或发电机）。有时也可能按其它原则选择，例如，为提高计算的收敛性，可以选择出线数多或者靠近电网中心的发电厂母线作平衡节点。 以上三种节点四个运行参数P、Q、V、θ中已知量都是两个，待求量也是两个，只是类型不同而已。 二、 基本方程式 在潮流问题中，任何复杂的电力系统都可以归结为以下元件（参数）组成： （1）发电机（住入电流或功率） （2）负荷（负的注入电流或功率） （3）输电线支路（电阻、电抗） （4）变压器支路（电阻、电抗、变比） （5）母线上的对地支路（阻抗和导纳） （6）线路上的对地支路（一般为线路充电电容导纳） 集中了以上各种类型的元件的简单网络如图4－1(a)所示。
132 必须指出，如果仅研究稳态情况下的潮流而不涉及暂态过程的计算则不需要发电机和负荷的阻抗参数，只需要给出发电机和负荷的注入功率或电流，并且规定发电机和负荷的注入功率或电流取正，而负荷取负。 将图4—1（a）中的发电机和负荷节点用无阻抗线从网络中抽出（为不失一般性，将既非发电机又非负荷的浮动节点当作零注入功率的母线抽出网络之外），剩下的部分即由接地和不接地支路组成一个无源线性网络（图4－1（b）） 对于这个无源线性网络可用相应的

导纳矩阵（或阻抗矩阵）来描述，采用导纳矩阵焊时，节点注入电流和节点电压构成以下线性方程组 YVI Pg1+jQg1 IPg2+jQg2 2I 0+j0 3I 4I -PR4-jQR4 5I -PR5-jQR5 6I -PR6-jQR6 无源线性网络 （可用导纳 矩阵或阻抗 矩阵表示） 图4－1（b）潮流计算等值网络 ③ ④ ⑤ ① ⑥ ② 图4－1（a）潮流计算用的电网结构图
133 其中 321IIII 321VVVV 可展开为如下形式： n , , 2 , 1iY1jnjijiVI （4－1） 若 ZIV 可展开如下形式： n , , 2 , 1iZ1njjijiIV （4－2） 式中n为网络节点数 由于实际电网中测量的节点注入量一般不是电流而功率，因此必须将式（4－1）中的注入电流I用节点注入功率来表示。 根据电工理论，节点功率与节点电流之间的关系为： iiiiiIVjQPS （4－3） 式中 LDiGiiPPP LDiGiiQQQ 因此用导纳矩阵（4－1）式时，PQ节点可以表示为 iiiiiiVVSIjQ-P （4－4） 把帝个关系式代入式（4－1）中，得 n, , ,2 1ijQP1njjijiiiVYV （4－5） 比较式（4－1）和（4－5），由于功率代替电流的结果，使式（4－1）电流电压的线性方程组变量为功率和电压的非线性方程组，这
134 个非线性方程组就是潮流计算的基本方程。 式（4－5）是一组共有n个非线性方程组成的复数方程式，如果把实部和虚部分开便得到2n个实数方程，因此由该方程组可解出2n个运行参数。但是我们知道每一个节点都有4 个运行变量，即：节点的功率Pi、Qi，以及节点电压的幅值和相位（或对应于某一选定参考直角坐标的实部和虚部）。记作（Pi、Qi、Vi、θi）或（Pi、Qi、ei、fi），当节点数为n时，则共有4n个运行参数。 由2n个方程式要求出4n个运行参数是不可能的，只能求出2n个运行参数，而其余2n个应作为原始数据事先给定。这就得根据节点的分类，对每个节点的4个运行参数中的两个作为原始数据，而另外两个则作为待求量。 第三节 牛顿－拉夫逊法潮流计算 一、 牛顿－拉夫逊法概要 首先对一般的牛顿－拉夫逊法作一简单说明。已知一个变量X的函数 0Xf （4－6） 解此方程式时，由适当的近似值X（0）出发，根据 , ,2 1nXXXX'1nnnnff （4－7） 反复进行计算，当X（n）满足适当的收敛判定条件时就是（4－6）式的根。这样的方法就是所谓的牛顿－拉夫逊法。 式（4－7）就是取第n次近似解X（n）在曲线Xfy上的点nnXf , X处的切线与X轴的交点作下一次X（n+1）值的方法。参考图4－2（a）。在这一方法中为了能收敛于真解，

初值X（0）的选取及函数f（X）必须满足适当的条件，如图4－2（b）所示的那种情况就不能收敛或收敛到别的根上去。 这一方法还可以做下面的解释，设第n次迭代得到的解与真值之
135 差，即)(nX的误差为时，则 0)()(nXf （4—8） 把)()(nXf在)(nX附近对用泰勒级数展开 0)(!2)()()()(''2)(')()(nnnnXfXfXfXf （4—9） 上式略去2以下的项 0)()()(')(nnXfXf （4—10） )(nX的误差可近似由上式计算出来 图4-2 )()()(')(nnXfXf （4—11） 比较式（4—7）和（4—11），可以看出牛顿—拉夫逊法的修正量和)(nx的误差的一次项相等。 用同样的方法考虑，给出对n个变量n21X, ,X ,X的n个方程式0X , X, Xf 0X , X, Xf0X , X, Xfn 21nn 212n 211 （4－12） 对其近似解'n'2'1X , , X, X的修正量n 21X , ,X , X，可以解下面的方程式来确定 x o y y=f(x) X(t+1) X(t) x o y )(tx f(x(t)) )1(nx y=f(x) )()1(nxf
136 n21nn2n1nn22212n12111'n'2'1'n'2'12'n'2'11XXX xf xfxfxf xfxfxf xfxf X , X, XfnX , X, XfX , X, Xf （4－13） 式（4－13）的右边的矩阵的jixf等都是对于'n'2'1X , , X, X的值。这一矩阵称为雅可比（Jacobi）矩阵。 按上述得到的修正量n 21X , ,X , X后，得到如下关系： n'n''n2'2''21'1''1XXX , , XX ,XXΔXX 这比'n'2'1X , , X, X进一步接近于真值。这一步骤在收敛到希望的值以前重复进行。一般要反复计算到满足 ε , X , max121n2111nnnnnnnXXXXX 时为止。ε为预先规定的小正数，此处niX是第n次迭代Xi的近似值。 一、 牛顿－拉夫逊法潮流计算 把牛顿法用于潮流计算，要求将潮流方程改写成形如方程式（4－12）所示的形式。为此，首先应将潮流方程（4－5）的变形式jnjijiVYV1iijQPi的右端展开，并且分开实部和虚部。采用直角坐标时，节点电压可表示为： iijfeVi 节点导纳矩阵元素则表示为： ijijijjBGY 将上述表示式代入jnjijiVYV1iijQPi的右端，展开并分出实部和虚部，便得：
137 n1jn1jjijjijijijjijiin1jn1jjijjijijijjijiieBfGefBeGfQeBfGffBeGeP （4－14） 按照上节的分类，PQ节点的有功功率和无功功率给定的，第I个节点的给这功率设为Pis和Qis。假定系统中的第1，2，………m号节点为PQ节点，对其中每一个节点可列n1jn1jjijjijiiijjijiisIisin1jn1jjijjijiiijjijiisiisi0eBfGefBeGf-QQQQ0eBfGffBeGe-PP-PP （i=1，2，…………，m） （4－15） PV节点的有功功率和节点电压幅值是给定的。假定系统中的第m+1，m+2，………n-1号节点为PV节点，则对其中每一节点可以列写方程： 0feVVVV0eBfGffBeGe-PP-PP2i2i2is2i2is2in1jn1jijjijijijjijiisiisij（4－16） 1-n ,

, 2m , 1mi 第n号节点为平衡节点，其电压nnjfenV是给定的，故不参加迭代。 式（4－15）和（4－16）总共包含了2（n-1）个方程，待求的变量有1-n1-n2 211f , c , , f , e , f ,e也是2（n-1）个。我们还可以看到，方程式（4－15）和（4－16）已经具备方程组（4－12）的形式： V-JW （4－16）’ 式中
138 1-n21-n1m21mmm11VPVPQPQPW 1-n1-n1m1mmm11fefefefeV 1-n211-n211-m211-m2m2m21211211-n1-n1-n1-n1-m1-n1-m1-nm1-nm1-n11-n11-n1-n21-n2-1m21-m2m2m212121-n1m1-n1m1-m1m1-m1mm1mm1m11m11m1-nm1-nm1-mm1-mmmmmm1m1m1-nm1-nm1-mm1-mmmmmm1m1m1-n11-n11m11m1m1m111111-n11-n11m11m1m1m11111fVeVfVeVfVeVfVeVfPePfPePfPePfPePfVeVfVeVfVeVfVeVfPePfPePfPePfPePfQeQfQeQfQeQfQeQfPePfPePfPePfPePfQeQfQeQfQeQfQeQfPePfPePfPePfPeP1-n1-n1n1m1m1m1m1m1m1m1mnnnnnJ 上述方程中雅可比矩阵的各元素，可以对（4－15）和（4－16）式求偏导数获得。 当ij时，对角元素是
139 iinjinjinjinjif2fVe2eVfBeGfBeGfQfGeBeBfGeQfGeBeBfGfPfBeGfBeGePi2ii2i1iiiiijijjijii1iiiiijijiijii1iiiiiiijjijii1iiiiijijjijii （4—17） 当ij时，矩阵中非对角元素是 iijfGiijjijiiijiijjijieBfVfPfBeGfVeP （4－18） 022jijifVeV 由以上表达式不难看出，雅可比矩阵有以下特点： （1） 雅可比矩阵中的诸元素都是节点电压的函数，因此在迭代过程中，它们将随着各节点电压的变化而不断地改变； （2） 矩阵是不对称的； （3） 由式（4－18）可以看出，当导纳矩阵中的非对角元素Yij为零时，雅可比矩阵中相对应的元素也是零，即矩阵是非常稀疏的。因此，修正方程的求解同样可以应用稀疏矩阵的求解技巧。正是由于这一点才使N－R法获得广泛的应用。 三、牛顿法的框图及求解过程 程序框图如下：
140 用牛顿法计算潮流时，有以下的步骤： （1）给这各节点电压初始值00f , e （2） 将以上电压初始值代入式（4－15）和（4－16），求出修正方程式的常数项向量0200V ,Q , P； （3） 将电压初始值再代入式（4－17）和（4－18），求出修正对于PQ结点，按（4－15）计算QP,对于PV结节,，按（4－16）计算iiUP, 启 动 置0 给定电压初值00,iife 用eee1fff1修正和结点电压 形成导纳矩阵 输出出 输入原始数据 以1 用（4－17）、（4－18）式计算雅可比矩阵各元素 以ffee1,1 解修正方程（4－16），求fe, 按系统的潮流分析分布计算 平衡结点功率及线路功率 是否QP, 否 是 图4－3 牛顿法程序框图
141 方程式中系数矩阵（雅可比矩阵）的各元素； （4） 解修正方程式（4－16），求出修正量00f , e； （5） 修正各节点电压 001001fffeee （6）将11f , e再代入（4－15），（4－16）式求1211V , Q , P； （7）校验是否收敛，即 εQ , maxkikiP （8

）如果收敛，迭代到此结束，进一步计算各线路潮流和平衡节点功率，并打印输出结果。如果不收敛，转回第（2）步进行下一次迭代计算，直到收敛为止。 四、实例及程序清单 例4－1 试用牛顿－拉夫逊法计算图4－4所示电力系统的潮流分布。 P4=5 j0.015 ⑤ U4=1.05 1:1.05 ③ 0.08+j0.30 -j0.25 -j0.25 ② 1.05:1 ①20.03 U1=1.05 δ1=00 3.7+j1.3 1.6+j0.8 2+j1 ④ j0.25 j0.25 0.04+j0.25 0.1+j0.35 图4－4
142 解 ㈠ 需要输入的数据 ⑴ NN－节点数、NL－线路数、ISB－平衡母线节点号（此例为节点1）、PR－误差精度,如果要输入则输入eps即可. ⑵ 请输入由支路参数形成的矩阵B1 矩阵B1的每行是由下列参数构成的： ○1某支路的首端号P。 ○2末端号Q;且P143 可得到第一次迭代时的修正方程式。 （五）方程式，求11,iife。 （六）各节点电压，即得出第一次迭代后各节点的电压值。 （七）计算步骤迭代下去，当收敛精度取410（即PR=0.0001）时，需要进行四次迭代。求出了各节点电压后，即可求各支路的潮流分布。见图4－5 程序清单及打印结果： %本程序的功能是用牛顿－拉夫逊法进行潮流计算 n=input('请输入节点数:n='); nl=input('请输入支路数:nl='); isb=input('请输入平衡母线节点号:isb='); pr=input('请输入误差精度:pr='); B1=input('请输入由支路参数形成的矩阵:B1='); B2=input('请输入各节点参数形成的矩阵:B2='); 1.0500/21.84750 5.0000-j1.8134 ⑤ ③ 5.0000+ j1.4285 1.6-j0.8 0.8821\-4.77940+ ④ 2-j1 1.4157+j0.2442 -1.2775-j0.2033 ② 1.0634/-4.28200 ① 3.7-j1.3 2.5795- j1.9746 2.5795-j2.2995 1.0500/0.00000 图4－5 例4—1结果
144 Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n); O=zeros(1,n);S1=zeros(nl); for i=1:nl if B1(i,6)==0 p=B1(i,1); q=B1(i,2); else p=B1(i,2); q=B1(i,1); end Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); Y(q,p)=Y(p,q); Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; end %求导纳矩阵 disp('导纳矩阵Y='); disp(Y); G=rea

l(Y);B=imag(Y); for i=1:n e(i)=real(B2(i,3)); f(i)=imag(B2(i,3)); V(i)=B2(i,4); end for i=1:n S(i)=B2(i,1)-B2(i,2); B(i,i)=B(i,i)+B2(i,5); end P=real(S);Q=imag(S); ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; while IT2~=0 IT2=0;a=a+1; for i=1:n
145 if i~=isb C(i)=0; D(i)=0; for j1=1:n C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1); D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1); end P1=C(i)*e(i)+f(i)*D(i); Q1=f(i)*C(i)-D(i)*e(i); %求'P,Q' V2=e(i)^2+f(i)^2; if B2(i,6)~=3 DP=P(i)-P1; DQ=Q(i)-Q1; for j1=1:n if j1~=isb & j1~=i X1=-G(i,j1)*e(i)-B(i,j1)*f(i); X2=B(i,j1)*e(i)-G(i,j1)*f(i); X3=X2; X4=-X1; p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2; elseif j1==i & j1~=isb X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i); p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2; end
146 end else DP=P(i)-P1; DV=V(i)^2-V2; for j1=1:n if j1~=isb & j1~=i X1=-G(i,j1)*e(i)-B(i,j1)*f(i); X2=B(i,j1)*e(i)-G(i,j1)*f(i); X5=0; X6=0; p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2; elseif j1==i & j1~=isb X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X5=-2*e(i); X6=-2*f(i); p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2; end end end end end %求雅可比矩阵 for k=3:N0 k1=k+1;N1=N; for k2=k1:N1 J(k,k2)=J(k,k2)./J(k,k); end
147 J(k,k)=1; if k~=3 k4=k-1; for k3=3:k4 for k2=k1:N1 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2); end J(k3,k)=0; end if k==N0,break;end for k3=k1:N0 for k2=k1:N1 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2); end J(k3,k)=0; end else for k3=k1:N0 for k2=k1:N1 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2); end J(k3,k)=0; end end end for k=3:2:N0-1 L=(k+1)./2; e(L)=e(L)-J(k,N); k1=k+1;
148 f(L)=f(L)-J(k1,N); end for k=3:N0

DET=abs(J(k,N)); if DET>=pr IT2=IT2+1; end end ICT2(a)=IT2; ICT1=ICT1+1; end %用高斯消云法解"w=-J*V" disp('迭代次数'); disp(ICT1); disp('没有达到精度要求的个数'); disp(ICT2); for k=1:n V(k)=sqrt(e(k)^2+f(k)^2); O(k)=atan(f(k)./e(k))*180./pi; E(k)=e(k)+f(k)*j ; end disp('各节点的实部电压标么值E为(节点号从小到大排):'); disp(E); disp('各节点的电压大小V为(节点号从小到大排):'); disp(V); disp('各节点的电压相角O为(节点号从小到大排):'); disp(O); for p=1:n C(p)=0; for q=1:n
149 C(p)=C(p)+conj(Y(p,q))*conj(E(q)); end S(p)=E(p)*C(p); end disp('各节点的功率S为(节点号从小到大排):'); disp(S); disp('各条支路的首端功率Si为(顺序同您输入B1时一样):'); for i=1:nl if B1(i,6)==0 p=B1(i,1); q=B1(i,2); else p=B1(i,2); q=B1(i,1); end Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5)))); disp(Si(p,q)); end disp('各条支路的末端功率Sj为(顺序同您输入B1时一样):'); for i=1:nl if B1(i,6)==0 p=B1(i,1); q=B1(i,2); else p=B1(i,2); q=B1(i,1); end Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5)))); disp(Sj(q,p)); end disp('各条支路的功率损耗DS为(顺序同您输入B1时一样):'); for i=1:nl if B1(i,6)==0
150 p=B1(i,1); q=B1(i,2); else p=B1(i,2); q=B1(i,1); end DS(i)=Si(p,q)+Sj(q,p); disp(DS(i)); end 输入数据如下： 请输入节点数:n=5 请输入支路数:nl=5 请输入平衡母线节点号:isb=1 请输入误差精度:pr=0.00001 请输入由支路参数形成的矩阵:B1=[1 2 0.03i 0 1.05 0;2 3 0.08+0.3i 0.5i 1 0;2 4 0.1+0.35i 0 1 0;3 4 0.04+0.25i 0.5i 1 0;3 5 0.015i 0 1.05 1] 请输入各节点参数形成的矩阵:B2=[0 0 1.05 1.05 0 1;0 3.7+1.3i 1 0 0 2;0 2+1i 1 0 0 2;0 1.6+0.8i 1 0 0 2;5 0 1.05 1.05 0 3] 结果： 导纳矩阵Y= Columns 1 through 4 0 -33.3333i 0 +31.7460i 0 0 0 +31.7460i 1.5846 -35.7379i -0.8299 + 3.1120i -0.7547 + 2.6415i 0 -0.8299 + 3.1120i 1.4539 -66.9808i -0.6240 + 3.9002i 0 -0.7547 + 2.6415i -0.6240 + 3.9002i 1.3787 - 6.2917i 0 0 0 +63.4921i 0 Column 5 0 0 0 +63.4921i
151 0 0 -66.6667i 迭代次数 5 没有达到精度要求的个数 7 8 8 6 0 各节点的实部电压标么值E为(节点号从小到大排): Columns 1 through 4 1.0500 1.0335 - 0.0774i 1.0260 + 0.3305i 0.8592 - 0.0718i Column 5 0.9746 + 0.3907i 各节点的电压大小V为(节点号从小到大排): 1.0500 1.0364 1.0779 0.8622 1.0500 各节点的电压相角O为(节点号从小到大排): 0 -4.2819 17.8535 -4.7785 21.8433 各节点的功率S为(节点号从小到大排): Columns 1 through

4 2.5794 + 2.2994i -3.7000 - 1.3000i -2.0000 - 1.0000i -1.6000 - 0.8000i Column 5 5.0000 + 1.8131i 各条支路的首端功率Si为(顺序同您输入B1时一样): 2.5794 + 2.2994i -1.2774 + 0.2032i 0.1568 + 0.4713i 1.5845 + 0.6726i 5.0000 + 1.8131i 各条支路的末端功率Sj为(顺序同您输入B1时一样):
152 -2.5794 - 1.9745i 1.4155 - 0.2443i -0.1338 - 0.3909i -1.4662 - 0.4091i -5.0000 - 1.4282i 各条支路的功率损耗DS为(顺序同您输入B1时一样): -0.0000 + 0.3249 0.1381 - 0.0412i 0.0230 + 0.0804i 0.1184 + 0.2635i 0.0000 + 0.3849i 第四节 PQ分解法潮流计算 一 、PQ分解法的基本方程式 60年代以来N—R法曾经是潮流计算中应用比较普遍的方法，但随着网络规模的扩大（从计算几十个节点增加到几百个甚至上千个节点）以及计算机从离线计算向在线计算的发展，N—R法在内存需要量及计算速度方面越来越不 适应要求。70年代中期出现的快速分解法比较成功的解决了上述问题，使潮流计算在N—R法的基础上向前迈进了一大步，成为取代N—R法的算法之一。 快速分解法（又称P—Q分解法）是从简化牛顿法极坐标形式计算潮流程序的基础上提出来的。它的基本思想是根据电力系统实际运行特点：通常网络上的电抗远大于电阻值 ，则系统母线电压副值的微小变化V对母线有功功率的改变P影响很小。同样，母线电压相角的少许改变，也不会引起母线无功功率的明显改变Q。因此，节点功率方程在用极坐标形式表示时，它的修正方程式可简化为：
153 VVLHQP/00 （4—19） 这就是把2（n—1）阶的线性方程组变成了两个n—1阶的线性方程组，将P和Q分开来进行迭代计算，因而大大地减少了计算工作量。但是，H，L在迭代过程中仍然在不断的变化，而且又都是不对称的矩阵。对牛顿法的进一步简化（也是最关键的一步），即把（4—19）中的系数矩阵简化为在迭代过程中不变的对称矩阵。 在一般情况下，线路两端电压的相角ij是不大的（不超过10○～20○）。因此，可以认为： ijijijijBGsin1cos （4—20） 此外，与系统各节点无功功率相应的导纳BLDi远远小于该节点自导纳的虚部，即 iiiiLDiBVQB2 因而 iiiiBVQ2 （4—21） 考虑到以上关系，式（4—19）的系数矩阵中的各元素可表示为： ijjiijBVVH （i,j=1,2,………，n-1） （4—22） ijjiijBVVL （i,j=1,2,……………，m） (4—23) 而系数矩阵H和L则可以分别写成： 11,1122,1111,1111,222222121211,1121211111nnnnnnnnnnnnVBVVBVVBVVBVVBVVBVVBVVBVVBVH =1211,12,11,11,222211,11211121nnnnnnnnVVVBBBBBBBBBVVV

154 =11DDBVV （4—24） mmmmmmmmmmmmmVBVVBVVBVVBVVBVVBVVBVVBVVBVL22122222212121121211111 =mmmmmmmmVVVBBBBBBBBBVVV2121222211121121 =22''DDVBV （4—25） 将（4—24）和（4—25）式代入（4—19）中，得到 11DDVBVP VBVQD''2 用11DV和12DV分别左乘以上两式便得： 111'DDVBPV （4—26） VBQVD''12 （4—27） 这就是简化了的修正方程式，它们也可展开写成： 1122111,12,11,11,222211,11211112211nnnnnnnnnnVVVBBBBBBBBBVPVPVP （4—28）
155 mmmmmmmmmVVVBBBBBBBBBVQVQVQ212122221112112211 （4—29） 在这两个修正方程式中系数矩阵元素就是系统导纳矩阵的虚部，因而系数矩阵是对称矩阵，且在迭代过程中保持不变。这就大大减少了计算工作量。 用极坐标表示的节点功率增量为： njijijijijjiisinjijijijijjiisiBGVVQQBGVVPP110)cossin(0)sincos( （4—30） 式（4—28）、（4—29）和（4—30）构成了P—Q分解法迭代过程的基本方程式。 二、计算步骤和程序框图 （1） 给定各节点电压的初始值)0()0(,iiV； （2） 代入式（4—30）计算各节点有功功率iP，并求出iiVQ/； （3） 解修正方程式（4—28），得出各节点电压相角修正量i； （4） 修正各节点电压的相角i )()()1(kikiki （5）式（4—30）求得各节点无功功率误差iQ，并求出iiVQ/ （6）求解修正方程式（4—29），得出各节点电压幅值的修正量iV； （7）修正各节点电压的幅值iV， )()()1(kikikiVVV （8） 回（2）进行迭代，直到各节点功率误差iP及iQ都满足收敛条件P—Q分解法程序框图： 一、 程序清单及打印结果
156 是 置KQ=0 KP=0? 否 是 是 置KP=1 K+1K 计算平衡机节点功率及全部线路功率
出 Max?}{)(QkiQ 解修正方程（4—29）求)(kiV )()()1(kikikiVVV 否 否 解修正方程（4—28）求)(ki 用公式（4—30）计算不平衡功率)(kiP，计算)()(/kikiVP 输入原始数据 形成矩阵B’和B’’并进行三角分解 设PQ节点电压初值，各节点电压相角初值 置迭代计数k=0 Kp=1,kQ=1 ?}max{')(pkiP )()()1(kikiki置1QK 用公式（4—30）计算不平衡功率)(kiQ，计算)()(/kikiVQ 置Kp=0 KQ=0否 是
157 %本程序的功能是用P-Q分解法进行潮流计算 n=input('请输入节点数:n='); nl=input('请输入支路数:nl='); isb=input('请输入平衡母线节点号:isb='); pr=input('请输入误差精度:pr='); B1=input('请输入由支路参数形成的矩阵:B1='); B2=input('请输入由节点参数形成的矩阵:B2='); na=input('请输入PQ节点数na='); Y=zeros(n);YI=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);O=zeros(1,n); for i=1:nl if B1(i,6)==0 p=B1(i,1); q=B1(i,2); e

lse p=B1(i,2); q=B1(i,1); end Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); YI(p,q)=YI(p,q)-1./B1(i,3); Y(q,p)=Y(p,q); YI(q,p)=YI(p,q); Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; YI(q,q)=YI(q,q)+1./B1(i,3); Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; YI(p,p)=YI(p,p)+1./B1(i,3); end %求导纳矩阵 G=real(Y);B=imag(YI);BI=imag(Y); for i=1:n S(i)=B2(i,1)-B2(i,2); BI(i,i)=BI(i,i)+B2(i,5);
158 end P=real(S);Q=imag(S); for i=1:n e(i)=real(B2(i,3)); f(i)=imag(B2(i,3)); V(i)=B2(i,4); end for i=1:n if B2(i,6)==2 V(i)=sqrt(e(i)^2+f(i)^2); O(i)=atan(f(i)./e(i)); end end for i=2:n if i==n B(i,i)=1./B(i,i); else IC1=i+1; for j1=IC1:n B(i,j1)=B(i,j1)./B(i,i); end B(i,i)=1./B(i,i); for k=i+1:n for j1=i+1:n B(k,j1)=B(k,j1)-B(k,i)*B(i,j1); end end end end p=0;q=0;
159 for i=1:n if B2(i,6)==2 p=p+1;k=0; for j1=1:n if B2(j1,6)==2 k=k+1; A(p,k)=BI(i,j1); end end end end for i=1:na if i==na A(i,i)=1./A(i,i); else k=i+1; for j1=k:na A(i,j1)=A(i,j1)./A(i,i); end A(i,i)=1./A(i,i); for k=i+1:na for j1=i+1:na A(k,j1)=A(k,j1)-A(k,i)*A(i,j1); end end end end ICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1; while ICT2~=0|ICT3~=0 ICT2=0;ICT3=0;
160 for i=1:n if i~=isb C(i)=0; for k=1:n C(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+BI(i,k)*sin(O(i)-O(k))); end DP1(i)=P(i)-V(i)*C(i); DP(i)=DP1(i)./V(i); DET=abs(DP1(i)); if DET>=pr ICT2=ICT2+1; end end end Np(K)=ICT2; if ICT2~=0 for i=2:n DP(i)=B(i,i)*DP(i); if i~=n IC1=i+1; for k=IC1:n DP(k)=DP(k)-B(k,i)*DP(i); end else for LZ=3:i L=i+3-LZ; IC4=L-1; for MZ=2:IC4 I=IC4+2-MZ;
161 DP(I)=DP(I)-B(I,L)*DP(L); end end end end for i=2:n O(i)=O(i)-DP(i); end kq=1;L=0; for i=1:n if B2(i,6)==2 C(i)=0;L=L+1; for k=1:n C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k))); end DQ1(i)=Q(i)-V(i)*C(i); DQ(L)=DQ1(i)./V(i); DET=abs(DQ1(i)); if DET>=pr ICT3=ICT3+1; end end end else kp=0; if kq~=0; L=0; for i=1:n if B2(i,6)

==2 C(i)=0;L=L+1;
162 for k=1:n C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k))); end DQ1(i)=Q(i)-V(i)*C(i); DQ(L)=DQ1(i)./V(i); DET=abs(DQ1(i)); end end end end Nq(K)=ICT3; if ICT3~=0 L=0; for i=1:na DQ(i)=A(i,i)*DQ(i); if i==na for LZ=2:i L=i+2-LZ; IC4=L-1; for MZ=1:IC4 I=IC4+1-MZ; DQ(I)=DQ(I)-A(I,L)*DQ(L); end end else IC1=i+1; for k=IC1:na DQ(k)=DQ(k)-A(k,i)*DQ(i); end
163 end end L=0; for i=1:n if B2(i,6)==2 L=L+1; V(i)=V(i)-DQ(L); end end kp=1; K=K+1; else kq=0; if kp~=0 K=K+1; end end end disp('迭代次数'); disp(K); disp('每次没有达到精度要求的有功功率个数为'); disp(Np); disp('每次没有达到精度要求的无功功率个数为'); disp(Nq); for k=1:n E(k)=V(k)*cos(O(k))+V(k)*sin(O(k))*j ; O(k)=O(k)*180./pi; end disp('各节点的电压标么值E为(节点号从小到大排):');
164 disp(E); disp('各节点的电压大小V为(节点号从小到大排):'); disp(V); disp('各节点的电压相角O为(节点号从小到大排):'); disp(O); for p=1:n C(p)=0; for q=1:n C(p)=C(p)+conj(Y(p,q))*conj(E(q)); end S(p)=E(p)*C(p); end disp('各节点的功率S为(节点号从小到大排):'); disp(S); disp('各条支路的首端功率Si为(顺序同您输入B1时一样):'); for i=1:nl if B1(i,6)==0 p=B1(i,1); q=B1(i,2); else p=B1(i,2); q=B1(i,1); end Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5)))); disp(Si(p,q)); end disp('各条支路的末端功率Sj为(顺序同您输入B1时一样):'); for i=1:nl if B1(i,6)==0 p=B1(i,1); q=B1(i,2); else p=B1(i,2); q=B1(i,1);
165 end Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5)))); disp(Sj(q,p)); end disp('各条支路的功率损耗DS为(顺序同您输入B1时一样):'); for i=1:nl if B1(i,6)==0 p=B1(i,1); q=B1(i,2); else p=B1(i,2); q=B1(i,1); end DS(i)=Si(p,q)+Sj(q,p); disp(DS(i)); end 输入数据为： 请输入节点数:n=5 请输入支路数:nl=5 请输入平衡母线节点号:isb=1 请输入误差精度:pr=0.0001 请输入由支路参数形成的矩阵:B1=[1 2 0.03i 0 1.05 0;2 3 0.08+0.3i 0.5i 1 0;2 4 0.1+0.35i 0 1 0;3 4 0.04+0.25i 0.5i 1 0;3 5 0.015i 0 1.05 1]

请输入由节点参数形成的矩阵:B2=[0 0 1.05 1.05 0 1;0 3.7+1.3i 1.05 0 0 2;0 2+1i 1.05 0 0 2;0 1.6+0.8i 1.05 0 0 2;5 0 1.05 1.05 0 3] 请输入PQ节点数na=3 结果： 迭代次数 8 每次没有达到精度要求的有功功率个数为
166 4 4 3 4 4 4 3 0 每次没有达到精度要求的无功功率个数为 3 3 3 3 3 3 1 0 各节点的电压标么值E为(节点号从小到大排): Columns 1 through 4 1.0500 1.0335 - 0.0774i 1.0260 + 0.3305i 0.8592 - 0.0718i Column 5 0.9746 + 0.3907i 各节点的电压大小V为(节点号从小到大排): 1.0500 1.0364 1.0779 0.8622 1.0500 各节点的电压相角O为(节点号从小到大排): 0 -4.2818 17.8529 -4.7781 21.8427 各节点的功率S为(节点号从小到大排): Columns 1 through 4 2.5794 + 2.2993i -3.6999 - 1.3000i -2.0001 - 1.0000i -1.6000 - 0.8000i Column 5 5.0000 + 1.8130i 各条支路的首端功率Si为(顺序同您输入B1时一样): 2.5794 + 2.2993i -1.2773 + 0.2031i 0.1568 + 0.4713i 1.5845 + 0.6725i 5.0000 + 1.8130i 各条支路的末端功率Sj为(顺序同您输入B1时一样): -2.5794 - 1.9744i 1.4154 - 0.2443i -0.1338 - 0.3909i
167 -1.4662 - 0.4091i -5.000 – 1.4281i 各条支路的功率损耗DS为(顺序同您输入B1时一样): 0.0000 + 0.3249i 0.1381 - 0.0412i 0.0230 + 0.0804i 0.1184 + 0.2634i -0.0000 + 0.3849i 习 题 4—1．电力系统如图4—1所示。 各元件阻抗标么值为z10=－j30 ,z20=－j34,z30=－j29, z12=0.08+j0.40,z23=0.10+j0.40, z13=0.12+j0.50,z34=j0.30。 设节点1平衡节点，节点4 为PV节点，节点2、3为 PQ节点。若给定：V1=1.05， S2=0.55+j0.13S3=0.3+j0.18，P4=0.5，V4=1.10， 试求：用牛顿—拉夫逊法计算 的潮流分布。（采用计算机编 程方法） 4—2．用P—Q分解法对题4—1 电力系统进行潮流计算。（采用 计算机编程方法） 4—3．如图4—2所示,已知参数 的标么值如下：R12=0.02， R23=0.04，R34=0.04，R14=0.01， Z12 Z13 4 3 2 1 Z23 Z30 Z34 S2 S3 P4 Z20 Z10 0.9：1 图4—1 4 1 2 R14 R34 3 S2 S3 S1 图4—2 R12 R23
168 S1=0.3，S2=－0.2，S3=0.15， 节点4为平衡节点，V4=1.0， 试用牛顿—拉夫逊法作潮流 计算（采用计算机编程方法） 4—4．简单电力系统如图4—3所示，已知 各段线路阻抗和节点功率为：Z12=10+j16， Z13=13.5+j21，Z23=24+j22， SLD2=20+j15MVA，SLD3=25+j18MVA， 节点1为平衡节点，V1=1150KV， 试用牛顿—拉夫逊法计算潮流。 （采用计算机编程方法） 1 3 2 SLD2 110KV 图4—3 SLD3