基于Simulink与FlightGear联合建模并仿真多旋翼无人机在平衡态的动态控制

基于Simulink与FlightGear联合建模并仿真多旋翼无人机在平衡态的动态控制

本次实验中主要关注的是对无人机的动态控制进行建模并仿真，即探究计算主板输出PWM波到整个无人机系统，该系统输出变量的变化（例如：位置、速度、位姿、角速度等）。为什么要选择对平衡态的无人机进行建模呢？因为当俯仰角和滚转角过大时，计算主板对俯仰角和滚转角的控制模型从线性控制器变为非线性控制，控制难度增大，解析计算难度增大。当然并不是无法对无人机的高速高机动建模，比如SE（3）微分平坦算法就可以解决该类问题。但非线性控制不是本文主题（我太菜，所以只能在平衡态建模）。

贴个图说明一下当姿态角幅度过大时，线性控制器已经无法良好跟踪信号。

一、建模

这部分内容参考了北航可靠飞行组的实验代码，我从头梳理一遍之后，自己按照物理模型一步一步画simulink框图。这样自己走一遍代码更有体会。

贴一下资料网站

代码与实验课件：
http://rfly.buaa.edu.cn/course.html#ZhCourse

控制平台FlightGear与PX4工具包：https://rflysim.com/zh/3_Using/ControllerDesignAndSimulationPlatform.html#id2

首先我们打开动力系统模拟网站http://flyeval.com/paper/，选择无人机的元器件型号，随后计算生成该系统的物理参数

根据悬停时力平衡的关系
$$mg=4c_T{\omega }^2$$

可以解出悬停时电机的转速，我们将它作为电机的初始转速。

将这些物理参数作为初始化值，写成Matlab脚本:

clear;

%电机相关参数
% 油门——电机稳定转速关系式：wss = Cr * σ + wb
MotorInit_RPM = 557.1420;   % 给定电机初始转速
Motor_Cr = 737.24;          % 油门——电机稳定转速曲线斜率
Motor_wb = 161.87;          % 油门——电机稳定转速曲线截距
Motor_Tm = 0.0895;          % 电机惯性时间常数
Motor_Jm = 2.13e-4;         % 电机螺旋桨转动惯量

%螺旋桨相关参数
Rotor_CT = 1.051e-5;         % 螺旋桨推力系数
Rotor_CM = 1.789e-7;         % 螺旋桨力矩系数 

%机体气动参数
uav_Cd =  0.228;             % 阻力系数
uav_mass = 1.5;              % 多旋翼总体质量
uav_R = 0.225;               % 多旋翼机架半径
uav_Jxx = 1.663e-2;          % x轴转动惯量（单位： kg.m^2）
uav_Jyy = 1.663e-2 ;         % y轴转动惯量（单位： kg.m^2）
uav_Jzz = 3.033e-2;          % z轴转动惯量（单位： kg.m^2）
uav_J = [uav_Jxx,0,0;0,uav_Jyy,0;0,0,uav_Jzz];
uav_GravityAcc = 9.8;        % 机体在重力场所受加速度

%计算悬停时所需油门
%公式: mg = 4*CT*ωss^2
%     ωss = CR*σ + ωb
THR = (sqrt(uav_mass*uav_GravityAcc/4/Rotor_CT)-Motor_wb)/Motor_Cr;

得到机体自身的物理参数后，我们打开Simulink环境，开始搭建动态控制模型。按照如下的控制流程框图，我更偏向于按照从底层到顶层的顺序来搭建模型。

首先我们先搭建动力单元模型。主板输出的PWM波范围为[1000,2000]，我们需要提前对PWM波进行归一化处理，转换成油门值。

接下来我们对动力单元进行建模，根据动力单元公式：

$${\omega }_{ss}=C_r\sigma +{\omega }_b$$

$${\omega }_i=\frac{1}{T_{m}s+1}{\omega }_{ss}$$

先对第一条公式进行建模，输入为油门 输出为目标稳定转速

我们可以很明显的发现第二个公式实际上是一个一阶动态系统，为了控制的方便，我们将其化为状态空间方程式的形式。

$$x\prime =Ax+Bu$$

$$y=Cx+D$$

设ωi作为输出变量，Tmω为状态变量，套进表达形式里推一推可以得到

$$x=T_m{\omega }_{ss}$$

$$x\prime =\frac{T_m}{T_{m}s+1}{\omega }_{ss}\prime =\frac{T_{m}s}{T_{m}s+1}{\omega }_{ss}$$

T m s T m s + 1 ω s s = ( A T m T m s + 1 + B ) ω s s \frac{T_ms}{T_ms+1} \omega_{ss} =( \frac{AT_m}{T_ms+1} +B)\omega_{ss}