外彈道設(shè)計文檔

南京理工大學(xué)外彈道設(shè)計理論與方法課程作業(yè)學(xué)生姓名：學(xué)號：學(xué)院：課程：外彈道設(shè)計理論與方法

2015

年

1月

作業(yè)：完成小口徑火炮用榴彈的外彈道優(yōu)化設(shè)計數(shù)學(xué)模型的程序設(shè)計。要求：(1)模型：，，分別代表無量綱化后的變量；(2)設(shè)計變量取m，ln，lb，lc；(3)約束條件R域為：(4)已知條件：，，，m，，，，。備注：用0.4kg無量綱化，用d＝0.035米無量綱化。初始點

1、目標(biāo)函數(shù)，，分別代表無量綱化后的變量；二、優(yōu)化設(shè)計變量選取彈重m、彈頭部長ln、圓柱部長lc、彈尾部長lb作為設(shè)計變量。三、約束關(guān)系約束條件R域為：4、優(yōu)化設(shè)計過程4.1建立優(yōu)化模型要求：1)彈丸飛行時間少2)彈丸氣動外形好（阻力?。?)彈丸飛行穩(wěn)定；4)初速與彈重組合較佳；5)達到目標(biāo)時的存速較大?？梢越⑷缦聝?yōu)化模型，并化為標(biāo)準(zhǔn)形式（1）。目標(biāo)函數(shù):（1）優(yōu)化模型中的輸入?yún)⒘繛?，輸出的?yōu)化目標(biāo)為飛行時間最小。優(yōu)化模型可以整理為圖1所示的流程。圖1優(yōu)化流程圖輸入初始參數(shù)后，優(yōu)化程序調(diào)用優(yōu)化算法和外彈道子程序，并進行循環(huán)計算判定優(yōu)化結(jié)果是否達到所需的精度。優(yōu)化程序的兩個子模塊是優(yōu)化過程的關(guān)鍵：優(yōu)化算法和外彈道解算。前者決定著初始點按照何種規(guī)則向最優(yōu)解逼近，后者用于計算特定輸入?yún)?shù)對應(yīng)的彈道飛行時間。4.2外彈道解算模塊本文模型針對的是小口徑榴彈的飛行時間解算。一方面外彈道優(yōu)化要有很快的計算速度；另一方面彈丸在良好的飛行狀態(tài)下，質(zhì)點彈道基本上反映了彈丸的實際飛行情況，因而在外彈道優(yōu)化設(shè)計中采用質(zhì)點彈道即可。外彈道模塊的解算過程為：外彈道函數(shù)的輸入為優(yōu)化函數(shù)的自由變量矩陣X0，分別對應(yīng)著彈丸質(zhì)量、彈頭部長度、彈底部長度和彈尾部長度，并將4個自由變量分別進行了無量綱化處理，輸出為彈丸在空中的飛行時間t。4.3約束坐標(biāo)輪換法模塊圖2坐標(biāo)輪換法示意圖外彈道優(yōu)化設(shè)計問題一般歸結(jié)為有約束非線性規(guī)劃問題，有著自身的復(fù)雜性。對于有約束的非線性規(guī)劃問題，在最優(yōu)設(shè)計中有一個基于懲罰函數(shù)和障礙函數(shù)的序列無約束最小化方法，把求解一個有約束問題轉(zhuǎn)化為求解無序約束的最優(yōu)化問題。由于懲罰函數(shù)具有簡單、易行等特點，在外彈道優(yōu)化設(shè)計中，研究人員常采用懲罰函數(shù)法將有約束最優(yōu)化問題轉(zhuǎn)化為無約束最優(yōu)化問題，然后采用優(yōu)化理論中的直接方法，如模式搜索法、powell法等求解。本文模型的函數(shù)為質(zhì)點外彈道方程組，函數(shù)形式較為復(fù)雜，故采用簡單有效的坐標(biāo)輪換法。坐標(biāo)輪換法的基本思想如圖2所示，即，對于多變量輸入，可以先選擇一維方向（將其他變量視為常值），按照一定的搜索方向從初始點開始搜索對應(yīng)于目標(biāo)的最優(yōu)點。然后，用得到一維方向最優(yōu)解去替代變量X中的對應(yīng)值，再選擇第二個方向進行搜索最優(yōu)解，依次類推。在搜索完所有維度后對搜索后的最優(yōu)解和搜索前的初值經(jīng)行比較，如果滿足收斂的精度要求則停止搜索，如不滿足則進行下一輪的搜索。五、優(yōu)化模型的程序?qū)崿F(xiàn)選定好優(yōu)化方法后，采用matlab軟件進行相應(yīng)的程序編制，以完成所需達到的計算要求。5.1程序的圖形用戶界面設(shè)計本文設(shè)計的界面如圖所示：5.2程序的運行打開程序后，先核對程序的默認(rèn)參數(shù)，然后點擊“開始優(yōu)化”按鈕即可進入外彈道的優(yōu)化過程，如上圖所示。6、優(yōu)化結(jié)果運行編制好程序，經(jīng)過一定時間的計算得到如下結(jié)果：附錄1.外彈道計算程序(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*d^2/4;V=sqrt(2*Energy/m);x=X(:,1);y=X(:,2);theata=deg2rad(theata);Cx0n=[1.1500

0.38301.2500

0.38201.3500

0.37501.4500

0.36601.5500

0.35601.6500

0.34601.7500

0.33701.8500

0.32801.9500

0.32002.0500

0.31302.5000

0.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*H^2-0.0267*H^3;Cx=i*interp1(Cx0n(:,1),Cx0n(:,2),Ma);dV_ds=-0.5*rou*V*S*Cx

/m-g*sin(theata)/V;

dtheata_ds=-g*cos(theata)/V^2;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)

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;

end---------------------------------------------------------functionvalue=pi_y(y)value=1;Rd=287.05;if(y<=9300)

value=(1-2.1904e-5*y)^5.4;

endif(y>9300&&y<12000)

value=0.2922575*exp(-2.1206426*(atan((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));end---------------------------------------------------------functionstep=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)function[xxx,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(1)-X(k,1,1))^2+(x0(2)-X(k,1,2))^2+(x0(3)-X(k,1,3))^2+(x0(4)-X(k,1,4))^2)>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;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,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;mintime=min;---------------------------------------------------------function[xxx,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(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(1)+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(1)+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);whilea