基于MATLAB程序PID控制/参数自动优化+飞机的小扰动线性化建模

基于MATLAB程序PID控制/参数自动优化+飞机的小扰动线性化建模

说明：纯m代码根据教科书上的参考数据，完成小扰动线性化建模，计算飞机横纵向两个状态方程；
以纵向状态方程为例，使用pid控制律，进行姿态控制，同时自动计算控制器的性能指标，并设计程序自动优化控制参数

以下文字及示例代码仅供参考

• ✅ 飞机纵向与横向运动的小扰动线性化建模
• ✅ 基于状态空间模型的PID控制器设计
• ✅ 系统响应性能指标计算（如上升时间、超调量等）
• ✅ 使用 fmincon 或 ga 实现 PID参数自动优化
• 📌 全部使用纯 .m 文件编写，便于调试和扩展

📌 资料说明

✅ 功能模块

1. 飞机动力学建模

   • 纵向与横向运动方程推导
   • 小扰动线性化处理
   • 构造状态空间模型 A, B, C, D

2. PID 控制器设计

   • 基于反馈误差设计控制律
   • 支持姿态角跟踪（俯仰角、高度等）

3. 系统性能评估

   • 计算阶跃响应的性能指标：上升时间、调节时间、超调量、稳态误差

4. 参数自动优化

   • 使用 MATLAB 优化工具箱（如 fmincon, ga）进行参数整定

5. 仿真验证

   • 使用 Simulink 或 m 文件进行闭环仿真
   • 绘图展示系统响应

🧠 核心原理简述

1️⃣ 小扰动线性化建模（参考《Aircraft Dynamics and Automatic Control》）

纵向运动方程（以速度、迎角、俯仰角速率、俯仰角为状态）：

$$
\begin{bmatrix}
\dot{u} \
\dot{\alpha} \
\dot{q} \
\dot{\theta}
\end{bmatrix}

\begin{bmatrix}
X_u & X_\alpha & 0 & -g \cos(\theta_0) \
Z_u/U_0 & Z_\alpha/U_0 & (U_0 + Z_q)/U_0 & -g \sin(\theta_0)/U_0 \
M_u & M_\alpha & M_q & M_\theta \
0 & 0 & 1 & 0
\end{bmatrix}
\cdot
\begin{bmatrix}
u \
\alpha \
q \
\theta
\end{bmatrix}
+
\begin{bmatrix}
X_{\delta_e}/U_0 \
Z_{\delta_e}/U_0 \
M_{\delta_e} \
0
\end{bmatrix}
\cdot \delta_e
$$

其中：

• $ u $: 速度扰动
• $ \alpha $: 迎角
• $ q $: 俯仰角速率
• $ \theta $: 俯仰角
• $ \delta_e $: 升降舵偏转
• $ U_0 $: 平衡飞行速度
• $ g $: 重力加速度
• 系数如 $ X_u, Z_\alpha, M_q $ 等来自气动导数（可查手册或文献）

📁 文件结构（建议）

Aircraft_PID_Control/
│
├── build_state_space.m     # 构造状态空间矩阵 A/B/C/D
├── pid_controller.m        # PID 控制律函数
├── performance_index.m     # 性能指标计算函数
├── optimize_pid.m          # PID参数自动优化主程序
├── simulate_pid.m          # 主仿真脚本
└── utils.m                 # 工具函数（如绘图、归一化等）

🧪 示例代码（简化版）

1️⃣ build_state_space.m —— 构造纵向状态空间模型

function [A, B] = build_longitudinal_model()
% 参数设置（单位统一为 SI 制）
U0 = 200;       % 平衡速度 (m/s)
g = 9.81;       % 重力加速度
theta0 = 0;     % 平飞

% 气动导数（示例值，可根据实际飞机数据替换）
Xu = -0.05;
Xa = 0.1;
Zu = -0.2;
Za = -3.0;
Zq = -1.0;
Me = -0.5;
Mu = 0.01;
Ma = 0.3;
Mq = -2.0;

% 状态方程系数矩阵 A
A = [
    Xu,           Xa,             0,              -g*cos(theta0);
    Zu/U0,        Za/U0,          (U0 + Zq)/U0,   -g*sin(theta0)/U0;
    Mu,           Ma,             Mq,             0;
    0,            0,              1,              0
];

% 输入矩阵 B
Be = [
    Xa/U0;
    Za/U0;
    Me;
    0
];
B = Be;
end

2️⃣ pid_controller.m —— PID 控制律函数

function u = pid_controller(e, edot, eint, Kp, Ki, Kd)
% PID 控制律
u = Kp * e + Ki * eint + Kd * edot;
end

3️⃣ performance_index.m —— 性能指标计算函数

function J = performance_index(y, t, y_ref)
% 计算性能指标：ITAE + 超调量惩罚
e = abs(y - y_ref);

% ITAE 指标
J_itae = sum(t' .* e') * mean(diff(t));

% 超调量
overshoot = max(y) - y_ref(end);
if overshoot < 0, overshoot = 0; end

% 总目标函数
J = J_itae + 10 * overshoot;
end

4️⃣ simulate_pid.m —— 主仿真脚本（含闭环控制）

function [t, y, u] = simulate_pid(Kp, Ki, Kd)
% 设置采样时间和仿真时间
dt = 0.1;
t_sim = 0:dt:30;

% 获取状态空间模型
[A, B] = build_longitudinal_model();
C = eye(4); D = zeros(4,1);
sys = ss(A, B, C, D, dt);

% 初始状态
x0 = [0; 0; 0; 0];  % u=0, alpha=0, q=0, theta=0

% 参考信号（俯仰角变化）
y_ref = 5 * pi/180 * ones(size(t_sim));  % 5度

% 初始化变量
y = [];
u = [];
eint = 0;
prev_error = 0;

for k = 1:length(t_sim)
    % 当前输出
    x = lsim(sys, u(end+1,:), t_sim(k), x0)';
    y(end+1,:) = x(4);  % 输出俯仰角

    % 计算误差
    error = y_ref(k) - y(end);
    eint = eint + error * dt;
    edot = (error - prev_error) / dt;

    % PID 控制
    u(end+1,:) = pid_controller(error, edot, eint, Kp, Ki, Kd);

    % 更新状态
    x0 = x;
    prev_error = error;
end

% 返回结果
t = t_sim;
end

5️⃣ optimize_pid.m —— PID 参数自动优化主程序

function best_params = optimize_pid()
% 目标函数
obj_func = @(params) objective_function(params);

% 初始猜测
K0 = [1, 0.1, 0.5];

% 边界约束
lb = [0.1, 0, 0];
ub = [10, 2, 5];

% 使用遗传算法优化
options = optimoptions('ga', 'Display', 'iter', 'PopulationSize', 50, ...
                       'MaxGenerations', 100);
[best_params, fval] = ga(obj_func, 3, [], [], [], [], lb, ub, [], options);

% 显示最优参数
disp(['Best PID parameters: Kp=', num2str(best_params(1)), ...
      ', Ki=', num2str(best_params(2)), ', Kd=', num2str(best_params(3))]);
end

% 目标函数定义
function J = objective_function(params)
Kp = params(1);
Ki = params(2);
Kd = params(3);

% 仿真获取响应
[t, y, u] = simulate_pid(Kp, Ki, Kd);
y_ref = 5 * pi/180 * ones(size(t));

% 计算性能指标
J = performance_index(y, t, y_ref);
end

🧩 拓展功能建议

• ✅ 添加横向通道控制（滚转、侧滑）
• ✅ 使用非线性模型验证控制器鲁棒性
• ✅ 加入传感器噪声和执行机构延迟
• ✅ 使用 Simulink 实现可视化仿真
• ✅ 引入在线辨识技术更新模型参数

📦 数据集 & 文献参考

气动导数可参考如下资料：

• Brian L. Stevens, Frank L. Lewis - Aircraft Control and Simulation
• NASA Technical Reports
• AIAA Journal Papers on Flight Dynamics

📝 总结

✅ 结语

本项目实现了从飞机建模到 PID 控制再到参数自动优化的完整流程，适用于教学演示、科研仿真、课程设计等场景。所有代码均为 .m 文件，无需编译即可运行。