The equations below are from the last post, where we see we matched our Bouabdallah paper‘s equations for rotational motion.

We know some propeller input is going to be a part of the applied torque. There is going to be a set of equations for linear motion too, but let’s clarify what is going on here regarding torques about our 3 body axes. Let’s look at these equations and…

Review Rigid body rotation gyroscopic effect

Remember Newton’s Second Law ? This tells us that the sum of the forces on a mass in the i’th direction must equal the mass of the body times the acceleration in that direction.

If we take one of the equations above and rewrite it we can recognize this familiar form, albeit in rotational nomenclature. The rightmost term is the gyroscopic effect of the entire rigid body. It tells us that if there is rotational rate about y- and z-axes a torque is induced about the X-axis. If we neglect this term you see we’d be left with   in simple second-law form.

It’s a little confusing, but all we’ve done is already offer one item to the, “sum-of-torques”: the gyroscopic effect resulting from the rigid body rotation. Notice when we moved the term in parenthesis to the other side it looks like an added torque now but we swapped the subtraction in parenthesis so the math is the same.

Gyroscopic effect of Propellers rotation.

“tail-dragger” airplane example

The video below to describes a phenomenon familiar to pilots, particularly pilots of powerful, “tail draggers” like the P-51 Mustang. A propeller is a large spinning disk. When we tip the plane of rotation of the propeller: pitch the airplane, a torque on the body results due to this gyroscopic effect.

Imagine you are in the cockpit of the P-51, speeding down the runway (lucky you). That monster prop has a healthy angular momentum vector pointing perpendicular to the prop. Your tail is dragging until you get enough lift from your wings (due to streamline curvature, if you recall a recent post!). At some point your tail will lift and you’ll be parallel to the ground before your front wheels leave the ground. This looks like a quick, “pitch” of say 15-degrees based on the photo above. Let’s say it took a second for the tail to rise as we’re nearing take-off, so the pitch rate is 15-degrees per-second as the tail wheel leaves the ground.

That same pitch is pitching the plane of the spinning propeller forward, or more to the vertical. When you watch the video below and inspect the equations below you will see that this, “pitch forward” of the propeller at the above angular rate of the airplane going from dragging-to-level results in a gyroscopic torque about the airplane’s vertical axis: you feel the plane want to yaw (twist) so you’ll need to correct with some rudder.

The effect can be negligible

Our propellers are horizontal in the body-plane, but there is a torque induced on the body that results from propeller rotation and rotation of the propeller plane about the x- and y-axis of our system too…at least in theory. It may be negligible depending on propeller moment of inertia relative to other effects, but let’s model it anyway.

We can choose to include it or neglect it later when we perform a sensitivity analysis. This is where we can look at rough magnitude of various effects and decide to include them in our model or cancel the terms in the equations: ignore them as too small to make a difference. We’ll do this systematically and not guess at this point. Besides, this is a neat phenomenon for us to appreciate in detail.

Video Explaining Equation Derivation

Gyroscopic effect of our propellers

The following PDF derivations arrive at the same terms as the Bouabdallah paper (equation 8). There’s a sign convention difference between the paper and these derivations, but the result is the same:

The first equation tells us that there will be a torque about our body X-axis when there is angular rate-of-change to our pitch:  (not the dot, implying pitch rate). This means our quadcopter will feel a  body roll-inducing torque when it pitches. This is due to the gyroscopic effects of the 4 propellers.

The second equation tells us that there will be a torque about our body Y-axis when there is angular rate-of-change to our roll: . This means our quadcopter will feel a  body pitch-inducing torque when it rolls.  Again, this is due to the gyroscopic effects of the 4 propellers.

The hand calculations below arrive at the above equations following the process in the video. The video above explains how to go about the hand calculations for a single propeller.

Hand Calculation PDF Notes

RB 2018-07-11 17.27.01

Conclusion

If you expand these equations you’ll see I just performed the steps in the above video Eight times and combined the equations. We need eight steps because we have 4 propellers and we need to account for the plane-of-rotation rate-of-change about both the x- and y-axis. Rotation about the body z-axis contributes no gyroscopic effect from the propellers because it is just like twisting the axis of the propeller. The plane-of-rotation of the propeller does not change in this case so no torque is induced on the body.

We’ve now accounted for gyroscopic effect of the rigid body rotations and the propeller rotations.

Next we’ll add the propeller thrust components. These will be our inputs to counteract the above effects and unknown disturbances (wind) in order to maintain the attitude of our quadcoptor and move it where we want it to go.

This post explains how we determine propeller thrust and drag factors for our quadcopter project.

The last couple posts have been working-out the sum-of-torques on our quad-copter. We covered gyroscopic effect of the total air-frame in the, “equations of motion” post a few weeks ago.

Next we looked at the torque induced by the gyroscopic effects of the spinning propellers in the last post.

Now we’ll consider the propeller drive inputs.

Propeller Drive Control Input

We started the, “equations of motion” with the sketch below, so we’ll use it again here as we model the command inputs from the propellers: .

ref_frames

Torque About Body X-axis

Assume the length of the arm from the center of the quadcopter to the propeller is, ‘l’ then…

The, ‘F’ terms are thrust forces. We know from an earlier post that this thrust from a propeller is going to be a function of propeller speed squared: :

Where, ‘b’ is a, “thrust factor” we need to determine. After you watch the video below you’ll understand how we get…

where…

our torque about the X-axis equation above becomes…

Torque About Body Y-axis

By the same logic above for the X-axis torque, this Y-axis torque becomes…

Torque About Body Z-axis

The torque about the Z-Axis is proportional to the squared angular velocities of the propellers. We subtract the squared counter-clockwise props from the clockwise props to lump them. The proportionality constant (the multiplier) is going to be a, “drag factor”: d.

See the video for an explanation on how we get ‘d’ from the experimentally determined, “power coefficient” for a candidate propeller.

Deriving the, ‘b’ and, ‘d’ terms above

‘b’ is our thrust factor and, ‘d’ is our drag factor. They are arrived at through experimentally determined thrust and power coefficients for our propeller. We encountered the drag term earlier earlier (here) but we didn’t have insight into how it’s arrived at.

An excellent resource is offered by the University of Illinois and Urbana-Champaign: the UIUC Propeller Data Site. There are wind-tunnel data for a variety of hobby-craft propellers. The following video explains how we could arrive at our drag and thrust factors through measurement.

Propeller Thrust Coefficient

Let’s estimate the thrust coefficient from the, “static” data because a quadcoptor is not going to travel very fast along the axis of the propellers. A preliminary estimate for the static thrust coefficient after looking at a number of experimental data sets is…

From this we can calculate, ‘b’ as described above.

Propeller Drag Coefficient

Recall the motor+gearbox+propeller post. There we look at now to map the torque load of the propeller back through the gearbox to the motor shaft to get all the parameters correct in modelling the DC motor drive. In that post we have

“Aerodynamic torque for rotor R:        (page 9, Gafvert) for what looks like a comparable hobby-sized propeller lab set-up”.

We didn’t see how this term was derived though.

Confusion…

It’s a challenge (for me at least) to find a propeller, “drag coefficient” explicitly stated in many references. 

Then I read, “While an airfoil can be characterized by relations between angle of attack, lift coefficient and drag coefficient, a propeller can be described in terms of advance ratio, thrust coefficient, and power coefficient” here. That statement cleared some fog.

Next I found an online question from a person facing the same confusion I was experiencing. The question married the above quote to the same UIUC data I was looking at. The dialog there ultimately gets to a drag term resolved from the UIUC, “power coefficient”.

From there I made sense of it and as expressed in the video above. I derive my, ‘d’ term from the experimental, “power coefficient” as explained in the video. As the video explains, I’m estimating…

For the APC Carbon Fiber 7.8-inch propeller, if we plug-in numbers we calculate, ‘d’ as…

This is close to the number from page 9 of Gafvert above. We don’t know what specific propeller he is referencing, but we know it is in the family of these hobby propellers. We can be satisfied with this method for arriving at our drag factor from the UIUC data.

Conclusion…

We now have a good understanding of our thrust and drag factors for our propellers.

When we get into the design iterations we might try other propellers with different coefficients. We can revisit the UIUC data and compute a range of, ‘b’ and, ‘d’ terms as we did above. This is how we can assess propeller selection impact on performance of our quadrotor.

[latexpage]The last post was our introduction to the Linear Quadratic Regulator (LQR). We saw there that as we started with initial conditions or introduced a disturbance the LQR will drive the states to zero. In the simulations we saw the graphic of the copter converge on the zero state: zero roll, pitch, yaw, and respective rates all to zero.

The, “Control Law” (feedback gain K) was obtained through a solution of the matrix Ricatti equation buried in Matlab. We’ll dig into that math at some other time.

This same solution is relevant for the, “tracking” problem or servo case: when we desire the plant to be controlled to a particular set of non-zero state values.

Basic Idea

We can think of it as the reverse of our last simulated cases where we started with non-zero state initial conditions (some angles for roll, pitch, and yaw) and observed the platform converge on the zero-state attitude.

Think of this new problem as starting with zero initial conditions (a level quadcopter) and desiring the attitude, “step response” to a non-zero input.

This is a bit artificial for the quadcopter because a non-zero attitude will mean thrust in a particular lateral direction which we are ignoring at present. Our controller design will still be relevant, as this is the attitude control loop.

Ultimate Goal

Later we’ll introduce an outer position control loop that supplies the reference inputs to the attitude controller: the reference input that is the topic here. In flight it will not be a set value but a continually changing attitude reference input based on the outer position control loop’s desire to move the platform by altering the attitude and the resultant thrust vector.

Before we can illustrate how an outer position controller will command attitude to this inner attitude controller we need to understand how to introduce the reference for roll, pitch, and yaw.

The respective rates are states also, but our goal is to regulate these to zero as the attitude tracks the body angle input reference. Our reference vector will supply zero for the rate states.

Theory

Feedback Control of Dynamic Systems by Franklin, Paul, and Emami-Naeini introduces the topic via the material below. The following page copies are from a 1994 edition of their textbook.

Quadrotor Implementation

We are dealing with more states and a multi-input, multi-output (MIMO) problem. The equation to solve from the reference above is shown here with bold matrix elements representing matrices themselves.

Our ‘A‘ and ‘C‘ Matrices are 6×6. The ‘B‘ and ‘D‘ Matrices are 6×4. This results in the matrix needing an inverse being 12×10, which is non-invertible. On page 510 Stengel writes we can employ the right pseudoinverse in this case.

\begin{equation}
\begin{bmatrix}
\mathbf{N_{x}} \\
\mathbf{N_{u}} \\
\end{bmatrix}
=
\begin{bmatrix}
\mathbf{A} & \mathbf{B} \\
\mathbf{C} & \mathbf{D} \\
\end{bmatrix}^{-1}
\begin{bmatrix}
\mathbf{0} \\
\mathbf{1} \\
\end{bmatrix}
\end{equation}

We can then solve for the reference gains. The resulting N_x is 6×6. The resulting N_u is 6×4. We now have the matrices we need to operate on our state reference.

N_x and N_u are simplified into a composite gain in the second block diagram according to the reference textbook above:

\begin{equation}
\bar{N}=N_u+KN_x
\end{equation}

We’re now ready to see how the quadrotor LQ controller will track a reference input. The controller gain matrix K is from our LQR solution, so it’s the same controller.

Matlab Simulation

This Matlab script is a generalized version of the script in the last post covering the LQR simulation. In this script you will see the reference gain N is established and applied to a reference input. Yaw-axis sinusoidal reference tracking is illustrated in the following video generated by running the script.

You’ll need the same three additional files from the last post (see the matlab script section near the end).

Download the script, find yourself a Matlab seat, and explore!

Conclusion

We’ve covered how to solve for the reference input gain matrix and modified our LQR Matlab simulation script accordingly. The simple simulation shared in the video above indicates we’re stable and generally tracking the yaw sinusoid.

There’s more to do in this area: examine disturbance rejection, tracking performance, and general analysis of our LQR controller. Through this process we’d apply some performance assessment metrics and iteratively adjust our cost function weights.

We’re not yet to a real design and iteration process yet. We’re still building the tools and the understanding of how this MIMO Linear Quadratic quadcopter math works.

Happy Learning!

[latexpage]

In a previous post we got the Ardupilot simulator stable running the linear-quadratic attitude controller, but we had a hack in-place for normalizing our body torque commands from Newton-meters to -1 to +1 for input to the Ardupilot motor sub-system.

Purpose

This post uses a test stand to measure maximum available thrust over a range of supply (battery) voltages in order to gain a normalizing constant we need for the output commands of our linear quadratic controller. With normalized commands we drive -1 to +1 commands into native Ardupilot motor code which in-turn drives PWM outputs to electronic speed controllers (ESCs).

Equipment

After all the modelling and equation solving for the motor-prop assembly in previous posts I’ve been wanting to take some measurements. I found a reasonably priced test stand from Tyto Robotics: their model 1520. I purchased the add-on optical speed sensor. All-in with shipping was near $200.

Here’s the 1520 on my bench with the AIR2216/880Kv motor and T1045 prop that comes with the Hexsoon EDU-450 kit. On my setup here you can see a green LED on the RPM sensor board as well as the reflective tape on the motor rotor.

Problem

In the Linear-Quadratic controller implementation within the Ardupilot Attitude Controller we’re calculating roll, pitch, yaw command outputs as Torques. Challenge is, Ardupilot motor code desires normalized roll, pitch, and yaw commands scaled -1 to +1.

Focusing on a single axis, say roll, once we have a roll torque to command on each iteration we need to cram a normalized roll (-1 to +1) into the Ardupilot motors object. From there it will result in PWM output to electronic speed controllers (ESCs) to drive the motors.

Roll torque is differential thrust from a cross-body pair of propellers multiplied by the arm length. If we divide our controllers command torque by arm length we get this differential thrust command but…

How to normalize it for the Ardupilot motors and PWM stack of code?

• We need a maximum thrust
  • with added twist that it’s voltage-dependent as battery voltage drops during flight.

Solution

Characterize the motor-prop assembly maximum thrust as a function of battery voltage. Call it, Thrust_max.

Given Thrust_max and also run-time battery voltage sampled in real-time, scale Thrust_max down to a Thrust_available.

Then at run-time we’ll normalize that differential thrust for roll and pitch with a divide-by Thrust_available (first clamping the differential command to +/- Thrust_available ). This will give us a normalized -1 to +1 roll and pitch to command the resident motors object in Ardupilot.

Yaw is different as it’s driven by propeller drag force and resultant moment induced about the craft’s center. We’ll look at that later.

Max thrust derating by battery voltage

The first plot below illustrates what to expect for thrust as our Lithium Polymer (LiPo) battery voltage drops during flight. A single LiPo charged cell is nominally 3.7 volts to 4.2V. We’ll likely use a ‘4S’ configuration (4 cells in series) for nearly 15 Volts at max charge.

I don’t yet know what a low voltage cuttoff will be. I read maybe 3.2V per cell which means our minimum operating voltage might be say 10-12V. In any case I tested with my bench supply as follows in order to produce data seen here.

Procedure

1. Set bench supply at 15V.
   • This is source voltage to the ESC driving the motor on the test stand shown above.
2. Step ESC from 1000 to 2100 with intervals of 100.
   • This results in a 1-2ms pulse to the ESC in 100 microsecond steps.
3. Hold each step.
4. Sample motor speed and torque mid-step.
   • When the step’s steady-state is established.
5. On pause after last step.
   1. lower battery voltage.
6. Jump to step 2 above.

Measured Thrust over ESC range by Battery Voltage

We’ll get our data off the 2-D plot below. These 3-axis plots illustrate the test runs and reveal the relationship between battery voltage, ESC command, and measured motor speed and torque.

plot: Thrust output by ESC command input over range of Battery voltage.

Thrust derates with voltage because the maximum RPM drops with voltage. You can see this in the following plot where we’re plotting against measured RPM instead of command. We don’t actually drive thrust, we drive motor speed, which in turn provides thrust proportional to RPM squared as described in an earlier post.

In the perspective view here you can see it’s actually RPM that is dropping with lower voltage, and consequently thrust. Thrust it our output of interest though.

plot: Thrust output vs RPM output from ESC command over range of Battery voltage.

What we care about is the maximum thrust possible as a function of battery voltage. Note from the plot below how no-load battery voltage droops up to max thrust (RPM) for each step. Let’s take our voltage difference from the loaded maximum drive output for the first and last step. We’ll take the thrust difference and get our slope. It’s nicely linear.

\begin{equation} 
 \Delta V = 13.75 – 8.25 = 5.5V
\end{equation}

\begin{equation} 
 \Delta Thrust = 10.5 – 4.5 = 6N
\end{equation}

\begin{equation} 
 m = \dfrac {\Delta Thrust}{\Delta V} =\dfrac{6N}{5.5V}\approx 1\dfrac{N}{V}
\end{equation}

Plot: Maximum Thrust at Battery Voltage

Run-time equation for our normalizing value: Thrust_available

The equation here will replace the normalization constants described in an earlier diagram. Our interest then was to get stable with the simulator, so we hacked-in a scheme. Here we’ll have an acurate divisor for the -1 to +1 motors input range.

\begin{equation} 
 Thrust_{available} = Thrust_{cal} + m\cdot(V_{batt,current} – V_{batt,cal})
\end{equation}

Conclusion

We’ll confirm these results with a battery stack instead of a bench supply to finalize our numbers, but it looks like we have a tidy linear derating to apply here. We’ll calibrate somewhere near a fully-charged battery pack voltage and store these two parameters for Ardupilot:

• Thrust_cal (Newtons)
• ThrustV_cal (Volts)
• ThrustSlope (Newtons per Volt)
  • This is unity above but we’ll store it anyway

The above run-time equation and application of these parameters will cover cases where run-time voltage is higher than the calibration voltage.

In this post we’ll dissect and simulate basic concepts embodied in early gyro instrumentation applications in combat gun sights. We’ll build a deeper understanding of gyro mechanics through this application of their first principle: output torque resulting from precession. We’ll learn some interesting history too.

The basic problem to solve is this: when we fire a ballistic projectile how can the sight compensate to have the projectile intercept a moving target down range? We desire to hold the sight on the target and not guess how much to lead it by.

The British did early work in the Mk IID sight that was also produced in the US by the Sperry Gyroscope Company as the K-14 that saw use in US fighters during WWII. Charles Stark Draper started working in the same area in 1940 out of his Instrumentation lab at MIT.

Charles Stark Draper’s Gunsight Designs

Draper’s first efforts in naval gunnery stemmed from a prototype he built in 1940 illustrated here. If some independent sensor could measure the crossing speed of the target paper, “Doc” is targeting and simultaneously know the line-of-sight of the barrel axis this independent actor could command the site in the “shoebox”.

The complicating factor is that the shoebox is the actor measuring the target crossing rate (drop rate too) by way of integrated gyroscopes. It is the, “computer” that assumes its own rates represent the target’s. This will present a problem we’ll cover below: tracking instability.

From this early, “shoe box” effort Draper designed anti-aircraft fire control systems for WWII ships. This became the Mark 14 sight manufactured by the Sperry Gyroscope Company. There’s some interesting information here.

John H. Blakelock’s, Automatic Control of Aircraft and Missiles, Second Edition Chapter 9 and Appendix J present Draper’s later work on the A-1 air-to-air sight. Blakelock writes on page 324 that Dr. Draper, “was probably the first person to completely understand the tracking instability phenomenon and was thus able to solve this troublesome problem”. This might be generous towards, “Doc” as he was called, given the amount of varied work in this field. Click the image below for book information.

There’s a hint to this, “tracking instability” problem in training material for the Naval sight Mark 14 highlighted below. Our simulator will reveal why the Navy gave this advice. When a gunner increases slew rate of gun to catch-up to a target the pipper actually moves backwards in a manner that can confuse settling to smooth tracking. However, as you slow from ahead of a target this effect is less dramatic and easier to track to zero. We’ll make sense of this problem and identify how Draper and others mitigated the problem.

Draper’s A-1 Gunsight

Blakelock’s Appendix J introduces Draper’s later A-1 and A-4 air-to-air gun/rocket sight design. They were too late for WWII, but achieved some fame in the F86 Sabre vs. MIG battles over Korea in that later conflict.

The A-1 sight incorporated radar ranging. This eliminates the manual stadiametric ranging involving throttle grip twisting to close a concentric ring in the HUD on the wingspan of a target after having manually clicked the sight to specify target size. The promise of radar ranging eliminated this tedious process, but not without some growing pains as described in the article below, extracted from Air Power History (SUMMER 2009 – Volume 56, Number 2).

If you don’t read the entire article, scroll to the end and read the last paragraphs discussing the concepts of a, “heterogeneous engineer” and, “innovative departure”. These are relevant notions when attempting to gain traction when innovating in any industry or area.

Next Up…

With the above history a guide, the next post will present a simplified dynamic model we can use to investigate the, “tracking instability” problem in particular. There will be some matlab code to share so you’ll be able to play with the model yourself.