Learning Objectives

After completing this lab you will be able to:

• derive the equations of motion for longitudinal drone flight dynamics
• formulate the explicit first order ordinary differential equations for the system
• translate ordinary differential equations into a computer function that evaluates the equations at any given point in time
• develop a function that evaluates state dependent input functions
• develop a function that evaluates state and input dependent output variables
• numerically integrate ordinary differential equations with Octave/Matlab's ode45
• create complete and legible plots of the resulting input, state, and output trajectories
• create a report with textual explanations and plots of the simulation
• design a PID controller that stabilizes the drone in pitch and tracks a path

Introduction

In this lab, you will investigate the longitudinal (altitude and pitch) flight dynamics of a quadcopter traveling at a constant forward speed. This vehicle relies completely on the independently controlled thrust generated by pairs of rotors at the front and rear of the vehicle to change altitude and its pitch orientation. The vehicle is very much unstable and requires controllers to maintain a desired pitch angle and to control altitude. The following video shows an autonomous quadcopter that can avoid obstacles using a combination of vision detection to determine relative locations of the obstacles and a feedback controller to avoid the obstacles given the computed desired path of motion. You will develop the feedback controller portion of a system like the one shown.

Here is another similar system using a fixed wing autonomous aircraft at higher speeds:

System Description

You will model the longitudinal dynamics of a drone flying at a constant forward speed \(U\). The vehicle can be modeled as a single rigid body with pitch centroidal moment of inertia \(I\) and mass \(m\). Use an inertial coordinate system with coordinate \(x\) to measure the horizontal travel and coordinate \(z\) for the vertical travel. Use the standard of positive \(z\) pointing downwards. The front rotors are separated from the rear rotors by distance \(a + b\) and the mass center lies at distance \(a\) from the front rotor and distance \(b\) from the rear rotor.

The mass center will have a velocity \(U\) in the \(x\) direction which is constant and a variable velocity \(W\) in the vertical direction. The rigid body will have a pitch angle \(\theta\) and a pitch angular velocity \(q\).

The vehicle will have four forces acting on it. The two pairs of rotors produce the thrust force \(F_f\) at the front of the vehicle and the rotors at the rear produce thrust force \(F_r\). The weight of the vehicle \(mg\) is the third force and the fourth force is the drag force acting perpendicular to the quadcopter's body \(F_D\) at the center of pressure (equal to aerodynamic center in this case) which is equidistant from the to rotors. The drag force should be modeled by the equation:

This drag force will not affect motion in the \(x\) direction since \(U\) is specified to be constant but you will need to determine its contribution to resisting vertical motion. Make sure that the force is always in the opposite direction of the component of \(\bar{v} = U \hat{x} + W \hat{z}\) that is perpendicular to the vehicle (same as \(F_D\)) using the signum function, sign() in Matlab/Octave. Neglect any lift that may be produced due to the vehicles angle of attack or body characteristics.

Figure 1: Schematics of the longitudinal quadcopter dynamics. The vectors are all drawn with a common sign convention, positive in the direction of the arrow.

h = gca;  % gets a handle to the most recent axis
set(h, 'YDir', 'reverse');  % reverses the ordinate

Initial Conditions

For the simulations, set the initial conditions all to zero.

See Integrating the Equations for how to set up the initial condition vector. Make sure that your initial conditions are arranged in the same order as your state variables.

Time Steps

Simulate the system for 20 seconds with time steps of 1/20th of a second for all simulations.

Animation

The following function animate_drone.m function will create an animation of your drone flight simulation. You can use this to visualize the simulation and assess the controller's performance.

function animate_drone(ts, ys, p)
% ANIMATE_DRONE - Creates an animation of the drone flight given the output
% trajectories.
%
% Syntax: animate_drone(ts, ys, p)
%
% Inputs:
%   ts - Vector of time values, size 1x1000.
%   ys - Matrix of output values , size 1000x10
%        The outputs are, in order:
%        1: theta, pitch angle [rad]
%        2: z, altitude [m]
%        3: q, pitch rate [rad/s]
%        4: w, vertical velocity [m/s]
%        5: thetac, cumlative pitch error [rad]
%        6: zc, cumlative altitude error [rad]
%        7: Ff, Front rotor thrust [N]
%        8: Fr, Rear rotor thrust [N]
%        9: x, longitudinal position [m]
%        10: zd, desired altitude [m]
%   p - Constant parameter structure that includes at least a and b.

figure()
hold on
box on

extents = 2.0;  % meters

delt = ts(2) - ts(1);

thetas = ys(:, 1);
zs = ys(:, 2);
xs = ys(:, 9);
zds = ys(:, 10);

desired_path = plot(xs(1:2), zds(1:2), 'r');
actual_path = plot(xs(1:2), zs(1:2), 'b');

for i = 2:length(ts)
    theta = thetas(i);
    z = zs(i);
    x = xs(i);
    % update the desired path
    set(desired_path, 'XData', xs(1:i), 'YData', zds(1:i));
    % update the acutal path
    set(actual_path, 'XData', xs(1:i), 'YData', zs(1:i));
    % plot the "wings"
    plot([x - p.b*cos(theta), x + p.a*cos(theta)], ...
         [z + p.b*sin(theta), z - p.a*sin(theta)], 'k', 'linewidth', 2)
    % plot the mass center
    plot(x, z, 'ob', 'MarkerFaceColor', 'b')
    h = gca;
    set(h, 'YDir', 'reverse');
    xlim([x - extents, x + extents])
    ylim([z - extents, z + extents])
    xlabel('x [m]')
    ylabel('z [m]')
    t_str = sprintf('Time = %1.2f s', ts(i));
    title(t_str);
    pause(delt);
end

hold off

end

Deliverables

In your lab report, show your work for creating and evaluating the simulation model. Include any calculations you had to do, for example those for state equations, initial conditions, input equations, time parameters, and any other parameters. Additionally, provide the indicated plots and answer the questions below. Append a copy of your Matlab/Octave code to the end of the report. The report should follow the report template and guidelines.

Submit a report as a single PDF file to Canvas by the due date that addresses the following items:

1. Create a function defined in an m-file that evaluates the right hand side of the ODEs, i.e. evaluates the state derivatives. See Defining the State Derivative Function for an explanation.
2. Create one function defined in an m-file that calculates the control input. See Time Varying Inputs for an explanation.
3. Create a function defined in an m-file that calculates the requested outputs. See Outputs Other Than the States and Outputs Involving State Derivatives for an explanation.
4. Create a script in an m-file that utilizes the above functions to simulate system for the final path tracking simulation. This should setup the constants, integrate the dynamics equations, and plot each state, and output versus time. See Integrating the Equations for an explanation.
5. Derive the equations of motion of the system. Include your derivation in the report and the resulting equations.
6. Develop two PID feedback controllers using manual tuning. You should try to minimize steady state error, oscillations, the time constant, and overshoot in that order of importance for both the altitude tracking and the pitch stabilization. Show time history plots of the step responses for pitch and altitude using your final gain selection.
7. Present a single controlled simulation of the vehicle and explain the behaviors you observe in each of the ten output variables using knowledge and principles you have learned in the class.

Assessment Rubric