clc;
clear;
close all;
global setup
% 基础参数设置
setup.n = [30]; % LG点个数
setup.ns = [5]; % 状态变量个数
setup.nu = [1]; % 控制变量个数
Re = 6378145; % 地球半径
rho0 = 1.225; % 海平面大气密度
H = 7200; % 大气范围
g0 = 9.80665; % 重力加速度
s = 15; % 参考面积
mach = 340;
setup.Re = Re;
setup.rho0 = rho0;
setup.s = s;
setup.H = H;
setup.g0 = g0;
setup.mach = mach;
mI = 656080; % 运载器起飞质量
lI = 0.9; % 初始推重比
F0 = mI * lI * g0;
F0_ = 1.1 * F0;
setup.F0_ = F0_;
% 过程约束
qmax = 50e3;
nmax = 5;
setup.qmax = qmax;
setup.nmax = nmax;
m0 = 136080; % 初始质量
mf = 60400; % 末端质量
t0 = 0; % 初始时间
t1 = 200;
r0 = 2000 + Re; % 初始地心距
v0 = 200; % 初始速度
gamma0 = 30/57.3; % 初始弹道倾角
rf = 65000 + Re; % 终端高度
vf = 5 * mach; % 终端速度
gammaf = 20/57.3; % 终端弹道倾角
alpha0 = 20/57.3; % 起始攻角
alphaf = 0/57.3; % 终端攻角
% 无量纲化参数
tc = (Re/g0)^0.5; % 时间尺度
rc = Re; % 长度尺度
vc = (Re*g0)^0.5; % 速度尺度
Pc = m0 * g0; % 力尺度
mc = m0; % 质量尺度
setup.tc = tc;
setup.rc = rc;
setup.vc = vc;
setup.Pc = Pc;
setup.mc = mc;
% 保存原始设置用于蒙特卡洛
orig_setup = setup;
orig_m0 = m0;
orig_r0 = r0;
orig_v0 = v0;
% 蒙特卡洛参数
num_runs = 50; % 运行次数
initial_velocities = zeros(num_runs, 1); % 存储初始速度
final_velocities = zeros(num_runs, 1); % 存储终端速度
% 初始速度随机范围
v0_min = 150; % 最小初始速度 (m/s)
v0_max = 250; % 最大初始速度 (m/s)
% 优化选项设置
options = optimset('TolX',1e-8, 'TolFun',1e-6, 'TolCon',1e-6, ...
'MaxFunEvals',50000, 'Algorithm','sqp', ...
'MaxIter',50, 'Display','off');
% 主循环 - 蒙特卡洛模拟
for run = 1:num_runs
fprintf('蒙特卡洛运行 %d/%d...\n', run, num_runs);
% 重置全局设置
setup = orig_setup;
% 随机生成初始速度 (在指定范围内均匀分布)
v0_perturbed = v0_min + (v0_max - v0_min) * rand();
initial_velocities(run) = v0_perturbed;
% 其他初始状态保持不变
r0_perturbed = orig_r0;
m0_perturbed = orig_m0;
% 重新计算无量纲初始状态
r0dot = r0_perturbed / rc;
v0dot = v0_perturbed / vc;
gamma0dot = gamma0;
m0dot = m0_perturbed / mc;
alpha0dot = alpha0;
X0 = [r0dot; v0dot; gamma0dot; m0dot; alpha0dot];
setup.X0 = X0;
setup.r0dot = r0dot;
setup.v0dot = v0dot;
setup.m0dot = m0dot;
setup.gamma0dot = gamma0dot;
setup.alpha0dot = alpha0dot;
% 重新初始化阶段参数
phase = struct();
iphase = 1;
phase.t0{iphase} = t0/tc;
phase.tf{iphase} = t1/tc;
phase.m0{iphase} = m0_perturbed/mc;
phase.mf{iphase} = mf/mc;
% 生成Legendre-Gauss点
[phase.tau{iphase}, phase.w{iphase}, phase.D{iphase}] = LG_Points(setup.n(iphase));
% 初始猜测
phase.X{iphase} = [r0dot v0dot gamma0dot phase.m0{iphase} alpha0dot;
rf/rc vf/vc gammafdot phase.mf{iphase} alphafdot];
phase.Xguess{iphase} = interp1([-1;1], phase.X{iphase}, [-1; phase.tau{iphase}; 1], 'spline');
phase.Uguess{iphase} = (1/57.3) * ones(setup.n(iphase), 1);
phase.x0guess{iphase} = [reshape(phase.Xguess{iphase}, [], 1); reshape(phase.Uguess{iphase}, [], 1)];
% 边界约束
rmindot = Re / rc;
rmaxdot = 2 * Re / rc;
vmindot = -10000 / vc;
vmaxdot = -vmindot;
gammamindot = 0;
gammamaxdot = 80/57.3;
alphamindot = -10/57.3;
alphamaxdot = 30/57.3;
dalphamin = -5/57.3;
dalphamax = 5/57.3;
phase.XU{iphase} = [rmaxdot*ones(setup.n(iphase)+2,1) vmaxdot*ones(setup.n(iphase)+2,1) ...
gammamaxdot*ones(setup.n(iphase)+2,1) phase.m0{iphase}*ones(setup.n(iphase)+2,1) ...
alphamaxdot*ones(setup.n(iphase)+2,1)];
phase.UU{iphase} = dalphamax * ones(setup.n(iphase), 1);
phase.xU{iphase} = [reshape(phase.XU{iphase}, [], 1); reshape(phase.UU{iphase}, [], 1)];
phase.XL{iphase} = [rmindot*ones(setup.n(iphase)+2,1) vmindot*ones(setup.n(iphase)+2,1) ...
gammamindot*ones(setup.n(iphase)+2,1) phase.mf{iphase}*ones(setup.n(iphase)+2,1) ...
alphamindot*ones(setup.n(iphase)+2,1)];
phase.UL{iphase} = dalphamin * ones(setup.n(iphase), 1);
phase.xL{iphase} = [reshape(phase.XL{iphase}, [], 1); reshape(phase.UL{iphase}, [], 1)];
setup.phase = phase;
x0 = [phase.x0guess{1}; phase.tf{1}];
xU = [phase.xU{1}; phase.tf{1}];
xL = [phase.xL{1}; phase.t0{1}];
% 设置终端约束
setup.oe = [rf/rc, vf/vc]; % 无量纲高度和速度
% 运行优化
try
[x, ~] = fmincon(@Objective, x0, [], [], [], [], xL, xU, @NonlconFun, options);
% 提取终端速度
n_points = setup.n(1) + 2;
v_final_index = n_points + 1; % 第二个状态变量的终端值索引
v_final_dot = x(v_final_index);
final_velocities(run) = v_final_dot * vc;
fprintf('运行 %d 完成: 初始速度 = %.2f m/s, 终端速度 = %.2f m/s\n', ...
run, v0_perturbed, final_velocities(run));
catch
final_velocities(run) = NaN;
fprintf('运行 %d 失败\n', run);
end
end
% 过滤失败的结果
valid_runs = ~isnan(final_velocities);
initial_velocities = initial_velocities(valid_runs);
final_velocities = final_velocities(valid_runs);
% 生成散点图
if ~isempty(final_velocities)
figure('Color', 'white', 'Position', [100, 100, 900, 650]);
% 创建散点图
scatter(initial_velocities, final_velocities, 80, 'filled', ...
'MarkerFaceColor', [0.1, 0.5, 0.8], ...
'MarkerEdgeColor', 'k', ...
'MarkerFaceAlpha', 0.7, ...
'LineWidth', 1.2);
% 添加网格和美化
grid on;
box on;
set(gca, 'FontSize', 12, 'LineWidth', 1.2, 'GridAlpha', 0.3);
% 添加标题和标签
title('终端速度 vs 初始速度', 'FontSize', 18, 'FontWeight', 'bold');
xlabel('初始速度 (m/s)', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('终端速度 (m/s)', 'FontSize', 14, 'FontWeight', 'bold');
% 添加参考线 (目标终端速度)
hold on;
target_vf = vf;
yline(target_vf, 'r--', 'LineWidth', 2.5, 'Alpha', 0.8);
% 添加理论最大终端速度线
theoretical_max = 1.05 * max(final_velocities);
yline(theoretical_max, 'g-.', 'LineWidth', 2, 'Alpha', 0.7);
% 添加趋势线
p = polyfit(initial_velocities, final_velocities, 1);
x_fit = linspace(min(initial_velocities), max(initial_velocities), 100);
y_fit = polyval(p, x_fit);
plot(x_fit, y_fit, 'b-', 'LineWidth', 2.5);
% 添加统计信息
corr_coef = corr(initial_velocities, final_velocities);
mean_final = mean(final_velocities);
std_final = std(final_velocities);
text(min(initial_velocities)+5, max(final_velocities)-100, ...
{sprintf('相关系数: %.4f', corr_coef), ...
sprintf('终端速度均值: %.2f m/s', mean_final), ...
sprintf('终端速度标准差: %.2f m/s', std_final)}, ...
'FontSize', 12, 'BackgroundColor', 'white', 'EdgeColor', 'k');
% 添加图例
legend('优化结果', '目标终端速度', '理论最大速度', '线性趋势线', ...
'Location', 'southeast', 'FontSize', 10);
% 添加颜色条表示速度增量
delta_v = final_velocities - initial_velocities;
colorbar_label = '速度增量 (m/s)';
% 保存结果
save('velocity_analysis.mat', 'initial_velocities', 'final_velocities');
saveas(gcf, 'terminal_vs_initial_velocity.png');
disp('分析完成,结果已保存。');
else
disp('所有运行均失败,无法生成结果。');
end
% 辅助函数:生成Legendre-Gauss点
function [tau, w, D] = LG_Points(n)
[tau, w] = legendre_points(n);
D = collocation_matrix(tau);
end
function [x, w] = legendre_points(n)
% 返回n阶Legendre多项式的根(LG点)和权重
beta = 0.5./sqrt(1-(2*(1:n-1)).^(-2)); % 三项递推系数
T = diag(beta,1) + diag(beta,-1); % Jacobi矩阵
[V, D] = eig(T); % 特征值分解
x = diag(D); % LG点
w = 2*V(1,:).^2; % 权重
% 排序
[x, idx] = sort(x);
w = w(idx)';
end
function D = collocation_matrix(tau)
% 计算微分矩阵
n = length(tau);
D = zeros(n);
for j = 1:n
for k = 1:n
if j ~= k
D(j,k) = lagrange_derivative(tau, j, k);
end
end
end
for j = 1:n
D(j,j) = -sum(D(j,:));
end
end
function d = lagrange_derivative(tau, j, k)
% 计算Lagrange基函数在节点处的导数
n = length(tau);
prod_val = 1;
for m = 1:n
if m ~= j && m ~= k
prod_val = prod_val * (tau(k) - tau(m))/(tau(j) - tau(m));
end
end
d = prod_val / (tau(j) - tau(k));
end
% 目标函数
function J=Objective(x)
global setup
J=x(end)*setup.tc;
end
% 非线性约束函数
function [c, ceq] = NonlconFun(x)
global setup
% 提取参数
rho0 = setup.rho0;
Re = setup.Re;
H = setup.H;
s = setup.s;
g0 = setup.g0;
qmax = setup.qmax;
nmax = setup.nmax;
rc = setup.rc;
vc = setup.vc;
tc = setup.tc;
mc = setup.mc;
n = setup.n;
ns = setup.ns;
nu = setup.nu;
np = length(n);
% 初始化约束
ceq = [];
c = [];
% 处理优化变量
nx = n(1)*(ns(1)+nu(1)) + 2*ns(1);
xi{1} = x(1:nx);
% 初始化存储
Xi = cell(1, np);
Ui = cell(1, np);
% 提取终端时间
tf{1} = x(end);
% 处理每个阶段
for i1 = 1:np
% 提取状态变量
X = zeros(n(i1)+2, ns(i1));
for i2 = 1:ns(i1)
start_idx = (i2-1)*(n(i1)+2) + 1;
end_idx = i2*(n(i1)+2);
X(:, i2) = xi{i1}(start_idx:end_idx);
end
Xi{i1} = X;
% 提取控制变量
U = zeros(n(i1), nu(i1));
start_control = ns(i1)*(n(i1)+2);
for i3 = 1:nu(i1)
start_idx = start_control + (i3-1)*n(i1) + 1;
end_idx = start_control + i3*n(i1);
U(:, i3) = xi{i1}(start_idx:end_idx);
end
Ui{i1} = U;
% 计算动力学方程
F = MyDeo(X, U, i1);
% 动力学约束
t0 = setup.phase.t0{i1};
D = setup.phase.D{i1};
ceqd = D*X(1:end-1, :) - (tf{i1}-t0)/2*F;
ceqdr = reshape(ceqd, [], 1);
% 终端约束
w = setup.phase.w{i1};
ceqe = X(1,:) + (tf{i1}-t0)/2*w'*F - X(end,:);
ceqe = ceqe(:);
% 路径约束 - 动压约束
v = X(2:end-1, 2); % 无量纲速度
r = X(2:end-1, 1); % 无量纲地心距
r_dim = r * rc; % 实际地心距
rho = rho0 * exp(-(r_dim - Re)/H);
q = 0.5 * rho .* (v * vc).^2; % 实际动压
cq = q - qmax; % 动压约束
% 路径约束 - 法向过载
Cl = 0.5; % 升力系数(简化假设)
L = 0.5 * rho .* (v * vc).^2 * s * Cl; % 升力
n_load = L ./ (X(2:end-1, 4) * mc * g0); % 法向过载
cn = n_load - nmax; % 过载约束
% 组合约束
ceq = [ceq; ceqdr; ceqe];
c = [c; cq; cn];
end
% 初始状态约束
ceqs = Xi{1}(1, :) - setup.X0';
ceq = [ceq; ceqs(:)];
% 终端约束
rf = Xi{1}(end, 1);
Vf = Xi{1}(end, 2);
oe = setup.oe;
ceqo = [rf; Vf] - oe(:);
ceq = [ceq; ceqo];
end
% 动力学方程函数
function F = MyDeo(X, U, iphase)
global setup
% 提取参数
Re = setup.Re;
rho0 = setup.rho0;
H = setup.H;
g0 = setup.g0;
s = setup.s;
F0_ = setup.F0_;
rc = setup.rc;
vc = setup.vc;
tc = setup.tc;
mc = setup.mc;
% 状态变量:r, V, gamma, m, alpha
r = X(2:end-1, 1); % 无量纲地心距
V = X(2:end-1, 2); % 无量纲速度
gamma = X(2:end-1, 3); % 弹道倾角(弧度)
m = X(2:end-1, 4); % 无量纲质量
alpha = U(:, 1); % 攻角(控制变量,弧度)
% 计算有量纲值
r_dim = r * rc; % 地心距(m)
V_dim = V * vc; % 速度(m/s)
m_dim = m * mc; % 质量(kg)
% 大气密度
height = r_dim - Re; % 高度
rho = rho0 * exp(-height / H);
% 空气动力学系数(简化模型)
Cd = 0.5 + 0.1*alpha.^2; % 阻力系数
Cl = 2*pi*alpha; % 升力系数
% 气动力
q = 0.5 * rho .* V_dim.^2; % 动压
D = q * s .* Cd; % 阻力(N)
L = q * s .* Cl; % 升力(N)
% 推力(假设恒定)
Thrust = F0_; % 固定推力(N)
% 重力
g = g0 * (Re ./ r_dim).^2; % 重力加速度(m/s^2)
% 动力学方程(有量纲)
drdt = V_dim .* sin(gamma); % 高度变化率(m/s)
dVdt = (Thrust .* cos(alpha) - D) ./ m_dim - g .* sin(gamma); % 速度变化率(m/s^2)
dgammadt = (Thrust .* sin(alpha) + L) ./ (m_dim .* V_dim) - (g ./ V_dim - V_dim ./ r_dim) .* cos(gamma); % 弹道倾角变化率(rad/s)
dmdt = -Thrust / (300 * g0); % 质量变化率(kg/s),假设比冲为300秒
% 攻角变化率(控制量)
dalphadt = U(:, 1); % 控制变量就是攻角变化率(rad/s)
% 将状态变化率无量纲化
drdt_dot = drdt / (rc/tc); % 无量纲位置变化率
dVdt_dot = dVdt / (vc/tc); % 无量纲速度变化率
dgammadt_dot = dgammadt * tc; % 无量纲弹道倾角变化率
dmdt_dot = dmdt / (mc/tc); % 无量纲质量变化率
dalphadt_dot = dalphadt * tc; % 无量纲攻角变化率
% 组装F
F = [drdt_dot, dVdt_dot, dgammadt_dot, dmdt_dot, dalphadt_dot];
end出现以下报错:函数或变量 'gammafdot' 无法识别。
出错 Monta_v (第 129 行)
rf/rc vf/vc gammafdot phase.mf{iphase} alphafdot];