Lab 3:A Disturbance in the Force

The questions below are due on Friday October 30, 2026; 10:00:00 AM.

You are not logged in.
Please Log In for full access to the web site.
Note that this link will take you to an external site (https://shimmer.mit.edu) to authenticate, and then you will be redirected back to this page.

You are viewing the PRELIMINARY 6.3100 Fall 2026 web site. Please note that information on this page may change throughout the semester, please check back frequently.

Quick Reference

NOTE

When the sketch is uploaded to the Teensy, it sets K_p, K_d, and K_i to 0.
To get non-zero values, you must set them via the serial plotter send window.

Serial Plotter Send Window Indices and Offsets

Task Summary

Week A

Add a second motor to your arm, download the new sketch, re-calibrate the angle sensor offset, and level the arm (Checkoff 1).
Set K_i = 4.0, K_p = 0.5, and find the smallest value of K_d for which the system is stable (the step response settles to its steady state after about seven visible oscillations).
Measure the frequency of the oscillation in the step response, and the magnitude and phase of sinusoidal steady-state responses for several input frequencies (Checkoff 2).
Set the input to zero (6 0 in send window), then disturb the arm by fanning it, and find the "fanning resonance frequency" that maximizes the disturbance response.
Consider how the following are related:
- the frequency of the oscillations in the step response,
- the frequency that maximizes the sinusoidal steady-state response,
- and the frequency that maximizes the disturbance response.
Set K_i=0, K_p = 2, and for this case find a new smallest K_d for which the step response settles after about seven visible oscillations. Turn off the input (6 0), and find the fanning resonance frequency (Checkoff 3).
Download the IPython file, ctPropArm.py, to a directory, start ipython in that directory, and at the python prompt type %load ctPropArm.py

Double the fun

Add a second motor to the right side of the arm, and reclamp the arm to the lab table. As shown in the pictures below:

Be sure that the propellers on BOTH motors are smooth side up!
Be sure to use a SHORT axle to prevent over-rotation, long axles can get tangled in the wires.
Be sure to separate the propeller connector from arm connector by three brick lengths on both the left and right sides.
Securely clamp your arm to the side of the lab table. The clamp should be close to the PCB and the propeller arm should be at least two inches away from the lab table (to minimize air resistance).
Check that the wiring allows the arm to rotate from -180 degrees to +180 degrees, or from -\pi to \pi radians (you will need to remove the axle to check this, be sure to reattach it afterwards).

Power up the Propeller Arm

Connect both power supplies (the 5V USB power supply AND the 5V 10 amp power supply) to the controller board and connect your laptop to the Teensy. Download the sketch linked at the top of this page, upload it to the Teensy, and start the serial plotter. Note that the offsets for K_p, K_d, and K_i are all zero. Below you will change these gain values to get enough feedback for the arm to level itself.

Hold the arm level and ensure the plotter shows a zero angle. If not, adjust the AngleOffset in the Teensy sketch (you may be able to use the same AngleOffset from the week 2 sketch). Once the sensor is calibrated, stop holding the arm, and try to level the arm by adjusting the black potentiometer (in the upper right hand corner of the board). The potentiometer can be used to adjust the relative thrust of the right and left propellers, so that the arm holds itself level (set K_p = 0.1 first, otherwise the arm will just drift).

Tip: If your motors are not turning, here are some troubleshooting steps:

Is the high-current power supply (round connector) plugged in?
Did you download the sketch for Lab 3 and upload it to your Teensy?
Check the connection between each motor and its extension wires. It is very easy to bend the pins in the extension wire connector, THIS IS THE MOST COMMON ISSUE!!! Get new extension wires!
Check the connection of the extension leads to the screw terminals on the circuit board. There should be about 3mm of bare silver wire for the connector to grab, and it should be invisible outside the screw terminal. If uninsulated silver wire is visible outside the screw terminal, that wire could "short" and damage the board.
Use a multimeter (there are about a dozen of them distributed around the lab, ask for help) to check the potential difference across the two terminals at the controller board. Touch the multimeter voltage probes to the screws and check if there is a non-zero voltage; if not, feel free to ask staff for help.

Tip: Make sure that both propellers are blowing air downwards. If they are not, there are two ways to fix this: swapping wires or editing the code. EDIT THE CODE, it is much easier. Find lines 278 and 280 and swap the order of B1L and B2L or A2R and A1R, depending on which motor needs to be changed; note that these ports are labeled on the PCB.

The figure below shows serial plotter output from running this lab's sketch on a propeller arm (from a prior term of 6.310!!!) with a sinusoidal input (in red) and the resulting measured angle (in blue). Notice the two intervals measured on the plot, one denoted with a black double arrow and the other with a green double arrow. How does the ratio of those two intervals relate to the phase difference between the two plots?

Checkoff 1:

Show that the powered arm is level and the plotted arm rotation is zero.
How does the ratio of the black to green intervals in the figure above relate to the phase difference between the desired arm rotation (the input) sinusoid (red) and the measured arm rotation sinusoid (blue)?
Estimate the phase shift in the figure above. Is the phase shift of the measured arm rotation positive or negative (assuming the desired arm rotation has zero phase shift)?

Steps, Sinusoids, and Disturbances

Step response

Set K_p=0.5; K_d=0.3; K_i=4.0 and then type "6 0.1" and "5 0.05" to generate an input that is high for 10 seconds, then low for 10 seconds, and then repeats (note that the plot window is 10 seconds wide, and the plot points are 20 milliseconds apart). Then try a few different values for K_d. You should be able to find a K_d between 0.1 and 0.3 that produces a QUALITATIVELY similar step response (its step response should oscillates as much) to the step response shown in the following video (an earlier version of the propeller arm was used in this video, do not try to compare numerical values!!). Notice that in the video, the step response visibly oscillates about seven times before settling; the step response of your arm with your selected K_d should oscillate similarly (seven visible oscillations).

Enter the frequency of oscillation (in Hz) seen in your barely-stable step response:

Measure Frequency Response

In theory, measuring the frequency response of a STABLE linear system is easy:

Pick a frequency and apply a sinusoid of that frequency to your system's input.
Wait for the system's natural responses to die out (they must die out eventually, since the system is stable). What remains will be sinusoidal at the same frequency as the input.
Compare the input sinusoid to the output sinusoid and compute the amplitude ratio and the phase difference.

We can measure the frequency response for our two-propeller arm system using the above approach as long as we are careful about the three issues described below.

Monitor Motor Commands

While the two-propeller arm is intended to cancel the gravity-induced nonlinearity of the one-propeller arm, the cancellation is only perfect if the arm is perfectly symmetric (so balance the arm as best you can!). In addition, there is still a "clipping" nonlinearity that arises when the controller-generated motor commands exceed power supply limitations. This nonlinearity can usually be avoided by reducing the magnitude of the sinusoidal input, provided the clipping has been detected.

Monitor the motor command plot whenever measuring frequency response, to ensure the commands are not being "clipped" by power supply limits. The clipping effect is easy to spot, as shown in the figure below. The black arrows point to cases where the motor command (the green curve) flattens, an indication that the controller's commands are being clipped because of power supply limitations.

The clipping problem can usually be eliminated by reducing the input amplitude, as shown in the picture below. Avoid reducing the input amplitude too much though, because the microcontroller won't be able to represent the signal accurately. Use index 6 in the plotter window to set a reasonable sinusoidal amplitude (keep the magnitude of the amplitude at least 0.03, i.e. no closer to zero than -0.03; smaller amplitudes are not converted accurately by the controller board).

NOTE: USE NEGATIVE AMPLITUDES. THE SKETCH WILL GENERATE STEPS INSTEAD OF SINES IF YOU USE POSITIVE AMPLITUDES.

Wait for Steady-State

It can take several seconds to achieve sinusoidal steady-state; a typical case is shown in the picture below. Be sure to wait long enough (ten seconds or so) before pausing the plotter and taking measurements.

Spread Samples

In last week's lab, we sent an update to the serial plotter every 10 milliseconds. For this week, the sketch default is to send an update to the serial plotter every 20 milliseconds. At higher sinusoidal frequencies, it may be helpful to decrease the time between plotter updates to measure time more accurately. The control algorithm runs every millisecond, so the fastest serial plotter update is once per millisecond (zero skips). To set the fastest update rate, type "4 -19" into the plotter send window (as was done to generate the plot below). Note that the default offset for the number of skips is set to 19 in the code, and the value you send is added to that offset. As another example, to update the serial plotter every 10 milliseconds, we want to skip 9 samples, so type "4 -10" (to reduce 19 skips to 9) in the plotter send window.

Checkoff 2

For the second checkoff, please measure the frequency response of your propeller arm, with K_p=0.5, K_i=4.0 and the value for K_d you found above. You should notice that for a very narrow range of input sinusoid frequencies, the arm response increases by about five fold relative to its low-frequency response. We refer to this phenomenon as a resonance, and the center of the associated frequency range as the resonance frequency.

Remember that the default time between plot points is twenty milliseconds.
Measure the frequency of the oscillation in the step response.
Measure the amplitude ratio (measured amplitude / input amplitude) and phase shift of sinusoidal steady-state responses for several input frequencies. Try to find the lowest frequency for which the phase shift is about -45 degrees (anywhere between -30 and -60 degrees is fine), and the lowest frequency for which there is a -135-degree phase shift. We suggest you use index 6 in the plotter window to set a very small sinusoidal amplitude (-0.03), because, at just the right frequency, the arm response can become quite large.

Checkoff 2:

How does phase change with frequency?
How much do the amplitude and phase change at frequencies that are twenty percent above and below the frequency of the peak response?
What frequency gives you a phase shift of -45 degrees? -135 degrees?
What is the relationship between the frequency of the oscillations in the step response and the frequency that maximizes the sinusoidal steady-state response?

Once Again, with Feeling

Set the arm input to zero (6 0 in the send window), but do not change the feedback gains (K_i=4, K_p=0.5, and K_d is the just-barely-stable value you found in the last checkoff), and try to "fan" your system into resonance. You should be able to excite a resonance, like in the video below (an earlier version of the propeller arm was used in this video, so again, do not try to compare numerical values!!). Note that the resonance frequency for the arm in the video below is probably higher than yours, and you will try the higher resonance frequency case next. At higher resonance frequencies, it can be hard to "match" your fanning to the resonance, so please watch until (or skip to) the end to see what that looks like!!

Now try the higher resonance frequency case: set K_i=0, K_p = 2, and adjust K_d so that you get a "barely stable" step response like in the video below (an earlier version of the propeller arm was used in this video!!).

Checkoff 3

For this checkoff,

Set K_i=4, K_p=0.5, and K_d to the just-barely-stable value you found in the last checkoff.
TURN OFF THE INPUT, and try disturbing the arm by fanning it. Find the "fanning resonance frequency" that maximizes the response.
Set K_i=0, K_p=2.0, and K_d to the just-barely-stable value for this no-integral-feedback case.
Measure the frequency of the oscillation in the step response, and determine the input frequency that maximizes the sinusoidal steady-state response.
TURN OFF THE INPUT, and try disturbing the arm by fanning it. Find the "fanning resonance frequency" that maximizes the response.
For these two sets of controller gains, what did you observe about the relationship between the frequency of oscillations in the step response, the frequency that maximizes the sinusoidal steady-state response, and the fanning frequency that maximizes the disturbance response?

Finally, verify that you can run ipython with the control library by downloading our IPython file, ctPropArm.py, to a directory.

YOU DO NOT HAVE TO USE IPython if you prefer another python environment!!! ALSO, the file contains lines that say FIXME, do not worry about correcting those, we will address them in part B!

Open a terminal window in that directory (PowerShell on Windows, Terminal on Macs) and type ipython and return. NOT python but ipython for interactive python.

If ipython is not recognized, use pip to install ipython. Type pip install ipython followed by a return, then pip install control followed by a return. You may also need to install numpy and matplotlib using pip.

If python and pip are not installed on your laptop, but pip3 and python3 are, try the same commands with pip3 and ipython3. If that does not work, we suggest installing the Anaconda package here, but you can also install python from python.org.

Once you have installed ipython and started it, you should get a prompt. At the prompt type %load ctPropArm.py and then press return. You should see plots and a command window like in the screenshot below.

Checkoff 3:

Demonstrate for the two sets of controller gains that you can fan your arm into resonance, and explain what you observed about the relation between the oscillations in the arm's step response, the arm's frequency response, and the arm's response to fanning disturbances.
Show that IPython (or the environment of your choice) works on your laptop.

STOP HERE FOR WEEK A!!!!

Update the Setup and the Sketch!

Updates to the Teensy sketch

Please update the nominal command specified on line 15 in the Teensy sketch to 0.5 as in

#define nomCmd 0.5 // Motor command that nominally hold arm level. and please update the offsetExtraBack parameter specified on line 39 to zero as in

#define offsetExtraBack 0

Updates to the Arm Hardware

Please update the propeller arm by moving the propellers further apart and by adding weight to both the left and right propellers, as shown in the figure below.

Why the changes?

The above changes help stabilize the arm and make its behavior more linear, but more importantly, the increase in arm inertia results in a more compelling control problem (we also eliminated a poorly-justified approach to reducing propeller noise so you can replace it with a better-justified strategy).

Continuous Time (CT) Modeling

Last lab we derived an approximate discrete-time (DT) model for the innately continuous-time (CT) propeller arm. We chose discrete-time partly to match our innately discrete-time Teensy-based controller, and partly because discrete-time models, at least the ones formulated as systems of linear difference equations, have explicit solutions in the form of easily computed sums. For this and the next few labs, we will do the reverse. That is, we will derive CT models for innately-CT physical artifacts, and create approximate CT models for our innately-DT controllers.

We are switching to CT because the concepts we are exploring next, transfer functions and frequency response, have a simpler and more intuitive description in CT. And that simplicity will allow us to gain new insights into control system design.

Doubly-Driven Arm

The larger number of forces and torques for the doubly-driven arm are depicted using a more abstract picture (below) than we used for the singly-driven case. In this picture we show the motor commands (C_{nom} and c), the measured and desired arm angles (\theta_a and \theta_d), and the motor thrust (F_{nom} and f_{thrust}) and gravity (F_{gravity}) forces.

This picture also explicitly illustrates the arm torque (\Gamma_{thrust}) due to the difference in right and left motor thrusts. In the last lab, we did not model disturbances, and implicitly assumed the net torque on the arm was due to motor thrust alone. Under that assumption, \Gamma_{thrust} was an unnecessary intermediate in the relation between motor command and arm acceleration. This lab is focused on modeling disturbances, and since our example disturbances (fanning and dropping weights) generate torques, their effect on arm acceleration will be easier to model if we can add them to an explicit \Gamma_{thrust}.

Discrete-Time (DT) to Continuous Time (CT)

In the first propeller-arm lab, we related the change in arm angle, \theta_a, to the arm angular velocity, \omega_a,

\theta_a[n] \approx \theta_a[n-1] + \Delta T \cdot \omega_a[n-1], \;\;or\;\;equivalently,\;\; \frac{\theta_a[n] - \theta_a[n-1]}{\Delta T} \approx \omega_a[n-1],

an approximation to the CT differential equation,

\frac{d\theta_a(t)}{dt} = \omega_a(t),\;\;\;or\;\;\; \theta_a(t) = \int \omega_a(\tau)) d\tau.

Similarly, we related the change in arm angular velocity to the arm angular acceleration, \alpha_a,

\omega_a[n] = \omega_a[n-1] + \Delta T \cdot\alpha_a[n-1], \;\;or\;\;equivalently,\;\; \frac{\omega_a[n] - \omega_a[n-1]}{\Delta T} \approx \alpha_a[n-1],

an approximation to the CT differential equation,

\frac{d\omega_a(t)}{dt} = \alpha_a(t),\;\;\;or\;\;\; \omega_a(t) = \int \alpha_a(\tau)) d\tau.

Block Diagrams

Simple Model, Simple Control

In our simplest arm model, angular acceleration is assumed proportional to motor command, and in the simplest controller, the motor command is assumed proportional to angle error. In discrete time,

\alpha_a[n] = \gamma_a c[n] \;\;\; c[n] = K_p \left( \theta_d[n] - \theta_a[n] \right),

and in continuous time

\alpha_a(t) = \gamma_a c(t) \;\;\; c(t) = K_p \left( \theta_d(t) - \theta_a(t) \right).

where K_p is the gain of the proportional feedback.

We can combine the equations for the arm, from the previous section, and the equations for the motor and controller (above) into a CT block diagram, shown below.

Note that in the CT diagram, we combined blocks with integral operators rather than combining blocks with derivative operators, a choice that parallels our intuition. That is, the system's controller sets the acceleration, which integrates to velocity, which in turn integrates to position, which we then measure and feed back.

Simple Model, PID control

We know that proportional feedback alone can not stabilize the propeller arm, so derivative feedback is essential, and we also know that integral feedback improves steady-state tracking. Including all three (P,I,D) in a CT block diagram is easy, because the integral and derivative relations can be represented explicitly as integral and derivative blocks (see below), rather than by approximations based on differences.

Note that even though we are representing the arm control system with a CT block diagram, the derivative and integral controller terms will still be approximated on the microcontroller using scaled angle error differences and running sums. Nevertheless, we will continue to approximate the DT behavior with a CT representation, because we can model the system with simple block diagrams and develop intuition and insight.

Non-instant Motor Model

As we saw in the first propeller-arm lab, a more accurate model for the relationship between motor command and arm acceleration is

\alpha_a[n] \approx \alpha_a[n-1] + \Delta T \beta \left(\gamma_{a} c[n-1] - \alpha_a[n-1] \right) ,

where \beta and \gamma are both positive. Reorganizing the above difference equation

\frac{\alpha_a[n] - \alpha_a[n-1]}{\Delta T} = \beta \left(\gamma_{a} c[n-1] - \alpha_a[n-1] \right) ,

leads to an analogous differential equation,

\frac{d \alpha_a(t)}{dt} = \beta \left(\gamma_{a} c(\tau) - \alpha_a(\tau) \right),

or equivalently, the integral equation

\alpha_a(t) = \int \beta \left(\gamma_{a} c(\tau) - \alpha_a(\tau) \right) d \tau.

Adding the motor model to the CT block diagram leads to the figure below.

Transfer Functions (Background Section you can Skip)

The above block diagram, often referred to as an integrator-adder-gain block diagram (for obvious reasons), offers very limited insight into system behavior. Each of the blocks indicates a relationship between an input and output variable, e.g. the arm rotational velocity is the integral of arm acceleration, but there is no simple way to compose the blocks. We might imagine a double-integral block to synthesize a relation between arm acceleration and arm angle, but how will we compose the blocks and the feedback into a relation between desired arm angle, \theta_d, and measured arm angle \theta.

In this lab we are transitioning from DT difference equations and matrix representations to CT differential equations and transfer functions, because transfer functions can be easily composed. Since transfer functions are rational functions, that is, ratios of polynomials in s (where we have seen that s can represent the differentiation operator, or can be a complex number, depending on how we are interpreting the transfer function), when we compose transfer functions using addition, subtraction, multiplication or division, the result will always be a ratio of polynomials in s, that is a rational function of s, so the result will be another transfer function.

Transfer Function Notation

As discussed in lecture, regardless of interpretation, we replace differentiation with multiplication by s, as in

\frac{dy(t)}{dt} \rightarrow s Y

and integration with multiplication by \frac{1}{s} as in

\int y(\tau)d\tau \rightarrow \frac{1}{s}Y,

where the switch to upper case is intended to indicate either that the operator s is acting on an entire time function, or to indicate a scale factor for an assumed time dependence e^{st} as in

\frac{dy(t)}{dt} = \frac{d Ye^{st}}{dt} = sYe^{st}.

We denote the transfer function for the plant (for our arm, from motor command to propeller arm angle) as the rational function H(s), and denote the numerator and denominator polynomials as

H(s) = {\sum_{\tilde{l}=0}^{\tilde{L}_H} b^H_{\tilde{l}} s^{\tilde{l}} \over \sum_{l=0}^{L_H} a^H_l s^l} \equiv \frac{n_H(s)}{d_H(s)}.

In the above formula, \tilde{L}_H and L_H are the orders of n_H and d_H, the numerator and denominator polynomials of the H(s), respectively.

Similarly, the transfer function for the controller, K(s), is a rational function whose details are denoted as

K(s) = {\sum_{\tilde{l}=0}^{\tilde{L}_K} b^K_{\tilde{l}} s^{\tilde{l}} \over \sum_{l=0}^{L_K} a^K_l s^l} = \frac{n_K(s)}{d_K(s)},

where \tilde{L}_K and L_K are the orders of n_K and d_K, the numerator and denominator polynomials of the K(s), respectively.

We denote the transfer function that relates the measured arm angle to the desired arm angle as G(s), which can be determined from K(s), H(s) and the feedback configuration, by applying Black's formula. For the case of the arm controller, in which the measured arm angle is fed back and compared to a desired arm angle, G(s) is given by

G(s) = \frac{K(s)H(s)}{1 + K(s)H(s)} = {\sum_{\tilde{l}=0}^{\tilde{L}_G} b^G_{\tilde{l}} s^{\tilde{l}} \over \sum_{l=0}^{L_G} a^G_l s^l} = \frac{n_G(s)}{d_G(s)},

where we know that G(s) is a rational function (ratios of polynomials) because addition, subtraction, multiplication and division of rational functions result in rational functions.

Transfer Function Properties

Once we have determined

G(s) = \frac{n_G(s)}{d_G(s)},

we can compute natural frequencies (to verify stability), and calculate frequency response (to assess performance).

The system described by G(s) is stable if its natural frequencies have negative real parts. As we saw in lecture, if the complex number \lambda_i is one the system's natural frequencies, then d_G(\lambda_i) = 0, and therefore

G(\lambda_i) = \infty.

We can refer to the \lambda_i's as the nature frequencies of the system described by G(s), or as the poles of G(s). Regardless, for the system described by G(s) to be stable,

\mbox{Real}\big(\mbox{poles}(G(s)\big) < 0,

where the poles of G(s) are the roots of it denominator polynomial d_G(s).

If the poles of G(s) are stable, then the response to a sinusoidal input will EVENTUALLY be a sinusoid (sinusoidal steady state). Sines and cosines are weighted combinations of complex exponentials, and since exponentials have analytically simpler derivatives, we use in any mathematical analysis. That is,

e^{j\omega t} = \cos{\omega t} + j\sin{\omega t},

and the sinusoidal-steady-state (eventual) response is given by

G(j \omega) e^{j\omega t}

where G(j \omega) = G(s)|_{s = j \omega} and \omega = 2 \pi f is the frequency in radians/second (and f is the frequency in Hertz).

End Of Skippable Background Section

Frequency Response

So far, we have designed discrete-time controllers by trying to minimize steady-state error by maximizing K_p, within the constraints associated with stability (recall that in discrete time, natural frequencies are stable if are all less than one in magnitude). We could follow the same path in continuous-time, minimizing steady-state errors by maximizing K_p while maintaining stability,

\mbox{Real}\big(\mbox{poles}(G(s)\big) < 0.

But, then switching to CT would be about trading the stability region for DT, the interior of the unit disk, for the stability region for CT, the open left-half plane. And that does not seem worth the effort.

Whether in DT or CT, we can use frequency response to change our approach to controller design. As we will start to see in this lab, and explore further in the postlab, we can design good controllers if we know the actuator's frequency response (its sinusoidal-steady-state response to a sinusoidal inputs as a function of frequency).

Another reason for switching to CT before examining frequency-response-based controller design is that we can capitalize on our shared experience with the physical world, experience that gives us a lot of intuition about CT frequency response (e.g. the bass response of headphones, the fastest human drum roll, rate of leg pumping on a playground swing).

Arm Block Diagram using Transfer Functions

Replacing the integrals in the CT block diagram with \frac{1}{s} is the first step in using transfer functions, as in the diagram below.

REMINDER: USE CAPITAL LETTERS when using transfer functions. That is, use capital letters when substituting, for example, arm rotational velocity

\omega_a(t) = \Omega e^{st}.

As noted in class and in the prelab, blocks described by tranfer functions can be composed by multiplication, and internal feedback loops can be simplified using Black's formula, leading to the following diagram below.

What is the command to acceleration transfer function, H_{c2a}(s), and what is the acceleration to \theta_a transfer function H_{a2t}(s) (where the "t" stands for theta)? See below if you need a hint.

We can simplify the diagram further by denoting H(s) = H_{c2a}(s) H_{a2t}(s), simplifying the arm block diagram into two transfer functions, K(s), the controller, and H(s), the actuator, where each is a ratio of polynomials in s, as shown below.

Numerical Tools for CT and Frequency Response.

We recommend you add plt.close('all') to the python code right above the line s = ct.tf('s'). This will close all the plots.

We will use numerical tools in ipython to help simulate our system. After opening ipython, be sure you load our file ctPropArm.py (linked at the top of the page). To load the file, at the ipython prompt type

%load ctPropArm.py

After you hit return, you should see output like the image below. If you only see one plot, the plots may be on top of each other, try moving the plot to see if there are others beneath it.

You can use ipython to manipulate continuous-time transfer functions. First, you have to tell python to use s as the CT transform variable, by typing (this is already in the ctPropArm.p file).

>> s = ct.tf('s')

After that, you can define a CT transfer function. For example, suppose your system is described by a set of differential equations

\frac{d v(t)}{dt} = x(t), \;\;\; \frac{d w(t)}{dt} = v(t), \;\;\; \frac{d y(t)}{dt} = -0.1y(t) + w(t),

then

Y(s) = \frac{1}{s+0.1}\frac{1}{s}\frac{1}{s} X(s) = H(s) X(s),

where H(s) is the CT transfer function given by

H(s) = \frac{1}{(s+0.1)s^2}.

In ipython, if you type

>> H = (1/(s+0.1))*(1/s)*(1/s)

followed by

>> print(H)

ipython will respond with

H =
 
        1
  -------------
  s^3 + 0.1 s^2

(If you just type H instead of print(H), you will get the output in a hard to decipher format). Then if you type

>> pole(H)

ipython will respond with something like

[0, 0, -0.1000]

as you would expect since in CT, two of your poles are at zero. Note that the is no specific meaning to the order of the poles, so you might see,

[-0.1000, 0, 0],

or the since the poles can be complex in general, you might see

[-0.1+0.j  0. +0.j  0. +0.j],

all of which are equivalent.

Many (feedback, step, etc), and we will also make extensive use of frequency plots. You can use the classic bode command

bode(H)

which plots magnitude in decibles versus frequency in radians/second. We have created a version that plots magnitude versus frequency in Hertz (cycles per second), which is easier to compare to experiments. To use our non-classical Bode type

bodeF(H)

Our example file should have produced frequency response outputs for both of the two commands above, see the image above.

For example, consider a system with input x(t) and output y(t) described by a first-order differential equation,

\frac{d y(t)}{dt} = -y(t) + 1000 x(t),

where y is the output and x is the input. Suppose we know that x(t) = \cos(\omega t) = \cos(2\pi f t) where \omega is the frequency in radians per second and f is the frequency in cycles per second. Approximately what is y(t) in the limit as t \rightarrow \infty ?

The transfer function is easy to compute using the techniques described in lecture and above, and is H(s) = \frac{1000}{s+1}. If x(t) = u(t), the unit step, then x(t)=0 for t < 0, and x(t) = 1 for t \ge 0. To plot the response to the unit step using ipython, in the command window type (be sure you have loaded our ctPropArm.py file),

H = 1000/(s+1)

to define the transfer function to ipython, and then

step(H)

to plot y(t) given x(t) = u(t). Note that y(t) rises to 1000 with a time constant of one second (about two-thirds of the rise in one second).

The frequency response of this first-order system can be plotted using

bode(H)

bodeF(H)

Our ipython file uses the command feedback to apply Black's formula to K(s)H(s), but the feedback command applies a more general form of Black's formula than we have seen so far. The feedback command takes two transfer function arguments, H_1 and H_2 to compute an input-output transfer function. In particular, if

Y = H_1(s) (Y_d - H_2(s)Y),

then

G = feedback(H_1, H_2)

will compute the input-output transfer function

Y = G(s)Y_d,

where

G(s) = \frac{H_1(s)}{1+H_1(s)H_2(s)}.

If we use unity-gain feedback with this first order system, then

x(t) = y_{desired}(t) - y(t)

and by using H_1(s) = K(s)H(s), H_2(s) = 1 in the generalized Black's formula above, the closed-loop system function, G(s), is given by

G(s) = \frac{K(s)H(s)}{1 + K(s)H(s)} = \frac{1000}{s+1001}

which ipython can determine for us using the feedback command

G = feedback(1000/(s+1),1)

To plot the step response of the feedback system,

step(G)

The frequency response also shows the impact of including the first-order system in a closed loop, as can be seen in the plot generated by

bode(G)

bodeF(G)

What is different about the open-loop to the closed-loop step and frequency responses?

If we want to determine y(t) given x(t) = \cos{\omega t} in the limit of large t (what we call sinusoidal steady state), we can use the frequency response. That is,

y(t) = |G(j\omega)|\cos{(\omega t + \angle G(j\omega))}

The magnitude and angle of G(j\omega) can be determined from the Bode plot, as as shown in the figure below (the figure below uses plots from another application, matlab, so your plots will not be identical).

Determining response at a particular frequency from the Bode plot.

Implement the model

Determine an expression for Hc2a.

Enter your expression for Hc2a, as a ratio of polynomials in s, in the box below, using Python syntax (e.g., ** for exponentiation) and use, gamma_a for \gamma_a, beta for \beta, and common constants as e, j, pi, etc.

Determine an expression for Ha2t.

Enter your expression for Ha2t, as a ratio of polynomials in s, in the box below, using Python syntax (e.g., ** for exponentiation) with Greek letters written in English (e.g., alpha for \alpha) and common constants as e, j, pi, etc.

Implement your expressions for Hc2a and Ha2t in the python ctPropArm.py linked at the top of this page.

Running the script after you have corrected the transfer functions should produce plots like (but not exactly like) the ones below. On the left is our rescaled Bode plot in which magnitude is plotted directly, and the logrithmically-placed frequencies are in Hz (cycles per second). In the middle is a classic bode plot for which the magnitude is plotted using the logrithmic dB scale and frequency in radians/second. And the plot on the right shows the associated step response

Calibrate your model

In order to find values of \beta and \gamma_a so that your model predicts your measured data as you vary controller gains (K_p, K_d, and K_i), you could try to match model and measured step responses for a particular set of gains (just like in lab 2). Instead, we are going to use frequency response (aka sinusoidal steady-state response) to try to find \beta and \gamma_a by pick controller gains so that we see modest peaking in the arm's sinusoidal steady state response, and then match it to the frequency response of our python model. We strongly recommend using the values for gains suggested in the steps below, as it is very easy to pick values that lead to nonlinear effects.

Set K_p=1.6, K_i=0.
On your propeller arm, find a derivative gain, K_d, so that the maximum (over frequency) response to a 0.01 amplitude input sinusoid is three times larger, (has amplitude 0.03). Record the maximizing frequency.
In your python model, set K_p, K_d, and K_i to the same gains used on your propeller arm, and adjust the \beta and \gamma in your model so that its frequency response matches the behavior of your arm (same amplitude ratio, same maximizing frequency).Neither \beta or \gamma should be greater than one hundred or less than fifteen.
Check that your calibrated model predicts the sinsoidal steady state at twice the maximizing frequency.
Validate your model by repeating steps one and two for the gains K_p=0.8 and K_i=0. Compare to the response predicted by your calibrated python model.

As an example, below is the sinusoidal-steady-state response of an arm in which we used four-brick weights one each propellor motors (YOURS SHOULD ONLY HAVE THREE). The response is at the output-maximizing input frequency. Note that the output amplitude is roughly four times the input amplitude and the frequency is approximately 2.1 hertz.

The frequency response for a well-calibrated python model is below, and as the zoomed in version makes clear, the peak frequency is about 2.1 hertz and the magnitude is four.

Checkoff 4:

Show that you have updated to the propeller arm with the wider spread of propeller motors and the pair of three-brick weights.
Show your calibration results. Show your plot of the frequency response of the model, and compare results at the maximizing frequencies to those in your experiment. Did you get a good match between the experiment and the python model? How did you check?
In the description of the python feedback command, we gave the equation for a generalized feedback system, Y = H_1(s)(Y_d - H_2(s)Y). Draw a block diagram for this generalized feedback system showing the input Y_d, the output Y , a summer, and the placement of the transfer function blocks for H_1 and H_2.

Disturbing Problems

When we dropped a weight on the arm or fanned it, as shown in the images above, we were testing our arm controller's ability to reject those disturbances (that is, maintain a desired arm angle). We also observed that if the derivative gain, (K_d,) is high, then motor commands are noisy, presumably because any noise in the angle measurement is magnified by the differencing used to estimate derivatives. Weight-dropping and fanning are examples of torque disturbances, and in this case, they primarily effect arm angle. The angle measurement noise is an example of measurement disturbance, and this case,it primarily effects the motor command.

If we want to design controllers that minimize angle displacement due to torque disturbance or minimize command noise due to measurement disturbance, the first step would be to determine the transfer functions that relate disturbance to output of interest, and then try to design controllers that minimize those transfer functions. Disturbance transfer functions are typically much harder to measure than input/output transfer functions, so typically we calibrate a model using input/output behavior (like in checkoff 4), and then analyze the model to determine the disturbance transfer functions.

Torque disturbances

It is easiest to add disturbance torques (generated by dropping or fanning) to the CT block diagram by thinking of them as disturbance accelerations (which is just a rescaling), as shown below.

Notice that we introduced \alpha_{thrust}, the acceleration due to the thrust torque generated by the motor-propeller pairs, and added a disturbance acceleration, \alpha_{disturb} (as might be generated by "dropping" or "fanning"). Then the sum of two accelerations becomes the arm angular acceleration.

What is the transfer function from \alpha_{disturb} to \theta_a, which we will denote as G_{disturb2\theta_a}(s), or more compactly, G_{d2t}. To determine G_{d2t}, we suggest you convert the above diagram using transfer functions (e.g. K(s), H_1(s), and H_2(s) as we did above). Then, rather than trying to use Black's formula directly, try writing two equations and solving. That is, given

\Theta_a = H_{a2t}(s) \left( A_{disturb} + A_{thrust} \right),

write an additional equation for A_{thrust} in terms of \Theta_d and \Theta_a, and then solve the two equations to determine a relation between A_{disturb} and \Theta_a (assuming \Theta_d = 0).

Measurement Disturbances

To determine the impact of measurement noise on motor command, we can start by adding a measurement disturbance to our CT block diagram, \theta_{noise}, as shown in the diagram below.

We can use the above block diagram to help us determine the transfer function from measurement disturbance, \theta_{noise}, to motor command, c, using an approach similar to the above approach used to determine G_{disturb2\theta_a}(s). We will refer to this transfer function as G_{noise2c} (s).

Checkoff 5:

What are G_{disturb2\theta_a}(s) and G_{noise2c}(s).
Modify the python script to plot the frequence responses associated with these two disturbance transfer functions, but please use only the Bode function, NOT BodeF. Disturbance transfer functions often vary by orders of magnitude over frequencies of interest, so are best respresented in using log magnitude.
Show that your model can predict the most disturbing frequency by measuring the most disturbing fanning frequency for K_p = 3 and K_d set low enough that the step response oscillates about seven times, and then compare that frequency to your model's prediction.
What does G_{noise2c}(j\omega) tell you about minimizing command noise? Which causes more command noise, low frequency or high frequency measurement noise? Note that the propeller arm system samples at 1000 times per second, so what do you think is the highest noise frequency of interest?
What can you say about the disturbance frequency response when K_i > 0? What does that tell you about steady-state errors?

Checkoff 6 - Design a less disturb-able controller

Now that you have a calibrated numerical model for your system, and can use the model to plot G_{d2t}(j\omega) (to within a scale factor), and you can pick K_p, K_i and K_d to minimize the maximum disturbance response. Please do so, and show that your new controller is less "disturb-able". How does command noise change with changing K_d? What did G_{noise2c}(j\omega) for your less disturbable controller tell you about command noise versus K_d? Does that match what you see when you plot the motor commands?

Rather than setting the three gains using the serial plotter send window, you can "preset" the gains in the Teensy file. For example, if you want to preset K_d = 0.2, K_i=0, and K_p = 2.0, you would set offsetKi, offsetKp, and offsetKd as shown in the figure below.

Checkoff 6:

Use your calibrated model to plot the disturbance frequency response for your disturbance rejecting controller and compare it to the disturbance frequency response for the K_i = 0, K_p = 3.0 case from checkoff 5. (use the value you found for K_d in checkoff five).
Demonstrate your disturbance-rejecting controller on the arm by showing limited arm motion due to fanning disturbances at a few different fanning frequencies.
If you double both K_p and K_d for your disturbance-rejecting controller, how do your measured and modeled disturbance frequency response change? What about the command noise?
Try tapping your arm with a stick for both your disturbance-rejecting controller and the controller with double K_p and K_d. What differences do you notice?

Postlab - Bury the Lead

As we have seen above, we have accurate tracking and good disturbance rejection at frequencies for which K(j\omega)H(j\omega) >> 1. If we are only concerned with very low frequencies, or steady-state, this can usually be achieved by setting K_i > 0 in a PID controller. For this postlab, however, we will focus on more dynamic tracking and disturbance rejection (as is common in high performance systems). Our propeller arm test case input will be a 0.5 hertz square wave with amplitude 0.05 radians, and our "disturbing" test will be removing two bricks from ONE SIDE of the arm.

In the Teensy sketch linked below, we have used the presets to set up the above test case using a PD controller with K_p = 1 and K_d=0.2. You can download the sketch, compile it, run it on your propeller arm (after updating the angle offset and motor directions), and verify that the controller performs well. We ask that you then "disturb" the propeller arm by removing two bricks from one side of the arm, and examine the angle error caused by this disturbance.

You will need a better controller to reduce the angle error, and designing that controller is the topic of this postlab. In particular, you will start by trying to use a PD controller by increasing K_p to 5 and increasing K_d just enough to make the arm response reasonably stable. Seems simple enough, but there will be trouble, as you will discover when you plot the motor command and see all the noise.

Lead controllers can reduce command noise IF their maximum high frequency gain is low enough. For our arm, limiting that gain to 150 is certain to reduce command noise, but adhering to that constraint will also make it difficult to match the stability of the PD controller. But please try! You can then try designing a lead-squared controller, to see if you can get better stability while still adhering to the high frequency gain limit of 150.

We have provided a python script that will help you design your lead and lead squared controllers using phase margin and to help you determine how to map your CT lead controller into a DT controller, implement that on the Teensy microcontroller, and test its performance. Please be sure you understand the last two questions on prelab 3, particularly ( compensation ), BEFORE TRYING THIS PART OF THE LAB.

Lead Controller

If we use a PD controller, then

K(s) = K_p + s K_d,

with the associated frequency response

K(j\omega) = K_p + j\omega K_d,

where \omega is frequency in radians/second, and K_p and K_d are the to-be-designed proportional and derivative gains respectively.

Alternatively, if we are using a lead controller then,

K(s) = K_p \frac{s_p}{s_z} \left( \frac{s - s_z}{s- s_p} \right),

where s_z < 0 , s_p < 0 and K_p > 0 are to be designed. As we will see below, some of our analyses will simplify if we substitute pzratio*s_z for s_p, leading to

K(s) = K_p pzratio \left( \frac{s - s_z}{s- pzratio*s_z} \right),

with associated frequency response

K(j\omega) = K_p pzratio \left( \frac{j\omega - s_z}{j\omega- pzratio*s_z} \right),

where again s_z < 0 , and pzratio > 0 (usually between 10 and 100), and K_p > 0 are the to-be-designed zero, pole-zero ratio, and gain of the lead.

A lead squared controller would be

K_p \left[pzratio_2 \left( \frac{j\omega - s_{z2}}{j\omega- pzratio_2*s_{z2}} \right)\right]^2,

where we have used pzratio_2 and s_{z2} to clarify that the best lead squared and best lead controllers will likely have different zeros and ratios.

Finally, please determine the high frequency limit of the magnitude of the transfer function from measurement noise to motor command for the lead and lead squared controllers. To ensure low motor command noise,

limit_{\omega \rightarrow \infty} |G_{noise2c}(j\omega)| \lt 150!!

Obviously (is it?) the above limit is violated when we use a PD controller. If our controller is a lead controller, the above constraint implies that Kp*pzratio < 150, do you see why? What if our controller is a lead-squared controller?

DT Mapping

Our controller determines the Teensy command, c(t), from measurements of the separation distance error, e(t). If we use a PD controller, then

c(t) = K_p e(t) + K_d \frac{d e(t)}{dt},

where K_p and K_d are the proportional and derivative gains respectively. But since the Teensy operates on a clock, and only samples y(t) and updates c(t) once every \Delta T, we can approximate the continuous-time controller with a difference equation,

c(n \Delta T) = c[n] \approx K_p e[n] + K_d \frac{e[n] - e[n-1]}{\Delta T},

where \Delta T is the Teensy update period (1 millisecond for the propeller arm).

Alternatively, if we are using a lead controller then

\frac{d c(t) }{dt} - pzratio*s_z*c(t) = K_0*pzratio*\left( \frac{d e(t)}{dt} - s_z e(t) \right)

where s_z < 0 , s_p < 0 , and K_0 > 0 are the to-be-designed zero, pole, and steady-state gain of the lead.

You can determine the difference equation for c[n] in the case of a lead controller, using the same approximation applied to deriving the difference equation for the PD controller. The python script calculates those coefficients for you, and prints them out. And, if you uncomment the last four lines of that python script, you can get the difference equation coefficients calculated with a more accurate method (and you can use the same function for your lead squared controller). PLEASE BE SURE YOU COPY AND PASTE AT LEAST 5 DIGITS PAST THE DECIMAL POINT! There is detailed cancellation in these numbers!!!

Note also that the more accurate coefficients are printed in the z-transform format, and then in a format that gives you more digits of accuracy. The translation to a difference equation is NOT a deep idea worth mulling over, if it is unclear to you, just follow the suggestions in the comments and look it up on-line.

Implementation

Line 139 in the Teensy sketch below,

dCmd= 0.0*pastDCmd + 0.0*Err - 0.0*pastErr;  // FIX THIS!!!

computes dCmd (which is the Teensy name for c[n]) from pastDCmd (which is the Teensy name for c[n-1], Err (which is the Teensy name for e[n]), and pastErr (which is the Teensy name for e[n-1]). To change the lead controller, you need only replace the three coefficients with values determined from mapping your CT lead controller to DT (printed out by the python script below).

For the lead squared controller, you will update line 141, using coefficients generated by the python script.

Checkoff 7

PLEASE USE THE TEENSY AND PYTHON FILES FOR THIS POSTLAB:

POSTLAB Links

Sketch for the POSTLAB is here.
IPython Arm Model File for the POSTLAB here.

In the above python script and in the above Teensy sketch, we have implemented almost everything you need to use phase margin ideas to design PD and lead controllers. If you run the ctPropArmClass.py file in ipython, it produces three plots (though they may be on top of each other): frequency responses for open-loop controller transfer functions, step responses, and closed loop transfer functions (see the image below).

If you download and run the 6310_PROP_26B (AFTER you UPDATE the angleoffset, line 13, and FIX the motor directions, lines 153 and 155) and open the Teensy serial plotter, you should see arm step responses from a PD controller with K_p = 1.

Lines that need to be changed in the python script and Teensy sketch are labeled "FIX THIS".

The postlab, step-by-step.

Please

Download the lead postlab Teensy sketch using the above link, fix its angle offset and motor directions, upload it to your arm, and verify the controller is working well. You should see plots like in the images below (first with the commands plotted and then without).

Remove two bricks from ONE side of your propeller arm, like in the image below.

Examine the behavior of your disturbed propeller arm with the default PD controller (K_p = 1.) You should see one of the two plots below, depending on whether you removed the two bricks from the right or left side.

Set K_p=5 and K_d = 0.6 in the Teensy presets (lines 34 and 33), and reupload your sketch. Adjust K_d if your arm oscillates too much. You should see a result like the one below. Although the arm angle and desired value do not match perfectly, raising K_p from one to five substantially improved the match.

Update the python model with your arm's Gamma and Beta, and the value you found for Kd. Try running the python script, and compare its predicted step response for the PD-controlled arm to the behavior of your undisturbed PD-controlled arm, i.e., with three legos on each arm.
Optimize sz in the python script to get the best phase margin for a lead controller (the phase margin is printed at the top of the bode plot), but DO NOT INCREASE the maximum value of the magnitude of the measurement noise to command transfer function (so do not change pzratio)! BE SURE TO ADD THE LEAD CONTROLLED SYSTEM TO THE STEP RESPONSE PLOTS!
Try your lead controller by updating line 139 in the Teensy sketch with the coefficients printed by the python script, then upload the sketch to your arm, and type 7 1 in the send window to run the lead controller (type 7 0 in the send window to switch back to PD controller). When running the lead controller, be sure to examine the command noise, is it lower than the command noise when running the PD controller? Is your lead controller as stable as your PD controller? Is it possible to make your lead controller as stable as the PD controller given the limit on pzratio?
Uncomment the very last four lines of the python script, rerun the script, and examine the coefficients that the script prints out (check the ct.c2d documentation online). They should the same order of magnitude as the difference equation coefficients computed using the difference approximation and are printed multiple ways to make sure that you can see all the digits of these numbers. Try using them to replace the coefficients on line 139 in your Teensy sketch. BE SURE THAT YOU COPY THE COEFFICIENTS WITH AT LEAST 5 DIGITS PAST THE DECIMAL POINT
Design a lead-squared controller using the python script. Square the lead transfer function, but adjust the pzRatio so you do not increase the maximum magnitude of the transfer function relating measurement noise to command (maximum magnitude should be less than 150!!). Be sure to add your lead squared results to the step response and frequency response plots in the script.
Implement the lead squared controller using the coefficients computed by ct.c2d by editing line 141 in the Teensy sketch. Upload the sketch, then type 7 2 in the send window if you want to try the lead squared controller.
Compare the controllers. Our lead squared controller results are below. Note the improved match between measured and desired angle, and the reduction in command noise.

Checkoff 7:

YOU SHOULD BE ABLE TO DEMONSTRATE YOUR WORKING LEAD SQUARED CONTROLLER AND BE ABLE TO ANSWER ALL THE CHECKOFF QUESTIONS IN ORDER TO GET THE CHECKOFF!!!

What were your values for Kd, sz, and pzRatio for your PD, lead and lead squared controllers when Kp=5?
How did you use the bode plots to design your controllers?
Was the lead squared controller better? Why?