简介本资源是一套面向自动化、电子信息与应用数学等专业高年级本科生及飞行控制初学者的MATLAB教学仿真平台聚焦小型旋翼无人机六自由度非线性动力学建模与动态响应分析。资源提供R2014a、R2019b、R2024b三版本兼容环境涵盖113个文件以34个参数配置.mat文件、22个核心算法.m脚本如mavDynamics.m、FandM.m、8个Simulink仿真模型.slx及17个.xml工程配置文件为主总大小仅368KB轻量易部署。已有75人学习下载适用于控制系统设计、飞行器建模课程实验及运动规划算法验证场景。用户可直接运行runfg_small_UAV.bat启动仿真通过修改质量特性、惯性张量与气动系数等配置参数快速开展构型对比与动态特性扫描分析所有关键模块均含中文注释模型严格基于牛顿-欧拉方程构建融合气动、推进与重力效应具备工程级建模精度与教学级可读性。基于Matlab的小型无人机6自由度非线性建模与仿真分析以前总觉得无人机建模仿真是件加分项的事真机飞得好就行。直到前段时间帮朋友调一台自制四旋翼姿态控制参数在悬停时明明很稳一做大角度翻转机动姿态就直接发散光靠真机试飞反复试错烧掉好几组桨我才意识到一个完整的6自由度非线性模型在Matlab里有多重要。后来我把整个飞行器的运动方程原原本本搬进Simulink用非线性模型做故障复现和控制参数预调问题才真正定位清楚。这篇东西就把我这次基于Matlab的小型无人机6自由度非线性建模与仿真分析过程完整梳理一遍从坐标系定义讲到方程推导再到Simulink实现和模型校核适合正在做飞控算法验证、无人机课程设计或者想入门飞行器建模的读者参考。1. 非线性建模为什么绕不开从一次悬停实验说起1.1 线性化模型的适用边界到底在哪大多数飞控教材在讲控制器设计时都会先做一步小扰动线性化假设无人机悬停在某个平衡点附近姿态角和角速度变化都很小于是把含有三角函数、速度叉乘项的非线性方程直接用泰勒展开近似成线性方程。这个思路在理论分析阶段完全没问题线性化模型最大的价值在于能直接用特征值、根轨迹、频域响应这些工具来分析系统的稳定性和控制器参数设计思路清晰运算成本低。但问题藏在小扰动这三个字里面。所谓的小究竟是多小从我的工程经验看当滚转角或俯仰角的偏差超过15度、角速度超过0.5 rad/s、线速度超过5 m/s时线性化模型和真实系统之间的误差就会迅速拉大。如果仿真任务只是验证悬停点的控制器稳定性线性模型是够用的可一旦涉及快速机动、抗风分析、轨迹跟踪这类大幅运动场景线性模型的假设就撑不住了。更关键的是很多控制器的参数在悬停点整定得再好在大姿态角下也会失效因为动力学本身的状态依赖关系已经变了。我朋友那台四旋翼就是这个典型情况悬停状态调试得很稳PID参数放上去的姿态阶跃响应曲线看起来也漂亮但真机做翻滚动作时角速度一上去机身突然剧烈抖动甚至翻转。用回放日志里的电机转速和姿态数据去驱动线性模型发现模型输出和真机姿态完全是两回事。这时候才意识到不是控制算法的问题是模型根本没法描述大机动工况下的真实动力学。1.2 非线性模型里到底装进了哪些非线性那么所谓的非线性模型到底比线性模型多了哪些东西其实本质上就是没有省略运动方程里的完整物理项。对一个6自由度刚体飞行器来说非线性来源主要集中在四个方面。第一是姿态运动学里的三角函数。欧拉角速率和机体坐标系下的角速度之间有一个变换矩阵里面含有正弦、余弦和正切项这些函数天然是高度非线性的。姿态角一大变换关系就不是简单的比例关系了。第二是动力学方程里的惯性叉乘项。角动量在机体坐标系下求导时会出现矢量叉乘形式的耦合项比如角速度向量与惯量矩阵的乘积再叉乘角速度向量。这项描述的是旋转坐标系带来的表观力飞行器在高速滚转时这项的作用非常明显也是大机动时线性化模型最容易丢掉的部分。第三是执行机构的非线性。无刷电机加螺旋桨产生的推力理想情况下近似和转速的平方成正比转速响应本身也有延迟。这个平方关系在油门大范围变化时绝不是一个小增量线性关系能够拟合的。第四是气动力的非线性。低速飞行时可以近似简化但速度稍微提高气动阻力和速度的平方项就占了主导线性阻力的误差可以大到几倍。理解这些非线性来源最大的好处是能帮你在仿真发散或结果异常时快速判断到底是模型哪里没建对还是算法本身有问题。2. 建模前的物理基础坐标系、欧拉角与刚体运动学2.1 NED坐标系与机体坐标系的选择逻辑在正式开始写方程之前坐标系约定必须先定清楚这一步如果含糊后面所有推导都是一笔糊涂账。我用的惯性系是NED也就是北东地坐标系X轴指向正北Y轴指向正东Z轴垂直指向地面。机体坐标系的原点取在飞行器质心上X轴指向机头方向Y轴指向右翼Z轴按右手定则指向机腹方向。为什么选NED而不是大多数人习惯的ENU东-北-天坐标系原因其实很实际地面站、飞控、惯性测量单元输出的数据流方向都是按NED约定的姿态解算出来的四元数和欧拉角也遵循这个约定。如果在仿真里用ENU导航数据的符号和方向就得在代码里来回翻转非常容易搞错。相信我真实项目中因为坐标系方向写反导致仿真的水平位置方向和飞控日志反着的案例比想象中多得多。姿态角用Z-Y-X欧拉顺序定义先绕惯性系Z轴偏航得到偏航角ψ再绕新的Y轴俯仰得到俯仰角θ最后绕新的X轴滚转得到滚转角φ。这个顺序是航空航天里的通用约定三个旋转复合起来就构成了从机体坐标系到惯性坐标系的旋转矩阵。旋转矩阵R的具体形式是教科书级的公式但它本质上的作用可以理解成一个翻译机把机体坐标系下表示的向量比如机体速度、推力方向翻译到惯性坐标系下。翻译的规则取决于飞行器当前的头朝哪、机身倾了多少。如果姿态角很小旋转矩阵可以近似成一个单位矩阵加一个小反对称阵这就是线性化模型的坐标关系而姿态角一大三个三角函数互相耦合只有完整的旋转矩阵才能准确描述。2.2 运动学方程位置和姿态的积分关系运动学方程只关心几何关系不关心力是怎么产生的。平动运动学很简单惯性系下的位置变化率等于机体速度向量通过旋转矩阵变换到惯性系下的结果。写成方程就是[ \begin{bmatrix} \dot{x} \ \dot{y} \ \dot{z} \end{bmatrix} R \begin{bmatrix} u \ v \ w \end{bmatrix} ]这里u、v、w是机体坐标系下的三个线速度分量。这个方程几乎没有争议只是旋转矩阵的代入问题。姿态运动学稍微复杂一点。很多人容易在这里犯迷糊机体坐标系下的角速度p、q、r和欧拉角速率ψ̇、θ̇、φ̇之间并不是简单的相等关系。原因在于欧拉角本身是绕不同旋转轴定义的而这套旋转轴是逐次转出来的不是机体坐标系的固定轴。把它们投影到机体坐标系上就得到了一个含三角函数的变换矩阵[ \begin{bmatrix} p \ q \ r \end{bmatrix} \begin{bmatrix} 1 0 -\sin\theta \ 0 \cos\phi \sin\phi\cos\theta \ 0 -\sin\phi \cos\phi\cos\theta \end{bmatrix} \begin{bmatrix} \dot{\phi} \ \dot{\theta} \ \dot{\psi} \end{bmatrix} ]反解欧拉角速率时这个矩阵的逆矩阵里会出现(\tan\theta)项。当俯仰角θ接近±90度时正切项趋向无穷大数值上就会爆炸——这就是欧拉角表示的奇异性问题。做悬停和普通飞行仿真时θ一般不会接近90度所以用欧拉角没问题但如果你要对倒飞、筋斗这类机动做仿真就一定要换成四元数表示姿态否则仿真会在接近奇异点时直接崩溃。2.3 动力学方程力、力矩与加速度的关系动力学方程处理的是为什么物体会动的问题。用的是牛顿-欧拉方程但要注意方程必须写在机体坐标系下才最方便因为转动惯量、推力、力矩都是相对于机体轴给出的。平动动力学方程在机体坐标系下写为[ m \dot{\vec{v}} \vec{\omega} \times (m \vec{v}) \vec{F}_{ext} m R^T \vec{g}_0 ]其中(\vec{\omega})是机体角速度向量(\vec{F}_{ext})是作用在机身上的合力推力、阻力等(\vec{g}_0)是惯性系下的重力加速度向量在NED系里是([0,0,g]^T)乘上旋转矩阵转置把它变换到机体坐标系里。这里那个叉乘项(\vec{\omega} \times (m \vec{v}))就是惯性坐标系求导和机体坐标系求导之间的转换项。原理可以这样理解惯性系下动量对时间求导等于机体系下的变化率加上坐标旋转带来的表观变化。如果飞行器在高速旋转哪怕机体系下的速度不变惯性系下动量方向其实也在变这一项正是描述这个现象的。很多简化模型会把这一项丢掉低速飞行时可能感觉不到但大角速度机动时它是不可忽略的。角运动方程也一样[ I \dot{\vec{\omega}} \vec{\omega} \times (I \vec{\omega}) \vec{M} ]这里I是3x3转动惯量矩阵在机体坐标系下可以近似为对角阵。(\vec{\omega} \times (I \vec{\omega}))是角动量转移项描述的是转动坐标系下的角动量变化。一般四旋翼的Ixx和Iyy比较接近Izz略大这一项的耦合效应在滚转俯仰同时发生时特别明显。把运动学方程和动力学方程合起来六个自由度对应12个状态变量位置x、y、z姿态角φ、θ、ψ线速度u、v、w角速度p、q、r。后续所有仿真实现都围绕这12个状态展开。3. 完整的六自由度非线性方程组与参数来源3.1 作用在机身上的力与力矩清单要把方程写完整关键是把所有作用在机身上的力和力矩列全。很多初版模型要么少算了一项要么把方向搞反结果仿真出来的行为非常诡异。完整的受力清单如下重力。大小是质量乘以重力加速度方向在NED系下恒为Z轴正方向向下。在机体坐标系下表示时需要把重力向量乘上旋转矩阵的转置所以姿态角会直接影响重力在机体轴上的分量。这个投影关系看似简单实际上是大姿态角模型的核心特征之一。旋翼推力。每个旋翼产生的推力沿机体Z轴负方向也就是向上方向大小与电机转速的平方近似成正比(F_i k_F \omega_i^2)其中kF是推力系数。四个旋翼合力之和就是[ F_{thrust} k_F (\omega_1^2 \omega_2^2 \omega_3^2 \omega_4^2) ]反扭矩。每个旋翼转动时空气会给机身一个反方向力矩大小为(M_i k_M \omega_i^2)方向与旋翼转向相反。四旋翼相邻电机转向相反正是为了在悬停时让反扭矩相互抵消。这个力矩差是偏航控制的核心。陀螺力矩。螺旋桨本身是一个高速旋转的转子当机身姿态变化时转子角动量方向跟着改变会产生陀螺效应力矩。数学上是(M_{gyro} J_p \cdot \vec{\omega}_b \times \vec{\Omega})其中Jp是单个旋翼转子的转动惯量(\vec{\Omega})是电机转子的转速矢量沿机体Z轴。这个力矩在快速滚转俯仰时会体现为姿态交叉耦合做精准控制仿真时必须包含。气动阻力。完整的多旋翼气动阻力建模非常复杂工程上常采用简化模型低速时用线性阻尼(F_d -k_d \vec{v})中高速时用二次阻力(F_d -\frac{1}{2}\rho C_d A |v| \vec{v})。在小型无人机仿真里线性阻尼项通常足够描述低速悬停附近的特性但如果做高速前飞或抗风仿真就得用二次项才比较贴合实际。3.2 参数测定与估算一个典型小型四旋翼的数值模型参数对不对直接决定仿真结果可不可信。参数不能拍脑袋随意估至少要有一个来源依据。以我这次搭建的轴距450mm小型四旋翼为例关键参数如下表参数符号数值来源机身质量m1.5 kg电子秤直接称重重力加速度g9.81 m/s²常数机臂长度l0.25 m游标卡尺测量转动惯量X轴Ixx0.021 kg·m²三线摆实测转动惯量Y轴Iyy0.024 kg·m²三线摆实测转动惯量Z轴Izz0.037 kg·m²三线摆实测推力系数kF1.63e-6 N/(rad/s)²静拉力台标定反扭矩系数kM3.6e-8 N·m/(rad/s)²悬停转速反推转子转动惯量Jp1.2e-5 kg·m²厂家参数估算线性阻尼系数kd0.03 N·s/m真机日志拟合这里转动惯量我强烈建议实测不要只靠SolidWorks估算。空机CAD模型算出来的转动惯量往往偏低因为电调、电池扎带、减震球这类装配件在3D模型里不够精确而且电池的实际布局位置对惯量影响很大。三线摆的测量方法其实不复杂做一个圆形托盘用三根等长细线悬挂起来把无人机放在托盘上让系统做微小扭转振动测周期就可以算出总惯量再减去托盘自身的惯量就得到无人机惯量。这个方法精度足够工程使用成本也低。推力系数需要在静拉力台上标定。把电机和桨正装在一个力传感器上给不同油门信号测量对应推力。转速可以用转速计或者电调遥测数据读出来然后拟合(F k_F \omega^2)。反扭矩系数更难直接测通常的做法是让四旋翼在真机悬停此时偏航力矩平衡通过飞控日志读出四个电机的转速差反推出kM的大致范围。还有一个更省事的办法用厂家提供的螺旋桨推力扭矩系数数据表来查。没有条件实测时至少也要保证参数的量级合理否则仿真悬停需要的转速和真机差太多那模型就没有意义了。3.3 完整状态方程的组装所有力和力矩都齐了之后就可以组装成12维非线性状态方程了。为了后续在Matlab里写代码方便我习惯用一个结构体数组来组织状态向量[ x [x,y,z,\phi,\theta,\psi,u,v,w,p,q,r]^T ]然后每个状态的变化率写成显式方程。这里全部写出来位置方程 [ \dot{x} (\cos\theta\cos\psi)u (\sin\phi\sin\theta\cos\psi - \cos\phi\sin\psi)v (\cos\phi\sin\theta\cos\psi \sin\phi\sin\psi)w ] [ \dot{y} (\cos\theta\sin\psi)u (\sin\phi\sin\theta\sin\psi \cos\phi\cos\psi)v (\cos\phi\sin\theta\sin\psi - \sin\phi\cos\psi)w ] [ \dot{z} (-\sin\theta)u (\sin\phi\cos\theta)v (\cos\phi\cos\theta)w ]欧拉角速率方程把姿态运动学矩阵求逆 [ \dot{\phi} p q\sin\phi\tan\theta r\cos\phi\tan\theta ] [ \dot{\theta} q\cos\phi - r\sin\phi ] [ \dot{\psi} q\frac{\sin\phi}{\cos\theta} r\frac{\cos\phi}{\cos\theta} ]线速度方程 [ \dot{u} rv - qw - g\sin\theta ] [ \dot{v} pw - ru g\sin\phi\cos\theta ] [ \dot{w} qu - pv g\cos\phi\cos\theta - \frac{k_F}{m}(\omega_1^2\omega_2^2\omega_3^2\omega_4^2) ]角速度方程 [ \dot{p} \frac{I_{yy}-I_{zz}}{I_{xx}}qr - \frac{J_p}{I_{xx}}q\Omega_{sum} \frac{l k_F}{I_{xx}}(\omega_2^2 - \omega_4^2) ] [ \dot{q} \frac{I_{zz}-I_{xx}}{I_{yy}}pr \frac{J_p}{I_{yy}}p\Omega_{sum} \frac{l k_F}{I_{yy}}(\omega_1^2 - \omega_3^2) ] [ \dot{r} \frac{I_{xx}-I_{yy}}{I_{zz}}pq \frac{k_M}{I_{zz}}(\omega_1^2 - \omega_2^2 \omega_3^2 - \omega_4^2) ]其中(\Omega_{sum} \omega_1 - \omega_2 \omega_3 - \omega_4)注意这里是带符号的转速和和电机转向有关。这个方程组的推导过程虽然长了点但每一步都有明确的物理来源后面在Matlab里实现时就变成了一个几乎可以翻译成代码的公式清单。4. Simulink模型搭建从方程到仿真模块4.1 总体架构设计为什么我要拆成三个子系统在Simulink里搭模型我强烈建议不要把所有方程塞进一个巨大的框图里而是按物理逻辑拆成三个子系统力与力矩计算模块、刚体动力学模块、运动学更新模块。之所以这样拆最直接的好处是复用性强。力与力矩模块的输入是四个电机转速命令输出是合力和合力矩刚体动力学模块的输入是合力和合力矩输出是机体线加速度和角加速度运动学更新模块的输入是线速度和角速度输出是位置和姿态角。这样拆完之后每个模块都可以单独测试。更实际的好处是后面要接入控制器时你只需要在力与力矩模块前加一个控制信号接口完全不用改动动力学核心要做故障注入比如某个电机失效也只需要在力与力矩模块里单独处理对应电机的出力就行。顶层接口设计如下模型最外层输入是四个电机转速命令(\omega_1)到(\omega_4)输出是12维状态向量也可以直接把位置和姿态角单独引出来给Scope或者Record模块。第二个好处是方便调试。仿真结果不对时可以先把每个子系统的中间输出拉出来看。比如机体角速度是否符合物理直觉推力计算有没有在某个转速区间出现明显突变这些中间层信号的可见性能帮你省下大量排查时间。4.2 MATLAB Function实现核心方程实现方式上我选择用MATLAB Function模块写状态方程而不是用一堆积分器、增益、求和模块去搭框图。原因很简单对12维非线性方程来说框图的连线复杂度太高容易连错且不易检查而MATLAB Function里的代码几乎可以直接对照第三节的公式逐行翻译可读性和可维护性都强得多。下面这段代码就是核心状态方程的完整实现。我把它写在MATLAB Function模块里输入是时间t实际未用、状态向量x和参数结构体params输出是状态导数x_dot。function x_dot quadrotor_dynamics(x, params) % 状态向量解包 % x [x,y,z,phi,theta,psi,u,v,w,p,q,r] phi x(4); theta x(5); psi x(6); u x(7); v x(8); w x(9); p x(10); q x(11); r x(12); % 电机转速输入由控制器或参数配置给出 w1 params.w1; w2 params.w2; w3 params.w3; w4 params.w4; % 三角函数预计算避免重复计算 sp sin(phi); cp cos(phi); st sin(theta); ct cos(theta); ss sin(psi); cs cos(psi); % 旋转矩阵R机体坐标系-NED惯性系 R11 ct*cs; R12 sp*st*cs - cp*ss; R13 cp*st*cs sp*ss; R21 ct*ss; R22 sp*st*ss cp*cs; R23 cp*st*ss - sp*cs; R31 -st; R32 sp*ct; R33 cp*ct; % 位置导数 x_dot zeros(12,1); x_dot(1) R11*u R12*v R13*w; x_dot(2) R21*u R22*v R23*w; x_dot(3) R31*u R32*v R33*w; % 姿态角导数欧拉角速率方程 if abs(ct) 1e-4 % 接近俯仰±90度欧拉角奇异这里做保护 x_dot(4) 0; x_dot(5) 0; x_dot(6) 0; else x_dot(4) p q*sp*tan(theta) r*cp*tan(theta); x_dot(5) q*cp - r*sp; x_dot(6) (q*sp r*cp)/ct; end % 合力与合力矩 Fz -params.kF * (w1^2 w2^2 w3^2 w4^2); % 机体Z轴负方向推力 tau_phi params.l * params.kF * (w2^2 - w4^2); tau_theta params.l * params.kF * (w1^2 - w3^2); tau_psi params.kM * (w1^2 - w2^2 w3^2 - w4^2); % 陀螺力矩Omega_sum为带符号转速和 Omega_sum w1 - w2 w3 - w4; gyro_p -params.Jp * q * Omega_sum / params.Ixx; gyro_q params.Jp * p * Omega_sum / params.Iyy; % 线速度导数 x_dot(7) r*v - q*w - params.g * st; x_dot(8) p*w - r*u params.g * sp * ct; x_dot(9) q*u - p*v params.g * cp * ct Fz/params.m; % 角速度导数 x_dot(10) ((params.Iyy - params.Izz)/params.Ixx)*q*r gyro_p tau_phi/params.Ixx; x_dot(11) ((params.Izz - params.Ixx)/params.Iyy)*p*r gyro_q tau_theta/params.Iyy; x_dot(12) ((params.Ixx - params.Iyy)/params.Izz)*p*q tau_psi/params.Izz; end这段代码里有两个细节值得特别说明。第一三角函数项我全部预先算好再复用避免在每个方程里重复计算虽然对单个仿真影响不大但做蒙特卡洛参数扫描时能明显省时间。第二姿态角导数里加了俯仰角接近90度时的保护实际项目中欧拉角奇异是仿真崩溃最常见的原因之一宁可在这里做保护也不能让积分器直接算出一个NaN。线速度方程里的Fz是四个推力求和方向沿机体Z轴负方向也就是向上。很多新手容易在这里把符号搞反结果仿真一开始无人机就往下掉怎么调参数都救不回来。4.3 初始化脚本、配平初值与求解器设置Simulink模型里的参数不应该硬编码在模块里而是通过一个初始化脚本统一管理。我习惯用一个.m文件定义所有的参数并打包成结构体% quad_initialize.m params.m 1.5; % kg params.g 9.81; % m/s^2 params.l 0.25; % m params.Ixx 0.021; % kg*m^2 params.Iyy 0.024; % kg*m^2 params.Izz 0.037; % kg*m^2 params.kF 1.63e-6; % N/(rad/s)^2 params.kM 3.6e-8; % N*m/(rad/s)^2 params.Jp 1.2e-5; % kg*m^2 params.kd 0.03; % N*s/m % 悬停配平转速总推力重力 % m*g kF * 4 * w_hover^2 w_hover sqrt(params.m * params.g / (4 * params.kF)); params.w1 w_hover; params.w2 w_hover; params.w3 w_hover; params.w4 w_hover;配平初值这一步非常关键。直接给一个所有状态都为0的初值电机转速也给一个随意值那仿真一开始的加速度就会非常大飞机会迅速坠落曲线很难看也看不出真实动力学特性。正确的做法是先把悬停转速算出来让初始时刻合力矩为零、推力等于重力然后在这个配平点附近做扰动实验。这样模型在时间原点附近的行为才可控便于分析。求解器设置方面如果只是做一般性的开环动力学仿真固定步长ode4四阶龙格库塔加1毫秒步长就够了。四旋翼的刚体动力学带宽一般在10Hz到50Hz这个量级1ms的采样足够覆盖。但如果后面把电机动态响应也建模进去电机转速环的响应时间可能只有几十毫秒这时候建议把步长缩到0.1ms或0.2ms否则高频动态会失真。做硬件在环仿真时步长还必须要和飞控的控制频率对齐一般250Hz或500Hz这个要单独确认。5. 仿真分析与模型校核模型建得对不对5.1 开环悬停响应什么样的结果算模型正确模型搭完第一步要做的不是急着接控制器而是先跑开环悬停扰动响应。具体做法是用配平初值给一个微小初始扰动比如让滚转角初始值设为2度然后观察姿态角的响应曲线。这里有个反直觉的经验一个没接控制器的四旋翼开环模型姿态响应应该是发散的或者等幅振荡的而不是稳定的。原因很简单四旋翼本身是一个静不稳定的系统没有姿态反馈任何扰动都会让它继续偏离。如果你发现开环悬停仿真里姿态角自己慢慢回到了0那大概率模型里有问题可能是陀螺力矩符号搞反了或者某个阻尼项加得过大掩盖了真实动力学特性。以我搭建的模型为例给2度滚转角初始偏差仿真10秒滚转角会呈振荡发散趋势同时偏航角也会因为反扭矩耦合出现漂移。这才是健康模型的输出。如果曲线一开始就冲出去变成NaN那通常是数值问题优先检查姿态角导数里的奇异保护有没有生效以及积分器步长是否合适。5.2 非线性模型与线性化模型的对比实验模型通过开环合理性检查之后接下来做线性化对比实验。这一步主要是为了验证非线性模型的正确性同时也能量化线性化模型的误差边界。在Simulink里可以把模型在悬停配平点做线性化分析得到状态矩阵A再画出线性模型的阶跃响应与非线性模型在同样输入下的响应做对比。小扰动情况下比如滚转阶跃指令只有5度两个模型的响应曲线几乎重合把阶跃指令加大到30度非线性模型和线性模型的峰值角速度、超调量就开始有肉眼可见的差异指令到60度时两个模型的响应行为就完全不同了线性化模型可能预测出一个稳定收敛的响应而非线性模型则表现出明显的耦合与发散。这个对比实验的价值在于第一它验证了非线性模型在平衡点附近和线性模型的一致性说明方程推导没有大问题第二它给出了线性化模型的适用范围为后续控制器设计提供参考边界。做论文或者技术报告时这组对比数据也是非常有力的素材。5.3 与真实飞行日志校核参数修正的闭环纯仿真层面的验证还不够模型最终要能用必须和真机数据校核。方法不复杂把一次真实试飞的飞控日志导出来取出电机转速指令、姿态角、角速度这些信号。然后把电机转速作为模型输入在Matlab里重新跑一遍仿真对比模型输出的姿态角和真机记录的姿态角。第一次对比通常会有偏差原因就在气动阻尼系数、推力系数这些参数和真实值有差距。这里有一个很实用的调参思路先观察悬停段的误差如果模型姿态振荡频率比真机高说明转动惯量偏小了如果模型响应比真机迟钝说明惯量偏大或者阻尼系数给大了。然后做一段机动段对比重点看大姿态角下误差走向如果误差在大姿态时快速增大优先检查陀螺力矩项和气动阻力模型是否需要从线性升级成二次阻尼。这个过程本质上是一个以模型输出和真机输出之间的误差最小化为目标的参数辨识不需要什么复杂算法手动根据试飞数据逐项修正几轮精度就能达到工程可用的水平。以我的经验经过三四轮修正后模型姿态输出和真机日志的平均误差能控制在5度以内这足够支撑控制器参数预调了。5.4 仿真发散原因排查清单以我踩过的坑和帮别人排查的经验把最常见的问题列成一张表方便仿真异常时快速定位现象可能原因排查方法仿真几秒内出现NaN单位混用度/弧度或初值不符合物理范围检查所有角度是否已转为弧度检查欧拉角奇异保护是否生效无人机一开始就猛坠推力系数错误或悬停转速未配平核对悬停转速公式检查推力方程里符号方向姿态振荡发散速度快得不正常转动惯量过小或陀螺力矩缺失核对三线摆实测值检查角速度方程里Jp项是否加入悬停时姿态角无扰动却自己漂移重心与机体坐标系原点不重合检查Z轴力矩是否为零电机转向和安装方向是否一致小扰动响应和线性化模型完全一致正常现象说明模型基本正确继续做大角度扰动对比大角度机动时仿真崩溃欧拉角奇异或步长过大改用四元数姿态表示或缩小固定步长到0.1ms这里面最容易被忽略的是单位问题。我见过一个同学在Simulink里把角度初值设置成了30运行几秒姿态就完全乱掉排查到最后发现其实他把30度直接当成30弧度用了。这在仿真里是灾难性的错误Matlab本身不会提示因为数值上积分器照样能跑只是结果完全失真。排查的方法论其实就一句话不要先怀疑求解器和积分器先相信你的模型方程和参数。大多数发散问题都出在建模阶段——符号方向、单位、参数量级这三类错误占据了九成以上的排查时间。搞完模型校核之后这个非线性6自由度模型就可以干很多事了接PID控制器做姿态轨迹仿真、做电机失效故障注入、做参数敏感性分析、甚至配合硬件在环平台做半实物仿真。我个人的体会是模型的价值不在于复杂本身而在于你遇到问题时能快速定位到底是模型假设错误还是控制器设计错误。从那次帮朋友调试四旋翼之后我基本养成了习惯任何飞行器项目起步先在Matlab里搭完整非线性模型哪怕真机还没装配完成模型也能把很多问题提前暴露出来。这套流程对一个刚开始做无人机建模的人而言一开始可能会觉得方程推导枯燥、参数标定繁琐但咬牙把模型搭完并且和真机数据校核通过之后后续的控制算法开发会很轻松。本文还有配套的精品资源点击获取
