發(fā)電機(jī)最優(yōu)勵磁控制系統(tǒng)設(shè)計(jì)

1、發(fā)電機(jī)最優(yōu)勵磁控制系統(tǒng)設(shè)計(jì)1 線性最優(yōu)控制理論1.1 線性最優(yōu)控制原理線性最優(yōu)控制理論所研究的核心問題是選取最優(yōu)的控制規(guī)律，使控制系統(tǒng)在特定指標(biāo)下的性能為最優(yōu)?？刂葡到y(tǒng)框如圖1所示。圖1 最優(yōu)控制框圖線性定常系統(tǒng)狀態(tài)空間方程的一般形式為： （1）式中，為狀態(tài)系數(shù)矩陣；為控制系數(shù)矩陣；為n維狀態(tài)向量；為r維控制向量。如圖1所示，如果要改善系統(tǒng)性能，可以引入狀態(tài)反饋構(gòu)成閉環(huán)系統(tǒng)。反饋系統(tǒng)的狀態(tài)向量為： （2）式中，為狀態(tài)反饋增益矩陣。將式（1）和式（2）合并，可以得到： （3）可見，最優(yōu)控制的本質(zhì)就是如何選取反饋矩陣，以使它在給定控制規(guī)律下達(dá)到特定條件下的最優(yōu)。1.2 二次型性能指標(biāo)假定為系統(tǒng)的實(shí)

2、際響應(yīng)，為系統(tǒng)預(yù)期的響應(yīng)。最優(yōu)控制性能指標(biāo)是使兩者的偏差最小，即 （4）式中，是隨函數(shù)而變化的一個(gè)泛函數(shù)，是在時(shí)間內(nèi)求取偏差平方的積分，稱為二次型指標(biāo)。如果以為實(shí)際狀態(tài)向量，以為預(yù)期狀態(tài)向量，則要求狀態(tài)向量偏差最小的二次型性能指標(biāo)為： （5）在二次型性能指標(biāo)中，也需要引入對控制量的約束。如果沒有這個(gè)約束，所設(shè)計(jì)的控制器中，控制量的變化范圍可能會很大，難以實(shí)現(xiàn)。引入控制量約束的二次型指標(biāo)為： （6）式中和表示狀態(tài)向量和控制向量的權(quán)矩陣。為了便于分析，通常把系統(tǒng)平衡點(diǎn)置于狀態(tài)空間的原點(diǎn)，即，則式（6）可以變換為： （7）以上式作為性能指標(biāo)設(shè)計(jì)最優(yōu)控制系統(tǒng)，可以證明這個(gè)最優(yōu)控制規(guī)律是存在且唯一的，表

3、達(dá)式為： （8）式中，為最優(yōu)反饋增益矩陣。為維對稱常矩陣，矩陣是黎卡梯方程：的解。2 單機(jī)無窮大系統(tǒng)數(shù)學(xué)模型與無窮大系統(tǒng)并聯(lián)運(yùn)行的同步發(fā)電機(jī)的基本數(shù)學(xué)模型包括4個(gè)部分：與系統(tǒng)并聯(lián)的同步發(fā)電機(jī)基本方程組，發(fā)電機(jī)電勢、定子電壓和電流之間的聯(lián)系方程，功率方程以及轉(zhuǎn)子繞組動態(tài)方程。2.1與系統(tǒng)并聯(lián)的同步發(fā)電機(jī)基本方程組與系統(tǒng)并聯(lián)的同步發(fā)電機(jī)基本方程組決定了運(yùn)行中的發(fā)電機(jī)的穩(wěn)態(tài)與過渡過程的全部動力學(xué)特性。轉(zhuǎn)子運(yùn)動方程的形式如下： （9）式中，用秒，用弧度，用秒，、均為標(biāo)幺值。在實(shí)用計(jì)算中，一般采用，這里的是機(jī)組的阻尼系數(shù)，用標(biāo)幺值，單位用弧度。2.2發(fā)電機(jī)電勢、定子電壓和電流之間的聯(lián)系方程將觀察發(fā)電機(jī)的

4、坐標(biāo)系從abc坐標(biāo)系轉(zhuǎn)換到dq0坐標(biāo)系，可以極大地簡化發(fā)電機(jī)的數(shù)學(xué)模型。這是因?yàn)閍bc坐標(biāo)系固定在三相對稱的定子上，由于轉(zhuǎn)子的旋轉(zhuǎn)，對于每一坐標(biāo)軸方向上的磁路磁導(dǎo)率都是時(shí)間t的周期函數(shù)，不便于計(jì)算和分析。而dq0坐標(biāo)系是固定在轉(zhuǎn)子上，坐標(biāo)軸與磁路相對靜止，這樣磁導(dǎo)率就是一常數(shù)，因而惦記的各個(gè)參數(shù)也會成為與時(shí)間t無關(guān)的常數(shù)。由發(fā)電機(jī)電勢與電流矢量圖可以直接得到如下關(guān)系式組： （10）由上述的發(fā)電機(jī)定子電勢、電壓和電流的關(guān)系式組，可以解出定子電流為： （11） （12） （13）由矢量圖又可以得到： （14） （15） （16）代入式（11）和（12），可以解得： （17）2.3 功率方程與無窮大

5、系統(tǒng)并聯(lián)的發(fā)電機(jī)有功功率的表達(dá)式為： （18）以和為自變量的電磁功率表達(dá)式如下： （19）2.4 轉(zhuǎn)子繞組動態(tài)方程在列寫轉(zhuǎn)子繞組的實(shí)用動態(tài)方程時(shí)，忽略阻尼繞組，只列寫勵磁繞組動態(tài)方程。另外假設(shè)定子繞組中磁鏈?zhǔn)娇梢酝蛔兊?，勵磁繞組磁鏈?zhǔn)睾?，這是因?yàn)槎ㄗ永@組中暫態(tài)電流的非周期分量比勵磁繞組非周期分量衰減快得多。帶負(fù)載的發(fā)電機(jī)在轉(zhuǎn)子d軸上的總磁鏈為： （20）式中，是有用磁通；是轉(zhuǎn)子漏磁通；是電樞反應(yīng)磁通；是勵磁繞組與定子d軸繞組的互感；是勵磁繞組的漏感。假定過渡過程中轉(zhuǎn)速標(biāo)幺值為，式（20）可以寫為： （21）可以求得： （22）此處，為勵磁繞組的自感抗。假設(shè)定子d軸總磁鏈與端電壓q軸分量標(biāo)幺值相

6、等。類似于轉(zhuǎn)子繞組，可以得到： （23）這里為發(fā)電機(jī)d軸同步電抗，其中是定子漏電抗；為d軸暫態(tài)電抗；稱為d軸暫態(tài)電抗后電勢。由勵磁繞組等值電路可得回路方程式： （24）以乘以上式，可得在空載情況下勵磁繞組的動態(tài)方程： （25）上式中，為定子繞組開路時(shí)勵磁繞組的時(shí)間常數(shù)。若以乘以式（24），可得： （26）以上就是得到了主發(fā)電機(jī)勵磁繞組的動態(tài)方程。2.5 基本方程組的偏差化與線性化由式（9）可以得到偏差方程： （27）根據(jù)式（18）和式（19），將電磁功率方程進(jìn)行偏差化和線性化，可得： （28） （29）式中，；。將式（28）和（29）代入式（27），得到線性化的轉(zhuǎn)子運(yùn)動方程： （30） （31

7、）由，可得到運(yùn)行角偏差和機(jī)端電壓偏差表示的電磁功率偏差： （32）其中，。于是轉(zhuǎn)子的運(yùn)動偏差方程有： （33）勵磁繞組的動態(tài)偏差方程為： （34）3 最優(yōu)勵磁控制系統(tǒng)設(shè)計(jì) 由式（33）和式（34）可以得到系統(tǒng)狀態(tài)方程為： （35） 由于不便于測量，所以采用機(jī)端電壓來代替，即 （36）又由于在實(shí)際系統(tǒng)中采用的是、和的一組狀態(tài)量，通過狀態(tài)變量置換，再結(jié)合式（35），可以得到、和所表示的狀態(tài)方程： （37）如果發(fā)電機(jī)采用可控硅勵磁系統(tǒng)，則可以忽略勵磁時(shí)間常數(shù)，那么控制量就是勵磁繞組電壓，系統(tǒng)狀態(tài)方程可以表示為形式的三階形式： （38）最優(yōu)控制量4 最優(yōu)控制設(shè)計(jì)算例發(fā)電機(jī)電氣參數(shù)：隱極機(jī)縱軸電抗，縱軸

8、暫態(tài)電抗，縱軸開環(huán)時(shí)間常數(shù)，綜合阻尼系數(shù)；變壓器參數(shù)：；線路參數(shù)：；系統(tǒng)參數(shù)：。所研究的發(fā)輸電系統(tǒng)的示意圖見圖2。圖2 系統(tǒng)示意圖設(shè)計(jì)最優(yōu)運(yùn)行點(diǎn)為，時(shí)的最優(yōu)控制器。4.1 計(jì)算最優(yōu)控制參數(shù)通過式（18）和式（19），根據(jù)最優(yōu)運(yùn)行點(diǎn)和，計(jì)算得到最優(yōu)運(yùn)行點(diǎn)的相關(guān)參數(shù)如下：由式（18）可得：由式（19）可得：由式（17）可得：數(shù)據(jù)代入到相應(yīng)公式，可以求得；選取和作為權(quán)矩陣，進(jìn)而可以得到狀態(tài)反饋矩陣：可知最優(yōu)控制量。上述計(jì)算過程使用MATLAB完成，見報(bào)告附錄中的程序1。4.2 驗(yàn)證最優(yōu)勵磁控制系統(tǒng)性能由于勵磁控制系統(tǒng)參數(shù)是采用線性化狀態(tài)方程而得到的，而實(shí)際系統(tǒng)的狀態(tài)方程是非線性的，所以使用MATLA

9、B的ode45算法，對單機(jī)系統(tǒng)中發(fā)電機(jī)的運(yùn)動微分方程進(jìn)行計(jì)算，通過觀察發(fā)電機(jī)的轉(zhuǎn)子角和輸出電磁功率的變化情況，得到最優(yōu)勵磁控制系統(tǒng)的控制性能。單機(jī)系統(tǒng)的發(fā)電機(jī)運(yùn)動方程如下：此運(yùn)動方程以轉(zhuǎn)子角、角速度和q軸暫態(tài)電勢作為三個(gè)狀態(tài)量。微分方程中，；作為輸入量，引入最優(yōu)狀態(tài)反饋矩陣，即用三個(gè)狀態(tài)量進(jìn)行表示，為其中， 和作為與和有關(guān)的變量，分別用式（19）和式（17）表示。4.2.1 靜態(tài)穩(wěn)定性系統(tǒng)在運(yùn)行至5s時(shí)刻，加入一個(gè)時(shí)間長達(dá)1秒的電磁功率的擾動，模擬負(fù)荷的突然變化，即5s到6s之間從原來的0.55增加到0.6，6s之后回到最初設(shè)計(jì)的穩(wěn)態(tài)運(yùn)行點(diǎn)。進(jìn)行小擾動的仿真。仿真程序見報(bào)告附錄中的程序27和1

10、3。仿真的功角和電磁功率波形圖如下：圖3 靜態(tài)穩(wěn)定功角仿真圖從波形圖可以看出，在發(fā)生小擾動后，采用最優(yōu)控制理論所設(shè)計(jì)的勵磁系統(tǒng)，能夠快速地平息系統(tǒng)的振蕩，轉(zhuǎn)子角和電磁功率均能夠迅速地回到最初設(shè)計(jì)的最優(yōu)運(yùn)行點(diǎn)，提高了電力系統(tǒng)的靜態(tài)穩(wěn)定性。可見，發(fā)電機(jī)加入最優(yōu)勵磁控制器可以有效地改善系統(tǒng)穩(wěn)定性。4.2.2 暫態(tài)穩(wěn)定性系統(tǒng)運(yùn)行至5s的時(shí)刻，在其中一條線路首端即圖2中的d處，發(fā)生三相短路故障，故障0.1s后保護(hù)動作切除故障線路。程序見報(bào)告附錄中的程序2、3和812。該過程中的功角變化波形圖如下：圖4 暫態(tài)穩(wěn)定功角仿真圖由圖4可以看出，系統(tǒng)在5s時(shí)刻發(fā)生三相短路后功角和電磁功率急劇上升，5.1s切除故障

11、線路后系統(tǒng)迅速地恢復(fù)穩(wěn)定，整個(gè)過程不超過5秒鐘。因此，最優(yōu)勵磁控制器也可以很好地維持電力系統(tǒng)的暫態(tài)穩(wěn)定性。5 結(jié)論 本報(bào)告中先介紹了最優(yōu)控制理論的相關(guān)原理，然后通過對與電網(wǎng)并聯(lián)的發(fā)電機(jī)基本方程組的研究從而建立了發(fā)電機(jī)與勵磁系統(tǒng)的線性化狀態(tài)方程，采用MATLAB計(jì)算最優(yōu)狀態(tài)反饋參數(shù)并代入到系統(tǒng)的非線性微分方程組，然后對系統(tǒng)的非線性微分方程組進(jìn)行了時(shí)域仿真。仿真結(jié)果表明，最優(yōu)勵磁控制器可以很好地提高系統(tǒng)的靜態(tài)穩(wěn)定性和暫態(tài)穩(wěn)定性。可見，最優(yōu)勵磁控制器有效地增加了發(fā)電機(jī)勵磁系統(tǒng)的阻尼，抑制了系統(tǒng)受到擾動時(shí)的振蕩，加快系統(tǒng)恢復(fù)穩(wěn)定的速度。附 錄程序1：%輸入發(fā)電機(jī)與系統(tǒng)的相關(guān)參數(shù)Xd=2.5;Xdz=0

12、.5;XT=0.01;Xq=Xd;XL=0.5*1.5;Td0=10;Vs=1;D=5;H=8;Pe=0.55;ang=65*pi/180;f0=50;%中間計(jì)算過程Xe=XT+XL;Xds=Xd+XT+XL;Xdzs=Xdz+XT+XL;Xqs=Xq+XT+XL;Eq=Pe*Xds/(Vs*sin(ang)Eqz=(Pe-0.5*Vs2*(Xdzs-Xqs)/(Xdzs*Xqs)*sin(2*ang)*Xdzs/(Vs*sin(ang)Td=Td0*Xdzs/Xds;SE=Eq*Vs*cos(ang)/Xds+Vs*Vs*(Xds-Xqs)/(Xds*Xqs)*cos(2*ang);SEz=

13、Eqz*Vs*cos(ang)/Xdzs+Vs*Vs*(Xdzs-Xqs)/(Xdzs*Xqs)*cos(2*ang);RE=Vs*sin(ang)/Xds;REz=Vs*sin(ang)/Xdzs;Vt=sqrt(Eq2*Xe2+(Vs*cos(ang)*Xd)2+2*Xe*Xd*Eq*Vs*cos(ang)/Xds2+(Vs*sin(ang)*Xq/Xqs)2)dVtEq=(Eq*Xe2+Xe*Xd*Vs*cos(ang)/(Xds2)*(Eq2*Xe2+Vs2*(cos(ang)2*Xd2+2*Xe*Xd*Eq*Vs*cos(ang)/(Xds2)+(Vs*sin(ang)*Xq/Xqs

14、)2)(-0.5);dVtang=0.5*(Vs2*Xq2*sin(2*ang)/Xqs2-(Vs2*Xd2*sin(2*ang)+2*Xe*Xd*Eq*Vs*sin(ang)/Xds2)*(Eq2*Xe2+(Vs*cos(ang)*Xd)2+2*Xd*Xe*Eq*Vs*cos(ang)/Xds2+(Vs*sin(ang)*Xq/Xqs)2)(-0.5);RV=RE/dVtEq;SV=SE-RV*dVtang;w0=2*pi*f0;%計(jì)算得到A,B,K矩陣Q=1 0 0;0 100 0;0 0 5000;R=1;A=(SE-SV)/(Td*SV) SEz -RV*SE/(Td*SV);-w0/

15、H -D/H 0;(SE-SV)/(Td*RV*SV) (SEz-SV)/RV -SE/(Td*SV)B=REz/Td0;0;REz/(Td0*RV)K=lqr(A,B,Q,R)程序2：%故障前或擾動前的系統(tǒng)狀態(tài)function dy=gz1(t,y);Xd=2.5;Xdz=0.5XT=0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;Eq0=1.9784;Pe=y(3)*Vs*sin(y(1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1);Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds

16、)/(Xdzs*Xds)*cos(y(1);Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1);Td0=10;Td=Td0*Xdzs/Xds;dy(1)=y(2)-100*pi;dy(2)=2*pi*50*0.55/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td0*(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy(2);d

17、y(3);程序3：%計(jì)算故障前或擾動前的電磁功率function Pe1=dcgl1(y)Pe1=zeros(length(y(:,1),1);Xd=2.5;Xdz=0.5;XT=0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;for i=1:length(y(:,1);Pe1(i)=y(i,3)*Vs*sin(y(i,1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1);End程序4：%某線路首端三相短路故障中的系統(tǒng)狀態(tài)function dy=gz2(t,y);Xd=2.5;Xdz=0.5XT=

18、0.01;Xe=XT;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;Pe=0;Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds)/(Xdzs*Xds)*cos(y(1);Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1);Eq0=1.9784;Td0=10;Td=Td0*Xdzs/Xds;dy(1)=y(2)-100*pi;dy(2)=2*pi*50*0.55/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-6

19、8.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td0*(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy(2);dy(3);程序5：%計(jì)算故障中的電磁功率function Pe4=dcgl4(y)Pe4=zeros(length(y(:,1),1);for i=1:length(y(:,1);Pe4(i)=0;End程序6：%切除故障線路后的系統(tǒng)狀態(tài)function dy=gz3(t,y);Xd=2.5;Xdz=0.5XT=0.01;XL=1.5;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;Pe=y(3)*Vs*

20、sin(y(1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1);Eq=Pe*Xds/(Vs*sin(y(1);Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1);Eq0=1.9784;Td0=10;Td=Td0*Xdzs/Xds;dy(1)=y(2)-100*pi;dy(2)=2*pi*50*0.55/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487

21、)+Eq0)/Td0-y(3)/Td+1/Td0*(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy(2);dy(3);程序7：%計(jì)算故障后的電磁功率function Pe5=dcgl5(y)Pe5=zeros(length(y(:,1),1);Xd=2.5;Xdz=0.5;XT=0.01;XL=1.5;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;for i=1:length(y(:,1);Pe5(i)=y(i,3)*Vs*sin(y(i,1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1);e

22、nd程序8：%發(fā)生擾動的系統(tǒng)狀態(tài)function dy=rd(t,y);Xd=2.5;Xdz=0.5;XT=0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;Eq0=1.9784;Pe=y(3)*Vs*sin(y(1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1);Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds)/(Xdzs*Xds)*cos(y(1);Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1);Td0=10;Td=Td0*X

23、dzs/Xds;dy(1)=y(2)-100*pi;dy(2)=2*pi*50*0.6/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;dy(3)=(-46.537*(Pe-0.6)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td0*(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy(2);dy(3);程序9：%計(jì)算擾動中的電磁功率function Pe2=dcgl2(y)Pe2=zeros(length(y(:,1),1);Xd=2.5;Xdz=0.5;XT=0.01;X

24、L=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;for i=1:length(y(:,1);Pe2(i)=y(i,3)*Vs*sin(y(i,1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1);end程序10：%擾動后的系統(tǒng)狀態(tài)function dy=rdh(t,y);Xd=2.5;Xdz=0.5;XT=0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;Eq0=1.9784;Pe=y(3)*Vs*sin(y(1)/Xdzs+0.5*Vs2*(Xdzs-Xds)

25、/(Xdzs*Xds)*sin(2*y(1);Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds)/(Xdzs*Xds)*cos(y(1);Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1);Td0=10;Td=Td0*Xdzs/Xds;dy(1)=y(2)-100*pi;dy(2)=2*pi*50*0.55/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3

26、)/Td+1/Td0*(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy(2);dy(3);程序11：%計(jì)算擾動后的電磁功率function Pe3=dcgl3(y)Pe3=zeros(length(y(:,1),1);Xd=2.5;Xdz=0.5;XT=0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;for i=1:length(y(:,1);Pe3(i)=y(i,3)*Vs*sin(y(i,1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1);end程序12：%設(shè)定仿真

27、時(shí)間為40秒,5秒時(shí)電磁功率受到擾動,6秒時(shí)刻擾動結(jié)束t0=0;t1=5;t2=6;t3=40;%擾動前，以選定的最優(yōu)運(yùn)行點(diǎn)為初值T1,Y1=ode45(gz1,t0 t1,65*pi/180 100*pi 1.0487);Pe1=dcgl1(Y1);ac=Y1(length(Y1),1);wc=Y1(length(Y1),2);ec=Y1(length(Y1),3);%擾動中，以擾動前狀態(tài)量的最后一組解為擾動中的初值T2,Y2=ode45(rd,t1 t2,ac wc ec);Pe2=dcgl2(Y2);ac2=Y2(length(Y2),1);wc2=Y2(length(Y2),2);ec2=Y2(length(Y2)