
简介本资源是一套面向控制工程与无人机方向初学者及课程设计者的Simulink建模仿真实践材料聚焦四旋翼无人机的运动学与动力学建模、PD串级控制器设计及闭环轨迹跟踪验证。资源包含51个文件以36个MAT数据文件存储仿真参数与结果、2个SLX主模型文件分别对应简易建模与复杂建模双版本、7个L文件可能为S-Function或自定义模块、3个M脚本用于初始化或后处理为核心辅以PNG示意图与少量ZIP嵌套包整体压缩包仅647KB轻量易部署。已有2025人学习下载适合高校自动化/航空航天专业学生开展课程实验、毕设建模或控制算法入门实践。用户可直接运行两个层级的Simulink模型基础版支持手动PID参数调试与稳定性验证进阶版实现70秒内从起飞、三目标点连续跟踪到返航的完整轨迹闭环附带清晰的控制器分层结构与角速度推导逻辑便于理解建模原理与控制策略落地细节。1. 四旋翼无人机 Simulink 建模不是“搭积木”而是对刚体动力学边界的精确刻画很多初学者打开 Simulink 后直接拖 PID 模块、Scope 和 Step 信号源以为“连上线就能飞”——结果仿真跑起来姿态角发散、位置漂移、甚至出现 NaN 报错。这不是模型没调好而是建模起点错了四旋翼本质是强耦合、非线性、欠驱动的六自由度刚体系统其运动由四个旋翼推力与力矩的微分关系决定。本文提供的两套 Simulink 模型恰恰对应工程实践中两种不可替代的建模范式简易建模用于控制律快速验证与参数敏感性分析复杂建模用于闭环轨迹跟踪性能评估与物理约束嵌入。前者用代数方程一阶惯性近似替代真实动力学后者严格基于牛顿-欧拉方程推导出 12 维状态空间3 位置 3 速度 3 姿态角 3 角速度并显式引入电机响应延迟、气流扰动、陀螺仪零偏等真实效应。适合控制算法工程师做控制器选型也适合飞控系统集成人员做硬件在环HIL前的全链路验证。2. 从牛顿-欧拉方程出发构建四旋翼动力学核心模块2.1 为什么必须用牛顿-欧拉而非拉格朗日——刚体建模的坐标系选择逻辑四旋翼属于典型的空间刚体运动其质心平动与绕质心转动存在强耦合。若采用拉格朗日方法需定义广义坐标并求解复杂的动能/势能表达式且难以直观嵌入电机推力分配矩阵而牛顿-欧拉法直接在惯性系I与机体坐标系B双框架下分别列写平动与转动方程物理意义清晰便于 Simulink 中用坐标变换模块如Rotation Angles to Direction Cosine Matrix实现帧间映射。本项目复杂建模中状态变量定义为x [px; py; pz; vx; vy; vz; phi; theta; psi; p; q; r]; % px/py/pz: 惯性系下位置 (m) % vx/vy/vz: 惯性系下线速度 (m/s) % phi/theta/psi: 滚转/俯仰/偏航角 (rad) % p/q/r: 机体坐标系下角速度 (rad/s)提示Simulink 中所有角度单位必须统一为弧度rad否则sin/cos计算将严重失真。若输入为度数务必通过deg2rad()转换切勿依赖 Scope 自动标尺。2.2 推力与力矩生成从电机 PWM 到机体总力的映射链四旋翼的控制输入是四个电机的转速或等效 PWM 占空比输出是作用于质心的合力与合力矩。该映射由两步构成2.2.1 电机模型平方关系与饱和限制单个电机产生的推力近似满足 $F_i k_f \omega_i^2$其中 $k_f$ 为推力系数N·s²/rad²。Simulink 中用Math Function模块设为square接Gain模块实现但必须叠加Saturation模块限制 $\omega_i$ 范围如 0–600 rad/s否则高速下推力爆炸式增长导致数值不稳定。2.2.2 推力分配矩阵解耦刚体六自由度的关键四个旋翼推力 $[F_1,F_2,F_3,F_4]^T$ 通过固定几何构型X 型或型合成总力 $[F_x,F_y,F_z]^T$ 与总力矩 $[\tau_\phi,\tau_\theta,\tau_\psi]^T$。以 X 型为例分配矩阵为$$ \begin{bmatrix} F_x \ F_y \ F_z \ \tau_\phi \ \tau_\theta \ \tau_\psi \end{bmatrix}\begin{bmatrix} 0 0 0 0 \ 0 0 0 0 \ 1 1 1 1 \ 0 l 0 -l \ -l 0 l 0 \ -b b -b b \end{bmatrix} \begin{bmatrix} F_1 \ F_2 \ F_3 \ F_4 \end{bmatrix} $$其中 $l$ 为旋翼到质心距离m$b$ 为扭矩系数N·m·s²/rad²。该矩阵在 Simulink 中用Matrix Multiply模块实现输入为 4×1 向量输出为 6×1 向量。注意矩阵元素符号必须与实际旋翼旋转方向严格匹配顺时针/逆时针否则偏航方向相反。参数符号典型值物理含义Simulink 设置位置推力系数$k_f$2.5e-5电机推力与转速平方的比例系数Gain模块增益值扭矩系数$b$1.0e-6电机反扭矩与推力的比例系数分配矩阵最后一行系数臂长$l$0.25旋翼中心到无人机质心距离分配矩阵第4、5行系数最大转速$\omega_{max}$600电机物理极限转速Saturation上限2.3 坐标系转换欧拉角奇点规避与方向余弦矩阵实现姿态更新需将机体坐标系下的角速度 $[p,q,r]$ 积分得到欧拉角 $[\phi,\theta,\psi]$。直接积分会导致万向节锁死$\theta \pm \pi/2$ 时 $\dot{\psi}$ 无定义。本项目复杂模型采用方向余弦矩阵DCM 四元数更新双路径设计主路径用Quaternion Dynamics模块来自 Aerospace Blockset解算四元数 $q[q_0,q_1,q_2,q_3]$再经Quaternion to Rotation Angles转为欧拉角备用路径当 $\theta$ 接近 $\pm85^\circ$ 时自动切换至 DCM 积分Direction Cosine Matrix模块避免奇点。% 在 MATLAB Function 模块中实现四元数微分方程 function dq fcn(p, q) % q [q0;q1;q2;q3], p [p;q;r] Omega [0, -p, -q, -r; ... p, 0, r, -q; ... q, -r, 0, p; ... r, q, -p, 0]; dq 0.5 * Omega * q; end该函数封装在MATLAB Function模块中输入为角速度向量pqr输出为四元数导数dq再接入Integrator模块完成积分。关键参数积分器初始值设为[1;0;0;0]零姿态四元数采样时间设为0.0011 kHz过低会导致姿态抖动过高则积分误差累积。3. PD 串级控制器设计与 Simulink 实现细节3.1 为何选用 PD 而非 PID——四旋翼控制中的微分先行与抗积分饱和四旋翼位置控制存在显著滞后姿态变化 → 机体倾斜 → 水平加速度 → 位置变化形成二级惯性环节。若在位置环加入积分项易因传感器噪声或模型误差导致积分饱和引发持续振荡。本项目采用位置环 PD 姿态环 PD 的串级结构其优势在于位置环输出为期望姿态角$\phi_d,\theta_d$不直接输出力避免过载姿态环输出为期望力矩响应快且无稳态误差需求微分项提供阻尼抑制高频振荡且可通过Derivative模块的滤波时间常数如0.01s抑制噪声放大。3.2 位置控制器从参考轨迹到期望姿态的映射位置控制器接收三维参考轨迹 $ x_r,y_r,z_r $ 与当前状态 $[x,y,z,v_x,v_y,v_z]$输出期望滚转/俯仰角$$ \phi_d \frac{1}{g} \left( \ddot{x}r k{px}(x_r - x) k_{dx}(\dot{x}r - v_x) \right), \quad \theta_d \frac{1}{g} \left( \ddot{y}r k{py}(y_r - y) k{dy}(\dot{y}_r - v_y) \right) $$其中 $g9.81$ m/s²。Simulink 中用Second-Order Filter模块生成 $\ddot{x}_r,\ddot{y}_r$避免数值微分噪声Gain模块实现比例/微分增益Sum模块完成误差计算。典型增益配置适用于 1 kg 无人机环节参数数值效果说明位置环比例$k_{px},k_{py}$2.5提升位置跟踪响应速度过高导致超调位置环微分$k_{dx},k_{dy}$1.8抑制位置振荡过大会放大噪声高度环比例$k_{pz}$4.0高度控制需更强刚度因重力主导高度环微分$k_{dz}$2.2防止高度突变时电机功率骤增3.3 姿态控制器角速度闭环与执行器饱和处理姿态控制器接收期望姿态角 $[\phi_d,\theta_d,\psi_d]$ 与当前姿态 $[\phi,\theta,\psi]$先经PID Controller模块仅启用 P/D生成期望角速度 $[p_d,q_d,r_d]$再与当前角速度 $[p,q,r]$ 构成内环输出期望力矩$$ \tau_\phi k_{p\phi}(\phi_d - \phi) k_{d\phi}(p_d - p), \quad \tau_\theta k_{p\theta}(\theta_d - \theta) k_{d\theta}(q_d - q), \quad \tau_\psi k_{p\psi}(\psi_d - \psi) k_{d\psi}(r_d - r) $$注意PID Controller模块中Derivative gain必须设为非零值如0.15否则角速度环开环姿态无法稳定。同时在力矩输出端添加Saturation模块上下限 ±2.0 N·m防止电机过载烧毁。3.4 控制器参数整定Ziegler-Nichols 法在 Simulink 中的实操步骤关闭微分项将所有k_d设为 0仅保留比例增益逐步增大 $k_p$从 0.1 开始每次增加 0.5运行 10 秒仿真观察姿态角响应记录临界振荡增益 $k_u$当姿态角出现等幅振荡时记下此时 $k_p$ 值如 $k_u3.2$计算 PD 参数按 Z-N 公式 $k_p 0.45k_u$, $k_d 0.15k_u \cdot T_u$其中 $T_u$ 为振荡周期如 1.2 s验证与微调将计算值填入模型运行轨迹跟踪仿真若超调 15%则 $k_p$ 减 10%若调节时间 3 s则 $k_d$ 增 20%。4. 简易建模与复杂建模的仿真对比验证方法4.1 简易模型用代数方程替代微分方程的适用边界简易模型将动力学简化为位置动力学$\ddot{x} g \tan\phi$, $\ddot{y} g \tan\theta$, $\ddot{z} u_z - g$姿态动力学$\ddot{\phi} u_\phi$, $\ddot{\theta} u_\theta$, $\ddot{\psi} u_\psi$其中 $u_z,u_\phi,u_\theta,u_\psi$ 为控制器直接输出。该模型在 Simulink 中仅需Integrator×6 Trigonometric Function模块计算开销极低。但仅适用于小角度$|\phi|,|\theta|15^\circ$与低速$v2$ m/s场景。验证方法在Scope中同时显示简易模型与复杂模型的 $z$ 轴位置当高度阶跃响应超调量相差 5% 且调节时间差 0.3 s 时可判定简易模型在此工况下有效。4.2 复杂模型轨迹跟踪验证三段式参考路径的构造与评估复杂模型验证采用三段连续路径上升段0–20 s$z 0.5t$$xy0$水平移动段20–45 s$x 2\sin(0.1t), y 2\cos(0.1t), z 10$返回段45–70 s沿直线回归原点 $(0,0,0)$。在 Simulink 中用Repeating Sequence模块生成分段函数Mux模块合并三轴参考信号。评估指标需导出至 MATLAB 工作区% 在仿真配置中启用 Log data to workspace变量名设为 simout t simout.time; pos_err sqrt((simout.signals.values(:,1)-x_ref).^2 ... (simout.signals.values(:,2)-y_ref).^2 ... (simout.signals.values(:,3)-z_ref).^2); fprintf(最大位置误差: %.3f m\n, max(pos_err)); fprintf(平均位置误差: %.3f m\n, mean(pos_err));合格标准最大误差 0.3 m平均误差 0.15 m。若不达标优先检查推力分配矩阵符号、DCM 更新频率应 ≥500 Hz、以及Solver设置推荐ode4变步长相对误差1e-4。4.3 仿真发散的三大根源与定位命令当仿真出现NaN或剧烈发散时按以下顺序排查现象定位命令根本原因解决方案初始时刻即发散sim(model_name, TraceExecution, on)状态初值不满足静平衡条件将初始姿态角设为[0;0;0]初始角速度设为[0;0;0]初始推力设为mg/4运行中突然发散set_param(model_name,Solver,Fixed-step); set_param(model_name,FixedStep,0.0001)变步长求解器在刚性系统中步长过大改用ode1Euler或ode3Bogacki-Shampine步长固定为1e-4姿态角跳变plot(t, simout.signals.values(:,7:9)); grid on欧拉角奇点未处理替换Euler Angle积分为Quaternion积分或启用Direction Cosine Matrix模块5. 从 Simulink 模型到可部署代码C 代码生成的关键配置项5.1 模型配置参数确保生成代码符合实时性要求在Model Configuration Parameters中必须设置Solver→Type:Fixed-step,Solver:ode3Bogacki-ShampineFixed step size:0.001Code Generation→System target file:ert.tlcEmbedded CoderTarget hardware vendor:GenericOptimization→Default parameter behavior:InlinedSignal storage reuse:onCustom Code→Include directories: 添加$(MATLAB_ROOT)/extern/includeLibraries:libeng.lib libmx.lib。提示若目标平台为 STM32需在Hardware Implementation中指定Device vendor:STMicroelectronicsDevice type:STM32F407VG并勾选Enable support for processor-in-the-loop (PIL)。5.2 关键模块替换避免不支持代码生成的模块以下模块在生成 C 代码时会报错必须替换Simulink 模块替换方案替换后模块说明Derivative用一阶低通滤波器替代Transfer Fcn$1/(0.01s1)$避免数值微分噪声放大MATLAB Function重写为Stateflow或Atomic SubsystemStateflow Chart支持代码生成的确定性逻辑Scope删除或替换为To WorkspaceTo Workspace变量名logsout仅用于离线分析不参与实时控制5.3 生成代码并验证三步验证法确保功能等价编译验证运行slbuild(model_name)检查model_name_ert_rtw文件夹是否生成model_name.c/h数值一致性验证在 MATLAB 中运行sim(model_name)导出yout再用生成的model_name_initialize()model_name_step()函数在 C 环境中运行相同步数对比输出数组实时性验证在目标板上运行生成代码用示波器测量model_name_step()执行时间必须 ≤1 ms对应 1 kHz 控制频率。若超时需启用Inline parameters并关闭Array bounds checking。最终生成的model_name.h中控制输入接口为extern void model_name_initialize(void); extern void model_name_step(double x_ref, double y_ref, double z_ref, double phi_ref, double theta_ref, double psi_ref, double* F1, double* F2, double* F3, double* F4);该函数每调用一次即完成一个控制周期的计算输出四个电机推力值N可直接映射至 PWM 占空比。本文还有配套的精品资源点击获取