三機(jī)九節(jié)點(diǎn)電力系統(tǒng)仿真matlab

1、電力系統(tǒng)仿真作業(yè)-三機(jī)九節(jié)點(diǎn)電力系統(tǒng)暫態(tài)仿真學(xué)院：能源與動(dòng)力工程學(xué)院專(zhuān)業(yè)：電力系統(tǒng)及其自動(dòng)化學(xué)號(hào)：姓名：于永生導(dǎo)師：授課教師：目錄一、概述1二、課程主要任務(wù)11.系統(tǒng)數(shù)據(jù)12.潮流計(jì)算23.負(fù)荷等效和支路簡(jiǎn)化44.求解電磁功率55.求解運(yùn)動(dòng)方程56.程序清單7(1).主程序：7(2).極坐標(biāo)轉(zhuǎn)換成直角坐標(biāo)函數(shù)pol2rect(V,del)16(3).直角坐標(biāo)轉(zhuǎn)換成極坐標(biāo)函數(shù)rect2pol(Z)16(4).求解微分方程所用的得到微分量的函數(shù)Gen_fw(t,X,Y_Gen,E,Pm0,Tj)16三、課程總結(jié)及心得體會(huì)16四、參考文獻(xiàn)17于永生電力系統(tǒng)仿真作業(yè)一、 概述在動(dòng)態(tài)穩(wěn)定分析中，系統(tǒng)由線

2、性化的微分方程組和代數(shù)方程組描寫(xiě)，并用經(jīng)典的或現(xiàn)代的線性系統(tǒng)理論來(lái)進(jìn)行穩(wěn)定分析，分析可以在時(shí)域或頻域進(jìn)行。當(dāng)用計(jì)算機(jī)和現(xiàn)代線性系統(tǒng)理論分析時(shí)，常把系統(tǒng)線性化的微分方程組和代數(shù)方程組消去代數(shù)變量，化為狀態(tài)方程形式，并廣泛采用特征分析進(jìn)行穩(wěn)定分析。電力系統(tǒng)是由不同類(lèi)型的發(fā)電機(jī)組、多種電力負(fù)荷、不同電壓等級(jí)的電力網(wǎng)絡(luò)等組成的十分龐大復(fù)雜的動(dòng)力學(xué)系統(tǒng)。其暫態(tài)過(guò)渡過(guò)程不僅包括電磁方面的過(guò)渡過(guò)程，而且還有機(jī)電方面的過(guò)渡過(guò)程。由此可見(jiàn)，電力系統(tǒng)的數(shù)學(xué)模型是一個(gè)強(qiáng)非線性的高維狀態(tài)方程組。在動(dòng)態(tài)穩(wěn)定仿真中使用簡(jiǎn)單的電力系統(tǒng)模型，發(fā)電機(jī)用三階模型表示。二、 課程主要任務(wù)本次課程主要應(yīng)用P. M. Anderson

3、 and A. A. Fouad編寫(xiě)的Power System Control and Stability一書(shū)中所引用的Western System Coordinated Council (WSCC)三機(jī)九節(jié)點(diǎn)系統(tǒng)模型。1. 系統(tǒng)數(shù)據(jù)其中，節(jié)點(diǎn)數(shù)據(jù)如下：%節(jié)點(diǎn)數(shù)據(jù)% 節(jié)點(diǎn) 電壓 電壓 發(fā)電機(jī) 發(fā)電機(jī) 負(fù)荷 負(fù)荷 節(jié)點(diǎn)% 號(hào) 幅值 相角 有功 無(wú)功 有功 無(wú)功 類(lèi)型(1PQ 2PV 3平衡)N= 1 1.040 0.7164 0.2705 0 0 3 2 1.0250 1.63 0.0665 0 0 2 3 1.0250 0.85 -0.1086 0 0 2 4 1 0 0 0 0 0 1 5

4、 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;其中，支路數(shù)據(jù)如下：% 線路數(shù)據(jù)% 首端 末端 電阻 電抗 電納（1/2) 變壓器非標(biāo)準(zhǔn)變比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

5、 1 2 7 0 0.0625 0 1 3 9 0 0.0586 0 1;發(fā)電機(jī)數(shù)據(jù)如下：% 發(fā)電機(jī) 母線 Xd Xd Td0 Xq Xq Tq0 Tj XfGe= 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;系統(tǒng)電路結(jié)構(gòu)拓?fù)鋱D如下：圖1 WSCC 3機(jī)9節(jié)點(diǎn)系統(tǒng)（所有參數(shù)以100MVA為基準(zhǔn)值的標(biāo)幺值）2. 潮流計(jì)

6、算首先進(jìn)行潮流計(jì)算，采用牛頓拉夫遜迭代法，電力系統(tǒng)潮流計(jì)算是電力系統(tǒng)運(yùn)行和規(guī)劃中最基本和最經(jīng)常的計(jì)算，其任務(wù)是在已知某些運(yùn)行參數(shù)的情況下，計(jì)算出系統(tǒng)中全部的運(yùn)行參數(shù)，一般來(lái)說(shuō)，各個(gè)母線所供負(fù)荷的功率是已知的，各個(gè)節(jié)點(diǎn)電壓是未知的（平衡節(jié)點(diǎn)除外），可以根據(jù)網(wǎng)絡(luò)結(jié)構(gòu)形成節(jié)點(diǎn)導(dǎo)納矩陣，然后由節(jié)點(diǎn)導(dǎo)納矩陣和網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)列寫(xiě)功率方程，由于功率方程里功率是已知的，電壓的幅值和相角是未知的，這樣潮流計(jì)算的問(wèn)題就轉(zhuǎn)化為求解非線性方程組的問(wèn)題了。為了便于用迭代法解方程組，需要將上述功率方程改寫(xiě)成功率平衡方程，并對(duì)功率平衡方程求偏導(dǎo)，得出對(duì)應(yīng)的雅可比矩陣，給未知節(jié)點(diǎn)賦電壓初值，一般為額定電壓，將初值帶入功率平衡方

7、程，得到功率不平衡量，這樣由功率不平衡量、雅可比矩陣、節(jié)點(diǎn)電壓不平衡量（未知的）構(gòu)成了誤差方程，解誤差方程，得到節(jié)點(diǎn)電壓不平衡量，節(jié)點(diǎn)電壓加上節(jié)點(diǎn)電壓不平衡量構(gòu)成新的節(jié)點(diǎn)電壓初值，將新的初值帶入原來(lái)的功率平衡方程，并重新形成雅可比矩陣，然后計(jì)算新的電壓不平衡量，這樣不斷迭代，不斷修正，一般迭代三到五次就能收斂。牛頓拉夫遜算法修正方程W = -JV其中W是節(jié)點(diǎn)不平衡量向量，包括有功，無(wú)功，電壓；J是雅克比矩陣；V是節(jié)點(diǎn)電壓修正量。令，則直角坐標(biāo)形式的功率不平衡量方程為PQ節(jié)點(diǎn)：PV節(jié)點(diǎn)：極坐標(biāo)形式的功率不平衡量方程雅可比矩陣J各元素的表達(dá)式當(dāng)ji時(shí)：當(dāng)j=i時(shí)：其中，。進(jìn)行牛頓拉夫遜算法迭代后得

8、到電壓幅值V和相角。3. 負(fù)荷等效和支路簡(jiǎn)化然后求出支路電流，將發(fā)電機(jī)內(nèi)電抗X加入系統(tǒng)導(dǎo)納矩陣，求出發(fā)電機(jī)內(nèi)電勢(shì)E。加入發(fā)電機(jī)內(nèi)節(jié)點(diǎn)后，系統(tǒng)導(dǎo)納矩陣變成12*12階的矩陣，并將負(fù)荷等效成阻抗。然后將支路導(dǎo)納矩陣分塊，如下：其中，A是3*3的方陣，E是3*9的矩陣，C是9*3的矩陣，D是9*9的方陣。經(jīng)過(guò)網(wǎng)絡(luò)簡(jiǎn)化得到故障前的3*3簡(jiǎn)化導(dǎo)納矩陣 （1）其中“inv（D）”是MATLAB中D矩陣的求逆。故障中導(dǎo)納矩陣的第七行和第七列從矩陣中刪除，此時(shí)有此時(shí)，A是3*3的方陣，E是3*8的矩陣，C是8*3的矩陣，D是8*8的方陣。簡(jiǎn)化矩陣求法如同公式（1）。故障后的Y矩陣相對(duì)于故障前的Y矩陣只是第5個(gè)

9、節(jié)點(diǎn)和第7個(gè)節(jié)點(diǎn)有變化，反映到12*12的矩陣中即為第（10,8），（8,10），（8,8），（10,10）位置的元素有變化，其中（10,8），（8,10）位置的元素變?yōu)榱悖?8,8)，（10,10）節(jié)點(diǎn)在故障前的基礎(chǔ)上加上（8,10）位置元素的值。然后簡(jiǎn)化導(dǎo)納矩陣的求法同式（1）。4. 求解電磁功率得到故障前，故障中，故障后三個(gè)不同的導(dǎo)納矩陣后，就開(kāi)始計(jì)算電磁功率和機(jī)械功率，機(jī)械功率等于穩(wěn)態(tài)的電磁功率中的有功分量。所以可以有Pe=real(E*I) (2)公式（2）中，E為發(fā)電機(jī)內(nèi)電勢(shì)，I為從發(fā)電機(jī)流出的電流。但在參考文獻(xiàn)Ramnarayan Patel, T. S. Bhatti and

10、D. P. Kothari.MATLAB/Simulink-based transient stability analysis of a multimachine power system中給出的電磁功率計(jì)算公式為：本人在此次仿真中用的是公式（2）計(jì)算得到的電磁功率。穩(wěn)態(tài)情況下有，機(jī)械功率Pme=Pe5. 求解運(yùn)動(dòng)方程發(fā)電機(jī)的運(yùn)動(dòng)方程可以寫(xiě)成常微分方程組：其中Pmi為第i個(gè)機(jī)組故障前穩(wěn)態(tài)的電磁功率。在本次仿真中Dji為零，即阻尼為零。仿真開(kāi)始，t=0時(shí)引入故障，0.083s后切除故障。求解運(yùn)動(dòng)方程后得到和曲線如下：3號(hào)發(fā)電機(jī)2號(hào)發(fā)電機(jī)圖1 曲線1號(hào)發(fā)電機(jī)3號(hào)發(fā)電機(jī)2號(hào)發(fā)電機(jī)圖2 曲線6. 程

11、序清單以下是我編寫(xiě)的仿真程序，除主程序外還包含三個(gè)函數(shù)：極坐標(biāo)轉(zhuǎn)換成直角坐標(biāo)函數(shù)pol2rect(V,del),其中輸入?yún)?shù)V為極坐標(biāo)向量的幅值，del為極坐標(biāo)向量的角度，函數(shù)返回值為一個(gè)復(fù)數(shù)，即為極坐標(biāo)向量在直角坐標(biāo)中的復(fù)數(shù)值；直角坐標(biāo)轉(zhuǎn)換成極坐標(biāo)函數(shù)rect2pol(Z),其中輸入?yún)?shù)Z為一個(gè)復(fù)數(shù)，函數(shù)返回值為一個(gè)二維向量，向量的第一個(gè)數(shù)為幅值，第二個(gè)數(shù)為相角；求解微分方程所用的得到微分量的函數(shù)Gen_fw(t,X,Y_Gen,E,Pm0,Tj),其中輸入?yún)?shù)X為和的迭代初值，Y_Gen，為求解所用的導(dǎo)納矩陣，這里是三階的方陣，對(duì)應(yīng)于故障前，故障中，和故障后的三個(gè)Y矩陣，E為發(fā)電機(jī)內(nèi)電勢(shì)，

12、Pm0為機(jī)械功率，等于穩(wěn)態(tài)時(shí)的有功功率，Tj為運(yùn)動(dòng)方程中的2H。(1). 主程序：clc;clear;%節(jié)點(diǎn)數(shù)據(jù)% 節(jié)點(diǎn) 電壓 電壓 發(fā)電機(jī) 發(fā)電機(jī) 負(fù)荷 負(fù)荷 節(jié)點(diǎn)% 號(hào) 幅值 相角 有功 無(wú)功 有功 無(wú)功 類(lèi)型(1PQ 2PV 3平衡% N= 1 1.040 0.71640.2705 0 0 3% 2 1.0250 1.63 0.0665 0 0 2% 3 1.0250 0.85 -0.1086 0 0 2% 4 1 0 0 0 0 0 1% 5 1 0 0 0 1.250.5 1% 6 1 0 0 0 0.90.3 1% 7 1 0 0 0 0 0 1% 8 1 0 0 0 1 0.35

13、 1% 9 1 0 0 0 0 0 1;% 線路數(shù)據(jù)% 首端 末端 電阻 電抗 電納（1/2) 變壓器非標(biāo)準(zhǔn)變比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;%形成節(jié)點(diǎn)導(dǎo)納矩陣fb=L(:,1);%起始節(jié)點(diǎn)編號(hào)tb=L(:,2);%末節(jié)點(diǎn)編號(hào)r=

14、L(:,3);%電阻x=L(:,4);%電抗b=L(:,5);%電納a=L(:,6);%變壓器非標(biāo)準(zhǔn)變比z=r + i*x;%z矩陣，復(fù)數(shù)，9乘1的矩陣y=1./z;%求倒數(shù)，9乘1的導(dǎo)納矩陣b=i*b;%虛數(shù)nb=max(max(fb),max(tb);%nb=9nl=length(fb);%支路個(gè)數(shù)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方陣的元素，非對(duì)角元（4,5）（4，6）（5，7）（6,9）（7，8）（8,9）（1,4）（2,7）（3,9） Y(tb(k),fb(k)=Y

15、(fb(k),tb(k);%Y方陣的元素，上一句的對(duì)稱(chēng)部分（5,4）（6,4）（7，5）（9，6）（8，7）（9,8）（4,1）（7,2）（9,3）endfor m=1:nb for n=1:nl if fb(n)=m%支路中的第n個(gè)數(shù)的首節(jié)點(diǎn)編號(hào)等于m Y(m,m)=Y(m,m)+y(n)/(a(n)2)+b(n);%Y矩陣中的對(duì)角元 elseif tb(n)=m%支路中的第n個(gè)數(shù)的末節(jié)點(diǎn)編號(hào)等于m Y(m,m)=Y(m,m)+y(n)+b(n);%因?yàn)橹肥鞘锥嘶蚰┒说扔趍，不能是首端編號(hào)等于末端編號(hào) end endend%節(jié)點(diǎn)數(shù)據(jù)% 節(jié)點(diǎn) 電壓 電壓 發(fā)電機(jī) 發(fā)電機(jī) 負(fù)荷 負(fù)荷 節(jié)點(diǎn)%

16、號(hào) 幅值 相角 有功 無(wú)功 有功 無(wú)功 類(lèi)型(1PQ 2PV 3平衡N= 1 1.040 0.71640.2705 0 0 3 2 1.0250 1.63 0.0665 0 0 2 3 1.0250 0.85 -0.1086 0 0 2 4 1 0 0 0 0 0 1 5 1 0 0 0 1.250.5 1 6 1 0 0 0 0.90.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;%節(jié)點(diǎn)數(shù)bus=N(:,1);%節(jié)點(diǎn)編號(hào)type=N(:,8);%節(jié)點(diǎn)類(lèi)型V=N(:,2);%額定電壓del=N(:,3);%電壓相角P

17、g=N(:,4);%發(fā)電機(jī)有功功率Qg=N(:,5);%發(fā)電機(jī)無(wú)功功率Pl=N(:,6);%負(fù)荷有功功率Ql=N(:,7);%負(fù)荷無(wú)功功率P=Pg-Pl;%發(fā)出減消耗Q=Qg-Ql;%發(fā)出減消耗Psp=P;%額定有功功率Qsp=Q;%額定無(wú)功功率G=real(Y);B=imag(Y);pv=find(type=3|type=2);%pv節(jié)點(diǎn)編號(hào)和平衡節(jié)點(diǎn)編號(hào)pq=find(type=1);%pq節(jié)點(diǎn)編號(hào)npv=length(pv);%npv=3npq=length(pq);%npq=6Tol=1;Iter=1;while(Tol1e-7) P=zeros(nbus,1);%迭代初值P為零 Q=

18、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;%額定無(wú)功功率減去流出或流入的無(wú)功功率 k=1; dQ=zeros(npq,1);%pq節(jié)點(diǎn)的無(wú)功功率，pq乘1的矩陣

19、,6乘1 for i=1:nbus if type(i)=1%pq節(jié)點(diǎn) dQ(k,1)=dQa(i); k=k+1; end end dP=dPa(2:nbus); Z=dP;dQ; %下面計(jì)算雅克比矩陣 J1=zeros(nbus-1,nbus-1);%有功功率對(duì)功角求偏微分 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)=J

20、1(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);%有功功率對(duì)電壓求偏微分 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)+

21、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);%無(wú)功功率對(duì)相角求偏微分 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(

22、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);%無(wú)功功率對(duì)電壓求偏微分 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

23、(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=JZ; 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);endVm = pol2rect(V,del); % 極坐標(biāo)轉(zhuǎn)換成直角坐

24、標(biāo)Del = 180/pi*del; % 弧度轉(zhuǎn)化成角度Iij = zeros(nb,nb);%支路電流for m = 1:nl%每一條支路都過(guò)一遍 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);endIij = sparse(Iij); % 線路數(shù)據(jù) % 首端 末端 電阻 電抗 電納（1/2) 變壓器非標(biāo)準(zhǔn)變比 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.

25、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);endfor k=1:9 for m=1:9 Yzengjia(3+k,3+m)=Y(k,m); endendYzengjia(1,4)=-Yz

26、engjia(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);%求負(fù)荷等效導(dǎo)納Y5=N(5,6)/(V

27、(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); endendE=zeros(3,9);for m=1:3 for n=4:12 E(m,n-3)=Yzengj

28、ia(m,n); endendC=zeros(9,3);for m=4:12 for n=1:3 C(m-3,n)=Yzengjia(m,n); endendD=zeros(9,9);for m=4:12 for n=4:12 D(m-3,n-3)=Yzengjia(m,n); endendYpre=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

29、(m,n+3);%因?yàn)镋dur的第四到第九列是零，而所以這里為了方便，直接取E的前六列 endendCdur=zeros(8,3);for m=1:8 for n=1:3 Cdur(m,n)=Yzengjiadur(m+3,n);%因?yàn)镃dur的第四到第九行是零，而所以這里為了方便，直接取C的前六行 endendDdur=zeros(8,8);for m=1:8 for n=1:8 Ddur(m,n)=Yzengjiadur(m+3,n+3); endendYdur=Adur-Edur*(inv(Ddur)*Cdur;%故障中的3*3矩陣Yzengjiaaft=Yzengjia;temp=Yz

30、engjiaaft(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); endendEaft=zeros(3,9);for m=1:3 for n=4:12 Eaft(m,n-3)=Yzengjiaaft(m,n); endendCaft=zeros(9,3);for

31、 m=4:12 for n=1:3 Caft(m-3,n)=Yzengjiaaft(m,n); endendDaft=zeros(9,9);for m=1:9 for n=1:9 Daft(m,n)=Yzengjiaaft(m+3,n+3); endendYaft=Aaft-Eaft*(inv(Daft)*Caft;%故障后的3*3矩陣%求電磁功率Enei=ones(3,1);%發(fā)電機(jī)內(nèi)電勢(shì)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

32、);polEnei1=rect2pol(Enei(1);%直角坐標(biāo)轉(zhuǎn)換成極坐標(biāo)后得到角度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;%

33、機(jī)械功率=故障前的電磁功率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,

34、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). 極坐標(biāo)轉(zhuǎn)換成直角坐標(biāo)函數(shù)pol2rect(V,del)function rect = pol2rect(rho,theta)rect = rho.*cos(theta) + j*rho.*sin(theta);(3). 直角坐標(biāo)轉(zhuǎn)換成極坐標(biāo)函數(shù)rect2pol(Z)function pol=rect2pol(Z)t=real(Z);k=imag(Z);pol(1)=(t2+k2)(1/2);%模pol(2)=atan(k/t);(4). 求解微分方程所用的得到微分量的函數(shù)Gen_fw(t,X,Y_