电力系统分析实验指导书

电力系统分析

实验指导书

郑州科技学院电气工程及其自动化教研室编

目 录

实验一 生成节点导纳矩阵的matlab程序设计 ................................................................................. 1 实验二 支路追加法生成节点阻抗矩阵的matlab程序设计 ........................................................ 5 实验三 牛顿-拉夫逊法电力系统潮流计算的matlab程序设计 ........................................................ 13 实验四 简化模型的静态稳定计算的matlab程序设计 .......................................................... 24 实验五 短路电流计算程序的实现 ....................................................................................................... 31

实验一 生成节点导纳矩阵的matlab程序设计

一、实验目的：

本实验通过对生成节点导纳矩阵的编程与调试，获得复杂电力系统的节点导纳矩阵，为电力系统的潮流计算程序编制打下基础。通过实验教学加深学生对节点导纳矩阵计算方法的理解，掌握变压器为非标准变比时的修正方法，熟悉MATLAB软件的使用方法，提高编制调试计算机程序的能力。

二、实验器材：

计算机、软件（已安装应用软件MATLAB等）、移动存储设备（学生自备，软盘、U盘等）。 三、实验内容：

1.节点导纳矩阵的形成

节点导纳矩阵的对角元Yii(i?1,2,...n)称自导纳，自导纳Yii是节点i以外的所有节点都接地，

??1的电压）时，由节点由节点i向整个网络看到而得到的导纳。在节点i加上单位大小的电压（Uii流向网络的电流就等于节点i的自导纳

?IYii?i?Ui

??0,j?iUj通过计算不难发现，Yii就等于与节点i连接的所有支路导纳的和。

节点导纳矩阵的非对角元Yij(i?1,2,...n；j?1,2,...n；i?j)称互导纳，即把节点j以外的节点全接地，而在节点j加以单位电压时，由节点i流向j的电流加上负号的值。

?IYij?i?Uj

??0,k?jUk同样，Yij其实就是连接节点i和节点j支路的导纳值再加上负号。显然，Yij?Yji。而且，如节点i，j之间没有直接联系，也不计两支路之间，例如两相邻电力线路之间的互感时， Yij?Yji?0。互导纳的这些性质决定了节点导纳矩阵是一个对称的稀疏矩阵。而且，由于每个节点所连接的

支路数总有一定限度，随着网络中节点数的增加，非零元素数相对愈来愈少，节点导纳矩阵的稀疏度，即零元素数与总元素数的比值也就愈来愈高。

导纳矩阵的计算归结如下：

(1) 导纳矩阵的阶数等于电力系统网络的节点数。

(2) 导纳矩阵各行非对角元素中非零元素的个数等于对应节点所连的不接地支路数。 (3) 导纳矩阵的对角元素，即各节点的自导纳等于相应节点所连支路的导纳之和。

Yii??yij

j?i式中，yij为节点i与节点j间支路阻抗Zij倒数，符号j?i表示j属于i或与i连接的j，即?内只包括与节点i直接相连的节点j。当节点i有接地支路时，还应包括j?0的情况。

(4)导纳矩阵非对角元素Yij等于节点i与节点j之间的导纳的负数

Yij??1??yij zijL当i、j之间有多条并联支路时，求Yij时应求所有并联支路导纳的代数和的负数

(s) Yij???yijs?1式中，L表示i、j之间并联支路数的条数。

按照以上计算式，对于实际网络均可根据给定的支路参数和连接情况，直观而简单地求出导纳

1

矩阵。可以看出，用以上计算公式求得的导纳矩阵与根据定义得到的导纳矩阵是完全一致的。

2.变压器为非标准变比时的修正

无论采用有名制或标么制，凡涉及多电压级网络的计算，在精确计算时都必须将网络中所有参数和变量按实际变比归算到同一电压等级。实际上，在电力系统计算中总是有些变压器的实际变比不等于变压器两侧所选电压基准值之比，也就是不等于标准变比，而且变压器的变比在运行中是可以改变的。这将使每改变一次变比都要重新计算元件参数，很不方便。下面将介绍另一种可等值地体现变压器电压变换功能的模型。

等值变压器模型

Z??I11?ZT

1:KI22?Z??U1?U2

?图1 非标准变比时的修正电路

??I11KZTI22U1?KZTK?1K2ZT1?KU2

? 图2 以变压器阻抗表示

增加非零非对角元素为

Yij?Yji??节点i的自导纳，增加一个改变量为

YT KYii?K?11YT?YT?YT KK节点j的自导纳，也增加一个改变量为

Yjj?11?K1YT?Y?YT TKK2K23. 程序代码及说明

n=input('请输入节点数:n='); nl=input('请输入支路数:nl=');

%isb=input('请输入平衡母线节点号:isb=');

2

%pr=input('请输入误差精度:pr=');

B1=input('请输入由各支路参数形成的矩阵:B1='); 2=input('请输入各节点参数形成的矩阵:B2=');

X=input('请输入由节点号和接地支路参数形成的矩阵:X=');

Y=zeros(n); e=zeros(1,n);f=zeros(1,n);V=zeros(1,n); O=zeros(1,n);S1=zeros(nl); % zeros(m,n)产生m×n的零矩阵，zeros(n)产生n×n的全0方阵。 for i=1:n %节点数

if X(i,2)~=0; %matlab关系符号：>（大于），>=（大于等于），<（小于），<=（小于等于）， ==（等于）~=（不等于） p=X(i,1);

Y(p,p)=1./X(i,2); %接地支路, “./”点除代表矩阵对应元素相除 end end

for i=1:nl %支路数

if B1(i,6)==0 %折算到哪一侧的标志,0非标准变比在q侧,1非标准变比在在p侧 p=B1(i,1);q=B1(i,2); % B1(i,1), 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));

%计算节点p和q间的互导纳,B1(i,3)为支路阻抗,B1(i,5)为变压器支路变比k, B1(i,5)=1时不是变压器支路

Y(q,p)=Y(p,q);

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %对角元, 节点q的自导纳, B1(i,4)为支路对地容抗

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元,节点p的自导纳 end

%输出导纳矩阵 disp('导纳矩阵Y='); disp(Y);

四、实验数据：

电力系统接线图如图3所示，根据图中数据可得B1和X矩阵：

P5=5j0.0151:1.05③j0.252+j1j0.25j0.25④1.6+j0.80.04+j0.250.1+j0.350.8+j0.30j0.253.7+j1.3②1.05:1j0.03①U1=1.05δ1=0°⑤U5=1.05

3

图3 电力系统接线图 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]; X=[1 0;2 0;3 0;4 0;5 0];

五、实验过程

1.手算图3所示网络的导纳矩阵。

2.认真阅读源程序，建立matlab的M程序，输入节点导纳矩阵计算程序。 3.输入实验数据,用B1和X矩阵代替B1和X输入语句并调试程序。 B1矩阵

B1(i,1)为支路始端编号p；

B1(i,2)为支路末端编号q，且p

B1(i,4)为支路对地容抗，无填0；

B1(i,5)为变压器支路变比，非变压器支路填1；

B1(i,6)为变压器变比折算到哪一侧的标志,0非标准变比在q侧,1非标准变比在在p侧，非变压器支路填0； X矩阵

X(i,1)为节点编号；

X(i,2)为该节点的接地支路阻抗，如果节点没有接地支路填”0”；

4.正确输入节点数和支路数进行计算，将计算结果与手算结果进行比较。 六、要求

认真书写实验报告并上交实验数据

4

实验二 支路追加法生成节点阻抗矩阵的matlab程序设计

一、实验目的：

本实验通过对生成节点阻抗矩阵的编程与调试，获得复杂电力系统的节点阻抗矩阵，为采用阻抗法的电力系统的短路电流计算程序编制打下基础。通过实验教学加深学生对节点阻抗矩阵计算方法的理解，掌握支路追加法生成节点阻抗矩阵的方法，熟悉MATLAB软件的使用方法，提高编制调试计算机程序的能力。 二、实验器材：

计算机、软件（已安装应用软件MATLAB等）、移动存储设备（学生自备，软盘、U盘等）。 三、实验内容：

1.节点阻抗矩阵的形成

目前求节点阻抗矩阵的方法，归纳起来主要有两种：一种是导纳矩阵求逆，间接求出阻抗矩阵；另一种是用支路追加法，直接形成节点阻抗矩阵。

支路追加法形成节点阻抗矩阵可以从网络的接线图中直接形成，它从某一个与地相连的支路开始，以后每次追加一条支路，直至形成整个网络的节点阻抗矩阵。该方法矩阵形成的规律性很强，易于理解和记忆。用支路追加法形成节点阻抗矩阵时，所追加的支路可以分成4类，即追加接地树支，追加树支，追加接地连支，追加连支。现以图1所示的含有3个独立节点的网络为例，别就上述4种情况阐述节点阻抗的形成。先将节点的电压、电流关系表示如下：

???Z?U111?????U2???Z21?U???Z?3??31Z12Z22Z32??Z13??I1????Z23??I2? ???Z33???I3?

图1 原始网络

(1)追加接地支路（0，4）

如图2所示，由于网络增加了一个节点4，所以原网络矩阵增加一阶，但是由于新增节点4与原网络矩阵增加一阶，但由于新增节点4与原网络3个节点1、2、3无电气直接联系，所以原来3

??zI?,其中z是新增支路的阻抗。这样节点方程修改为式各方程不收影响，只是新增一个方程U44（1），即节点阻抗矩阵原有的各元素均不变，新增的行、列元素均为零，只有新增的对角元素为z。

图2 追加接地树支

5

???Z?U1????11?U2??Z21??U???Z31?3??0???U4???(2)追加树支（2，4）

Z12Z22Z320??Z130??I1???IZ230???2? （1） ??Z330??I3???0z??I??4??

图3 追加树支

?，从如图3所示。由于网络增加了一个节点4，所以矩阵增加一阶，设节点4的注入电流为I4??I?,而其他节点注入电流不变，而新增了一个方程：原网络看进去，节点2的注入电流变为I24??U??zI?,所以节点方程改为 U424??ZI????????U1111?Z12(I2?I4)?Z13I3?Z11I2?Z12I2?Z13I3?Z12I4??ZI??Z(I??I?)?ZI??ZI??ZI??ZI??ZI?U22112224233212222233224??ZI????????U3311?Z32(I2?I4)?Z33I3?Z31I2?Z32I2?Z33I3?Z32I4??U??zI??ZI??ZI??ZI??(Z?z)I?U424211223233224

写成矩阵形式（2）。

???Z?U111?????U2??Z21??U???Z31?3??Z???U4???21Z12Z22Z32Z22Z13Z23Z33Z23??Z12??I1????Z22??I2? (2) ??Z32??I3???Z22?z??I??4??由上式可以看出，追加树支时，矩阵增加一阶，新增的行、列元素分别等于树支所接的原网络

节点2所对应的行、列元素，新对角元素等于树支所接的节点2的对角元素加上新增支路的阻抗值。

电力网络中包含有许多变压器。在追加变压器支路时，一般用一个等值阻抗同一个理想变压器相串联的支路来表示，具体可分为追加树支和追加连支两种情况。下面先讨论追加树支的情况，如图4所示，追加变压器树支（2，4），网络增加了一个节点4，

6

图4 追加变压器树支

?，从原网络看，节点2的注入电流为（I??KI?）设其注入电流为I424，所以有

??ZI????U1111?Z12(I2?KI4)?Z13I3??ZI????U2211?Z22(I2?KI4)?Z23I3 ??ZI??Z(I??KI?)?ZI?U33113224333另外还有

??K(U??KzI?)U424

??????U4?K?Z21I1?Z22(I2?KI4)?Z23I3?KzI4????整理后，有 U4???Z?U111?????U2??Z21??U???Z31?3??KZ???U4???21Z12Z22Z32KZ22Z13Z23Z33KZ23??KZ12??I1???IKZ22???2? ??KZ32??I3???2K(Z22?z)??I??4??由此可见，追加变压器树支和追加普通树支支路相似，只是在新增行、列的元素分别乘以变比K，新对角元素乘以变比K。

(3)追加接地连支（0，2）

2

图5 追加接地连支

如图5所示，由于追加接地连支，网络的节点数没有变化，故矩阵的阶次不变。对原网络来说，

???I??I?，其他节点注入电流不变，则各节点电压方程变为 节点2的注入电流变为I227

??ZI????????U1111?Z12(I2?I)?Z13I3?Z11I1?Z12I2?Z13I3?Z12I??ZI????????U2211?Z22(I2?I)?Z23I3?Z21I1?Z22I2?Z23I3?Z22I??ZI????????U3311?Z32(I2?I)?Z33I3?Z31I1?Z32I2?Z33I3?Z32I??zI???ZI????0??U2211?Z22I2?Z33I3?(Z22?z)I写成矩阵形式为

???Z?U111?????U2???Z21?U???Z31?3????0????Z21Z12Z22Z32?Z22Z13Z23Z33?Z23???Z12??I1????I?Z22??2? (3) ???Z32??I3?????Z22?z??I??可见，矩阵可暂时增加一阶，原矩阵元素不变，新行、列元素分别等于该追加连支的非零节点所对应的行列元素的负值；新对角元素等于该点的自阻抗加上连支阻抗z。形成了暂时增加一阶的节点阻抗矩阵以后，用高斯消去法消去矩阵的暂增行与列，原矩阵的元素Zij为

??Zij?Zij(4) 追加连支（2，3）

Zi2Z2jZ44(i,j?1,2,3) (4)

图6 追加连支

如图6所示，由于追加连支，所以矩阵阶次不变，与采用追加接地连支同样的处理方法，矩阵

?，由节点3流向节点2，从原网络看进去，节点2的注入电流为可暂时增加一阶。设连支电流为I??I?)，节点3的注入电流变为(I??I?)。则节点电压方程为 (I23?????UZ11Z12Z13(Z12?Z13)??I1?1????????IZ21Z22Z23(Z22?Z32)?U2?????2?

?U???????IZZZ(Z?Z)31323332333?3?????????0???(Z21?Z31)(Z22?Z32)(Z23?Z33)(Z22?Z33?Z23?Z32?z)???I? (5)

可见，追加连支后，原矩阵的元素不变，暂时增加的行、列元素分别等于该追加连支的两个节点所对应的行元素之差和列元素之差；新增对角元素为这两个节点自阻抗之和减去相互间的互阻抗

8

之和加上该连支阻抗。

形成了暂时增加一阶的节点阻抗矩阵之后，用高斯消去法消去暂增行、列，即得追加连支的节点阻抗矩阵，消元公式同式（4）。

图7 追加变压器连支

?,则由变压器支路流入 如图7，追加变压器连支（2，3）。由节点3流入理想变压器的电流为I?，从原网络看节点2的注入电流变为（I??KI?）节点2的电流为KI，节点3的注入电流变为2??I?），则： （I3??ZI?????I?)U（I1111?Z12(I2?KI)?Z133??ZI??Z(I??KI?)?Z(I??I?) U2211222233??ZI?????U3311?Z32(I2?KI)?Z33(I3?I)??K(U??zKI?) 另有： U32??U??zK2I??0 即，KU23?，U?代入整理后得 把U2322??(KZ?Z)I???(KZ21?Z31)I122322?(KZ23?Z33)I3?(KZ22?KZ32?KZ23?Z33?Kz)I?0

所以得

??????UZ11Z12Z13KZ12?Z13??I11????????Z21Z22Z23KZ22?Z23?U2?????I2?

?U???????IZ31Z32Z33KZ32?Z333?3????22????KZ21?Z31)(KZ22?Z32)(KZ23?Z33)(KZ22?KZ33?KZ23?KZ32?Kz)???0??（?I??前面已经讨论的都是变压器的漏抗归算至低压侧。如果变压器的漏抗归算至高压侧，则只需将变比变为1/K即可。

2. 程序代码及说明 n=input('请输入节点数:n='); nl=input('请输入支路数:nl=');

B=input('请输入由各支路参数形成的矩阵:B=');

9

B=[0 1 2i 0 1 0;0 2 4i 0 1 0;1 2 2i 0 1 0;0 3 20i 0 1 0;2 3 8i 0 1 0;1 3 5i 0 1 0];

m=0;Z=zeros(n); % zeros(n)产生n×n的全0方阵。

for k1=1:nl %支路数

p=B(k1,1);q=B(k1,2); % B(k1,1), B(k1,2)为支路编号 if B(k1,6)==0 % k=1./B(k1,5); else

k=B(k1,5); end

if p==0 %

if q>m % 追加接地树支 Z(q,q)=B(k1,3);m=m+1; else

for i=1:m

Z(i,m+1)=-Z(i,q);Z(m+1,i)=-Z(q,i); end

Z(m+1,m+1)=Z(q,q)+B(k1,3); for i=1:m for j=1:m

Z(i,j)=Z(i,j)-Z(i,m+1)*Z(m+1,j)./Z(m+1,m+1); end

Z(i,m+1)=0; end

for i=1:m+1 Z(m+1,i)=0; end end

else if q>m %追加不接地树支 for i=1:m

Z(i,q)=Z(i,p)*k;Z(q,i)=Z(p,i)*k; end

Z(q,q)=k^2*Z(p,p)+k^2*B(k1,3); m=m+1; else

for i=1:m

Z(i,m+1)=Z(i,p)*k-Z(i,q); Z(m+1,i)=Z(p,i)*k-Z(q,i); end

Z(m+1,m+1)=k^2*Z(p,p)+Z(q,q)-2*k*Z(p,q)+k^2*B(k1,3); for i=1:m

10

for j=1:m

Z(i,j)=Z(i,j)-Z(i,m+1)*Z(m+1,j)./Z(m+1,m+1); end

Z(i,m+1)=0; end

for i=1:m+1 Z(m+1,i)=0; end end end end

%输出节点阻抗矩阵

disp('节点阻抗矩阵Z='); disp(Z)

四、实验数据：

电力系统接线图如图8所示，根据图中数据可得B矩阵：

j2①j5③j20j2②j4j8 图8 电力系统接线图 B=[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];

五、实验过程

1.手算图3所示网络的节点阻抗矩阵。

2.认真阅读源程序，建立matlab的M程序，输入节点阻抗矩阵计算程序。 3.输入实验数据,用B矩阵代替B输入语句并调试程序。 B矩阵

B (i,1)为支路始端编号p；

B (i,2)为支路末端编号q，且p

B (i,4)为支路对地电纳，无填0；

B (i,5)为变压器支路变比，非变压器支路填1；

B (i,6)为变压器变比折算到哪一侧的标志,0非标准变比在q侧,1非标准变比在在p侧，非变压器支路填0；

11

4.正确输入节点数和支路数进行计算，将计算结果与手算结果进行比较。 六、要求

认真书写实验报告并上交实验数据

12

实验三 牛顿-拉夫逊法电力系统潮流计算的matlab程序设计

一、实验目的：

潮流计算是在给定电力系统网络结构、参数和决定系统运行状态的边界条件的情况下确定系统稳态运行状态的一种基本方法,它是电力系统分析中最基本、 最重要的计算,是系统安全、经济分析和实时控制与调度的基础。常规潮流计算的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。潮流计算的结果是电力系统稳定计算和故障分析的基础。本实验通过对电力系统潮流计算的计算机程序的编制与调试，获得对复杂电力系统进行潮流计算的计算机程序，使系统潮流计算能够由计算机自行完成，即根据已知的电力网的数学模型（节点导纳矩阵）及各节点参数，由计算程序运行完成该电力系统的潮流计算。通过实验教学加深学生对复杂电力系统潮流计算计算方法的理解，学会运用电力系统的数学模型，掌握潮流计算的过程及其特点，熟悉MATLAB软件的使用方法，提高编制调试计算机程序的能力，学习如何将理论知识和实际工程问题结合起来。 二、实验器材：

计算机、软件（已安装应用软件MATLAB等）、移动存储设备（学生自备，软盘、U盘等）。 三、实验内容：

1.采用直角坐标的牛顿-拉夫逊法潮流计算原理

牛顿-拉夫逊（Newton-Raphson）法是数学上解非线性方程式的有效方法，有较好的收敛性，将N-R（又称牛顿法）法用于潮流计算是以导纳矩阵为基础的，由于利用了导纳矩阵的对称性、稀疏性及节点编号顺序优化等技巧，使N-R法在收敛性、占用内存、计算速度等方面都有较大优势，成为20世纪60年代末期以后普遍采用的方法。

根据电力系统的实际运行条件，按给定变量的不同，一般将节点分为以下三种类型。 (1)PQ节点

这一类节点，事先给定的是节点功率（P、Q），待求的未知量是节点电压向量（U、θ），所以叫“PQ节点”。通常变电所母线都是PQ节点，当某些发电机的输出功率P、Q给定时，也作为PQ节点。PQ节点上的发电机称之为PQ机（或PQ给定型发电机）。在潮流计算中系统大部分节点属于PQ节点。

网络中还有一类既不接发电机、又没有负荷的联络节点（亦称浮动节点），也可以当作PQ节点，其PQ给定值为零。

(2)PV节点

这类节点给出的参数是该节点的有功功率P及电压幅值U，待求量为该节点的无功功率Q及电压向量的相角θ。这类节点在运行中往往要有一定可调节的无功电源，用以维持给定的电压值。通常选择有一定无功功率储备的发电机母线或变电所有无功补偿设备的母线作PV节点处理。PV节点上的发电机称之为PV机（或PV给定型发电机）。在电力系统中，这一类节点的数目很少。

(3)平衡节点

在潮流计算中，这类节点一般只设一个。对该节点，给定其电压值，并在计算中取该节点电压向量的方向作为参考轴，相当于给定该点电压向量的角度为零。也就是说，对平衡节点给定的运行参数是U和?，因此又称为U?节点，而待求量是该节点的P、Q，整个系统的功率平衡由这一节点承担。关于平衡点的选择，一般选择系统中担任调频调压的某一发电厂（或发电机），有时也可能按其他原则选择，例如，为提高计算的收敛性，可以选择出线数多或者靠近电网中心的发电厂作母线平衡点。

节点电压和导纳可表示为

13

???Ui?ei?jfi ???Yij?Gij?jBij将上式代入Pi?jQi?Ui*?YUijj?1n?j的右端，展开并分出实部和虚部，便得

nn?Pi?ei?(Gijej?Bijfj)?fi?(Gijfj?Bijej)??j?1j?1 （1） ?nn?Qi?fi?(Gijej?Bijfj)?ei?(Gijfj?Bijej)?j?1j?1?按照上节的分类，PQ节点的有功功率和无功功率是给定的，第i个节点的给定功率设为Pis和

··· ，m节点为PQ节点，对其中每一个节点可列方程： Qis。假定系统中的第1，2，

nn??Pi?Pis?Pi?Pis?ei?(Gijej?Bijfj)?fi?(Gijfj?Bijej)?0??j?1j?1 ?nn??Qi?Qis?Qi?Qis?fi?(Gijej?Bijfj)?ei?(Gijfj?Bijej)?0?j?1j?1? (i=1，2，?，m) (2) PV节点的有功功率和节点电压幅值是给定的。假定系统中的第m+1，m+2，···，n-1号节点为PV节点，则对其中每一节点可以列方程：

nn???Pi?Pis?Pi?Pis?ei?(Gijej?Bijfj)?fi?(Gijfj?Bijej)?0 ?j?1j?1??U2?U2?U2?U2?(e2?f2)?0isiisii?i (i=m+1，m+2，···，n-1） (3)

第n号节点为平衡点，其电压Un?en?jfn是给定的，故不参加迭代。

式（2)和（3）总共包含了2(n-1)个方程，待求的变量有e1,f1,e2,f2,?,en?1,fn?1也是2（n-1）个。对其近似解的修正量，可以由方程式(4)来确定：

?W??J?U (4） 式中，

???P??e1?1???Q???f?1???1????????????P?em???m???Qm???fm??W??? ?U???

?P?e?m?1??m?1???U2???f?m?1???m?1?????????P???e?n?1???n?1?2????Un?1????fn?1??14

???P1??e1????Q1??e1??????Pm??e1???Qm??eJ????P1?m?1??e1?2???Um?1??e??1????Pn?1??e1?2???Un?1???e1??P1?f1??Q1?f1???Pm?f1??Qm?f1??Pm?1?f12??Um?1?f1???Pn?1?f12??Un?1?f1......?..................??P1?em??Q1?em???Pm?em??Qm?em??Pm?1?em2??Um?1?em???Pn?1?em2??Un?1?em??P1?fm??Q1?fm???Pm?fm??Qm?fm??Pm?1?fm2??Um?1?fm???Pn?1?fm2??Un?1?fm??P1?em?1??Q1?em?1???Pm?em?1??Qm?em?1??Pm?1?em?12??Um?1?em???Pn?1?em?12??Un?1?em??P1?fm?1??Q1?fm?1???Pm?fm?1??Qm?fm?1??Pm?1?fm?12??Um?1?fm?1???Pn?1?fm?12??Un?1?fm?1........................??P1?en?1??Q1?en?1???Pm?en?1??Qm?en?1??Pm?1?en?12??Um?1?en?1???Pn?1?en?1??U2n?1?en?1??P?1?fn?1????Q1??fn?1??????Pm??fn?1???Qm???fn?1???Pm?1? ?fn?1????U2m?1??fn?1??????Pn?1??fn?1????U2n?1??fn?1??上述方程中雅可比矩阵的各元素，可以对式（2)和（3）求偏导数获得。当j=i时，对角元素是

n???Pi??e???(Gijej?Bijfj)?Gijej?Biifij?1?in???Pi??(Gf?Be)?Be?Gf?ijjijjiiiiii??fj?1?in???Qi?(Gf?Be)?Be?Gf?ijjijjiiiiii??eij?1?n??Qi????(Gijej?Bijfj)?Giiei?Biifi??fij?1 （5） 2???Ui???2ei?e?i???Ui2??2fi???fi当j?i时，矩阵中非对角元素是

???P??Qii????(Gijei?Bijfi)??fj??ej???P??Q?ii??Bijei?Gijfi （6） ??f?ej?j???U2i??U2i???0?fj???ej由以上表达式不难看出，雅可比矩阵有以下特点：

（1）雅可比矩阵中的诸元素都是节点电压的函数，因此在迭代过程中，它们将随着各节点电压的

变化而不断地改变。 （2）矩阵是不对称的。

由式（6）可以看出，当导纳矩阵中的非对角元素Yij为零时，雅可比矩阵中相对应的元素也是零，即矩阵是非常稀疏的。因此，修正方程的求解同样可以应用稀疏矩阵的求解技巧。正是由于这

15

一点才使牛顿-拉夫逊法获得广泛的应用。

用牛顿法计算潮流时，步骤如下：

（1）给出各节点电压初始值ei、fi(0)(0)。,

(0)（2）将以上电压初始值代入式（2）和（3），求出修正方程式的常数项向量?Pi、?Qi(0)、

?U2i素。

(0)。

（3）将电压初始值代入式（5）和（6），求出修正方程式中系数矩阵（雅可比矩阵）的各元（4）解修正方程式（4），求出修正量?ei（5）修正各节点电压

(0)，?fi(0)。

ei（6）将ei、fi(1)(1)?ei(0)??ei，fi(1)?fi(0)??fi(0)

(1)(0)(1)代入式（2）和（3），求?Pi，?Qi(1)，?U2i。

(1)（7）校验是否收敛，即

如果收敛，迭代到此结束，进一步计算各线路潮流和平衡节点功率，并打印输出结果；如果不收敛，转回第（2）步进行下一次迭代计算，直到收敛为止。

2. 程序代码及说明 n=input('请输入节点数:n='); nl=input('请输入支路数:nl=');

isb=input('请输入平衡母线节点号:isb='); pr=input('请输入误差精度:pr=');

B1=input('请输入由各支路参数形成的矩阵:B1='); B2=input('请输入各节点参数形成的矩阵:B2=');

X=input('请输入由节点号和接地支路参数形成的矩阵:X=');

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:n %节点数 if X(i,2)~=0; p=X(i,1);

Y(p,p)=1./X(i,2); %接地支路, “./”点除代表矩阵对应元素相除 end end

for i=1:nl %支路数

if B1(i,6)==0 %折算到哪一侧的标志,0非标准变比在q侧,1非标准变比在在p侧 p=B1(i,1);q=B1(i,2); % B1(i,1), 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; %对角元, 节点q的自导纳, B1(i,4)为支路对地容抗

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元,节点p的自导纳 end

16

%输出导纳矩阵 disp('导纳矩阵Y='); disp(Y);

G=real(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); %PV节点电压给定模值 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 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); 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);

17

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 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

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;

18

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;

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);

shita(k)=atan(f(k)./e(k))*180/pi; E(k)=e(k)+f(k)*j; end

disp('各节点的实际电压标么值E为(节点号从小到大排列):'); disp(E);

19

disp('各节点的电压大小V为(节点号从小到大排列):'); disp(V);

disp('各节点的电压相角时shita为(节点号从小到大排列):'); disp(shita); 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); 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

四、实验数据：

电力系统接线图如图1所示，根据图中数据可得B1、B2和X矩阵：

20

P5=5j0.0151:1.05③j0.250.8+j0.30j0.250.04+j0.25②1.05:1j0.03①⑤U5=1.052+j13.7+j1.30.1+j0.35U1=1.05δ1=0°j0.25j0.25④1.6+j0.8

图1 电力系统接线图

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]; X=[1 0;2 0;3 0;4 0;5 0];

五、实验过程

1.认真阅读源程序，建立matlab的M程序，输入潮流计算程序。

2.输入实验数据,用B1、B2和X矩阵代替B1、B2和X输入语句并调试程序。 B2矩阵

B2(i,1)为节点所接发电机功率P； B2(i,2)为节点的负荷功率SL； B2(i,3)为节点电压的初值；

B2(i,4)为PV节点电压的给定值；

B2(i,5)为节点所接的无功补偿设备容量；

B2(i,6)为节点的分类号：1-平衡节点，2-PQ节点，3-PV节点。 3. 运行程序，输入节点数n=5、支路数n1=5、输入平衡母线节点号:isb=1、输入误差精度:pr=0.00001，进行计算。 六、要求

认真书写实验报告，画出网络接线图，在图上标出各节点电压和各支路功率，并用箭头标出各支路功率方向。

附件：计算参考结果 请输入节点数:n=5 请输入支路数:nl=5

请输入平衡母线节点号:isb=1 请输入误差精度:pr=0.00001 导纳矩阵Y= Column 1

0 -33.33333333333334i 0 +31.74603174603175i

21

0 0 0 Column 2

0 +31.74603174603175i 1.58459249980427 -35.73785857758467i -0.82987551867220 + 3.11203319502075i -0.75471698113208 + 2.64150943396226i 0 Column 3

0 -0.82987551867220 + 3.11203319502075i 1.45390047967064 -66.98082109846432i -0.62402496099844 + 3.90015600624025i 0 +63.49206349206349i Column 4

0 -0.75471698113208 + 2.64150943396226i -0.62402496099844 + 3.90015600624025i 1.37874194213052 - 6.29166544020251i 0 Column 5

0 0 0 +63.49206349206349i 0 0 -66.66666666666667i 迭代次数 5

没有达到精度要求的个数

7 8 8 6 0

各节点的实际电压标么值E为(节点号从小到大排列): Column 1

1.05000000000000 Column 2

1.03351793598970 - 0.07738281826869i Column 3

1.02600730614413 + 0.33047260236366i Column 4

0.85915374471705 - 0.07182065733577i Column 5

0.97461502463687 + 0.39067320582818i

各节点的电压大小V为(节点号从小到大排列): Columns 1 through 3

22

1.05000000000000 1.03641083773551 1.07791610674215 Columns 4 through 5

0.86215042996073 1.05000000000000

各节点的电压相角时shita为(节点号从小到大排列): Columns 1 through 3

0 -4.28193006920383 17.85352958577940 Columns 4 through 5

-4.77851050308856 21.84331901623337 各节点的功率S为(节点号从小到大排列): Column 1

2.57942727562302 + 2.29940213367660i Column 2

-3.70000000000000 - 1.30000000000000i Column 3

-2.00000000000000 - 1.00000000000000i Column 4

-1.59999999999999 - 0.79999999999996i Column 5

5.00000000000000 + 1.81308400259910i

各条支路的首端功率Si为(顺序同您输入B1时一样): 2.57942727562302 + 2.29940213367659i -1.27736037177247 + 0.20317048483099i 0.15678764739549 + 0.47131477149638i 1.58454630565503 + 0.67255630190799i 5.00000000000000 + 1.81308400259910i

各条支路的末端功率Sj为(顺序同您输入B1时一样): -2.57942727562302 - 1.97448525632736i 1.41545369434497 - 0.24433316462154i -0.13381873506914 - 0.39092357835416i -1.46618126493085 - 0.40907642164580i -5.00000000000000 - 1.42822313728644i

各条支路的功率损耗DS为(顺序同您输入B1时一样): -0.00000000000000 + 0.32491687734924i 0.13809332257250 - 0.04116267979055i 0.02296891232635 + 0.08039119314222i 0.11836504072417 + 0.26347988026219i 0.00000000000000 + 0.38486086531266i

23

实验四 简化模型的静态稳定计算的matlab程序设计

一、实验目的：

电力系统的静态稳定是指电力系统受到小干扰后，不发生自发振荡或非周期性失步，自动恢复到初始运行状态的能力。电力系统几乎时时刻刻都受到小的扰动。例如，个别电动机地接入和切除或加负荷和减负荷；又如架空输电线因风吹摆动引起的线间距离的微小变化；另外，发电机转子的旋转速度也不是绝对均匀，即功角也是有微小变化的。因此，电力系统的静态稳定问题实际上就是确定某个系统的某个运行稳态能否保持的问题。通过实验教学加深学生对电力系统的静态稳定概念的理解，掌握应用小扰动法分析简单电力系统静态稳定的方法。熟悉MATLAB软件的使用方法，提高编制调试计算机程序的能力，提高工程计算的能力，学习如何将理论知识和实际工程问题结合起来。

二、实验器材：

计算机、软件（已安装应用软件MATLAB等）、移动存储设备（学生自备，软盘、U盘等）。 三、实验内容：

1. 电力系统的静态稳定

简单电力系统，如图1所示。假设发电机为隐极机，则在小扰动下，发电机的功角特性方程为

PEq?EqUxx?sin? （1）

图1 简单电力系统的功角特性曲线

在图1中，设原动机的输入功率PT不变，并略去摩擦、风阻损耗和定子回路中的电阻损耗，并设此时发电机向系统输送的有功功率为PT与发电机输出电磁功率Eq(0)。则原动机输入机械功率P?Eq相等，即PT=PEq(0)。在功角性曲线上，满足功率平衡条件的运行有两个，即运行点a、b，与其

相对应的功率角分别为?a、?b。下面分析系统在两个点运行时的稳定情况。

在a点，系统保持稳定运行，发电机电动势E与无穷大系统母线电压U之间的相角为?a。若

q..此时系统中出现一个微小的、瞬时出现但又立即消失的扰动使功率角增加一个微小增量??时，此时的功角由?a变成了?a'，由图1可见，发电机的输出功率也相应增大了?P而变成了PEq(0)??P。但此时发电机输入的机械功率PT不变。因此，a点输出的电磁功率

'PE'qa将大于输入的机械功率

24

PEq(0)。当扰动消失后，在制动转矩（功率）作用下发电机将减速，功率角将减小，经过一系列衰

减的震荡后，最后稳定在a点。如图2（a）中实线所示。同理，当发动机受到的扰动使?a减小了??而变成了?a'时，运行点从a点运行到a点，因此a点输出的电磁功率PEq??将大于输入的机械功率

''''PT。当扰动消失后，在驱动转矩（功率）作用下发电机将被加速，功率角将增大，经过一系列衰减

的震荡后，最后又稳定在a点。如图2（a）中虚线所示。从以上分析可以看出，运行于a点的系统在受到微小扰动后能回复到原来的运行状态，所以说，a点是静态稳定运行点。但a点不是唯一的稳定运行点，其实在图1点以后的所有点，即0???90时，皆为静态稳定运行点。

00

图2 受扰动后功角随时间变化情况

在b点，当系统受到瞬时的微小扰动后，若功角?b有一个微小的增量??而变成了?b'，相应的运行点从b点转移到b点，但发电机输出的电磁功率减小了?P而变为PEq(0)??P，当扰动消失后，在加速转矩（功率）作用下机组将加速，功角将增大，当功角增大后，与之对应的输出的电磁功率将进一步减小，如此循环，发电机将不断被加速，功角不断增大，最终使发电机与系统之间失去同步，运行点不能再回到b点，如图2（b）中实线所示。同样，当b点受到微小扰动使功角减小一个微量??时，输出的电磁功率将增加到b点相对应的值PE''，而PE''>PT，当扰动消失后，在

qbqb'''制动转矩（功率）作用下发动机将减速，使功角进一步减小，如此循环，最后经过一系列震荡后，直至运行于稳定点a，如图2（b）中虚线所示，所以b点不是静态稳定运行点。即90???180都不是静态稳定运行点。

通过以上分析，得出结论，对上述简单电力系统，当0???90，时，电力系统可保持静态

0000稳定运行，在此范围内，

dPEqd??0；而900???1800时，电力系统不能保持静态稳定运行，在此

范围内，

dPEqd??0，由此得到电力系统静态稳定的实用判据为

dPEqd?SEq?

25

式中SEq亦称为比整步功率。

发电机在一定的运行条件下可发出最大的功率，即功角特性曲线图1中c点对应的功率（??90）称为稳定功率极限，用PM表示即

0PM?EqUxd

c点是发电机运行稳定与不稳定的临界点，实际运行时，要求发电机运行点与功率极限要有一定的距离，即保持有一定的稳定储备系数，以便系统有能力应付经常出现的一些干扰而不致丧失静态稳定。

静态稳定储备系数定义为

Kp?PM?PEqPEq?100%

我国现行的《电力系统安全稳定导则》中规定，电力系统正常运行方式和正常检修运行方式下，

Kp≥16%~20%；在事故后的运行方式和特殊PM?EqUxd的运行方式时Kp≥10%。

对电力系统的静态稳定性作较严格的计算时，可应用小扰动法。

小扰动法是根据李雅普诺夫稳定性理论，以线性化分析为基础的分析方法。当受扰动的线性化微分方程组的特征方程式的根的实部皆为负值时，该系统是稳定的，当手扰动的线性化微分方程组的特征方程式的根实部有正值时，该系统是不稳定的。

应用小扰动法分析简单电力系统静态稳定的步骤： （1）列出系统中描述各元件运动状态的微分方程组。

（2）将以上非线性方程线性化处理，得到近似的线性微分方程式组。 （3）根据近似方程式根的性质（根实部的正、负性）判断系统的稳定性。

发电机组的阻尼作用包括由轴承摩擦和发电机转子与气体摩擦所产生的机械性阻尼作用，机械阻尼作用与发电机的实际转速有关。我们只介绍不计发电机阻尼作用的情况。

1)列写运动方程（由转子运动方程列写）

EqU?d??0?0?0?[P?P]?[P?sin?]?[PTEqTT?PEq(?)]?TJxd?TJ?dtTJ ??d?????0??dt上面为简单电力系统的状态方程，由于PEq(?)中含有sin?(由式（1）可以看到)，所以方程是非线性的。

2)将非线性方程线性化

令?0是系统的一个平衡点（a点，在功角特性中对应的功率为P，并且当系统受到小扰动时，0）

???0???表示对平衡点施加微小扰动后的运动状态，将其代入式（1）中并在?0附近按泰勒级

数展开得：

26

PEq(?)??EqUxd?EqUxd?sin(?0???)?sin?0?dPEqd??0EqUxd?21dPEqsin?0?????d?02!d?2dPEq?0??2??

???P0??PE?PT??PE式中，?0?EqUxd?sin?0，??E?d?Ed???，所以，

?0PT?PEq(?)???PE??dPEd??0??

用小扰动法分析时状态变量为??和??，可表示为

???0??? ???0????1???

小扰动方程为

d?d(?0???)d????????0??? dtdtdt?dPd?d(?0???)d???0???[PT?PEq(?)]??0EdtdtdtTJTJd?写成矩阵形式为

???0??

0??????????1dPEq???????Td?????J?0????00??????????T?EqUsin?? (2) ??0??????0??Jxjd???J??3)根据状态方程系数矩阵的特征值判定系统的稳定性

根据李雅普诺夫稳定性理论：如果状态方程系数矩阵的所有特征根都为负实数或是具有负实部的负数，则系统是稳定的；若特征根中出现一个零根或实部为零的一对虚根，则系统处于稳定的边界；若特征根有一个正实数或一对具有正实部的虚根，则系统是不稳定的，其中，特征值仅有一个是正实数，系统将非周期性失去稳定；特征值为一对具有实部的负数时，系统将周期性增幅振荡而失去稳定。

求式(2)二阶微分方程组的特征根为

?0dPEq?1,2???TJd?分析式（3），当

???0 (3)

d?Eqd????0?0时，?1,2为一个正实根和一个负实根，即??和??有随时间不

断单调增加的趋势，发电机相对于无限大系统非周期性失去同步，故系统是不稳定的。

当

d?Eqd????0?0时，?1,2为一对虚根，理论上??和??作等幅振荡，即状态变量?和?没有

稳定值，系统同样不稳定。但在以上的分析中，都未必考虑阻尼因素的影响，实际上，系统中由于

27

阻尼作用，??和??将作衰减的振荡，最后都稳定在最初值，系统恢复同步。

由以上分析可见，用小扰动法对简单系统稳定性分析的结果和用物理概念分析的结果是一致的，得到同一个静态稳定判据，即

d?Eqd??0

当发电机与无限大系统之间发生振荡（??和??振荡）或失去同步时，在发电机的转子回路，特别是在阻尼绕组中将有感应电流儿产生阻尼转矩或异步转矩。总的阻尼功率可近似表示为

?D?D?? （4）

式中，D为阻尼功率系数。

记及阻尼功率后，发电机转子转动方程为

d??????0 （5） dt???0

dPd??1Eq??(D???dtTJd???) （6）

式（5）和（6）写成矩阵形式为

0??????????1dPEq???????????TJd?式（7）特征方程的特征根为

???0????D??? （7） ??????TJ???0??1,2??对隐极机，由PEq?D1?2TJ2TJdPEqd??D2?4?0TJEqUxx?cos?

dPEqd????0 （8）

EqUxx?sin?，得

EQ0?(U0?Q0xq?U0)2?(P0xq?U0)2

?0?arctg特征根?具有负实部的条件为

P0xq?/U0U0?Q0xq?/U0

D?0??Eq?S?dPEq?d??当D?0，且D?4?0TJ2???0?0

d?Eqd????0时，?为两个负实根，系统在受到小扰动后，

发电机的状态变量?和?将按指数函数规律衰减到最初值；当D?0,但D?4?0TJ28

2d?Eqd?时，?

为一对具有负实部的共轭复根，这时系统在受到小扰动后，发电机状态变量?和?将做衰减的振荡，最后稳定在初值。

当D?0时，特征方程式的根?1,2至少有一个是正实数或两个都为具有正实部的共轭复根，无论SEq为何值，系统都是不稳定的。因此，要使系统保持静态稳定，必须同时满足方程中的两个条件，否则，系统就是不稳定的。

2. 程序代码及说明

S0=input('请输入初始功率:S0=');

V0=input('请输入无限大系统母线电压:V0='); Xd=input('请输入系统直轴等值电抗:Xd='); w0=input('请输入同步电角速度:w0='); Tj=input('请输入惯性时间常数:Tj='); zn=input('请输入综合阻尼系数:D=');

Eq=sqrt((V0+imag(S0)*Xd./V0)^2+(real(S0)*Xd./V0)^2); dtj0=atan(real(S0)*Xd./(V0*(V0+imag(S0)*Xd./V0))); Ps1=Eq*V0./Xd; C=w0;

D=-1./Tj*cos(dtj0)*V0./Xd; E=w0*zn./(2*Tj);

Kp=(Ps1-real(S0))./real(S0); root1=-E+sqrt(E^2+C*D); root2=-E-sqrt(E^2+C*D); disp('系统的静态稳定极限:Ps1='); disp(Ps1);

disp('系统的静态储备系数:Kp='); disp(Kp);

disp('系统的特征根1:root1='); disp(root1);

disp('系统的特征根2:root2='); disp(root2); if zn>=0

if(root1==conj(root2))|((real(root1)<0) & (real(root2)<0));%conj求复数的共轭,与\，或\

disp('该系统是静态稳定的!'); else

disp('该系统是不稳定的!'); end

elseif zn<0

disp('该系统是不稳定的!'); end

29

四、实验数据：

电力系统接线图如图3所示， 变压器T2与无限大系统相连接，系统直轴等值电抗标幺值为： Xd=XGd+XT1+XL/2+XT2=1.8+0.197+0.35+0.142=2.489

T1L1T2G~L2 图3 电力系统接线图

五、实验过程

1.手算图1所示网络的等值电抗。

2.认真阅读源程序，建立matlab的M程序，输入静态稳定计算程序。 3.运行程序，输入实验数据： 输入初始功率:S0=0.583+0.361i 输入无限大系统母线电压:V0=1 输入系统直轴等值电抗:Xd=2.489 输入同步电角速度:w0=100*pi 输入惯性时间常数:Tj=10 输入综合阻尼系数:D=0

4.分析计算结果，检验程序的正确性。 六、要求

认真书写实验报告并上交实验数据

30

实验五 短路电流计算程序的实现

一、实验目的：

掌握电力系统短路电流分析方法。 二、实验器材：

计算机、软件（已安装应用软件MATLAB等）、移动存储设备（学生自备，软盘、U盘等）。 三、实验内容

计算短路电流周期分量，如I??（I?）时，实际上就是求解交流电路的稳态电流，其数学模型也就是网络的线性代数方程，一般选用节点电压方程。方程的系数矩阵是对称的。在短路电流计算中变化的量往往是方程的常数项，需要多次求解线性方程组。 1．等值网络

图1给出了不计负荷情况下计算短路电流I??的等值网络。在图1（a）中G代表发电机端电压

???和x??，D表示负荷节点，f点为直接短路点。应用叠加原节点，发电机等值电势和电抗分别为Ed理如图1所示。正常运行方式为空载运行，网络中各点电压均为1；在故障分量网络中。只需作故

障分量的计算。由图1的故障分量网络可见，这个网络与潮流计算的网络的差别在于发电机节点上

??。当然如果短路计算中可以忽略线路电阻和电纳，而且不计变压器的实际变比，多接了对地电抗xd则短路计算网络较潮流计算网络简化，而且网络本身是纯感性的。

G1?x?d1?x?d2G2D1D2G1?x?d1?x?d2G2D1D2G1?x?d1?1?x?d2G2D1D2f???E2f???E2?Uf0f???E2??Uf0???E1???E1???E1??1zf

图1 在不计负荷情况下计算短路电流I″的等值电路

2． 用节电阻抗矩阵计算短路电流

如果已经形成了故障分量网络的节点阻抗矩阵，则矩阵中的对角元素就是网络从f点看进去的等值阻抗，又称为f点的自阻抗。Zfi为f点与i点的互阻抗，均用大写Z表示。由节点方程中的

??ZI???第f个方程：Uff11???ZffIf???ZfnIn。Zff为其它节电电流为零时，节点f的电压

和电流之比，即网络对f点的等值阻抗。

根据故障分量网络，直接应用戴维南定理可求得直接短路电流（由故障点流出）为

??If?Uf0Zff?zf （1）

?为f点短路前的电压。 式中，zf为接地阻抗；Uf0如果短路点为直接短路，则zf=0，在实用计算中采用（2）式

31

??If?Uf0Zff?1 （2） Zff因此，一旦形成了节点阻抗矩阵，任一点的短路电流即可方便地求出，即等于该点自阻抗（该点对角元素）的倒数。

3．计算节点电压和支路电流

?，各节点电压的故障分量为 由故障分量网络可知，只有节点f有节点电流?If???U?1??Z11?Z1f????????????0????????????U?f?=?Zf1?Zff?Z?fn????Z1f??????????????????=?Z?ff???? ????I?f?0????If??U?n????Zn1?Znf?Znn????????????Znf??所以,，各节点短路故障后的电压为

U?1?U?10?ΔU?1?U?10?Ζ1fΙ?f???????U??f?U?f0?ΔU?f?0?? ??????U?n?U?n0?ΔU?n?U?n0?ΖnfΙ??f??任一支路i-j的电流为

I??U?I?U?jijz ij式中，zij为i-j支路的阻抗。

32

(3) (4) (5)

源程序：

NF=input('请输入短路点的数目：NF='); n=input('请输入节点数:n='); nl=input('请输入支路数:nl=');

B=input('请输入由各支路参数形成的矩阵:B=');

V0= input('请输入由各节点的初电压标么值形成的列矩阵: V0='); D= input('请输入由短路号，短路点阻抗组成的矩阵D='); m=0;Z=zeros(n); % zeros(n)产生n×n的全0方阵。 V=zeros(n); I=zeros(nl);

for k1=1:nl %支路数

p=B(k1,1);q=B(k1,2); % B(k1,1), B(k1,2)为支路编号 if B(k1,6)==0 % k=1./B(k1,5);

33

else

k=B(k1,5); end

if p==0 %

if q>m % 追加接地树支 Z(q,q)=B(k1,3);m=m+1; else

for i=1:m

Z(i,m+1)=-Z(i,q);Z(m+1,i)=-Z(q,i); end

Z(m+1,m+1)=Z(q,q)+B(k1,3); for i=1:m for j=1:m

Z(i,j)=Z(i,j)-Z(i,m+1)*Z(m+1,j)./Z(m+1,m+1); end

Z(i,m+1)=0; end

for i=1:m+1

Z(m+1,i)=0; end end

else if q>m %追加不接地树支 for i=1:m

Z(i,q)=Z(i,p)*k;Z(q,i)=Z(p,i)*k; end

Z(q,q)=k^2*Z(p,p)+k^2*B(k1,3); m=m+1; else

for i=1:m

Z(i,m+1)=Z(i,p)*k-Z(i,q); Z(m+1,i)=Z(p,i)*k-Z(q,i); end

Z(m+1,m+1)=k^2*Z(p,p)+Z(q,q)-2*k*Z(p,q)+k^2*B(k1,3); for i=1:m for j=1:m

Z(i,j)=Z(i,j)-Z(i,m+1)*Z(m+1,j)./Z(m+1,m+1); end

Z(i,m+1)=0; end

for i=1:m+1

Z(m+1,i)=0; end

34

end end end

for k=1:NF

I(D(k,1),D(k,1))=V0(D(k,1),1)./(Z(D(k,1), D(k,1))+D(k,2)); ft=num2str(D(k,1)); ts1=('点短路时');

ts2=('电流的标么值If='); dn=strcat(ft,ts1,ts2); disp(dn);

disp(I(D(k,1),D(k,1))); for i=1:n

V(i,i)=V0(i,1)-I(D(k,1),D(k,1))*Z(i,D(k,1));%求各支路的短路电流标么值 end for i=1:nl

if B(i,6)==0 k=B(i,5); else k=1./B(i,5); end

p=B(i,1);q=B(i,2); if p==0

e=0;b=B(i,3);

I(i,i)=(e-V(q,q)./k)./b; else

I(i,i)=(V(p,p)-V(q,q)./k)./B(i,3); end end

disp('各节点的电压标么值U为（节点号从小到大排）：'); for i=1:n disp(V(i,i)); end

disp('各之路短路电流的标么值I为（顺序同您输入B时一样）：'); for i=1:nl disp(I(i,i)); end end 四 实验数据

所示电力系统，负荷全部略去，简化后的各电抗标么值注于等值网络中。试计算f点三相短路时的短路电流及网络中的电流分布。

35

B=[0 1 0.2i 0 1 0;0 2 4i 0 1 0;1 3 0.51i 0 1 0;2 3 0.59i 0 1 0;3 4 1.43i 0 1 0]; V0=[1;1;1;1]; D=[4,0];

五、实验过程

1.手算图1所示网络的等值电抗。

2.认真阅读源程序，建立matlab的M程序，输入短路电流计算程序。 3.运行程序，输入实验数据：

矩阵B的每行是由下列参数构成的：

? 某支路的首端号P；

? 某支路末端号Q，且P

? 变比折算到哪一侧的标志（如果支路的首端P处于高压侧则请输入“1”， 否则请输入“0”）。 矩阵D：由短路点号，短路点阻抗

列矩阵V0：由各节点的初电压标么值形成 4.分析计算结果，检验程序的正确性。 六、要求

认真书写实验报告并上交实验数据

请输入短路点的数目：NF=1 请输入节点数:n=4 请输入支路数:nl=5

请输入由各支路参数形成的矩阵:B=[0 1 0.2i 0 1 0;0 2 4i 0 1 0;1 3 0.51i 0 1 0;2 3 0.59i 0 1 0;3 4 1.43i 0 1 0];

请输入由各节点的初电压标么值形成的列矩阵: V0=[1;1;1;1];

36

请输入由短路号，短路点阻抗组成的矩阵D=[4,0];

37

本文来源：https://www.bwwdw.com/article/w8m6.html