南京理工大学
外弹道设计理论与方法
课程作业
学生姓名：学号：
学院：
课程：外弹道设计理论与方法
2015 年 1 月
作业：完成小口径火炮用榴弹的外弹道优化设计数学模型的程序设计。要求：(1) 模型：
，
，
分别代表
无量纲化后的变量；(2) 设计变量取m，ln，lb，lc；(3) 约束条件R域为：
(4) 已知条件：
，
，
，
m，
，
，
，
。
备注：
用0.4kg无量纲化，
用d＝0.035米无量纲化。
初始点
1、目标函数
，
，
分别代表
无量纲化后的变量；
二、优化设计变量
选取弹重m、弹头部长ln、圆柱部长lc、弹尾部长lb作为设计变量。
三、约束关系
约束条件R域为：
4、优化设计过程
4.1 建立优化模型
要求：
1) 弹丸飞行时间少
2) 弹丸气动外形好（阻力小）；
3) 弹丸飞行稳定；
4) 初速与弹重组合较佳；
5) 达到目标时的存速较大。
可以建立如下优化模型，并化为标准形式（1）。
目标函数:
（1）
优化模型中的输入参量为
，输出的优化目标为飞行时间最小。
优化模型可以整理为图1所示的流程。
图 1 优化流程图
输入初始参数后，优化程序调用优化算法和外弹道子程序，并进行循环计算判定优化结果是否达到所需的精度。优化程序的两个子模块是优化过程的关键：优化算法和外弹道解算。前者决定着初始点按照何种规则向最优解逼近，后者用于计算特定输入参数对应的弹道飞行时间。
4.2 外弹道解算模块
本文模型针对的是小口径榴弹的飞行时间解算。一方面外弹道优化要有很快的计算速度；另一方面弹丸在良好的飞行状态下，质点弹道基本上反映了弹丸的实
际飞行情况，因而在外弹道优化设计中采用质点弹道即可。外弹道模块的解算过程为：外弹道函数的输入为优化函数的自由变量矩阵X0，分别对应着弹丸质量、弹头部长度、弹底部长度和弹尾部长度，并将4个自由变量分别进行了无量纲化处理，输出为弹丸在空中的飞行时间t。
4.3 约束坐标轮换法模块
图 2 坐标轮换法示意图
外弹道优化设计问题一般归结为有约束非线性规划问题，有着自身的复杂性。对于有约束的非线性规划问题，在最优设计中有一个基于惩罚函数和障碍函数的序列无约束最小化方法，把求解一个有约束问题转化为求解无序约束的最优化问题。由于惩罚函数具有简单、易行等特点，在外弹道优化设计中，研究人员常采用惩罚函数法将有约束最优化问题转化为无约束最优化问题，然后采用优化理论中的直接方法，如模式搜索法、powell法等求解。
本文模型的函数为质点外弹道方程组，函数形式较为复杂，故采用简单有效的坐标轮换法。坐标轮换法的基本思想如图2所示，即，对于多变量输入
，可以先选择一维方向（将其他变量视为常值），按照一定的搜索方向从初始点开始搜索对应于目标的最优点。然后，用得到一维方向最优解去替代变量X中的对应值，再选择第二个方向进行搜索最优解，依次类推。在搜索完所有维度后对搜索后的最优解和搜索前的初值经行比较，如果满足收敛的精度要求则停止搜索，如不满足则进行下一轮的搜索。
五、优化模型的程序实现
选定好优化方法后，采用matlab软件进行相应的程序编制，以完成所需达到的计算要求。
5.1 程序的图形用户界面设计
本文设计的界面如图所示：
5.2 程序的运行
打开程序后，先核对程序的默认参数，然后点击“开始优化”按钮即可进入外弹道的优化过程，如上图所示。
6、优化结果
运行编制好程序，经过一定时间的计算得到如下结果：
附录
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)];
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.6
X(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 + 1
n = n + 1
end
xxx = x0;
mintime = min;
--------------------------------------------------------- function [xxx,mintime,alpha] = oneMIN(a,b,epsi,i,x0)
%global alpha
if i == 1
ii(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;
elseif i == 2
ii(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;
elseif i == 3
ii(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;
elseif i == 4
ii(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;
end
a1 = 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);
while abs(b - a) > epsi
if (f1 <= f2)
b = a2;
a2 = a1;
f2 = f1;
a1 = 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);
else
a = a1;
a1 = a2;
f1 = f2;
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);
end
end
alpha = (a + b) / 2;
xxx = [ii(1,1) * x0(1) + ii(1,2) * alpha, ii(2,1) * x0(2) + ii(2,2) * alpha, ii(3,1) * x0(3) + ii(3,2) * alpha, ii(4,1) * x0(4) + ii(4,2) * alpha];
mintime = fwddt(xxx);