彈箭飛行控制原理試驗(yàn)報(bào)告-自動(dòng)駕駛儀程序報(bào)告

《彈箭飛行控制原理試驗(yàn)報(bào)告-自動(dòng)駕駛儀程序報(bào)告》由會(huì)員分享，可在線閱讀，更多相關(guān)《彈箭飛行控制原理試驗(yàn)報(bào)告-自動(dòng)駕駛儀程序報(bào)告（21頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

1、《彈箭飛行控制原理》 實(shí)驗(yàn)報(bào)告 姓 名: 學(xué) 號(hào): 專 業(yè)：武器系統(tǒng)與發(fā)射工程 一、實(shí)驗(yàn)?zāi)康? 1. 深入理解彈箭飛行控制自動(dòng)駕駛儀的功能； 2. 掌握彈箭線性化傳遞函數(shù)的一般計(jì)算方法； 3. 掌握各種典型自動(dòng)駕駛儀結(jié)構(gòu)，具有一定地自動(dòng)駕駛儀參數(shù)設(shè)計(jì)與分析能力； 4. 增強(qiáng)對(duì)制導(dǎo)彈箭飛行控制系統(tǒng)的認(rèn)識(shí)、理解、工程設(shè)計(jì)能力，鍛煉通過(guò)實(shí)驗(yàn)驗(yàn)證分析控制系

2、統(tǒng)動(dòng)態(tài)品質(zhì)的動(dòng)手能力。 二、實(shí)驗(yàn)內(nèi)容和項(xiàng)目 1. 實(shí)驗(yàn)內(nèi)容 1）制導(dǎo)彈箭自動(dòng)駕駛儀功能及工作過(guò)程演示實(shí)驗(yàn) 結(jié)合課堂教學(xué)內(nèi)容和彈箭自動(dòng)駕駛儀模擬系統(tǒng)，觀摩自動(dòng)駕駛儀的功能和工作過(guò)程，理解其在實(shí)際飛行控制系統(tǒng)中的作用。 2）彈箭線性化傳遞函數(shù)的計(jì)算實(shí)驗(yàn) 根據(jù)課堂教學(xué)要求和制導(dǎo)彈箭飛行彈道方案，計(jì)算各個(gè)特征點(diǎn)上的動(dòng)力學(xué)系數(shù)，進(jìn)而獲得彈體傳遞函數(shù)系數(shù)。在進(jìn)行實(shí)驗(yàn)課時(shí)，學(xué)生可以將自己的計(jì)算結(jié)果與軟件中內(nèi)置的彈體傳函參數(shù)進(jìn)行比較，檢查自己的計(jì)算結(jié)果。 3) 彈箭典型控制系統(tǒng)的設(shè)計(jì)及性能分析實(shí)驗(yàn) 根據(jù)課堂教學(xué)和彈體傳遞函數(shù)計(jì)算結(jié)果，在實(shí)驗(yàn)課中，通過(guò)操作“彈箭飛行控制與軟件設(shè)計(jì)系統(tǒng)”，在自動(dòng)駕

3、駛儀結(jié)構(gòu)圖上，按照從內(nèi)環(huán)到外環(huán)的順序，逐步進(jìn)行時(shí)域分析和頻域分析，了解各個(gè)參數(shù)的影響和作用，并確定控制系統(tǒng)的控制參數(shù)。然后進(jìn)行六自由度彈道仿真驗(yàn)證控制參數(shù)的正確性。 2. 實(shí)驗(yàn)要求 1）獨(dú)立完成實(shí)驗(yàn)，撰寫試驗(yàn)報(bào)告； 2）提交實(shí)驗(yàn)相關(guān)程序。 三、實(shí)驗(yàn)人員和實(shí)驗(yàn)時(shí)間 設(shè)計(jì)人員： 夏斌 設(shè)計(jì)時(shí)間： 6.20 四、實(shí)驗(yàn)設(shè)備 彈箭飛行控制設(shè)計(jì)系統(tǒng) 五、實(shí)驗(yàn)步驟及過(guò)程 1. 根據(jù)初始彈道條件和彈道方案，計(jì)算方案彈道和導(dǎo)引彈道； 2. 分別根據(jù)設(shè)計(jì)要求，選擇下滑段、等高巡航段、末制導(dǎo)段典型特征點(diǎn)，計(jì) 算動(dòng)力學(xué)系數(shù)和傳遞函數(shù)； 3. 設(shè)計(jì)俯仰、滾轉(zhuǎn)和偏航通道的控制器參數(shù)，并

4、進(jìn)行時(shí)域和頻域的動(dòng)態(tài)品質(zhì) 分析； 4. 進(jìn)行六自由度彈道仿真驗(yàn)證，并繪制彈道參數(shù)曲線。六自由度彈道輸出參、數(shù)為：時(shí)間、x、y、z、V、彈道傾角θ、彈道偏角、質(zhì)量m、馬赫數(shù)Ma、俯仰角?、偏航角ψ、滾轉(zhuǎn)角γ、攻角α、側(cè)滑角β、速度傾斜角，俯仰指令（分段改變），偏航指令（分段改變），階段標(biāo)志（1、初始階段；2、中制導(dǎo)；3、末制導(dǎo)）； 5. 撰寫實(shí)驗(yàn)報(bào)告。 六、 實(shí)驗(yàn)結(jié)果與分析 1、 建立彈道模型 建立彈道模型，起始點(diǎn)坐標(biāo)為（x0,y0）=（0,7100m），初速V0=260m/s，目標(biāo)位置(xt,yt)=（12000m，0）。 1）初始發(fā)射段方程組： 初始發(fā)射段采用俯仰角控制

5、方案，方案俯仰角 2）巡航段方程組： 巡航段采用高度控制方案，方案高度，俯仰角 3）末端攻擊段方程組： 末端攻擊段采用比例導(dǎo)引，導(dǎo)引律為 2、 特征點(diǎn)的動(dòng)力學(xué)系數(shù)和傳遞函數(shù) 按正常式布局導(dǎo)彈， 1） 列寫a22等動(dòng)力學(xué)系數(shù)表達(dá)式 2） 列寫： 3） 寫出傳遞函數(shù)表達(dá)式 俯仰通道傳遞函數(shù)為 其中：a1=-a22;a2=-a24;a3=-a25;a4=a34;a5=a35; 參數(shù) 特征點(diǎn) a1 a2 a3 a4 a5 第一轉(zhuǎn)折點(diǎn) -0

6、.0051 -699.4136 -540.9301 2.4876 0.4131 -1.5108 0.0378 0.0471 0.5119 第二轉(zhuǎn)折點(diǎn) -0.0049 -805.5900 -615.3198 2.8095 0.4722 -1.6737 0.0352 0.0496 0.4563 動(dòng)壓最大點(diǎn) -0.0113 -1547.6 -1220.4 6.2421 -0.2262 -3.8355 0.0254 0.0795 0.2056 動(dòng)壓最小點(diǎn) -0.0047 -589.6635 -469.9341 2.2679 0.398

7、1 -1.4987 0.0380 0.0470 0.5167 攻角最大點(diǎn) -0.0051 -699.4138 -540.9297 2.4876 0.4131 -1.4987 0.0380 0.0470 0.5167 過(guò)載最大點(diǎn) -0.0051 -699.4136 -540.9301 2.4876 0.4131 -1.5108 0.0378 0.0471 0.5119 3、 俯仰通道控制器參數(shù)設(shè)計(jì)并進(jìn)行時(shí)域和頻域的動(dòng)態(tài)品質(zhì)分析 姿態(tài)控制回路 選取需用過(guò)載最大點(diǎn)作為特征點(diǎn)進(jìn)行參數(shù)設(shè)計(jì)。 參數(shù)a1 參數(shù)a2 參數(shù)a3 參數(shù)a4

8、舵機(jī)比例參數(shù) 舵機(jī)時(shí)間參數(shù) -0.0051 -699.4136 -540.9301 2.4876 1 1 在設(shè)計(jì)中只是從整個(gè)的傳遞函數(shù)著手，系統(tǒng)的開(kāi)環(huán)傳遞函數(shù): 設(shè)計(jì)結(jié)果： 內(nèi)回路 Kwz 頻域設(shè)計(jì) 時(shí)域設(shè)計(jì) 0.080 幅值裕度(db) 相角裕度(°) 穿越頻率(Hz) 超調(diào)量(%) 上升時(shí)間 (s) 調(diào)節(jié)時(shí)間 (s) Inf 91.3 55.8 463 0.005 0.240 外回路 控制參數(shù) 頻域設(shè)計(jì) 時(shí)域設(shè)計(jì) Kp Ki 賦值裕度(db) 相角裕度(°) 穿越頻率(Hz) 超調(diào)量(%) 上升時(shí)

9、間(s) 調(diào)節(jié)時(shí)間(s) 25 180 Inf 78.6 24.5 9 0.250 1.800 設(shè)計(jì)結(jié)果分析： 1）俯仰回路姿態(tài)控制內(nèi)回路頻域響應(yīng)Bode圖 2）俯仰回路姿態(tài)控制內(nèi)回路時(shí)域階躍響應(yīng)圖 3）俯仰回路姿態(tài)控制外回路頻域響應(yīng)Bode圖 4）俯仰回路姿態(tài)控制外回路時(shí)域階躍響應(yīng)圖 4、 進(jìn)行六自由度彈道仿真驗(yàn)證，并繪制彈道參數(shù)曲線 彈道曲線： 設(shè)計(jì)的理論參數(shù)曲線： （1）y-x曲線 （2）、y-t曲線 （3）、攻角Apa—t曲線 （4）、V—t曲線 （5）、彈道傾角θ-t曲線 （6）、俯仰舵偏角δz-t

10、曲線 七、 實(shí)驗(yàn)結(jié)論 1、各回路的作用： 內(nèi)環(huán)回路為阻尼回路，作用是提高增穩(wěn)回路的穩(wěn)定性，保證整個(gè)飛行過(guò)程中有較好的動(dòng)態(tài)性能，使彈體更好的跟隨動(dòng)態(tài)制導(dǎo)指令的變化。 外環(huán)回路是穩(wěn)增回路，作用是增加回路的穩(wěn)定性，改善導(dǎo)彈的動(dòng)態(tài)性能。 2、姿態(tài)控制回路參數(shù)設(shè)計(jì)步驟： 1）計(jì)算： 內(nèi)回路的傳遞函數(shù)為： W=KWZKM(T1S+1)TM2S2+2TM?M+KWZKMT1S+KWZKM+1 根據(jù)期望的阻尼系數(shù)為0.707來(lái)計(jì)算相應(yīng)的控制參數(shù)KWZ，表達(dá)式如下： KWZ=(-BK+AK22+AK2-BKAK2+AK44)/KM 其中的AK=2TM?M/T1,BK=2TM?M/

11、T1 通過(guò)公式求解出KWZ為-0.080 2）求 在設(shè)計(jì)時(shí)，應(yīng)將所在的支路斷開(kāi)，此時(shí)的開(kāi)環(huán)傳遞函數(shù)為 W=KPKMKWZ(T1S+1)S(TM2S2+2TM?M+KWZKMT1S+KWZKM+1) 所在的回路為增穩(wěn)回路，增穩(wěn)回路的帶寬WH與導(dǎo)彈固有頻率為：WH≥K1WN,將W換算成頻域形式，令|W|=1，計(jì)算出為25。 3）求Ki Ki所在的回路為制導(dǎo)回路，此時(shí)系統(tǒng)的開(kāi)環(huán)傳遞函數(shù)為 W=KiKPKMKWZ(T1S+1)S2(TM2S2+2TM?M+KWZKMT1S+KWZKM+1) 制導(dǎo)回路帶寬為K1Wn≥K2WT，，將W轉(zhuǎn)換成頻域形式，此時(shí)令|W|=1，計(jì)算出Ki為180.

12、 九、問(wèn)題分析及解決方法 （1）程序首先開(kāi)始是根據(jù)彈道方程，通過(guò)龍格庫(kù)塔法計(jì)算出彈道方程，上個(gè)學(xué)期使用c語(yǔ)言編寫程序，這個(gè)學(xué)期改用MATLAB編寫，在編寫的過(guò)程中在程序的調(diào)用，飛行各個(gè)過(guò)程的劃分中遇到了一些問(wèn)題，通過(guò)了向之前用過(guò)MATLAB的同學(xué)請(qǐng)教最后計(jì)算出飛行軌跡。 （2）接下來(lái)是進(jìn)行彈道特征點(diǎn)的選擇，我沒(méi)有編寫專門的程序進(jìn)行編寫，而是通過(guò)查看彈道軌跡數(shù)據(jù)手動(dòng)選擇，將選出的特征點(diǎn)的V,H,M輸進(jìn)程序計(jì)算KM T1 TM ?M。 （3）接下來(lái)是進(jìn)行KWZ KP KI的選擇。 在沒(méi)有控制的情況下，系統(tǒng)非常不穩(wěn)定 在程序剛開(kāi)始編制的時(shí)候，我先是使用了多次嘗試的方法，畫出KW

13、Z KP KI在不同情況下的各種時(shí)域響應(yīng)圖然后根據(jù)圖像選擇最后的情況，之后上課的PPT上講了詳細(xì)的KWZ KP KI的計(jì)算方法，我又重新計(jì)算了結(jié)果得到了更好的制導(dǎo)結(jié)果。 十、附件程序說(shuō)明 （1）MATLAB文件包含6個(gè)M函數(shù)，分別為主程序mainzhuhanshu.m，微分方程函數(shù)dxl，特征點(diǎn)計(jì)算好h1.m,h2.m,h3.m,控制參數(shù)計(jì)算函數(shù)r1.m. （2）程序圖像文件包含軌跡圖像和時(shí)域圖，伯德圖。 clc; clear; close all; R0=1.225*(1.0-2.0323*10^(-5)*7100)^4.83;%?????ü?èμ????? S=0.45

14、6; q0=0.5*R0*260*260; tspan=0:0.05:60; nspan=[260 0 0 7100 225 0 q0 0 0]; [t,y]=ode45('dxl',tspan,nspan); figure(1), plot(t,y(:,1),'r-','LineWidth',2),grid on;%?????ù?è??ê±??±??ˉ?ú?? title('V-T±??ˉ?ú??'), xlabel('t/s'); ylabel('v/(m/s)') figure(2), plot(y(:,3),y(:,4),'r-','LineWidth',2),g

15、rid on;%???????è???àà?±??ˉ?ú?? title('Y-X±??ˉ?ú??'), xlabel('x/m'); ylabel('y/m'); figure(3), plot(t,y(:,2),'r-','LineWidth',2),grid on;%????μˉμà??????ê±??±??ˉ?ú?? title('|è-T±??ˉ?ú??'), xlabel('t/s'); ylabel('theta/rad') figure(4), plot(t,y(:,4),'r-','LineWidth',2),grid on;%???????è??ê±??

16、±??ˉ?ú?? title('Y-T±??ˉ?ú??'), xlabel('t/s'); ylabel('y/m'); figure(5), plot(t,y(:,6),'r-','LineWidth',2),grid on;%????1￥????ê±??±??ˉ?ú?? title('|á-T±??ˉ?ú??'), xlabel('t/s'); ylabel('alpha/?￡'); figure(6), plot(t,y(:,7),'r-','LineWidth',2),grid on;%?????ˉ?1??ê±??±??ˉ?ú?? title('q-T±??ˉ?ú

17、??'), xlabel('t/s'); ylabel('q/pa'); figure(7), plot(t,y(:,8),'r-','LineWidth',2),grid on;%????Dèó?1y????ê±??±??ˉ?ú?? title('N-T±??ˉ?ú??'), xlabel('t/s'); ylabel('n/(m/s^2)'); figure(8), plot(t,y(:,9),'r-','LineWidth',2),grid on;%????????????ê±??±??ˉ?ú?? title('|?-T±??ˉ?ú??'), xlabel('t/s

18、'); ylabel('dz/?￡'); %1y??×?′óì??÷μ????? clc; clear; close all; S=0.456;%???y p=3000; m=225; L=2.358; Jz=60.85; Jy=60.85; Jx=3.45; %×?′ó1y??2?êy???? Vnmax=283.68;Hnmax=6986.8;qnmax=23026;mnmax=220.8; %1y??×?′óμ??à1?2?êy???? manmax=Vnm

19、ax/(20.046*(288.34-5.86*0.001*Hnmax)^0.5); MA=[0.6,0.8,0.9,1.0,1.2];%MA?a?ío?êy CX0=[-0.2000,-0.2210,-0.2315,-0.2420,-0.2500]; CXAA=[-0.0005,-0.0005,-0.0005,-0.0005,-0.0005]; CYA=[0.2505,0.253,0.254,0.255,0.256]; CYDZ=[0.04,0.047,0.043,0.045,0.045]; MZA=[-0.0290,-0.0298,-0.0300,-0.0302,-0.03

20、02]; MZDZ=[-0.0220,-0.0240,-0.0232,-0.0230,-0.0234]; MZWZ=[-0.002,-0.0015,-0.0015,-0.0012,-0.0016]; MXWX=[-0.0002,-0.00045,-0.0005,-0.00052,-0.00052]; MXDZ=[-0.0028,-0.0027,-0.0028,-0.0028,-0.0028]; %?¨á￠???ˉá|ó????ˉá|???μêyêy×é￡??a??oóμ?2??μ×?×?±? gnmax=9.806*(1.0-2*Hnmax/6371000);%??á|?ó?ù?è

21、μ????? R=1.225*(1.0-2.0323*10^(-5).*Hnmax)^4.83;%?????ü?èμ????? cx0=spline(MA,CX0,manmax); cxaa=spline(MA,CXAA,manmax); cya=spline(MA,CYA,manmax); cydz=spline(MA,CYDZ,manmax); mza=spline(MA,MZA,manmax); mzdz=spline(MA,MZDZ,manmax);%à?ó?èy′??ùì?2??μμ?μ??à1????ˉá|?μêy mzwz=spline(MA,MZWZ,manma

22、x); mxwx=spline(MA,MXWX,manmax); mxdz=spline(MA,MXDZ,manmax); a22=mzwz*qnmax*S*L*L/(Jz*Vnmax);a24=57.3*mza*qnmax*S*L/Jz;a25=57.3*mzdz*qnmax*S*L/Jz; a34=(57.3*cya*qnmax*S+p)/(mnmax*Vnmax);a35=57.3*cydz*qnmax*S/(mnmax*Vnmax); %????μˉì?′?μY?μêy Kmm=(-a25*a34+a35*a24)/(a22*a34+a24); %Km Tmm=1/

23、(-a24-a22*a34)^0.5 ; %Tm Emm=(-a22+a34)/(-a24-a22*a34)^0.5/2; %|?m T1m=a25/(a25*a34-a35*a24); %T1 %ì??÷·?3ì?μêy % P1m=-A33+A34-A22-A11; % P2m=A31*A14-A31*A13+A22*A33-A22*A34-A24+A33*A11-A34*A11+A22*A11; % P3m=-A21*A14-A31*A22*A14+A22*A31

24、*A13+A24*A33-A22*A33*A11+A22*A34*A11+A24*A11; % P4m=A21*A33*A14-A13*A21*A34+A24*A31*A13-A24*A33*A11; % num=[Kmm*T1m];den=[Tmm^2,2*Tmm*Emm,1]; %′?μYoˉêyμ?·?×óoí·???μ??μêy %?ˉ?1×?D?ì??÷μ????? clc; clear; close all; S=0.456;%???y p=3000; m=225; L=2.358;

25、 Jz=60.85; Jy=60.85; Jx=3.45; %?ˉ?1×?D?2?êy???? Vnmax=260;Hnmax=7100;qnmax=19507;mnmax=225; %?ˉ?1×?D?μ??à1?2?êy???? manmax=Vnmax/(20.046*(288.34-5.86*0.001*Hnmax)^0.5); MA=[0.6,0.8,0.9,1.0,1.2];%MA?a?ío?êy CX0=[-0.2000,-0.2210,-0.2315,-0.2420,-0.2500]; CXAA=[-0.0005,-0.0005,-0.0

26、005,-0.0005,-0.0005]; CYA=[0.2505,0.253,0.254,0.255,0.256]; CYDZ=[0.04,0.047,0.043,0.045,0.045]; MZA=[-0.0290,-0.0298,-0.0300,-0.0302,-0.0302]; MZDZ=[-0.0220,-0.0240,-0.0232,-0.0230,-0.0234]; MZWZ=[-0.002,-0.0015,-0.0015,-0.0012,-0.0016]; MXWX=[-0.0002,-0.00045,-0.0005,-0.00052,-0.00052]; MXD

27、Z=[-0.0028,-0.0027,-0.0028,-0.0028,-0.0028]; %?¨á￠???ˉá|ó????ˉá|???μêyêy×é￡??a??oóμ?2??μ×?×?±? gnmax=9.806*(1.0-2*Hnmax/6371000);%??á|?ó?ù?èμ????? R=1.225*(1.0-2.0323*10^(-5).*Hnmax)^4.83;%?????ü?èμ????? cx0=spline(MA,CX0,manmax); cxaa=spline(MA,CXAA,manmax); cya=spline(MA,CYA,manmax); cydz=s

28、pline(MA,CYDZ,manmax); mza=spline(MA,MZA,manmax); mzdz=spline(MA,MZDZ,manmax);%à?ó?èy′??ùì?2??μμ?μ??à1????ˉá|?μêy mzwz=spline(MA,MZWZ,manmax); mxwx=spline(MA,MXWX,manmax); mxdz=spline(MA,MXDZ,manmax); a22=mzwz*qnmax*S*L*L/(Jz*Vnmax);a24=57.3*mza*qnmax*S*L/Jz;a25=57.3*mzdz*qnmax*S*L/Jz; a34=(5

29、7.3*cya*qnmax*S+p)/(mnmax*Vnmax);a35=57.3*cydz*qnmax*S/(mnmax*Vnmax); %????μˉì?′?μY?μêy Kmm=(-a25*a34+a35*a24)/(a22*a34+a24); %Km Tmm=1/(-a24-a22*a34)^0.5 ; %Tm Emm=(-a22+a34)/(-a24-a22*a34)^0.5/2; %|?m T1m=a25/(a25*a34-a35*a24); %T1

30、 %?ˉ?1×?′óì??÷μ????? clc; clear; close all; S=0.456;%???y p=3000; m=225; L=2.358; Jz=60.85; Jy=60.85; Jx=3.45; %?ˉ?1×?′ó2?êy???? Vnmax=289.39;Hnmax=0;qnmax=51114;mnmax=189.0851; %?ˉ?1×?′óμ??à1?2?êy???? manmax=Vnmax/(20.046*(288.34-5.86*0.0

31、01*Hnmax)^0.5); MA=[0.6,0.8,0.9,1.0,1.2];%MA?a?ío?êy CX0=[-0.2000,-0.2210,-0.2315,-0.2420,-0.2500]; CXAA=[-0.0005,-0.0005,-0.0005,-0.0005,-0.0005]; CYA=[0.2505,0.253,0.254,0.255,0.256]; CYDZ=[0.04,0.047,0.043,0.045,0.045]; MZA=[-0.0290,-0.0298,-0.0300,-0.0302,-0.0302]; MZDZ=[-0.0220,-0.0240

32、,-0.0232,-0.0230,-0.0234]; MZWZ=[-0.002,-0.0015,-0.0015,-0.0012,-0.0016]; MXWX=[-0.0002,-0.00045,-0.0005,-0.00052,-0.00052]; MXDZ=[-0.0028,-0.0027,-0.0028,-0.0028,-0.0028]; %?¨á￠???ˉá|ó????ˉá|???μêyêy×é￡??a??oóμ?2??μ×?×?±? gnmax=9.806*(1.0-2*Hnmax/6371000);%??á|?ó?ù?èμ????? R=1.225*(1.0-2.0323

33、*10^(-5).*Hnmax)^4.83;%?????ü?èμ????? cx0=spline(MA,CX0,manmax); cxaa=spline(MA,CXAA,manmax); cya=spline(MA,CYA,manmax); cydz=spline(MA,CYDZ,manmax); mza=spline(MA,MZA,manmax); mzdz=spline(MA,MZDZ,manmax);%à?ó?èy′??ùì?2??μμ?μ??à1????ˉá|?μêy mzwz=spline(MA,MZWZ,manmax); mxwx=spline(MA,MXWX,ma

34、nmax); mxdz=spline(MA,MXDZ,manmax); a22=mzwz*qnmax*S*L*L/(Jz*Vnmax);a24=57.3*mza*qnmax*S*L/Jz;a25=57.3*mzdz*qnmax*S*L/Jz; a34=(57.3*cya*qnmax*S+p)/(mnmax*Vnmax);a35=57.3*cydz*qnmax*S/(mnmax*Vnmax); %????μˉì?′?μY?μêy Kmm=(-a25*a34+a35*a24)/(a22*a34+a24); %Km Tmm=1/(-a24-a22*a34)^0.5 ;

35、 %Tm Emm=(-a22+a34)/(-a24-a22*a34)^0.5/2; %|?m T1m=a25/(a25*a34-a35*a24); %T1 %n×?′óμ? Km1 =-1.5108; Tm1 =0.0378; Em1 = 0.0471; T11 = 0.5119; % ×?ì????????·éè??μúò???2?êyKwz Kwz=-0.080; % ì??÷μ?1μ?ê±óò?ìó| num=[Kwz*Km1*T11 Kwz*Km1];

36、den=[Tm1^2 2*Tm1*Em1 1]; Gs1=tf(num,den); [Gm1,Pm1,Wcg1,Wcp1]=margin(Gs1); sys1=feedback(Gs1,1); t=0:0.01:6; y=step(sys1,t); figure(1),plot(t,y); grid on title('Unit-Step Response kwz=-0.080'),xlabel('Time(seconds)'); figure(2),margin(num,den),grid on % éè??2?êyKp￡?ki Kp=25;%ó°?ìμ÷??ê±?

37、? Ki=180;%ó°?ì·ù?μ?￡?èoí?à???￡?èò2?áó°?ìμ÷??ê±??￡¨?úkp2?±???￡?êêμ±??D?ki?à???￡?è??′ó·ù?μ?￡?èò2??′ó;?úki2?±???￡???′ókp?à???￡?è??D?￡?·ù?μ?￡?è??D?￡???D?kp?à???￡?è??′ó￡? % ì??÷μ?1??ó|óúkiμ??μóò?ìó|ò??°ía???·±??·ê±óò?ìó| num=[Kp*Kwz*Km1*T11 Kwz*Km1*(Kp+Ki*T11) Kwz*Km1*Ki]; den=[Tm1^2 2*Tm1*Em1+K

38、wz*Km1*T11 1+Kwz*Km1 0 0]; figure(3),margin(num,den),grid on Gs13=tf(num,den); sys13=feedback(Gs13,1); y=step(sys13,t); figure(4), plot(t,y),grid on title('Unit-Step Response Kp=25 Ki=180'),xlabel('Time(seconds)'); function dy=dxl(t,y) dy=zeros(9,1);%dy(1)?a?ù?è??ê±??μ?±??ˉ?ê

39、 %dy(2)?aμˉμà??????ê±??μ?±??ˉ?ê %dy(3)?a?????àà???ê±??μ?±??ˉ?ê %dy(4)?a???è??ê±??μ?±??ˉ?ê %dy(5)?a?êá???ê±??μ?±??ˉ?ê %dy(6)?a1￥????ê±??±??ˉ?ê mc=0.7;p=3000;%??á÷á?ò??°·￠?ˉ?úí?á| MA=[0.6,0.8,0.9,1.0,1.2];%MA?a?ío?êy CX0=[-0.2000,-0.2

40、210,-0.2315,-0.2420,-0.2500]; CXAA=[-0.0005,-0.0005,-0.0005,-0.0005,-0.0005]; CYA=[0.2505,0.253,0.254,0.255,0.256]; CYDZ=[0.04,0.047,0.043,0.045,0.045]; MZA=[-0.0290,-0.0298,-0.0300,-0.0302,-0.0302]; MZDZ=[-0.0220,-0.0240,-0.0232,-0.0230,-0.0234]; %?¨á￠???ˉá|ó????ˉá|???μêyêy×é￡??a??oóμ?2??μ×?×

41、?±? C=20.046*(288.34-5.86*10^(-3).*y(4))^0.5;%ò??ùμ????? ma=y(1)./C; cx0=spline(MA,CX0,ma); cxaa=spline(MA,CXAA,ma); cya=spline(MA,CYA,ma); cydz=spline(MA,CYDZ,ma); mza=spline(MA,MZA,ma); mzdz=spline(MA,MZDZ,ma);%à?ó?èy′??ùì?2??μμ?μ??à1????ˉá|?μêy S=0.456; R=1.225*(1.0-2.0323*10^(-5)

42、*y(4))^4.83;%?????ü?èμ????? g=9.806*(1.0-2.*y(4)/6371000);%??á|?ó?ù?èμ????? YA=0.5*R.*y(1).^2*S.*cya;%éyá|1?óú1￥??μ?μ?êy cx=cx0+cxaa.*y(6).*y(6); dz=-mza/mzdz.*y(6); cy=cya.*y(6)+cydz.*dz; r=(y(3)-12000).^2+y(4).^2;%μˉ???àà? Dq=(y(1).*sin(y(2))*(y(3)-12000)-y(1).*cos(y(2)).*y(4))./r; X=0.

43、5.*R.*y(1).^2*S.*cx;%×èá|???? Y=0.5.*R.*y(1).^2*S.*cy;%éyá|???? %q=0.5*R.*y(1).*y(1); %?ˉ?1???? if y(4)>0 dy(1)=(p*cos(y(6)/57.3)+X-y(5)*g.*sin(y(2)))./y(5); dy(2)=(p*sin(y(6)/57.3)+Y-y(5)*g.*cos(y(2)))./(y(5).*y(1)).*(t<=6)+(0-y(2))*(t>6&t<20)+4*Dq.*(t>=20); dy(3)=y(1).*cos(y(2)); dy(4)=y(

44、1).*sin(y(2))*(t<=6)+(7000-y(4)).*(t>6&t<20)+y(1).*sin(y(2))*(t>=20); dy(5)=-mc; dy(6)=(5*sin(0.2*pi*t+pi)-y(2)*57.3).*(t<=6)+((y(5)*g)/(p/57.3+YA))*(t>6&t<20)+((y(5).*y(1).*4*Dq+y(5).*g.*cos(y(2)))/(p/57.3+YA))*(t>=20)-y(6); dy(7)=0.5*R*y(1)*y(1)-y(7); dy(8)=y(1)*dy(2)/g+cos(y(2))-y(8); dy(9)=0; end if y(4)<0 return end end