外彈道設(shè)計(jì)課件

1、南京理工大學(xué)外彈道設(shè)計(jì)理論與方法課程作業(yè)學(xué)生姓名：學(xué)號(hào)：學(xué)院：課程：外彈道設(shè)計(jì)理論與方法2015年1月作業(yè)：完成小口徑火炮用榴彈的外彈道優(yōu)化設(shè)計(jì)數(shù)學(xué)模型的程序設(shè)Sdds計(jì)。要求：模型：minf(X)=mint=minf，0vXWE41R,X:(Xi，X2，X3，X4)分別代表m,ln，lb，lc無量綱化后的變量；設(shè)計(jì)變量取m,ln,lb,lc；(3)約束條件R域?yàn)椋篻1(X)=m-m0_0g4(X)=ln-0g2(X)=l/d-1.6-0g5(X)=lb.0g3(X)=5.5(nbc).02_12一一(4)已知條件：g=9.8m/s2,E0=-m0V2=202500J,m0=0.5kg,0=2

2、000m,S=1二d2,Cx-i*Cxon,i=2.9-1.373H0.32H2-0.0267H3,4H=nb-0.3。備注：m用0.4kg無量綱化，ln,lb,lc用d=0.035米無量綱化。初始點(diǎn)m0=0.74kg5=3.0，b0=0.65%=1.7一、目標(biāo)函數(shù)“，、，、.Sddsminf(X)=mint=min0vXWE41R,X:(Xi,X2,X3,X4)分別代表m,ln,lb,lc無量綱化后的變量；二、優(yōu)化設(shè)計(jì)變量選取彈重m彈頭部長(zhǎng)ln、圓柱部長(zhǎng)lc、彈尾部長(zhǎng)lb作為設(shè)計(jì)變量。三、約束關(guān)系約束條件R域?yàn)椋篻1(X)=m-m；-0g4(X)=ln一0g2(X)=lc/d-1.6.0g5

3、(X)=lb0g3(X)=5.5-(：bc)-0四、優(yōu)化設(shè)計(jì)過程4.1 建立優(yōu)化模型要求：1)彈丸飛行時(shí)間少2)彈丸氣動(dòng)外形好(阻力小)；3)彈丸飛行穩(wěn)定；4)初速與彈重組合較佳；5)達(dá)到目標(biāo)時(shí)的存速較大?？梢越⑷缦聝?yōu)化模型，并化為標(biāo)準(zhǔn)形式(1)目標(biāo)函數(shù)：Sddsminf(X)=mint=min0vIm一-0.5(1)-九c<-1.6s.t.，n.，b.，c£5.5.-n0；-b<0優(yōu)化模型中的輸入?yún)⒘繛閄=m,h,九b$c,輸出的優(yōu)化目標(biāo)為飛行時(shí)間最小。優(yōu)化模型可以整理為圖1所示的流程。圖1優(yōu)化流程圖輸入初始參數(shù)后，優(yōu)化程序調(diào)用優(yōu)化算法和外彈道子程序，并進(jìn)行循環(huán)計(jì)算判

4、定優(yōu)化結(jié)果是否達(dá)到所需的精度。優(yōu)化程序的兩個(gè)子模塊是優(yōu)化過程的關(guān)鍵：優(yōu)化算法和外彈道解算。前者決定著初始點(diǎn)按照何種規(guī)則向最優(yōu)解逼近，后者用于計(jì)算特定輸入?yún)?shù)對(duì)應(yīng)的彈道飛行時(shí)間。4.2 外彈道解算模塊本文模型針對(duì)的是小口徑榴彈的飛行時(shí)間解算。一方面外彈道優(yōu)化要有很快的計(jì)算速度；另一方面彈丸在良好的飛行狀態(tài)下，質(zhì)點(diǎn)彈道基本上反映了彈丸的實(shí)際飛行情況，因而在外彈道優(yōu)化設(shè)計(jì)中采用質(zhì)點(diǎn)彈道即可。外彈道模塊的解算過程為：外彈道函數(shù)的輸入為優(yōu)化函數(shù)的自由變量矩陣X0,分別對(duì)應(yīng)著彈丸質(zhì)量、彈頭部長(zhǎng)度、彈底部長(zhǎng)度和彈尾部長(zhǎng)度，并將4個(gè)自由變量分別進(jìn)行了無量綱化處理，輸出為彈丸在空中的飛行時(shí)間to4.3 約束坐標(biāo)

5、輪換法模塊o圖2坐標(biāo)輪換法示意圖外彈道優(yōu)化設(shè)計(jì)問題一般歸結(jié)為有約束非線性規(guī)劃問題，有著自身的復(fù)雜性。對(duì)于有約束的非線性規(guī)劃問題，在最優(yōu)設(shè)計(jì)中有一個(gè)基于懲罰函數(shù)和障礙函數(shù)的序列無約束最小化方法，把求解一個(gè)有約束問題轉(zhuǎn)化為求解無序約束的最優(yōu)化問題。由于懲罰函數(shù)具有簡(jiǎn)單、易行等特點(diǎn)，在外彈道優(yōu)化設(shè)計(jì)中，研究人員常采用懲罰函數(shù)法將有約束最優(yōu)化問題轉(zhuǎn)化為無約束最優(yōu)化問題，然后采用優(yōu)化理論中的直接方法，如模式搜索法、powell法等求解。本文模型的函數(shù)為質(zhì)點(diǎn)外彈道方程組，函數(shù)形式較為復(fù)雜，故采用簡(jiǎn)單有效的坐標(biāo)輪換法。坐標(biāo)輪換法的基本思想如圖2所示，即，對(duì)于多變量輸入X=m,£n，*%,可以先選擇

6、一維方向（將其他變量視為常值），按照一定的搜索方向從初始點(diǎn)開始搜索對(duì)應(yīng)于目標(biāo)的最優(yōu)點(diǎn)。然后，用得到一維方向最優(yōu)解去替代變量X中的對(duì)應(yīng)值，再選擇第二個(gè)方向進(jìn)行搜索最優(yōu)解，依次類推。在搜索完所有維度后對(duì)搜索后的最優(yōu)解和搜索前的初值經(jīng)行比較，如果滿足收斂的精度要求則停止搜索，如不滿足則進(jìn)行下一輪的搜索。五、優(yōu)化模型的程序?qū)崿F(xiàn)選定好優(yōu)化方法后，采用matlab軟件進(jìn)行相應(yīng)的程序編制，以完成所需達(dá)到的計(jì)算要求。5.1程序的圖形用戶界面設(shè)計(jì)本文設(shè)計(jì)的界面如圖所示:舊untitled2外彈道優(yōu)化設(shè)計(jì)I-目標(biāo)詼政小U絲去愫H牽強(qiáng)的翅員-W史比設(shè)廿數(shù)字慢生二minr(x>n«it六與8t過餐后)

7、檢入?yún)⒒鵛-m,An.M>,Ac至關(guān)白r>小彈丸用g)0.S先也d-irrirnJ1351五度第呈小卷度AuCd>1.6和螭動(dòng)雁EU324口全酶最大粒度Lfd)5-.S飛行0蕓Sdcm)2DGD建頭紀(jì)5*嗎5CAU西大于零一步擊用工生內(nèi)嵬條正蚓丸要呈md.KgJ0.74三頭就投量MiIQ鼻M就檢度Ai>(d)也能Iri3CK¥Ac(d>11.7飛疔時(shí)同3.470S5.2程序的運(yùn)行打開程序后，先核對(duì)程序的默認(rèn)參數(shù)，然后點(diǎn)擊“開始優(yōu)化”按鈕即可進(jìn)入外彈道的優(yōu)化過程，如上圖所示。六、優(yōu)化結(jié)果運(yùn)行編制好程序，經(jīng)過一定時(shí)間的計(jì)算得到如下結(jié)果:附錄1.外彈道計(jì)算程序

8、(fwddt.m)functiontime=fwddt(var)m=var(1);lamnaN=var(2);lamnaB=var(3);lamnaC=var(4);d=0.035;X=0,0;Sd=2000;Energy=202500;theata=65;g=9.8;arc=0;h=0.02;n=0;iteraTime=0;tao0=288.9;rou0=1.2063;S=pi*dA2/4;V=sqrt(2*Energy/m);x=X(:,1);y=X(:,2);theata=deg2rad(theata);Cx0n=1.15000.38301.25000.38201.35000.37501

9、.45000.36601.55000.35601.65000.34601.75000.33701.85000.32801.95000.32002.05000.31302.50000.2880;while(arc<=Sd)temp=tao(y);rou=rou0*pi_y(y)*tao0/temp;Ma=V/(sqrt(temp)*20.047);H=lamnaN+lamnaB-0.3;i=2.9-1.373*H+0.32*W2-0.0267*W3;Cx=i*interp1(Cx0n(:,1),Cx0n(:,2),Ma);dV_ds=-0.5*rou*V*S*Cx/m-g*sin(thea

10、ta)/V;dtheata_ds=-g*cos(theata)/VA2;dx_ds=cos(theata);dy_ds=sin(theata);dt=h/V;x=x+RK(dx_ds,h);y=y+RK(dy_ds,h);V=V+RK(dV_ds,h);theata=theata+RK(dtheata_ds,h);arc=arc+h;n=n+1;iteraTime=iteraTime+dt;endtime=iteraTime;functiontemp=tao(y)temp=288.9;AA=230;BB=-6.328e-3;CC=1.172e-6;Rd=287.05;if(y<=9300

11、)temp=288.9-y*0.006328;endif(y>9300&&y<12000)temp=AA+(y-9300)*BB+pow(y-9300),2)*CC;endif(y>=12000&&y<30000)temp=221.5;endfunctionvalue=pi_y(y)value=1;Rd=287.05;if(y<=9300)value=(1-2.1904e-5*y)A5.4;endif(y>9300&&y<12000)value=0.2922575*exp(-2.1206426*(atan

12、(2.344*(y-9300)-6328)/32221.057)+0.19392520);endif(y>=12000&&y<30000)value=0.1937254*exp(-(y-12000)/6483.305);endif(y>30000)value=exp(-9.8/Rd*(y/221.5);endfunctionstep=RK(diff,h)k1=diff;k2=diff+h*k1/2;k3=diff+h*k2/2;k4=diff+h*k3;step=h*(k1+2*k2+2*k3+k4)/6;2.約束問題的坐標(biāo)輪換法函數(shù)(main.m)funct

13、ionxxx,mintime=main(epsi,x0)globalalphan=0;k=1;X(k,1,1)=0;X(k,1,2)=0;X(k,1,3)=0;X(k,1,4)=0;whilesqrt(x0-X(k,1,1)A2+(x0(2)-X(k,1,2)A2+(x0(3)-X(k,1,3)A2+(x0(4)-X(k,1,4)A2)>epsiX(k,1,1)=x0(1);X(k,1,2)=x0(2);X(k,1,3)=x0(3);X(k,1,4)=x0(4);alpha=0;xx,min,alpha=oneMIN(0.5,0.74,epsi,1,x0);X(k,2,1)=alpha;

14、X(k,2,2)=X(k,1,2);X(k,2,3)=X(k,1,3);X(k,2,4)=X(k,1,4);x0=X(k,2,1),X(k,2,2),X(k,2,3),X(k,2,4);alpha=0;xx,min,alpha=oneMIN(2.85,3.0,epsi,2,x0);X(k,3,1)=X(k,2,1);X(k,3,2)=alpha;X(k,3,3)=X(k,2,3);X(k,3,4)=X(k,2,4);x0=X(k,3,1),X(k,3,2),X(k,3,3),X(k,3,4);alpha=0;xx,min,alpha=oneMIN(0.5,0.8,epsi,3,x0);X(k,

15、4,1)=X(k,3,1);X(k,4,2)=X(k,3,2);X(k,4,3)=alpha;X(k,4,4)=X(k,3,4);x0=X(k,4,1),X(k,4,2),X(k,4,3),X(k,4,4);alpha=0;xx,min,alpha=oneMIN(1.6,1.7,epsi,4,x0);lamnaC>=1.6X(k,5,1)=X(k,4,1);X(k,5,2)=X(k,4,2);X(k,5,3)=X(k,4,3);X(k,5,4)=alpha;x0=X(k,5,1),X(k,5,2),X(k,5,3),X(k,5,4);%k=k+1n=n+1endxxx=x0;mintim

16、e=min;functionxxx,mintime,alpha=oneMIN(a,b,epsi,i,x0)%globalalphaifi=1ii(1,1)=0;ii(1,2)=1;ii(2,1)=1;ii(2,2)=0;ii(3,1)=1;ii(3,2)=0;ii(4,1)=1;ii(4,2)=0;elseifi=2ii(1,1)=1;ii(1,2)=0;ii(2,1)=0;ii(2,2)=1;ii(3,1)=1;ii(3,2)=0;ii(4,1)=1;ii(4,2)=0;elseifi=3ii(1,1)=1;ii(1,2)=0;ii(2,1)=1;ii(2,2)=0;ii(3,1)=0;ii

17、(3,2)=1;ii(4,1)=1;ii(4,2)=0;elseifi=4ii(1,1)=1;ii(1,2)=0;ii(2,1)=1;ii(2,2)=0;ii(3,1)=1;ii(3,2)=0;ii(4,1)=0;ii(4,2)=1;enda1=b-0.618*(b-a);xx1=ii(1,1)*x0+ii(1,2)*a1,ii(2,1)*x0(2)+ii(2,2)*a1,ii(3,1)*x0(3)+ii(3,2)*a1,ii(4,1)*x0(4)+ii(4,2)*a1;f1=fwddt(xx1);a2=a+0.618*(b-a);xx2=ii(1,1)*x0+ii(1,2)*a2,ii(2,1)*x0(2)+ii(2,2)*a2,ii(3,1)*x0(3)+ii(3,2)*a2,ii(4,1)*x0(4)+ii(4,2)*a2;f2=fwddt(xx2);whileabs(b-a)>epsiif(f1<=f2)b=a2;a2=a1;f2=f1;a1=b-0.618*(b-a);xx1=ii(1,1)*x0+ii(1,2)*a1,ii(2,1)*x0+ii(2,2)* a1,ii(3,1)*x0(3)+ii(3,2)*a1,ii(4,1)*x0(4)+ii(4,2)*a1;f1=fwddt(xx1);elsea=a1;a1=a2;f1=f2;a2