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

1、電力系統(tǒng)仿真作業(yè)-三機(jī)九節(jié)點(diǎn)電力系統(tǒng)暫態(tài)仿真學(xué)院：能源與動力工程學(xué)院專業(yè)：電力系統(tǒng)及其自動化學(xué)號：姓名：于永生導(dǎo)師：授課教師：目錄一、概述1二、課程主要任務(wù)11.系統(tǒng)數(shù)據(jù)12.潮流計(jì)算23.負(fù)荷等效和支路簡化44.求解電磁功率55.求解運(yùn)動方程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é)及心得體會16四、參考文獻(xiàn)17一、 概述在動態(tài)穩(wěn)定分析中，系統(tǒng)由線性化的微分方程組和代數(shù)

2、方程組描寫，并用經(jīng)典的或現(xiàn)代的線性系統(tǒng)理論來進(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)是由不同類型的發(fā)電機(jī)組、多種電力負(fù)荷、不同電壓等級的電力網(wǎng)絡(luò)等組成的十分龐大復(fù)雜的動力學(xué)系統(tǒng)。其暫態(tài)過渡過程不僅包括電磁方面的過渡過程，而且還有機(jī)電方面的過渡過程。由此可見，電力系統(tǒng)的數(shù)學(xué)模型是一個(gè)強(qiáng)非線性的高維狀態(tài)方程組。在動態(tài)穩(wěn)定仿真中使用簡單的電力系統(tǒng)模型，發(fā)電機(jī)用三階模型表示。二、 課程主要任務(wù)本次課程主要應(yīng)用P. M. Anderson and A. A.

3、Fouad編寫的Power System Control and Stability一書中所引用的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)% 號 幅值 相角 有功 無功 有功 無功 類型(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 1 0 0 0 1.

4、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 1 2 7 0 0.

5、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ì)算首先進(jìn)行潮流計(jì)算，采

6、用牛頓拉夫遜迭代法，電力系統(tǒng)潮流計(jì)算是電力系統(tǒng)運(yùn)行和規(guī)劃中最基本和最經(jīng)常的計(jì)算，其任務(wù)是在已知某些運(yùn)行參數(shù)的情況下，計(jì)算出系統(tǒng)中全部的運(yùn)行參數(shù)，一般來說，各個(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)列寫功率方程，由于功率方程里功率是已知的，電壓的幅值和相角是未知的，這樣潮流計(jì)算的問題就轉(zhuǎn)化為求解非線性方程組的問題了。為了便于用迭代法解方程組，需要將上述功率方程改寫成功率平衡方程，并對功率平衡方程求偏導(dǎo)，得出對應(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)電壓初值，將新的初值帶入原來的功率平衡方程，并重新形成雅可比矩陣，然后計(jì)算新的電壓不平衡量，這樣不斷迭代，不斷修正，一般迭代三到五次就能收斂。牛頓拉夫遜算法修正方程W = -JV其中W是節(jié)點(diǎn)不平衡量向量，包括有功，無功，電壓；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)行牛頓拉夫遜算法迭代后得到電壓幅值V和相角。3

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

9、，反映到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）位置元素的值。然后簡化導(dǎo)納矩陣的求法同式（1）。4. 求解電磁功率得到故障前，故障中，故障后三個(gè)不同的導(dǎo)納矩陣后，就開始計(jì)算電磁功率和機(jī)械功率，機(jī)械功率等于穩(wěn)態(tài)的電磁功率中的有功分量。所以可以有Pe=real(E*I) (2)公式（2）中，E為發(fā)電機(jī)內(nèi)電勢，I為從發(fā)電機(jī)流出的電流。但在參考文獻(xiàn)Ramnarayan Patel, T. S. Bhatti and D. P. Kotha

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

11、真程序，除主程序外還包含三個(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)納矩陣，這里是三階的方陣，對應(yīng)于故障前，故障中，和故障后的三個(gè)Y矩陣，E為發(fā)電機(jī)內(nèi)電勢，Pm0為機(jī)械功率，等于

12、穩(wěn)態(tài)時(shí)的有功功率，Tj為運(yùn)動方程中的2H。(1). 主程序：clc;clear;%節(jié)點(diǎn)數(shù)據(jù)% 節(jié)點(diǎn) 電壓 電壓 發(fā)電機(jī) 發(fā)電機(jī) 負(fù)荷 負(fù)荷 節(jié)點(diǎn)% 號 幅值 相角 有功 無功 有功 無功 類型(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

13、 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)編號tb=L(:,2);%末節(jié)點(diǎn)編號r=L(:,3);%電阻x

14、=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方陣的元素，非對角元（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

15、);%Y方陣的元素，上一句的對稱部分（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)編號等于m Y(m,m)=Y(m,m)+y(n)/(a(n)2)+b(n);%Y矩陣中的對角元 elseif tb(n)=m%支路中的第n個(gè)數(shù)的末節(jié)點(diǎn)編號等于m Y(m,m)=Y(m,m)+y(n)+b(n);%因?yàn)橹肥鞘锥嘶蚰┒说扔趍，不能是首端編號等于末端編號 end endend%節(jié)點(diǎn)數(shù)據(jù)% 節(jié)點(diǎn) 電壓 電壓 發(fā)電機(jī) 發(fā)電機(jī) 負(fù)荷 負(fù)荷 節(jié)點(diǎn)% 號 幅值 相角 有功

16、無功 有功 無功 類型(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)編號type=N(:,8);%節(jié)點(diǎn)類型V=N(:,2);%額定電壓del=N(:,3);%電壓相角Pg=N(:,4);%發(fā)

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

18、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節(jié)點(diǎn)的無功功率，pq乘1的矩陣,6乘1 for i=

19、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);%有功功率對功角求偏微分 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)

20、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)

21、; 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

22、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)*

23、(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)換成直角坐標(biāo)Del = 180/

24、pi*del; % 弧度轉(zhuǎn)化成角度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);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.153 1 6 9 0

25、.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)=-Yzengjia(1,1)

26、;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(5)2-j*(N(5

27、,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)=Yzengjia(m,n); en

28、dendC=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(m,n+3);%因?yàn)?/p>

29、Edur的第四到第九列是零，而所以這里為了方便，直接取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=Yzengjiaaft(8

30、,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 m=4:12 for

31、 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)電勢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=

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

33、功率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

34、(:,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