三机九节点电力系统仿真matlab

三机九节点电力系统

暂态仿真

学院：能源与动力工程学院 专业：电力系统及其自动化 学号：

姓名：*** 导师： 授课教师：

目录

一、 概述 ............................................................................................................................. 1 二、 课程主要任务 ................................................................................................. 1

1. 2. 3. 4. 5. 6.

系统数据 ........................................................................................................................ 1 潮流计算 ........................................................................................................................ 2 负荷等效和支路简化 ................................................................................................. 4 求解电磁功率 ............................................................................................................... 5 求解运动方程 ............................................................................................................... 5 程序清单 ........................................................................................................................ 7

(1). (2). (3). (4).

主程序： ............................................................................................................... 7 极坐标转换成直角坐标函数pol2rect(V,del) .................................................... 16 直角坐标转换成极坐标函数rect2pol(Z) .......................................................... 16 求解微分方程所用的得到微分量的函数Gen_fw(t,X,Y_Gen,E,Pm0,Tj) ....... 16

三、 课程总结及心得体会 ........................................................................... 16 四、 参考文献 ............................................................................................................. 17

于永生 电力系统仿真作业

一、 概述

在动态稳定分析中，系统由线性化的微分方程组和代数方程组描写，并用经典的或现代的线性系统理论来进行稳定分析，分析可以在时域或频域进行。当用计算机和现代线性系统理论分析时，常把系统线性化的微分方程组和代数方程组消去代数变量，化为状态方程形式，并广泛采用特征分析进行稳定分析。

电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞大复杂的动力学系统。其暂态过渡过程不仅包括电磁方面的过渡过程，而且还有机电方面的过渡过程。由此可见，电力系统的数学模型是一个强非线性的高维状态方程组。在动态稳定仿真中使用简单的电力系统模型，发电机用三阶模型表示。

二、 课程主要任务

本次课程主要应用P. M. Anderson and A. A. Fouad编写的《Power System Control and Stability》一书中所引用的Western System Coordinated Council (WSCC)三机九节点系统模型。

1. 系统数据

其中，节点数据如下：

%节点数据

% 节点 电压 电压 发电机 发电机 负荷 负荷 节点

% 号 幅值 相角 有功 无功 有功 无功 类型(1PQ 2PV 3平衡)

N=[ 1 1.04 0 0.7164 0.2705 0 0 3 2 1.025 0 1.63 0.0665 0 0 2 3 1.025 0 0.85 -0.1086 0 0 2 4 1 0 0 0 0 0 1 5 1 0 0 0 1.25 0.5 1 6 1 0 0 0 0.9 0.3 1 7 1 0 0 0 0 0 1 8 1 0 0 0 1 0.35 1 9 1 0 0 0 0 0 1]; 其中，支路数据如下： % 线路数据

% 首端 末端 电阻 电抗 电纳（1/2) 变压器非标准变比 L=[4 5 0.01 0.085 0.088 1 4 6 0.017 0.092 0.079 1 5 7 0.032 0.161 0.153 1 6 9 0.039 0.17 0.179 1 7 8 0.0085 0.072 0.0745 1

1

电力系统仿真作业 于永生

8 9 0.0119 0.1008 0.1054 1 1 4 0 0.0576 0 1 2 7 0 0.0625 0 1 3 9 0 0.0586 0 1]; 发电机数据如下：

% 发电机 母线 Xd Xd' Td0' Xq Xq' Tq0’ Tj Xf Ge=[ 1 1 0.1460 0.0608 8.96 0.0969 0.0969 0 47.28 0.0576 2 2 0.8958 0.1198 6.00 0.8645 0.1969 0.535 12.80 0.0625 3 3 1.3125 0.1813 8.59 1.2578 0.2500 0.600 6.02 0.0585]; 系统电路结构拓扑图如下：

图1 WSCC 3机9节点系统

（所有参数以100MVA为基准值的标幺值）

2. 潮流计算

首先进行潮流计算，采用牛顿拉夫逊迭代法，电力系统潮流计算是电力系统运行和规划中最基本和最经常的计算，其任务是在已知某些运行参数的情况下，计算出系统中全部

2

于永生 电力系统仿真作业

的运行参数，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点除外），可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵和网络拓扑结构列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。

牛顿拉夫逊算法修正方程

W = -JΔV

其中W是节点不平衡量向量，包括有功，无功，电压；J是雅克比矩阵；ΔV是节点电压修正量。

令

Vieijfi;YijGijjBij则直角坐标形式的功率不平衡量方程为 PQ节点：

，

PiPisei(GijejBijfj)fj(GijfjBijej)0j1j1nnnn

QiQisfi(GijejBijfj)ej(GijfjBijej)0j1j1

PV节点：

PiPisei(GijejBijfj)fj(GijfjBijej)0j1j1nn

Vi2Vis2Vi2Vis2(ei2fi2)极坐标形式的功率不平衡量方程

PiPisViVj(GijcosijBijsinij)0j1nn

QiQisViVj(GijsinijBijcosij)0j1

雅可比矩阵J各元素的表达式

HJMNL

当j≠i时：

3

电力系统仿真作业 于永生

HijPiBijeiGijfifj

NijPiGijeiBijfiejQiBijfiGijeifj

Mij

LijQiGijfiBijeiej

当j=i时：

HiiPiBiieiGiifibifi PiGiieiBiifiaiei QiGiieiBiifiaifi

NiiMiiLiiQiBiieiGiifibiei

n其中，

ai(GijejBijfj)i1nbi(GijfjBijej)i1。

进行牛顿拉夫逊算法迭代后得到电压幅值V和相角θ。

3. 负荷等效和支路简化

然后求出支路电流，将发电机内电抗X’加入系统导纳矩阵，求出发电机内电势E’。加入发电机内节点后，系统导纳矩阵变成12*12阶的矩阵，并将负荷等效成阻抗。然后将支路导纳矩阵分块，如下：

AEY

CD其中，A是3*3的方阵，E是3*9的矩阵，C是9*3的矩阵，D是9*9的方阵。 经过网络简化得到故障前的3*3简化导纳矩阵

4

于永生 电力系统仿真作业

Ypre=A-E*(inv(D))*C （1）

其中“inv（D）”是MATLAB中D矩阵的求逆。

故障中导纳矩阵的第七行和第七列从矩阵中删除，此时有

AEY

CD此时，A是3*3的方阵，E是3*8的矩阵，C是8*3的矩阵，D是8*8的方阵。简化

矩阵求法如同公式（1）。

故障后的Y矩阵相对于故障前的Y矩阵只是第5个节点和第7个节点有变化，反映到12*12的矩阵中即为第（10,8），（8,10），（8,8），（10,10）位置的元素有变化，其中（10,8），（8,10）位置的元素变为零，(8,8)，（10,10）节点在故障前的基础上加上（8,10）位置元素的值。然后简化导纳矩阵的求法同式（1）。

4. 求解电磁功率

得到故障前，故障中，故障后三个不同的导纳矩阵后，就开始计算电磁功率和机械功率，机械功率等于稳态的电磁功率中的有功分量。所以可以有

Pe=real(E*I) (2)

公式（2）中，E为发电机内电势，I为从发电机流出的电流。

但在参考文献Ramnarayan Patel, T. S. Bhatti and D. P. Kothari.MATLAB/Simulink-based transient stability analysis of a multimachine power system中给出的电磁功率计算公式为：

PeiEiGiiEiEjYijcos(ijij)

2j1jin本人在此次仿真中用的是公式（2）计算得到的电磁功率。 稳态情况下有，机械功率

Pme=Pe

5. 求解运动方程

发电机的运动方程可以写成常微分方程组：

n2HdiiE2GEEYcos(DP)jjmiiiiijijijijRdtj1ji

diRdt其中Pmi为第i个机组故障前稳态的电磁功率。在本次仿真中Djωi为零，即阻尼为零。仿真开始，t=0时引入故障，0.083s后切除故障。 求解运动方程后得到ω和δ曲线如下：

5

电力系统仿真作业 于永生

2号发电机

3号发电机

图1 ω曲线

2号发电机

3号发电机

1号发电机

图2 δ曲线

6

于永生 电力系统仿真作业

6. 程序清单

以下是我编写的仿真程序，除主程序外还包含三个函数：①极坐标转换成直角坐标函数pol2rect(V,del),其中输入参数V为极坐标向量的幅值，del为极坐标向量的角度，函数返回值为一个复数，即为极坐标向量在直角坐标中的复数值；②直角坐标转换成极坐标函数rect2pol(Z),其中输入参数Z为一个复数，函数返回值为一个二维向量，向量的第一个数为幅值，第二个数为相角；③求解微分方程所用的得到微分量的函数Gen_fw(t,X,Y_Gen,E,Pm0,Tj),其中输入参数X为ω和δ的迭代初值，Y_Gen，为求解所用的导纳矩阵，这里是三阶的方阵，对应于故障前，故障中，和故障后的三个Y矩阵，E为发电机内电势，Pm0为机械功率，等于稳态时的有功功率，Tj为运动方程中的2H。

(1). 主程序：

clc; clear; %节点数据

% 节点 电压 电压 发电机 发电机 负荷 负荷 节点

% 号 幅值 相角 有功 无功 有功 无功 类型(1PQ 2PV 3平衡 % N=[ 1 1.04 % 2 1.025 % 3 1.025

0 0

0.7164 0.2705 1.63

0 0

3

2 2

0.0665 -0.1086

0 0 0 0

0 0.85

% 4 1 0 % 5 1 0 % 6 1 0 % 7 1 0 % 8 1 0 % 9 1 0 % 线路数据

0 0 0 0 0 0 1.25

1

0.5 1

1 1];

0 0 0.9 0.3 1 0 0 0 0 0 0 0 0 0 0 1 0.35 1

% 首端 末端 电阻 电抗 电纳（1/2) 变压器非标准变比 L=[4 5 0.01 0.085 0.088 1 4 6 0.017 0.092 0.079 1 5 7 0.032 0.161 0.153 1 6 9 0.039 0.17 0.179 1 7 8 0.0085 0.072 0.0745 1 8 9 0.0119 0.1008 0.1054 1 1 4 0 0.0576 0 1 2 7 0 0.0625 0 1 3 9 0 0.0586 0 1]; %形成节点导纳矩阵 fb=L(:,1);%起始节点编号 tb=L(:,2);%末节点编号 r=L(:,3);%电阻 x=L(:,4);%电抗 b=L(:,5);%电纳

a=L(:,6);%变压器非标准变比

7

电力系统仿真作业

z=r + i*x;%z矩阵，复数，9乘1的矩阵 y=1./z;%求倒数，9乘1的导纳矩阵 b=i*b;%虚数

nb=max(max(fb),max(tb));%nb=9 nl=length(fb);%支路个数

Y=zeros(nb,nb);%Y是9乘9的支路矩阵

for k=1:nl

于永生

Y(fb(k),tb(k))=Y(fb(k),tb(k))-y(k)/a(k);%Y方阵的元素，非对角元（4,5）（4，6）（5，7）（6,9）（7，8）（8,9）（1,4）（2,7）（3,9）

Y(tb(k),fb(k))=Y(fb(k),tb(k));%Y方阵的元素，上一句的对称部分（5,4）（6,4）（7，5）（9，6）（8，7）（9,8）（4,1）（7,2）（9,3）

end

for m=1:nb for n=1:nl

if fb(n)==m%支路中的第n个数的首节点编号等于m Y(m,m)=Y(m,m)+y(n)/(a(n)^2)+b(n);%Y矩阵中的对角元 elseif tb(n)==m%支路中的第n个数的末节点编号等于m

Y(m,m)=Y(m,m)+y(n)+b(n);%因为支路是首端或末端等于m，不能是首端编号等于末端编号

end end end

%节点数据

% 节点 电压 电压 发电机 发电机 负荷 负荷 节点

% 号 幅值 相角 有功 无功 有功 无功 类型(1PQ 2PV 3平衡 N=[ 1 2 4 5 6 7 8 9

1.04 1.025

0 0

0.7164 0.2705 1.63

0 0

3

2 2

0.0665 -0.1086

0 0 0 0

3 1.025 0 0.85

1 0 1 0 1 0 1 0 1 0 1 0

0 0 0 0 0 0 1.25

1

0.5 1

1 1];

0 0 0.9 0.3 1 0 0 0 0 0 0 0 0 0 0 1 0.35 1

nbus=9;%节点数 bus=N(:,1);%节点编号 type=N(:,8);%节点类型 V=N(:,2);%额定电压 del=N(:,3);%电压相角 Pg=N(:,4);%发电机有功功率 8

于永生

电力系统仿真作业

Qg=N(:,5);%发电机无功功率 Pl=N(:,6);%负荷有功功率 Ql=N(:,7);%负荷无功功率 P=Pg-Pl;%发出减消耗 Q=Qg-Ql;%发出减消耗 Psp=P;%额定有功功率 Qsp=Q;%额定无功功率 G=real(Y); B=imag(Y);

pv=find(type==3|type==2);%pv节点编号和平衡节点编号 pq=find(type==1);%pq节点编号 npv=length(pv);%npv=3 npq=length(pq);%npq=6 Tol=1; Iter=1; while(Tol>1e-7)

P=zeros(nbus,1);%迭代初值P为零 Q=zeros(nbus,1);%迭代初值Q为零 for i=1:nbus for k=1:nbus

P(i)=P(i)+V(i)*V(k)*(G(i,k)*cos(del(i)-del(k))+B(i,k)*sin(del(i)-del(k))); Q(i)=Q(i)+V(i)*V(k)*(G(i,k)*sin(del(i)-del(k))-B(i,k)*cos(del(i)-del(k))); end end

dPa=Psp-P;%额定有功功率减去流出或流入的有功功率 dQa=Qsp-Q;%额定无功功率减去流出或流入的无功功率 k=1;

dQ=zeros(npq,1);%pq节点的无功功率，pq乘1的矩阵,6乘1 for i=1:nbus

if type(i)==1%pq节点 dQ(k,1)=dQa(i); k=k+1; end end

dP=dPa(2:nbus); Z=[dP;dQ];

%下面计算雅克比矩阵

J1=zeros(nbus-1,nbus-1);%有功功率对功角求偏微分 for i=1:(nbus-1) m=i+1; for k=1:(nbus-1)

9

电力系统仿真作业

n=k+1; if n==m for n=1:nbus

于永生

J1(i,k)=J1(i,k)+V(m)*V(n)*(-G(m,n)*sin(del(m)-del(n))+B(m,n)*cos(del(m)-del(n)));

end

J1(i,k)=J1(i,k)-V(m)^2*B(m,m); else

J1(i,k)=V(m)*V(n)*(G(m,n)*sin(del(m)-del(n))-B(m,n)*cos(del(m)-del(n))); end end end

J2=zeros(nbus-1,npq);%有功功率对电压求偏微分 for i=1:(nbus-1) m=i+1; for k=1:npq n=pq(k); if n==m for n=1:nbus

J2(i,k)=J2(i,k)+V(n)*(G(m,n)*cos(del(m)-del(n))+B(m,n)*sin(del(m)-del(n))); end

J2(i,k)=J2(i,k)+V(m)*G(m,m); else

J2(i,k)=V(m)*(G(m,n)*cos(del(m)-del(n))+B(m,n)*sin(del(m)-del(n))); end end end

J3=zeros(npq,nbus-1);%无功功率对相角求偏微分 for i=1:npq m=pq(i); for k=1:(nbus-1) n=k+1; if n==m for n=1:nbus

J3(i,k)=J3(i,k)+V(m)*V(n)*(G(m,n)*cos(del(m)-del(n))+B(m,n)*sin(del(m)-del(n)));

end

J3(i,k)=J3(i,k)-V(m)^2*G(m,m); else

J3(i,k)=V(m)*V(n)*(-G(m,n)*cos(del(m)-del(n))-B(m,n)*sin(del(m)-del(n))); end 10

于永生

电力系统仿真作业

end end

J4=zeros(npq,npq);%无功功率对电压求偏微分 for i=1:npq m=pq(i); for k=1:npq n=pq(k); if n==m for n=1:nbus

J4(i,k)=J4(i,k)+V(n)*(G(m,n)*sin(del(m)-del(n))-B(m,n)*cos(del(m)-del(n))); end

J4(i,k)=J4(i,k)-V(m)*B(m,m); else

J4(i,k)=V(m)*(G(m,n)*sin(del(m)-del(n))-B(m,n)*cos(del(m)-del(n))); end end end

J=[J1 J2;J3 J4];

%X=inv(J)*Z; X=J\\Z;

dTh=X(1:nbus-1); dV=X(nbus:end);

del(2:nbus)=dTh+del(2:nbus); k=1; for i=2:nbus if type(i)==1

V(i)=V(i)+dV(i-3); k=k+1; end end Iter=Iter+1; Tol=max(abs(Z)); end

Vm = pol2rect(V,del); % 极坐标转换成直角坐标 Del = 180/pi*del; % 弧度转化成角度

11

电力系统仿真作业

Iij = zeros(nb,nb);%支路电流

for m = 1:nl%每一条支路都过一遍 p = fb(m); q = tb(m);

于永生

Iij(p,q) = -(Vm(p) - Vm(q))*Y(p,q); % Y(m,n) = -y(m,n).. Iij(q,p) = -Iij(p,q); end

Iij = sparse(Iij);

% 线路数据

% 首端 末端 电阻 电抗 L=[4 5 0.01 0.085 4 6 0.017 0.092 5 7 0.032 0.161 6 9 0.039 0.17 7 8 0.0085 0.072 8 9 0.0119 0.1008 1 4 0 0.0576 2 7 0 0.0625 3 9 0 0.0586

Yzengjia=zeros(12,12);

Zci=[j*0.0608 j*0.1198 j*0.1813]; Yci=1./Zci; for k=1:3

Yzengjia(k,k)=Yci(k); end for k=1:9 for m=1:9

Yzengjia(3+k,3+m)=Y(k,m); end end

Yzengjia(1,4)=-Yzengjia(1,1); Yzengjia(4,1)=Yzengjia(1,4); Yzengjia(2,5)=-Yzengjia(2,2); Yzengjia(5,2)=Yzengjia(2,5); Yzengjia(3,6)=-Yzengjia(3,3); Yzengjia(6,3)=Yzengjia(3,6);

Yzengjia(4,4)=Yzengjia(4,4)+Yzengjia(1,1); Yzengjia(5,5)=Yzengjia(5,5)+Yzengjia(2,2); Yzengjia(6,6)=Yzengjia(6,6)+Yzengjia(3,3); 12

电纳（1/2) 变压器非标准变比 0.088 1 0.079 1 0.153 1 0.179 1 0.0745 1 0.1054 1 0 1 0 1 0 1]; 于永生

%求负荷等效导纳

Y5=N(5,6)/(V(5))^2-j*(N(5,7)/(V(5))^2); Y6=N(6,6)/(V(6))^2-j*(N(6,7)/(V(5))^2); Y8=N(8,6)/(V(8))^2-j*(N(8,7)/(V(5))^2); Yzengjia(8,8)=Yzengjia(8,8)+Y5; Yzengjia(9,9)=Yzengjia(9,9)+Y6; Yzengjia(11,11)=Yzengjia(11,11)+Y8;

A=zeros(3,3);

for m=1:3 for n=1:3

A(m,n)=Yzengjia(m,n); end end

E=zeros(3,9);

for m=1:3 for n=4:12

E(m,n-3)=Yzengjia(m,n); end end

C=zeros(9,3);

for m=4:12 for n=1:3

C(m-3,n)=Yzengjia(m,n); end end

D=zeros(9,9);

for m=4:12 for n=4:12

D(m-3,n-3)=Yzengjia(m,n); end end

Ypre=A-E*(inv(D))*C;%得到故障前的3*3矩阵

Yzengjiadur=Yzengjia; Yzengjiadur(:,[10])=[];

电力系统仿真作业

13

电力系统仿真作业

Yzengjiadur([10],:)=[]; Adur=A; Edur=zeros(3,8); for m=1:3 for n=1:8

于永生

Edur(m,n)=Yzengjiadur(m,n+3);%因为Edur的第四到第九列是零，而所以这里为了方便，直接取E的前六列

end end

Cdur=zeros(8,3); for m=1:8 for n=1:3

Cdur(m,n)=Yzengjiadur(m+3,n);%因为Cdur的第四到第九行是零，而所以这里为了方便，直接取C的前六行

end end

Ddur=zeros(8,8); for m=1:8 for n=1:8

Ddur(m,n)=Yzengjiadur(m+3,n+3); end end

Ydur=Adur-Edur*(inv(Ddur))*Cdur;%故障中的3*3矩阵

Yzengjiaaft=Yzengjia; temp=Yzengjiaaft(8,10); Yzengjiaaft(8,10)=0;

Yzengjiaaft(8,8)=Yzengjiaaft(8,8)+temp; Yzengjiaaft(10,8)=0;

Yzengjiaaft(10,10)=Yzengjiaaft(10,10)+temp; Aaft=zeros(3,3); for m=1:3 for n=1:3

Aaft(m,n)=Yzengjiaaft(m,n); end end

Eaft=zeros(3,9); for m=1:3 for n=4:12

Eaft(m,n-3)=Yzengjiaaft(m,n); end end

Caft=zeros(9,3); 14

于永生

for m=4:12 for n=1:3

Caft(m-3,n)=Yzengjiaaft(m,n); end end

Daft=zeros(9,9); for m=1:9 for n=1:9

Daft(m,n)=Yzengjiaaft(m+3,n+3); end end

Yaft=Aaft-Eaft*(inv(Daft))*Caft;%故障后的3*3矩阵

%求电磁功率

Enei=ones(3,1);%发电机内电势 Pe=ones(3,1);

Enei(2)=Vm(2)+Zci(2)*Iij(2,7); Enei(1)=Vm(1)+Zci(1)*(-Iij(4,1)); Enei(3)=Vm(3)+Zci(3)*Iij(3,9);

polEnei1=rect2pol(Enei(1));%直角坐标转换成极坐标后得到角度 polEnei2=rect2pol(Enei(2)); polEnei3=rect2pol(Enei(3));

Pe(1)=real(Enei(1)*((-Iij(4,1))'));%故障前的电气功率 Pe(2)=real(Enei(2)*((Iij(2,7))')); Pe(3)=real(Enei(3)*(Iij(3,9)));

delta=zeros(3,1); delta(1)=polEnei1(2); delta(2)=polEnei2(2); delta(3)=polEnei3(2);

Tj=[47.28;12.8;6.02];

Pm0=Pe;%机械功率=故障前的电磁功率 X0=[ones(3,1).*(2*pi*60),delta]; Y_Gen=Ydur;E=abs(Enei);

Time_s=[0,0.083];%0.083s后切除故障

[t1,Xout1]=ode45(@(t,X) Gen_fw(t,X,Y_Gen,E,Pm0,Tj),Time_s,X0); Time_s=[0.083,2];

Y_Gen=Yaft;X0=Xout1(end,:);

电力系统仿真作业

15

电力系统仿真作业 于永生

[t2,Xout2]=ode45(@(t,X) Gen_fw(t,X,Y_Gen,E,Pm0,Tj),Time_s,X0); XOUT=[Xout1;Xout2];T=[t1;t2]; XOUT=XOUT*180/pi;

figure,plot(T,XOUT(:,4),'r',T,XOUT(:,5),'g',T,XOUT(:,6),'y'); figure,plot(T,XOUT(:,5)-XOUT(:,4),'r',T,XOUT(:,6)-XOUT(:,4),'Y');

(2). 极坐标转换成直角坐标函数pol2rect(V,del)

function rect = pol2rect(rho,theta) rect = rho.*cos(theta) + j*rho.*sin(theta);

(3). 直角坐标转换成极坐标函数rect2pol(Z)

function pol=rect2pol(Z) t=real(Z); k=imag(Z);

pol(1)=(t^2+k^2)^(1/2);%模 pol(2)=atan(k/t);

(4). 求解微分方程所用的得到微分量的函数Gen_fw(t,X,Y_Gen,E,Pm0,Tj)

function [Xp] = Gen_fw(t,X,Y_Gen,E,Pm0,Tj) Wn=ones(3,1).*2*pi*60; D=zeros(3,1); dW=zeros(3,1); dDeta=zeros(3,1); W=X(1:3); Deta=X(4:6);

Gen_U=E.*exp(1i*Deta);

Pe=real(Gen_U.*conj(Y_Gen*Gen_U)); dW=(Pm0-Pe-D.*W)./Tj.*Wn; dDeta=W-Wn; Xp=[dW;dDeta]; End

三、 课程总结及心得体会

首先感谢任课老师对我仿真的指导与帮助，然后感谢我的导师对我的指导与关怀，还要感谢在程序编写过程中给过我无私帮助的其他同学。

虽然我在以前也编写过其他的程序，但这次课程的仿真设计实验让我学到了很多。首先，对于错误的查找，举个简单的例子，之前我的雅克比矩阵一直有问题，但找了好久没有找到，后来发现是我的节点类型与老师定义的节点类型不一样，将后面的自己的节点类型检索改过后雅克比矩阵就成功了，瞬时感觉豁然开朗。再次，对于某一个参数的计算可能有不同的方法，例如对于电磁功率如果按照文献Ramnarayan Patel, T. S. Bhatti and D. P. Kothari.MATLAB/Simulink-based transient stability analysis of a multimachine power system中给出的电磁功率计算公式来计算会变得很繁琐，但经过其他同学的帮助与指点，原来电磁

16

于永生 电力系统仿真作业

功率就是发电机输出的功率的实部，这样一来计算量就少了很多。最后，完成了这次作业，对于我本人也是一个极大的激励，我之前觉得这个作业是一个很大的挑战，经过这次完成作业我觉得我战胜了一个难题一样，心情无比愉悦。

但是，我的仿真还有很多不完善的地方，最明显的就是我的δ曲线的角度和老师的有很大误差，这是我今后要完善的地方。

四、 参考文献

[1]何仰赞.温增银. 电力系统分析（上）[M]. 武汉：华中科技大学出版社, 2002 [2]何仰赞.温增银. 电力系统分析（下）[M]. 武汉：华中科技大学出版社, 2002 [3]P. M. Anderson. A. A. Fouad. Power System Control and Stability [M]. (Iowa State University Press,Ames, IA, 1977)

[4]P. Kundur. Power System Stability and Control, EPRI Power System Engineering Series(M). (Mc Graw-Hill, New York, 1994)

17