四旋翼动力学和仿真翻译（Quadcopter Dynamics and Simulation）

本文翻译自Andrew Gibiansky的同名文章，该文献介绍了四旋翼的动力学模型和Matlab仿真的具体实现，对四旋翼入门非常有好处。原文如下
http://andrew.gibiansky.com/blog/physics/quadcopter-dynamics/

由于Neo已经于2014年对该文章前半部分进行了翻译（从Introduction到Physics），因此本文将从Simulation部分开始翻译。已经翻译部分如下

http://share.weiyun.com/9c1869b64a76fea67047adedc2f2eb14

本文的翻译工作由BUAA飞机总体设计HT11课程小组成员完成，具体分工如下:

孟德超 Page6~8，simulation和Control（PD control之前）
李志宇 page 8~11 PD control
刘倚铭 page11-14 PID control，
唐既尧 关智元 周震 page15 ~18 Automatic PID Tuning & conclusion

---

仿真

既然我们已经建立了描述此系统动力学的完整的方程组，因此我们可以搭建一个仿真环境，从而我们可以在仿真环境中实时观察到我们的变量和控制器的变化。尽管有很多更先进的方法可以选择，我们仍可以通过欧拉方程建立差分方程，进而模拟不同系统状态的发展情况。在Matlab中，这个仿真环境的部分代码如下

% Simulation times, in seconds.
start_time = 0;
end_time = 10;
dt = 0.005;
times = start_time:dt:end_time;

% Number of points in the simulation.
N = numel(times);

% Initial simulation state.
x = [0; 0; 10];
xdot = zeros(3, 1);
theta = zeros(3, 1);

% Simulate some disturbance in the angular velocity.
% The magnitude of the deviation is in radians / second.
deviation = 100;
thetadot = deg2rad(2 * deviation * rand(3, 1) - deviation);

% Step through the simulation, updating the state.
for t = times
    % Take input from our controller.
    i = input(t);

    omega = thetadot2omega(thetadot, theta);

    % Compute linear and angular accelerations.
    a = acceleration(i, theta, xdot, m, g, k, kd);
    omegadot = angular_acceleration(i, omega, I, L, b, k);

    omega = omega + dt * omegadot;
    thetadot = omega2thetadot(omega, theta); 
    theta = theta + dt * thetadot;
    xdot = xdot + dt * a;
    x = x + dt * xdot;
end

我们还需要计算力和力矩的函数

% Compute thrust given current inputs and thrust coefficient.
function T = thrust(inputs, k)
    % Inputs are values for ${\omega_i}^2$
    T = [0; 0; k * sum(inputs)];
end

% Compute torques, given current inputs, length, drag coefficient, and thrust coefficient.
function tau = torques(inputs, L, b, k)
    % Inputs are values for ${\omega_i}^2$
    tau = [
        L * k * (inputs(1) - inputs(3))
        L * k * (inputs(2) - inputs(4))
        b * (inputs(1) - inputs(2) + inputs(3) - inputs(4))
    ];
end

function a = acceleration(inputs, angles, xdot, m, g, k, kd)
    gravity = [0; 0; -g];
    R = rotation(angles);
    T = R * thrust(inputs, k);
    Fd = -kd * xdot;
    a = gravity + 1 / m * T + Fd;
end

function omegadot = angular_acceleration(inputs, omega, I, L, b, k)
    tau = torques(inputs, L, b, k);
    omegaddot = inv(I) * (tau - cross(omega, I * omega));
end

现在我们还需要的是一些物理常数，一个计算旋转矩阵的函数，以及一个把$\omega$转化为滚转，俯仰和偏航（roll,pitch,yaw）角的微分的函数以及它的反函数。这些函数没有写明。我们可以绘制一个三维坐标下的可视化的四旋翼并且让它在仿真运行时不停更新自己的位置。
注意：虽然原文没有写明该部分的实现，但是考虑到matlab仿真并非易事，此处给出一个译者实现的源代码链接：https://github.com/silverbulletmdc/quadcopter_simulator_in_matlab

四旋翼仿真。每个螺旋桨上的长条粗略的表示其相对升力大小。

控制

推导四旋翼的数学模型是为了协助我们为真实的四旋翼开发一个控制器。系统的输入包括四个电机的角速度，因为我们可以通过控制电压输出值直接控制角速度。注意在我们的简化模型里，我们仅仅采用角速度的平方${\omega_i}^2$而不使用角速度本身。为了计数简单，我们引入输入变量$\gamma_i={\omega_i}^2$。这是因为既然我们可以直接控置${\omega_i}$，那我们也可以直接控制$\gamma_i$。这样，我们就可以将我们的系统描述为状态空间中的一阶微分方程。令
$x_1$为四旋翼在空间中的位置，$x_2$为四旋翼的线速度，$x_3$为四旋翼的滚转、俯仰、偏航角，$x_4$为角速度。（注意这些变量都是三维向量）。有了这些变量，我们可以写出状态空间方程如下：

x1˙x2˙x3˙x4˙=x2=⎡⎣⎢00−g⎤⎦⎥+1mRTB+1mFD=⎡⎣⎢1000cϕ−sϕ−sθcθsϕcθcϕ⎤⎦⎥−1x4=⎡⎣⎢⎢τϕIxx−1τθIyy−1τψIzz−1⎤⎦⎥⎥−⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢Iyy−IzzIxxωyωz