A-A+

# 四旋翼动力学模型、仿真与控制（教程及代码）

2016年01月10日 阅读 1 次

# Introduction

A helicopter is a flying vehicle which uses rapidly spinning rotors to push air downwards, thus creating a thrust force keeping the helicopter aloft. Conventional helicopters have two rotors. These can be arranged as two coplanar rotors both providing upwards thrust, but spinning in opposite directions (in order to balance the torques exerted upon the body of the helicopter). The two rotors can also be arranged with one main rotor providing thrust and a smaller side rotor oriented laterally and counteracting the torque produced by the main rotor. However, these configurations require complicated machinery to control the direction of motion; a swashplate is used to change the angle of attack on the main rotors. In order to produce a torque the angle of attack is modulated by the location of each rotor in each stroke, such that more thrust is produced on one side of the rotor plane than the other. The complicated design of the rotor and swashplate mechanism presents some problems, increasing construction costs and design complexity.

A quadrotor helicopter (quadcopter) is a helicopter which has four equally spaced rotors, usually arranged at the corners of a square body. With four independent rotors, the need for a swashplate mechanism is alleviated. The swashplate mechanism was needed to allow the helicopter to utilize more degrees of freedom, but the same level of control can be obtained by adding two more rotors.

The development of quadcopters has stalled until very recently, because controlling four independent rotors has proven to be incredibly difficult and impossible without electronic assistance. The decreasing cost of modern microprocessors has made electronic and even completely autonomous control of quadcopters feasible for commercial, military, and even hobbyist purposes.

Quadcopter control is a fundamentally difficult and interesting problem. With six degrees of freedom (three translational and three rotational) and only four independent inputs (rotor speeds), quadcopters are severely underactuated. In order to achieve six degrees of freedom, rotational and translational motion are coupled. The resulting dynamics are highly nonlinear, especially after accounting for the complicated aerodynamic effects. Finally, unlike ground vehicles, helicopters have very little friction to prevent their motion, so they must provide their own damping in order to stop moving and remain stable. Together, these factors create a very interesting control problem. We will present a very simplified model of quadcopter dynamics and design controllers for our dynamics to follow a designated trajectory. We will then test our controllers with a numerical simulation of a quadcopter in flight.

We will start deriving quadcopter dynamics by introducing the two frames in which will operate. The inertial frame is defined by the ground, with gravity pointing in the negative zz direction. The body frame is defined by the orientation of the quadcopter, with the rotor axes pointing in the positive zz direction and the arms pointing in the xx and yy directions.

Quadcopter Body Frame and Inertial Frame

## Kinematics

Before delving into the physics of quadcopter motion, let us formalize the kinematics in the body and inertial frames. We define the position and velocity of the quadcopter in the inertial frame as x=(x,y,z)Tx=(x,y,z)T and x˙=(x˙,y˙,z˙)Tx˙=(x˙,y˙,z˙)T, respectively. Similarly, we define the roll, pitch, and yaw angles in the body frame as θ=(ϕ,θ,ψ)Tθ=(ϕ,θ,ψ)T, with corresponding angular velocities equal to θ˙=(ϕ˙,θ˙,ψ˙)Tθ˙=(ϕ˙,θ˙,ψ˙)T. However, note that the angular velocity vector ωθ˙ω≠θ˙. The angular velocity is a vector pointing along the axis of rotation, while θ˙θ˙ is just the time derivative of yaw, pitch, and roll. In order to convert these angular velocities into the angular velocity vector, we can use the following relation:

ω=⎡⎣⎢1000cϕsϕsθcθsϕcθcϕ⎤⎦⎥θ˙ω=[10−sθ0cϕcθsϕ0−sϕcθcϕ]θ˙

where ωω is the angular velocity vector in the body frame.We can relate the body and inertial frame by a rotation matrix RR which goes from the body frame to the inertial frame. This matrix is derived by using the ZYZ Euler angle conventions and successively “undoing” the yaw, pitch, and roll.

R=⎡⎣⎢cϕcψcθsϕsψcθcψsϕ+cϕsψsϕsθcψsϕcϕcθsψcϕcθcψsϕsψcϕsθsθsψcψsθcθ⎤⎦⎥R=

For a given vector v⃗ v→ in the body frame, the corresponding vector is given by Rv⃗ Rv→ in the inertial frame.

## Physics

In order to properly model the dynamics of the system, we need an understanding of the physical properties that govern it. We will begin with a description of the motors being used for our quadcopter, and then use energy considerations to derive the forces and thrusts that these motors produce on the entire quadcopter. All motors on the quadcopter are identical, so we can analyze a single one without loss of generality. Note that adjacent propellers, however, are oriented opposite each other; if a propeller is spinning “clockwise”, then the two adjacent ones will be spinning “counter-clockwise”, so that torques are balanced if all propellers are spinning at the same rate.

### Motors

Brushless motors are used for all quadcopter applications. For our electric motors, the torque produced is given by

τ=Kt(II0)τ=Kt(I−I0)

where ττ is the motor torque, II is the input current, I0I0 is the current when there is no load on the motor, and KtKt is the torque proportionality constant. The voltage across the motor is the sum of the back-EMF and some resistive loss:

V=IRm+KvωV=IRm+Kvω

where VV is the voltage drop across the motor, RmRm is the motor resistance, ωω is the angular velocity of the motor, and KvKv is a proportionality constant (indicating back-EMF generated per RPM). We can use this description of our motor to calculate the power it consumes. The power is

P=IV=(τ+KtI0)(KtI0Rm+τRm+KtKvω)Kt2P=IV=(τ+KtI0)(KtI0Rm+τRm+KtKvω)Kt2

For the purposes of our simple model, we will assume a negligible motor resistance. Then, the power becomes proportional to the angular velocity:

P(τ+KtI0)KvωKtP≈(τ+KtI0)KvωKt

Further simplifying our model, we assume that KtI0τKtI0≪τ. This is not altogether unreasonable, since I0I0 is the current when there is no load, and is thus rather small. In practice, this approximation holds well enough. Thus, we obtain our final, simplified equation for power:

PKvKtτω.P≈KvKtτω.

### Forces

The power is used to keep the quadcopter aloft. By conservation of energy, we know that the energy the motor expends in a given time period is equal to the force generated on the propeller times the distance that the air it displaces moves (Pdt=FdxP⋅dt=F⋅dx). Equivalently, the power is equal to the thrust times the air velocity (P=FdxdtP=Fdxdt).

P=TvhP=Tvh

We assume vehicle speeds are low, so vhvh is the air velocity when hovering. We also assume that the free stream velocity, vv∞, is zero (the air in the surrounding environment is stationary relative to the quadcopter). Momentum theory gives us the equation for hover velocity as a function of thrust,

vh=T2ρA−−−−√vh=T2ρA

where ρρ is the density of the surrounding air and AA is the area swept out by the rotor. Using our simplified equation for power, we can then write

P=KvKtτω=KvKτKtTω=T322ρA−−−−√.P=KvKtτω=KvKτKtTω=T322ρA.

Note that in the general case, τ=r⃗ ×F⃗ τ=r→×F→; in this case, the torque is proportional to the thrust TT by some constant ratio Kτ determined by the blade configuration and parameters. Solving for the thrust magnitude TT, we obtain that thrust is proportional to the square of angular velocity of the motor:

T=(KvKτ2ρA−−−−√Ktω)2=kω2T=(KvKτ2ρAKtω)2=kω2

where kk is some appropriately dimensioned constant. Summing over all the motors, we find that the total thrust on the quadcopter (in the body frame) is given by

TB=i=14Ti=k⎡⎣⎢00ωi2⎤⎦⎥.TB=∑i=14Ti=k[00∑ωi2].

In addition to the thrust force, we will model friction as a force proportional to the linear velocity in each direction. This is a highly simplified view of fluid friction, but will be sufficient for our modeling and simulation. Our global drag forces will be modeled by an additional force term

FD=⎡⎣⎢kdx˙kdy˙kdz˙⎤⎦⎥FD=[−kdx˙−kdy˙−kdz˙]

If additional precision is desired, the constant kdkd can be separated into three separate friction constants, one for each direction of motion. If we were to do this, we would want to model friction in the body frame rather than the inertial frame.

### Torques

Now that we have computed the forces on the quadcopter, we would also like to compute the torques. Each rotor contributes some torque about the body zz axis. This torque is the torque required to keep the propeller spinning and providing thrust; it creates the instantaneous angular acceleration and overcomes the frictional drag forces. The drag equation from fluid dynamics gives us the frictional force:

FD=12ρCDAv2.FD=12ρCDAv2.

where ρρ is the surrounding fluid density, AA is the reference area (propeller cross-section, not area swept out by the propeller), and CDCD is a dimensionless constant. This, while only accurate in some in some cases, is good enough for our purposes. This implies that the torque due to drag is given by

τD=12RρCDAv2=12RρCDA(ωR)2=bω2τD=12RρCDAv2=12RρCDA(ωR)2=bω2

where ωω is the angular velocity of the propeller, RR is the radius of the propeller, and bb is some appropriately dimensioned constant. Note that we’ve assumed that all the force is applied at the tip of the propeller, which is certainly inaccurate; however, the only result that matters for our purposes is that the drag torque is proportional to the square of the angular velocity. We can then write the complete torque about the zz axis for the iith motor:

τz=bω2+IMω˙τz=bω2+IMω˙

where IMIM is the moment of inertia about the motor zz axis, ω˙ω˙ is the angular acceleration of the propeller, and bb is our drag coefficient. Note that in steady state flight (i.e. not takeoff or landing), ω˙0ω˙≈0, since most of the time the propellers will be maintaining a constant (or almost constant) thrust and won’t be accelerating. Thus, we ignore this term, simplifying the entire expression to

τz=(1)i+1bωi2.τz=(−1)i+1bωi2.

where the (1)i+1(−1)i+1 term is positive for the iith propeller if the propeller is spinning clockwise and negative if it is spinning counterclockwise. The total torque about the zz axis is given by the sum of all the torques from each propeller:

τψ=b(ω12ω22+ω32ω42)τψ=b(ω12−ω22+ω32−ω42)

The roll and pitch torques are derived from standard mechanics. We can arbitrarily choose the i=1i=1 and i=3i=3 motors to be on the roll axis, so

τϕ=r×T=L(kω12kω32)=Lk(ω12ω32)τϕ=∑r×T=L(kω12−kω32)=Lk(ω12−ω32)

Correspondingly, the pitch torque is given by a similar expression

τθ=Lk(ω22ω42)τθ=Lk(ω22−ω42)

where LL is the distance from the center of the quadcopter to any of the propellers. All together, we find that the torques in the body frame are

τB=⎡⎣⎢⎢Lk(ω12ω32)Lk(ω22ω42)b(ω12ω22+ω32ω42)⎤⎦⎥⎥τB=[Lk(ω12−ω32)Lk(ω22−ω42)b(ω12−ω22+ω32−ω42)]

The model we’ve derived so far is highly simplified. We ignore a multitude of advanced effects that contribute to the highly nonlinear dynamics of a quadcopter. We ignore rotational drag forces (our rotational velocities are relatively low), blade flapping (deformation of propeller blades due to high velocities and flexible materials), surrounding fluid velocities (wind), etc. With that said, we now have all the parts necessary to write out the dynamics of our quadcopter.

### Equations of Motion

In the inertial frame, the acceleration of the quadcopter is due to thrust, gravity, and linear friction. We can obtain the thrust vector in the inertial frame by using our rotation matrix RR to map the thrust vector from the body frame to the inertial frame. Thus, the linear motion can be summarized as

mx¨=⎡⎣⎢00mg⎤⎦⎥+RTB+FDmx¨=[00−mg]+RTB+FD

where x⃗ x→ is the position of the quadcopter, gg is the acceleration due to gravity, FDFD is the drag force, and TBTB is the thrust vector in the body frame. While it is convenient to have the linear equations of motion in the inertial frame, the rotational equations of motion are useful to us in the body frame, so that we can express rotations about the center of the quadcopter instead of about our inertial center. We derive the rotational equations of motion from Euler’s equations for rigid body dynamics. Expressed in vector form, Euler’s equations are written as

Iω˙+ω×(Iω)=τIω˙+ω×(Iω)=τ

where ωω is the angular velocity vector, II is the inertia matrix, and ττ is a vector of external torques. We can rewrite this as

ω˙=⎡⎣⎢ω˙xω˙yω˙z⎤⎦⎥=I1(τω×(Iω)).ω˙=[ω˙xω˙yω˙z]=I−1(τ−ω×(Iω)).

We can model our quadcopter as two thin uniform rods crossed at the origin with a point mass (motor) at the end of each. With this in mind, it’s clear that the symmetries result in a diagonal inertia matrix of the form

I=⎡⎣⎢Ixx000Iyy000Izz⎤⎦⎥.I=[Ixx000Iyy000Izz].

Therefore, we obtain our final result for the body frame rotational equations of motion

ω˙=⎡⎣⎢⎢τϕIxx1τθIyy1τψIzz1⎤⎦⎥⎥⎡⎣⎢⎢⎢⎢IyyIzzIxxωyωzIzzIxxIyyωxωzIxxIyyIzzωxωy⎤⎦⎥⎥⎥⎥ω˙=[τϕIxx−1τθIyy−1τψIzz−1]−[Iyy−IzzIxxωyωzIzz−IxxIyyωxωzIxx−IyyIzzωxωy]

# Simulation

Now that we have complete equations of motion describing the dynamics of the system, we can create a simulation environment in which to test and view the results of various inputs and controllers. Although more advanced methods are available, we can quickly write a simulator which utilizes Euler’s method for solving differential equations to evolve the system state. In MATLAB, this simulator would be written as follows.

% 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;

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

% 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;
theta = theta + dt * thetadot;
xdot = xdot + dt * a;
x = x + dt * xdot;
end

We would then need functions to compute all of the physical forces and torques.

% 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

We would also need values for all of our physical constants, a function to compute the rotation matrix RR, and functions to convert from an angular velocity vector ωω to the derivatives of roll, pitch, and yaw and vice-versa. These are not shown. We can then draw the quadcopter in a three-dimensional visualization which is updated as the simulation is running.

Quadcopter Simulation. Bars above each propeller represent, roughly, relative thrust magnitudes.

# Control

The purpose of deriving a mathematical model of a quadcopter is to assist in developing controllers for physical quadcopters. The inputs to our system consist of the angular velocities of each rotor, since all we can control is the voltages across the motors. Note that in our simplified model, we only use the square of the angular velocities, ωi2ωi2, and never the angular velocity itself, ωiωi. For notational simplicity, let us introduce the inputs γi=ωi2γi=ωi2. Since we can set ωiωi, we can clearly set γiγi as well. With this, we can write our system as a first order differential equation in state space. Let x1x1 be the position in space of the quadcopter, x2x2 be the quadcopter linear velocity, x3x3 be the roll, pitch, and yaw angles, and x4x4 be the angular velocity vector. (Note that all of these are 3-vectors.) With these being our state, we can write the state space equations for the evolution of our state.

x1˙x2˙x3˙x4˙=x2=⎡⎣⎢00g⎤⎦⎥+1mRTB+1mFD=⎡⎣⎢1000cϕsϕsθcθsϕcθcϕ⎤⎦⎥1x4=⎡⎣⎢⎢τϕIxx1τθIyy1τψIzz1⎤⎦⎥⎥⎡⎣⎢⎢⎢⎢IyyIzzIxxωyωzIzzIxxIyyωxωzIxxIyyIzzωxωy⎤⎦⎥⎥⎥⎥x1˙=x2x2˙=[00−g]+1mRTB+1mFDx3˙=[10−sθ0cϕcθsϕ0−sϕcθcϕ]−1x4x4˙=[τϕIxx−1τθIyy−1τψIzz−1]−[Iyy−IzzIxxωyωzIzz−IxxIyyωxωzIxx−IyyIzzωxωy]

Note that our inputs are not used in these equations directly. However, as we will see, we can choose values for ττ and TT, and then solve for values of γiγi.

## PD Control

In order to control the quadcopter, we will use a PD control, with a component proportional to the error between our desired trajectory and the observed trajectory, and a component proportional to the derivative of the error. Our quadcopter will only have a gyro, so we will only be able to use the angle derivatives ϕ˙ϕ˙, θ˙θ˙, and ψ˙ψ˙ in our controller; these measured values will give us the derivative of our error, and their integral will provide us with the actual error. We would like to stabilize the helicopter in a horizontal position, so our desired velocities and angles will all be zero. Torques are related to our angular velocities by τ=Iθ¨τ=Iθ¨, so we would like to set the torques proportional to the output of our controller, with τ=Iu(t)τ=Iu(t). Thus,

⎡⎣⎢τϕτθτψ⎤⎦⎥=⎡⎣⎢⎢⎢⎢⎢Ixx(Kdϕ˙+KpT0ϕ˙dt)Iyy(Kdθ˙+KpT0θ˙dt)Izz(Kdψ˙+KpT0ψ˙dt)⎤⎦⎥⎥⎥⎥⎥[τϕτθτψ]=[−Ixx(Kdϕ˙+Kp∫0Tϕ˙dt)−Iyy(Kdθ˙+Kp∫0Tθ˙dt)−Izz(Kdψ˙+Kp∫0Tψ˙dt)]

We have previously derived the relationship between torque and our inputs, so we know that

τB=⎡⎣⎢Lk(γ1γ3)Lk(γ2γ4)b(γ1γ2+γ3γ4)⎤⎦⎥=⎡⎣⎢⎢⎢⎢⎢Ixx(Kdϕ˙+KpT0ϕ˙dt)Iyy(Kdθ˙+KpT0θ˙dt)Izz(Kdψ˙+KpT0ψ˙dt)⎤⎦⎥⎥⎥⎥⎥τB=[Lk(γ1−γ3)Lk(γ2−γ4)b(γ1−γ2+γ3−γ4)]=[−Ixx(Kdϕ˙+Kp∫0Tϕ˙dt)−Iyy(Kdθ˙+Kp∫0Tθ˙dt)−Izz(Kdψ˙+Kp∫0Tψ˙dt)]

This gives us a set of three equations with four unknowns. We can constrain this by enforcing the constraint that our inputs must keep the quadcopter aloft:

T=mg.T=mg.

Note that this equation ignores the fact that the thrust will not be pointed directly up. This will limit the applicability of our controller, but should not cause major problems for small deviations from stability. If we had a way of determining the current angle accurately, we could compensate. If our gyro is precise enough, we can integrate the values obtained from the gyro to get the angles θθ and ϕϕ. In this case, we can calculate the thrust necessary to keep the quadcopter aloft by projecting the thrust mgmg onto the inertial zz axis. We find that

Tproj=mgcosθcosϕTproj=mgcos⁡θcos⁡ϕ

Therefore, with a precise angle measurement, we can instead enforce the requirement that the thrust be equal to

T=mgcosθcosϕT=mgcos⁡θcos⁡ϕ

in which case the component of the thrust pointing along the positive zz axis will be equal to mgmg. We know that the thrust is proportional to a weighted sum of the inputs:

T=mgcosθcosϕ=kγiγi=mgkcosθcosϕT=mgcos⁡θcos⁡ϕ=k∑γi⟹∑γi=mgkcos⁡θcos⁡ϕ

With this extra constraint, we have a set of four linear equations with four unknowns γiγi. We can then solve for each γiγi, and obtain the following input values:

γ1γ2γ3γ4=mg4kcosθcosϕ2beϕIxx+eψIzzkL4bkL=mg4kcosθcosϕ+eψIzz4beθIyy2kL=mg4kcosθcosϕ2beϕIxx+eψIzzkL4bkL=mg4kcosθcosϕ+eψIzz4b+eθIyy2kLγ1=mg4kcos⁡θcos⁡ϕ−2beϕIxx+eψIzzkL4bkLγ2=mg4kcos⁡θcos⁡ϕ+eψIzz4b−eθIyy2kLγ3=mg4kcos⁡θcos⁡ϕ−−2beϕIxx+eψIzzkL4bkLγ4=mg4kcos⁡θcos⁡ϕ+eψIzz4b+eθIyy2kL

This is a complete specification for our PD controller. We can simulate this controller using our simulation environment. The controller drives the angular velocities and angles to zero.

Left: Angular velocities. Right: angular displacements. ϕϕ, θθ, ψψ are coded as red, green, and blue.

However, note that the angles are not completely driven to zero. The average steady state error (error after 10 seconds of simulation) is approximately 0.3. This is a common problem with using PD controllers for mechanical systems, and can be partially alleviated with a PID controller, as we will discuss in the next section.

In addition, note that since we are only controlling angular velocities, our positions and linear velocities do not converge to zero. However, the zz position will remain constant, because we have constrained the total vertical thrust to be such that it keeps the quadcopter perfectly aloft, without ascending or descending. However, this is really nothing more than a curiosity. With the limited sensing that we have available to us, there is nothing we can do to control the linear position and velocity of the quadcopter. While in theory we could compute the linear velocities and positions from the angular velocities, in practice the values will be so noisy as to be completely useless. Thus, we will restrict ourselves to just stabilizing the quadcopter angle and angular velocity. (Traditionally, navigation is done by a human, and stabilization is there simply to make control easier for the human operator.)

We have implemented this PD control for use in our simulation. The controller is implemented as a function which is given some state (corresponding to controller state, not system state) and the sensor inputs, and must compute the inputs γiγi and the updated state. The code for a PD control follows.

% Compute system inputs and updated state.
% Note that input = [$\gamma_1$, $\ldots$, $\gamma_4$]
function [input, state] = pd_controller(state, thetadot)
% Controller gains, tuned by hand and intuition.
Kd = 4;
Kp = 3;

% Initialize the integral if necessary.
if ~isfield(state, 'integral')
params.integral = zeros(3, 1);
end

% Compute total thrust
total = state.m * state.g / state.k / (cos(state.integral(1)) * cos(state.integral(2)));

% Compute errors
e = Kd * thetadot + Kp * params.integral;

% Solve for the inputs, $\gamma_i$
input = error2inputs(params, accels, total);

% Update the state
params.integral = params.integral + params.dt .* thetadot;
end

## PID Control

PD controllers are nice in their simplicity and ease of implementation, but they are often inadequate for controlling mechanical systems. Especially in the presence of noise and disturbances, PD controllers will often lead to steady state error. A PID control is a PD control with another term added, which is proportional to the integral of the process variable. Adding an integral term causes any remaining steady-state error to build up and enact a change, so a PID controller should be able to track our trajectory (and stabilize the quadcopter) with a significantly smaller steady-state error. The equations remain identical to the ones presented in the PD case, but with an additional term in the error:

eϕeθeψ=Kdϕ˙+KpT0ϕ˙dt+KiT0T0ϕ˙dtdt=Kdθ˙+KpT0θ˙dt+KiT0T0θ˙dtdt=Kdψ˙+KpT0ψ˙dt+KiT0T0ψ˙dtdteϕ=Kdϕ˙+Kp∫0Tϕ˙dt+Ki∫0T∫0Tϕ˙dtdteθ=Kdθ˙+Kp∫0Tθ˙dt+Ki∫0T∫0Tθ˙dtdteψ=Kdψ˙+Kp∫0Tψ˙dt+Ki∫0T∫0Tψ˙dtdt

However, PID controls come with their own shortcomings. One problem that commonly occurs with a PID control is known as integral windup.
In some cases, integral wind-up can cause lengthy oscillations instead of settling. In other cases, wind-up may actually cause the system to become unstable, instead of taking longer to reach a steady state.

If there is a large disturbance in the process variable, this large disturbance is integrated over time, becoming a still larger control signal (due to the integral term). However, even once the system stabilizes, the integral is still large, thus causing the controller to overshoot its target. It may then begin a series of dieing down oscillations, become unstable, or simply take an incredibly long time to reach a steady state. In order to avoid this, we disable the integral function until we reach something close to the steady state. Once we are in a controllable region near the desired steady state, we turn on the integral function, which pushes the system towards a low steady-state error.

With a properly implemented PID, we achieve an error of approximately 0.06 after 10 seconds.

We have implemented this PID control for use in simulation, in the same way as with the PD controller shown earlier. Note that there is an additional parameter to tune in a PID. The disturbances used for all the test cases are identical, shown to compare the controllers.

% Compute system inputs and updated state.
% Note that input = [$\gamma_1$, $\ldots$, $\gamma_4$]
function [input, state] = pid_controller(state, thetadot)
% Controller gains, tuned by hand and intuition.
Kd = 4;
Kp = 3;
Ki = 5.5;

% Initialize the integral if necessary.
if ~isfield(state, 'integral')
params.integral = zeros(3, 1);
params.integral2 = zeros(3, 1);
end

% Prevent wind-up
if max(abs(params.integral2)) > 0.01
params.integral2(:) = 0;
end

% Compute total thrust
total = state.m * state.g / state.k / (cos(state.integral(1)) * cos(state.integral(2)));

% Compute errors
e = Kd * thetadot + Kp * params.integral - Ki * params.integral2;

% Solve for the inputs, $\gamma_i$
input = error2inputs(params, accels, total);

% Update the state
params.integral = params.integral + params.dt .* thetadot;
params.integral2 = params.integral2 + params.dt .* params.integral;
end

## Automatic PID Tuning

Although PID control has the potential to perform very well, it turns out that the quality of the controller is highly dependent on the gain parameters. Tuning the parameters by hand may be quite difficult, as the ratios of the parameters is as important as the magnitudes of the parameters themselves; often, tuning parameters requires detailed knowledge of the system and an understanding of the conditions in which the PID control will be used. The parameters we chose previously were tuned by hand for good performance, simply by running simulations with many possibly disturbances and parameter values, and choosing something that worked reasonably well. This method is clearly suboptimal, not only because it can be very difficult and labor-intensive (and sometimes more or less impossible) but also because the resulting gains are not in any way guaranteed to be optimal or even close to optimal.

Ideally, we would be able to use an algorithm to analyze a system and output the “optimal” PID gains, for some reasonable definition of optimal. This problem has been studied in depth, and many methods have been proposed. Many of these methods require detailed knowledge of the system being modeled, and some rely on properties of the system, such as stability or linearity. The method we will use for choosing our PID parameters is a method known as extremum seeking.

Extremum seeking works exactly as the name implies. We define the “optimal” set of parameters as some vector θ⃗ =(Kp,Ki,Kd)θ→=(Kp,Ki,Kd) which minimizes some cost function J(θ⃗ )J(θ→). In our case, we would like to define a cost function that penalizes high error and error over extended durations of time. One candidate cost function is given by

J(θ⃗ )=1tftotft0e(t,θ⃗ )2dtJ(θ→)=1tf−to∫t0tfe(t,θ→)2dt

where e(t,θ⃗ )e(t,θ→) is the error in following some reference trajectory with some initial disturbance using a set of parameters θ⃗ θ→. Suppose we were able to somehow compute the gradient of this cost function, J(θ⃗ )∇J(θ→). In that case, we could iteratively improve our parameter vector by defining a parameter update rule

θ⃗ (k+1)=θ⃗ (k)αJ(θ⃗ )θ→(k+1)=θ→(k)−α∇J(θ→)

where θ⃗ (k)θ→(k) is the parameter vector after kk iterations and αα is some step size which dictates how much we adjust our parameter vector at each step of the iteration. As kk→∞, the cost function J(θ⃗ )J(θ→) will approach a local minimum in the space of PID parameters.The question remains as to how we can estimate J(θ⃗ )∇J(θ→). By definition,

J(θ⃗ )=(KpJ(θ⃗ ),KiJ(θ⃗ ),KdJ(θ⃗ )).∇J(θ→)=(∂∂KpJ(θ→),∂∂KiJ(θ→),∂∂KdJ(θ→)).

We know how to compute J(θ⃗ )J(θ→). Using this, we can approximate the derivative with respect to any of the gains numerically, simply by computing

KJ(θ⃗ )J(θ⃗ +δu^K)J(θ⃗ δu^K)2δ∂∂KJ(θ→)≈J(θ→+δ⋅u^K)−J(θ→−δ⋅u^K)2δ

where u^Ku^K is the unit vector in the KK direction. As δ0δ→0, this approximation becomes better. Using this approximation, we can minimize our cost function and achieve locally optimal PID parameters. We can start with randomly initialized positive weights, disturb the system in some set manner, evaluate J(θ⃗ )J(θ→) by simulating the system for different PID parameters, and then compute the gradient. Then, using the method of gradient descent, we can iteratively oprtimize our gains until we have some form of convergence.The gradient descent method does, however, have several problems. First of all, although it finds a local minimum, that minimum is only guaranteed to be a local minimum - there may be other minima which are better global minima. In order to avoid choosing suboptimal local minima in the cost function, we repeat our optimization several times, and choose the best result. We initialize our PID parameters randomly, so each time we run the optimization we will get a different result. In addition, instead of choosing disturbance and then optimizing the response to that disturbance, we use several random disturbances at each iteration and use the average response to compute costs and gradients. This ensures that our parameters are general and not optimized for a specific disturbance. In addition, we vary the step size and the number of disturbances to try per iteration, in order to increase the sensitivity of our results as our iteration continues. We stop iterations when we detect a steady state, which we do by computing a linear regression on the most recent costs and iterating until the slope is statistically indistinguishable from zero using a 99% confidence interval.

Using our quadcopter simulation, we can define a function that computes the cost for a given set of PID parameters.

function J = cost(theta)
% Create a controller using the given gains.
control = controller('pid', theta(1), theta(2), theta(3));

% Perform a simulation.
data = simulate(control);

% Compute the integral, $\frac{1}{t_f - t_0} \int_{t_0}^{t_f} e(t)^2 dt$
t0 = 0;
tf = 1;
J = 1/(tf - t0) * sum(data.theta(data.t >= t0 & data.t <= tf) .^ 2) * data.dt;
end

We can use this function to approximate a derivative with respect to a gain:

% Compute derivative with respect to first parameter.
delta = 0.01;
var = [1, 0, 0];
derivative = (cost(theta + delta * var) - cost(theta - delta * var)) / (2 * delta);

We can then use our gradient descent method (with all modifications described above) to minimize the cost function and obtain a good set of PID parameters. We can verify that our tuning method is working by visualizing the cost function versus the iteration number, and seeing that the cost function is indeed going down and stabilizing at a local minimum.

Cost function plotted as a function of iteration number, along with moving average. Tuning stops when the slope of the moving average becomes statistically indistinguishable from zero with a 99% confidence interval.

We can compare the manually-chosen PID parameters with those designed by the algorithm.

Top: Angular velocities and angular displacements, using manually tuned PID controller. Bottom: Angular velocities and angular displacements, using automatically tuned PID controller.

The automatically-chosen PID parameters do significantly better overall. They have significantly smaller swings in value, overshoot significantly less, and converge faster. However, the error in the angular displacement takes longer to converge to zero with the automatically tuned parameters than with the manually turned parameters, although the initial convergence is much better when the parameters are chosen via gradient descent. This is due to the fact that our cost function emphasizes squared error, and thus gives priority to minimizing overall error magnitude rather than long-term convergence. We could easily modify our cost function to give higher priority to long-term error, in which case the automatically-tuned parameters are likely to do better.

# Conclusion

We derived equations of motion for a quadcopter, starting with the voltage-torque relation for the brushless motors and working through the quadcopter kinematics and dynamics. We ignored aerodynamical effects such as blade-flapping and non-zero free stream velocity, but included air friction as a linear drag force in all directions. We used the equations of motion to create a simulator in which to test and visualize quadcopter control mechanisms.

We began with a simple PD controller. Although the PD controller worked, it left a significant steady-state error. In order to decrease the steady-state error, we added an integral term in order to create a PID controller. We tested the PID controller (with minor modifications to prevent integral wind-up) and found that it was better at preventing steady-state error than the PD controller when presented with the same disturbances and using the same proportional and derivative gains. We also found that tuning the PID controller was difficult, and would often lead to an unstable system for unknown reasons. In order to avoid the difficulty of PID tuning and find the optimal set of parameters, we used a gradient-descent based extremum seeking method in order to numerically estimate gradients of a cost function in PID-parameter space and iteratively choose a set of parameters to minimize the cost function. We found that the resulting controller was significantly better than the one using manually turned parameters.

## 2 条留言  访客：2 条  博主：0 条

1. opan08

你这样不标注一下人家文章的出处真的好吗？
原文：

• top

感谢指正！