空空导弹制导与控制系统设计

置顶樱桃小丸子123

已于 2024-10-22 10:40:34 修改

阅读量3k

点赞数 32

分类专栏：制导控制文章标签： matlab 算法

于 2024-10-21 12:07:11 首次发布

本文链接：https://blog.csdn.net/weixin_47115999/article/details/143100771

版权

制导控制专栏收录该内容

3 篇文章

订阅专栏

1、介绍

这篇文章主要通过构建空空导弹闭环控制系统，并结合线性二次调节器（LQR）与鲁棒控制，旨在优化导弹系统的响应性能，确保其在拦截目标时具备良好的稳定性和快速响应能力。文章将深入探讨空对空导弹对机动目标的拦截控制，具体目标如下：

重点分析鲁棒和自适应控制的。
结合制导律进行弹道仿真。

文章将详细阐述导弹的动态建模、控制器设计和系统分析，展示如何通过线性二次调节（LQR）方法实现高效的控制策略。此外，还通过案例仿真展示导弹在应对空中目标规避动作时的弹体响应控制及打击效果。

2、导弹弹道模型

AIM-120空空导弹长下面这样，我是基于这款空空导弹进行的研究：
在这里插入图片描述

2.1、导弹的气动计算公式

气动的计算直接使用如下公式进行计算，主要分为：

2.1.1、气动力计算

在这里插入图片描述

2.1.2 、气动力矩计算

在这里插入图片描述

2.2、导弹弹道方程的建模基础

对于刚体，可以应用牛顿定律来研究质心的移动，用动量矩定理研究刚体绕质心的转动。
在这里插入图片描述
导弹的空间运动可看成变质量物体的六自由度运动，由两个矢量方程描述。为研究方便起见，通常将矢量方程投影到坐标系上，写成三个描述导弹质心运动的动力学标量方程和三个导弹绕质心转动的动力学标量方程。

2.2.1、常用坐标系

导弹飞行力学中经常用到的坐标系有弹体坐标系、速度坐标系、地面坐标系和弹道坐标系，它们都是右手直角坐标系。

2.2.2、质心方程

导弹质心运动的动力学方程为：

在这里插入图片描述
导弹质心运动的运动学方程为：

2.2.3、绕质心方程

导弹绕质心的动力学方程为：
在这里插入图片描述
导弹绕质心运动学方程为：

2.2.4、导弹的弹道方程

综合前面所得的方程组即可得到导弹的弹道方程组为：
在这里插入图片描述

3、导弹的制导控制系统

导弹的制导控制系统是导弹飞行过程中用于实现对目标的追踪、飞行路径控制和命中目标的核心系统。它结合了传感器、计算机和执行机构，确保导弹能够在复杂动态环境下准确命中目标。以下是制导控制系统的主要组成部分和功能解释：
（1）制导系统
负责提供目标的位置信息，并根据导弹与目标的相对位置、速度和加速度来生成制导指令。常见的导引方法包括比例导引、追踪导引、预测导引等。
（2）控制系统
控制系统根据导引系统生成的制导指令，调整导弹的飞行姿态和运动轨迹。控制系统通过飞行控制面（如舵面）或推力矢量控制（TVC）来实现导弹的转向和加速。
控制算法通常采用经典的控制理论（如PID控制），也可能使用现代控制方法（如自适应控制、鲁棒控制、模糊控制等），以应对飞行过程中的非线性、不确定性和干扰。
（3）导航系统
导航系统用于实时确定导弹的飞行状态，包括位置、速度、加速度等参数。常用的导航方式包括惯性导航系统（INS）、全球导航卫星系统（GNSS）或者两者的组合导航（INS/GNSS）。陀螺仪和加速度计用于测量导弹的角速度和线加速度，经过积分运算得出导弹的速度和位置。
（4）飞行控制执行机构
执行机构接收控制系统的指令来调整导弹的舵面或者喷管，直接影响导弹的姿态和轨迹。执行机构通常包括伺服电机和液压装置等机械系统。

工作流程如下：

目标检测与跟踪：导弹在发射后，导引系统利用传感器（雷达、红外、激光等）获取目标的实时信息。
制导指令生成：导引系统计算导弹和目标的相对位置、速度，并生成相应的制导指令，比如通过比例导引律计算导弹所需的加速度。
姿态与轨迹控制：控制系统接收制导指令，生成控制信号，指挥执行机构调整导弹的姿态和飞行路径，确保导弹朝目标飞行。
飞行修正：导弹在飞行过程中，不断修正飞行轨迹和姿态，直到命中目标。

4、空对空导弹控制系统的鲁棒和自适应控制

建立一个空对空导弹拦截目标的闭环控制系统，并使用线性二次调节器（LQR）与鲁棒控制。
下面将从初始设置、系统模型、执行器模型、LQR控制设计、导弹控制系统的响应特性、舵面最大偏转角和速率、系统稳定性、增益和相位裕度、闭环系统的设计、状态空间模型、LQR控制器、观测器设计、灵敏度和灵敏度分析等模块逐步解析代码的框架和功能。

4.1、初始设置和变量定义

首先定义了仿真所需的基本参数，如频率范围、模拟时间、输入幅值等。w_rps 定义了频率范围，R2D 是弧度到角度的转换常数，time 设置了仿真时间。

clc;
clear all;
close all;

w_rps = logspace(-3,3,1000); % 频率范围 [10^(-3), 10^3] rad/s
R2D = (180/pi);              % 弧度转度的转换系数
time = 30;                   % 模拟时间
amp_g = 32.2 * 20;           % 20g的阶跃输入

4.2、系统模型

系统的开环状态矩阵 A 和输入矩阵 B 描述了飞行器的动态方程和控制输入的作用方式。代码计算系统的开环特征值并显示，以了解系统的动态特性。

A = [-1.57        0      0     1     0; 
         0     -.50    .17     0    -1; 
    -21.13  -2876.7  -2.10  -.14  -.05;
    -82.92   -11.22   -.01  -.57     0;
      -.19   -11.86   -.01     0  -.57];
  
B = [     0    -.10       0;
       -.07       0     .11;
    -1234.7  -30.49 -1803.2;
      -4.82 -119.65      -7;
      14.84     .27 -150.58];

eigs_ol = eig(A);
disp('Open Loop EigenValues');
disp(eigs_ol);

4.3、短周期动态模型

短周期动态模型简化了飞行器的动态方程，主要用于分析系统的短期响应。这里使用状态空间模型描述了短周期系统，并通过特征值分析其稳定性。

Asp = [-1.57   1.00;
      -82.92 -0.57];
Bsp = [-0.10;
         -119.65];
Csp = [Z_alpha 0; 1 0; 0 1];
Dsp = [Z_dele; 0; 0];

sys_sp_plant = ss(Asp,Bsp,Csp,Dsp);
Eigs_sp = eig(Asp);
disp('Short Period EigenValues');
disp(Eigs_sp);

4.4、执行器模型

执行器模型通过状态空间描述了执行器的动态特性。执行器的自然频率 wa_rps 和阻尼比 za 用于定义系统的响应。

wa_rps = 35*(2*pi);  % 自然频率，单位为rad/s
za = .71;            % 阻尼比
Aact = [0 1; -wa_rps^2  -2*za*wa_rps];
Bact = [0 ; wa_rps^2];
Cact = [1 0];
Dact = 0;
sys_actuator =  ss(Aact,Bact,Cact,Dact,'StateName', {'dele','dele_dot'},'InputName',{'dele_cmd'},'OutputName',{'dele'});

4.5、 LQR 控制设计

线性二次调节器（LQR）控制器用于设计最优控制律，代码中通过定义状态和控制权重矩阵 Q 和 R 来实现控制器设计。接着计算闭环系统的动态并分析其阶跃响应。

Q = zeros(3);  % 状态权重矩阵
R = 1;         % 控制输入权重

t = linspace(0, 2, 50); %time scale
q_eI  = logspace(-6,-3,100);
N = length(q_eI);

for ii = 1 : length(q_eI)
    Q(1,1) = q_eI(ii);
    [Kw,~,~] = lqr(Aw,Bw,Q,R); % LQR solver
    % Dynamic Controller Matrices
    Ac  = [0];
    Bc1 = [1 0 0];
    Bc2 = [-1];
    Cc  = -Kw(1);
    Dc1 = [0 -Kw(2:3)];
    Dc2 = [0];
    
    % Produce closed loop ss system
    sys_cl_act = ss(Acl_act,Bcl_act,Ccl_act,Dcl_act,'StateName',{'Az','alpha','q','dele','dele_dot'},'InputName',{'Az_cmd'},'OutputName',{'Az'});
    
    opt = stepDataOptions('StepAmplitude', amp_g); % Add in 20g
    [ydesign,tdesign,xdesign] = step(sys_cl_act,opt);
end

4.6、闭环系统动态分析

闭环系统结合了执行器和飞行器动态模型，并通过反馈控制器形成完整的控制系统。代码计算了系统的阶跃响应和相关的性能指标，如超调、上升时间等。

sys_cl_act = ss(Acl_act,Bcl_act,Ccl_act,Dcl_act,'StateName',{'Az','alpha','q','dele','dele_dot'},'InputName',{'Az_cmd'},'OutputName',{'Az'});
[yd,td,xd] = step(sys_cl_act,opt);
S = stepinfo(sys_cl_act);

4.6.1、闭环系统的响应特性

对闭环系统的阶跃响应进行分析来评估系统的性能。在控制系统中，阶跃响应能够直观反映系统对输入信号的动态反应，包括上升时间、超调、稳态误差等关键指标。

% 定义阶跃响应的输入选项，amp_g为阶跃信号的幅值
opt = stepDataOptions('StepAmplitude', amp_g); 
[ydesign,tdesign,xdesign] = step(sys_cl_act,opt); % 获取闭环系统的阶跃响应

% 获取系统响应的相关信息
S = stepinfo(sys_cl_act); 

% 计算最大舵角偏转和舵角速率
dmax(ii)  =  max(abs(xdesign(:,3)))*(180/pi); % 舵角偏转角度
ddmax(ii) =  max(abs(xdesign(:,4)))*(180/pi); % 舵角变化速率

4.6.2、闭环增益与稳定性裕度

为了确保导弹系统在各种不确定条件下的鲁棒性，计算了系统的环路增益和增益交叉频率，并分析了系统的增益裕度和相位裕度。通过计算增益裕度和相位裕度，代码能够评估系统的稳定性裕度，从而判断其是否在不确定性情况下保持稳定。增益裕度越大，系统能够承受的增益变化就越大，而不导致不稳定。相位裕度越大，系统对时滞的鲁棒性也越强。

% 定义增益矩阵
A_Lu   = [Ap zeros(4,1); Bc1*Cp Ac];
B_Lu   = [Bp; Bc1*Dp];
C_Lu   = -[Dc1*Cp Cc]; % 取负号以计算环路增益
D_Lu   = -[Dc1*Dp];

sys_Lu = ss(A_Lu,B_Lu,C_Lu,D_Lu); % 创建状态空间模型

% 获取增益裕度和相位裕度
S_Lu       = allmargin(sys_Lu); 
MGM        = max(S_Lu.GainMargin); % 最大增益裕度
GM(ii)     = 20.*log10(MGM); % 将增益裕度转化为分贝(dB)
PM(ii)     = S_Lu.PhaseMargin; % 相位裕度
w_rpsp(ii) = S_Lu.PMFrequency; % 相位裕度对应的频率
w_cg(ii)   = w_rpsp(ii)/2/pi;  % 交叉频率 (Hz)

4.6.3、敏感度和协同敏感度分析

敏感度函数和协同敏感度函数在控制系统的设计中扮演着关键角色。敏感度函数反映系统对输入扰动的敏感程度，而协同敏感度函数则衡量系统对外部噪声或不确定性变化的抵抗能力。

% 计算输入处的返回差异
sImAinv = ((w_rpsp(ii))*eye(max(size(A_Lu)))-A_Lu)\eye(max(size(A_Lu)));
Lu_eval = C_Lu*sImAinv*B_Lu + D_Lu; % 计算环路增益

% 返回差异
RD = svd(eye(max(size(Lu_eval))) + Lu_eval);
RDdesired(ii) = min(RD);

% 计算敏感度和协同敏感度
sensitivity           = svd(1/(1+Lu_eval)); 
sensitivity_dB(ii)    = abs(20*log10(sensitivity)); % 将敏感度转化为分贝
cosensitivity         = svd(Lu_eval/(1+Lu_eval));
cosensitivity_dB(ii)  = abs(20*log10(cosensitivity));

4.6.4、系统的极点与阻尼比

系统的稳定性不仅依赖于增益和相位裕度，还取决于其极点分布。通过分析极点位置和阻尼比，可以判断系统是否具有足够的阻尼以避免振荡，确保导弹在拦截目标时的平稳性。

% 获取闭环系统的阻尼比和极点
[Wn,Z,P]         = damp(sys_cl_act); 
damping_min(ii)  = min(Z); % 获取最小阻尼比
Eigs_Cl(:,ii)    = P;      % 保存闭环系统的极点

4.7、观测器设计

设计一个基于Leuenberger观测器的控制系统。观测器的主要作用是通过输出反馈来估计系统的状态，从而在不完全测量的情况下依然能够实现良好的控制效果。通过设计观测器，系统能够估计出一些难以直接测量的状态，从而提高控制精度。代码同时通过仿真展示了观测器与LQR控制器的组合效果，比较了估计状态与实际状态的差异。

% 使用 lqr() 设计观测器增益矩阵
%%  观测器设计
%% Leuenberger Observer
Cw_meas  = [1 0 0 ; 0 0 1];
nCw_meas = size(Cw_meas,1);
Qo   =  diag([1e4 1e4 1e5]); 
Ro   =  [1 0; 0 1];
ol   = [0;0];

%Define Observer Model
Awob = [0 C_reg ;ol Asp];  nAw = size(Awob,1);
Bwob = [D_reg;Bsp];
nBsp = size(Bsp,1);
Bcmd = [-1;0;0];
Cwob = [1 0 0];
Dwob = 0;

% Observer Integrators
Ai  = 0; nAi = size(Ai,1);
Bi1 = [1 0];
Bi2 = -1;
Ci  = [1;0]; %Error
Di1 = [0 0; 0 1];
Di2 = [0;0];

%solve dual ARE
% 使用 lqr() 函数设计观测器增益矩阵 L,Awob为观测器的系统矩阵,Cw_meas为测量输出矩阵
L = lqr(Awob',Cw_meas',Qo,Ro);
L = L';

% Building the state space model of controller,observer, integrator model
Ac  = [Ai 0 0 0; L*Ci (Awob - Bwob*(Kw) - L*Cw_meas)];
Bc1 = [Bi1;L*Di1];
Bc2 = [Bi2; L*Di2+Bcmd]; nBci2 = size(Bc2,2);
Cc  = [0 -Kw];
Dc1 = zeros(1,size(Bc1,2));
Dc2 = 0;
sys_obs_comp = ss(Ac,Bc1,Cc,Dc1);

%observer loop gain
A_Lu   = [Ap zeros(4,4); Bc1*Cp([1,3],:)  Ac];
B_Lu   = [Bp; Bc1*Dp([1,3],:) ];
C_Lu   = -[Dc1*Cp([1,3],:)  Cc]; % Change sign for loop gain
D_Lu   = -[Dc1*Dp([1,3],:) ]; % Change signe for loop gain
sys_Lu_ob = ss(A_Lu,B_Lu,C_Lu,D_Lu);

sImAinv_OB = ((w_rpsp(q_eI_datapoint))*eye(max(size(A_Lu)))-A_Lu)\eye(max(size(A_Lu)));
Lu_eval_OB  = C_Lu*sImAinv_OB*B_Lu + D_Lu;

% Closing the Loop
H = Dc1*Dp([1,3],:) ; nH = size(H,1);
IZ = inv(eye(nH) - H);
PL = Dp([1,3],:)*IZ*Dc1; nPL = size(PL,1);
IPL = eye(nPL) - PL;
A = [Ap+Bp*IZ*Dc1*Cp([1,3],:)    Bp*IZ*Cc;     Bc1*IPL*Cp([1,3],:)   Ac+Bc1*Dp([1,3],:)*IZ*Cc];
B = [Bp*IZ*Dc2 ; Bc2+Bc1*Dp([1,3],:)*IZ*Dc2];
C = [(IPL)*Cp([1,3],:)  Dp([1,3],:)*IZ*Cc];
D = Dp([1,3],:)*IZ*Dc2;

sys_combined = ss(A,B,C,D);

opt = stepDataOptions('StepAmplitude', amp_g); % Add in 20g
[yob,tob,xob] = step(sys_combined,opt);