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

电力系统仿真作业—--—---—————
三机九节点电力系统
暂态仿真
学院:能源与动力工程学院
专业：电力系统及其自动化
学号：
姓名：于永生
导师:
授课教师：

目录
一、 概述 1
二、 课程主要任务 1
1. 系统数据 1
2. 潮流计算 2
3. 负荷等效和支路简化 4
4. 求解电磁功率 5
5. 求解运动方程 5
6. 程序清单 7
(1). 主程序： 7
(2). 极坐标转换成直角坐标函数pol2rect(V,del) 16
(3). 直角坐标转换成极坐标函数rect2pol(Z) 16
(4). 求解微分方程所用的得到微分量的函数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
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. 潮流计算
首先进行潮流计算，采用牛顿拉夫逊迭代法，电力系统潮流计算是电力系统运行和规划中最基本和最经常的计算，其任务是在已知某些运行参数的情况下，计算出系统中全部的运行参数，一般来说，各个母线所供负荷的功率是已知的,各个节点电压是未知的（平衡节点除外),可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵和网络拓扑结构列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值,一般为额定电压，将初值带入功率平衡方程,得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值,将新的初值带入原来的功率平衡方程,并重新形成雅可比矩阵,然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛.
牛顿拉夫逊算法修正方程
W = -JΔV
其中W是节点不平衡量向量，包括有功，无功，电压；J是雅克比矩阵;ΔV是节点电压修正量。
令
,
则直角坐标形式的功率不平衡量方程为
PQ节点：
PV节点：
极坐标形式的功率不平衡量方程
雅可比矩阵J各元素的表达式
当j≠i时:
当j=i时:
其中，
。
进行牛顿拉夫逊算法迭代后得到电压幅值V和相角θ。
3. 负荷等效和支路简化
然后求出支路电流，将发电机内电抗X'加入系统导纳矩阵，求出发电机内电势E’。加入发电机内节点后,系统导纳矩阵变成12＊12阶的矩阵，并将负荷等效成阻抗。然后将支路导纳矩阵分块，如下：
其中,A是3＊3的方阵，E是3＊9的矩阵，C是9＊3的矩阵,D是9＊9的方阵。
经过网络简化得到故障前的3＊3简化导纳矩阵
（1）
其中“inv(D）”是MATLAB中D矩阵的求逆.
故障中导纳矩阵的第七行和第七列从矩阵中删除，此时有
此时,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中给出的电磁功率计算公式为：
本人在此次仿真中用的是公式（2）计算得到的电磁功率。
稳态情况下有，机械功率
Pme=Pe
5. 求解运动方程
发电机的运动方程可以写成常微分方程组:
其中Pmi为第i个机组故障前稳态的电磁功率.在本次仿真中Djωi为零，即阻尼为零。仿真开始，t=0时引入故障,0.083s后切除故障.
求解运动方程后得到ω和δ曲线如下：
3号发电机
2号发电机

图1 ω曲线
1号发电机
3号发电机
2号发电机

图2 δ曲线
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 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
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);%变压器非标准变比
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 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]；
nbus=9;％节点数
bus=N（：，1)；％节点编号
type=N（：，8)；％节点类型
V=N（:,2）;％额定电压
del=N（:，3)；%电压相角
Pg=N(：,4）；%发电机有功功率
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）
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
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； ％ 弧度转化成角度
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)；
% 线路数据
％ 首端 末端 电阻 电抗 电纳(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]；
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）;
％求负荷等效导纳
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］）=［］；
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）;
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，：)；
［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