外弹道优化设计问题一般归结为有约束非线性规划问题，有着自身的复杂性。对于有约束的非线性规划问题，在最优设计中有一个基于惩罚函数和障碍函数的序列无约束最小化方法，把求解一个有约束问题转化为求解无序约束的最优化问题。由于惩罚函数具有简单、易行等特点，在外弹道优化设计中，研究人员常采用惩罚函数法将有约束最优化问题转化为无约束最优化问题，然后采用优化理论中的直接方法，如模式搜索法、powell法等求解。
本文模型的函数为质点外弹道方程组，函数形式较为复杂，故采用简单有效的坐标轮换法。坐标轮换法的基本思想如图2所示，即，对于多变量输入X?[m,?n,?b,?c]，可以先选择一维方向（将其他变量视为常值），按照一定的搜索方向从初始点开始搜索对应于目标的最优点。然后，用得到一维方向最优解去替代变量X中的对应值，再选择第二个方向进行搜索最优解，依次类推。在搜索完所有维度后对搜索后的最优解和搜索前的初值经行比较，如果满足收敛的精度要求则停止搜索，如不满足则进行下一轮的搜索。
五、优化模型的程序实现
选定好优化方法后，采用matlab软件进行相应的程序编制，以完成所需达到的计算要求。
5.1 程序的图形用户界面设计
本文设计的界面如图所示：
5.2 程序的运行
打开程序后，先核对程序的默认参数，然后点击“开始优化”按钮即可进入外弹道的优化过程，如上图所示。
六、优化结果
运行编制好程序，经过一定时间的计算得到如下结果：
附录
1.外弹道计算程序(fwddt.m)
function time = 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.3830 1.2500 0.3820 1.3500 0.3750 1.4500 0.3660 1.5500 0.3560 1.6500 0.3460 1.7500 0.3370 1.8500 0.3280 1.9500 0.3200 2.0500 0.3130
2.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; end
time = iteraTime;
-----------------------------------------------------------------------------------------
function temp = 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; end
if (y > 9300 && y <12000)
temp = AA + (y - 9300) * BB + pow((y - 9300),2) * CC; end
if (y >= 12000 && y < 30000) temp = 221.5; end
--------------------------------------------------------- function value = pi_y(y)
value = 1; Rd = 287.05;
if (y <= 9300)
value = (1 - 2.1904e-5 * y)^5.4; end
if (y > 9300 && y <12000)
value = 0.2922575 * exp(-2.1206426 * (atan((2.344 * (y - 9300) - 6328) / 32221.057) + 0.19392520));
end
if (y >= 12000 && y < 30000)
value = 0.1937254 * exp(-(y - 12000)/6483.305); end
if (y > 30000)
value = exp(-9.8 / Rd * (y / 221.5)); end
--------------------------------------------------------- function step = 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.约束问题的坐标轮换法函数 (main.m)
function [xxx,mintime] = main(epsi,x0)
global alpha n = 0; k = 1;
X(k,1,1) = 0;X(k,1,2) = 0;X(k,1,3) = 0;X(k,1,4) = 0;
while sqrt((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) > epsi
X(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)];