Controls Engineering in FRC
Controls Engineering in FRC
Tyler Veness
Controls Engineering in the
FIRST Robotics Competition
Graduate-level control theory for high schoolers
Tyler Veness
Copyright © 2017-2020 Tyler Veness
HTTPS://GITHUB.COM/CALCMOGUL/CONTROLS-ENGINEERING-IN-FRC
Licensed under the Creative Commons Attribution-ShareAlike 4.0 International License (the “Li-
cense”). You may not use this file except in compliance with the License. You may obtain a copy
of the License at https://creativecommons.org/licenses/by-sa/4.0/. Unless required by applicable
law or agreed to in writing, software distributed under the License is distributed on an “AS IS” BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the
specific language governing permissions and limitations under the License.
Generated from commit 5132b8b made on April 23, 2020. The latest version can be downloaded
from https://file.tavsys.net/control/controls-engineering-in-frc.pdf.
Tree next to Oakes Circle at UCSC
Contents
Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
2 PID controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1 Proportional term 11
2.2 Derivative term 12
2.3 Integral term 14
2.4 PID controller definition 16
2.5 Response types 16
2.6 Manual tuning 17
2.7 Actuator saturation 17
2.8 Limitations 18
3 Transfer functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1 Laplace transform 19
3.2 Parts of a transfer function 19
3.2.1 Poles and zeroes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.2 Nonminimum phase zeroes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.3 Pole-zero cancellation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.3 Transfer functions in feedback 22
3.3.1 Root locus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4 Linear algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.1 Vectors 27
4.2 Linear combinations, span, and basis vectors 27
4.3 Linear transformations and matrices 27
4.4 Matrix multiplication as composition 27
4.5 The determinant 28
4.6 Inverse matrices, column space, and null space 28
4.7 Nonsquare matrices as transformations between dimensions 28
4.8 Eigenvectors and eigenvalues 28
4.9 Miscellaneous notation 28
5 State-space controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5.1 From PID control to model-based control 29
5.2 What is a dynamical system? 30
5.3 State-space notation 30
5.3.1 What is state-space? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.3.2 Benefits over classical control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.3.3 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.4 Controllability 32
5.5 Observability 32
5.6 Closed-loop controller 33
5.7 Pole placement 34
5.8 Linear-quadratic regulator 34
5.8.1 Bryson’s rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
5.8.2 Pole placement vs LQR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.9 Model augmentation 37
5.9.1 Plant augmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.9.2 Controller augmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.9.3 Observer augmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.9.4 Output augmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.9.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.10 Feedforward 39
5.10.1 Plant inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.10.2 Unmodeled dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.11 Integral control 42
5.11.1 Plant augmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.11.2 Input error estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
6 Digital control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
6.1 Phase loss 45
6.2 s-plane to z-plane 46
6.2.1 z-plane stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
6.2.2 z-plane behavior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
6.2.3 Nyquist frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
6.3 Discretization methods 48
6.4 Effects of discretization on controller performance 50
6.5 Matrix exponential 52
6.6 Taylor series 53
6.7 Zero-order hold for state-space 54
7 Nonlinear control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
7.1 Introduction 57
7.2 Linearization 57
7.3 Lyapunov stability 58
7.4 Affine systems 59
7.4.1 Feedback linearization for reference tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.5 Further reading 60
8 State-space applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
8.1 Elevator 61
8.1.1 Continuous state-space model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
8.1.2 Model augmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
8.1.3 Gravity feedforward . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
8.1.4 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
8.1.5 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
8.2 Flywheel 63
8.2.1 Continuous state-space model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
8.2.2 Model augmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
8.2.3 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
8.2.4 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
8.2.5 Flywheel model without encoder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
8.2.6 Voltage compensation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
8.3 Single-jointed arm 66
8.3.1 Continuous state-space model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
8.3.2 Model augmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
8.3.3 Gravity feedforward . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
8.3.4 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
8.3.5 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
8.4 Pendulum 69
8.4.1 State-space model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
8.5 Differential drive 70
8.5.1 Velocity subspace state-space model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
8.5.2 Linear time-varying model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
8.5.3 Improving model accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
8.5.4 Cross track error controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
8.5.5 Explicit time-varying control law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
8.5.6 Nonlinear observer design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
8.6 Ramsete unicycle controller 82
8.6.1 Velocity and turning rate command derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
8.6.2 Nonlinear controller equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
8.6.3 Linear reference tracker . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
8.7 Linear time-varying unicycle controller (cascaded) 87
8.7.1 Explicit time-varying control law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
IV System modeling
11 Calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
11.1 Derivatives 131
11.1.1 Power rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
11.1.2 Product rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
11.1.3 Chain rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
11.2 Integrals 132
11.2.1 Change of variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
11.3 Tables 133
11.3.1 Common derivatives and integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
12 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
12.1 Linear motion 135
12.2 Angular motion 135
12.3 Vectors 136
12.3.1 Basic vector operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
12.3.2 Cross product . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
12.4 Curvilinear motion 137
12.4.1 Differential drive . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
12.4.2 Mecanum drive . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
12.4.3 Swerve drive . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
V Motion planning
VI Appendices
H Derivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 215
H.1 Transfer function in feedback 215
H.2 Zero-order hold for state-space 216
H.3 Kalman filter as Luenberger observer 217
H.3.1 Luenberger observer with separate prediction and update . . . . . . . . . . . . . . . . . . . 218
H.4 Trapezoidal motion profile 218
Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
Online 221
Miscellaneous 222
Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223
Flower bush near Spencer’s Market in Santa Maria, CA
Preface
Motivation
I am the software mentor for a FIRST Robotics Competition (FRC) team. My responsibilities for
that include teaching programming, software engineering practices, and applications of control theory.
The curriculum I developed as of the spring of 2017 (located at https://csweb.frc3512.com/ci/) taught
rookies enough to be minimally competitive, but it provided no formal avenues of growth for veteran
students.
Also, out of a six week build season, the software team usually only got a few days with the completed
robot due to poor build schedule management. This led to two problems. First, two days is only
enough time to verify basic software functionality, not test and tune feedback controllers. Second,
this was also the first time the robot’s electromechanical systems have been tested after integration, so
any issues that arose consumed valuable software integration time while the team traced the problem
to a mechanical, electrical, or software cause.
This book expands my curriculum to cover control theory topics I learned in my graduate-level en-
gineering classes at University of California, Santa Cruz. It introduces state-space controllers and
serves as a practical guide for formulating and implementing them. Since state-space control utilizes
a system model, both problems mentioned earlier can be addressed. Basic software functionality
can be tested against it and feedback controllers can be tuned automatically based on system con-
straints. This allows software teams to test their work much earlier in the build season in a controlled
environment as well as save time during feedback controller design, implementation, and testing.
I originally started writing this book because I felt that the topic wasn’t difficult, but the information for
it wasn’t easily accessible. When I was a high school robotics team member, I had to learn everything
from scratch through various internet sources and eventually from college courses as part of my
bachelor’s degree. Control theory has a certain beauty to it, and I want more people to appreciate it
the way I do. Therefore, my goal is to make the learning process quicker and easier for future team
members by collating all the relevant information.
xii
Intended audience
This guide is intended to make an advanced engineering topic approachable so it can be applied by
those who aren’t experts in control theory. My intended audience is high school students who are
members of a FIRST Robotics Competition team. As such, they will likely have passing knowledge
of PID control and have basic proficiency in programming. This guide will expand their incomplete
understanding of control theory to include the fundamentals of classical control theory, enough linear
algebra to understand the notation and mechanics of modern control, and finally how to apply modern
control to design challenges they regularly see in their FRC robots from year to year.
Acknowledgements
I would like to thank my controls engineering instructors Dejan Milutinović and Gabriel Elkaim of
University of California, Santa Cruz. They taught their classes from a pragmatic perspective focused
on application and intuition that I appreciated. I would also like to thank Dejan Milutinović for
introducing me to the field of control theory and both of my instructors for teaching me what it
means to be a controls engineer.
Thanks to Austin Schuh from FRC team 971 for providing the final continuous state-space models
used in the examples section.
Trees by Baskin Engineering building at UCSC
0.1 Prerequisites
Knowledge of basic algebra and complex numbers is assumed. Some introductory physics and cal-
culus will be taught as necessary.
Part III, “Estimation and localization,” introduces the field of stochastic control theory. The Luen-
berger observer and the probability theory behind the Kalman filter is taught with several examples
of creative applications of Kalman filter theory.
Part IV, “System modeling,” introduces the basic calculus and physics concepts required to derive
the models used in the previous chapters. It walks through the derivations for several common FRC
subsystems. Then, methods for system identification are discussed for empirically measuring model
parameters.
Part V, “Motion planning,” covers planning how the robot will get from its current state to some
desired state in a manner achievable by its dynamics. It introduces motion profiles with one degree
of freedom for simple maneuvers. Trajectory optimization methods are presented for generating
profiles with higher degrees of freedom.
The appendices provide further enrichment that isn’t required for a passing understanding of the
material. This includes derivations for many of the results presented and used in the mainmatter of
the book.
The Python scripts used to generate the plots in the case studies double as reference implementa-
tions of the techniques discussed in their respective chapters. They are available in this book’s Git
repository. Its location is listed on the copyright page.
2 PID controllers . . . . . . . . . . . . . . . . . . . . . . 11
2.1 Proportional term
2.2 Derivative term
2.3 Integral term
2.4 PID controller definition
2.5 Response types
2.6 Manual tuning
2.7 Actuator saturation
2.8 Limitations
3 Transfer functions . . . . . . . . . . . . . . . . . . 19
3.1 Laplace transform
3.2 Parts of a transfer function
3.3 Transfer functions in feedback
This page intentionally left blank
Road near walking trail off of Rice Ranch Road in Santa Maria, CA
Control systems are all around us and we interact with them daily. A small list of ones you may have
seen includes heaters and air conditioners with thermostats, cruise control and the anti-lock braking
system (ABS) on automobiles, and fan speed modulation on modern laptops. Control systems monitor
or control the behavior of systems like these and may consist of humans controlling them directly
(manual control), or of only machines (automatic control).
How can we prove closed-loop controllers on an autonomous car, for example, will behave safely
and meet the desired performance specifications in the presence of uncertainty? Control theory is an
application of algebra and geometry used to analyze and predict the behavior of systems, make them
respond how we want them to, and make them robust to disturbances and uncertainty.
Controls engineering is, put simply, the engineering process applied to control theory. As such,
it’s more than just applied math. While control theory has some beautiful math behind it, controls
engineering is an engineering discipline like any other that is filled with trade-offs. The solutions
control theory gives should always be sanity checked and informed by our performance specifications.
We don’t need to be perfect; we just need to be good enough to meet our specifications (see section
0.4 for more on engineering).
1.1 Nomenclature
Most resources for advanced engineering topics assume a level of knowledge well above that which is
necessary. Part of the problem is the use of jargon. While it efficiently communicates ideas to those
within the field, new people who aren’t familiar with it are lost. See the glossary for a list of words
and phrases commonly used in control theory, their origins, and their meaning. Links to the glossary
are provided for certain words throughout the book and will use this color.
The system or collection of actuators being controlled by a control system is called the plant. A con-
troller is used to drive the plant from its current state to some desired state (the reference). Controllers
which don’t include information measured from the plant’s output are called open-loop controllers.
8 Chapter 1. Control system basics
Controllers which incorporate information fed back from the plant’s output are called closed-loop
controllers or feedback controllers. Figure 1.1 shows a plant in feedback with a controller.
Note that the input and output of a system are defined from the plant’s point of view. The negative
feedback controller shown is driving the difference between the reference and output, also known as
the error, to zero.
t K t
+
input open-loop output
∓
feedback
The open-loop gain is the total gain from the sum node at the input (the circle) to the output branch.
1.4 Why feedback control? 9
This would be the system’s gain if the feedback loop was disconnected. The feedback gain is the
total gain from the output back to the input sum node. A sum node’s output is the sum of its inputs.
Figure 1.4 is a block diagram with more formal notation in a feedback configuration.
+
X(s) P1 Y (s)
∓
P2
2. PID controllers
The PID controller is a commonly used feedback controller consisting of proportional, integral, and
derivative terms, hence the name. This chapter will build up the definition of a PID controller term
by term while trying to provide intuition for how each of them behaves.
First, we’ll get some nomenclature for PID controllers out of the way. The reference is called the
setpoint (the desired position) and the output is called the process variable (the measured position).
Below are some common variable naming conventions for relevant quantities.
where Kp is the proportional gain and e(t) is the error at the current time t.
Proportional gains act like “software-defined springs” that pull the system toward the desired posi-
tion. Recall from physics that we model springs as F = −kx where F is the force applied, k is
a proportional constant, and x is the displacement from the equilibrium point. This can be written
another way as F = k(0 − x) where 0 is the equilibrium point. If we let the equilibrium point be
our feedback controller’s setpoint, the equations have a one-to-one correspondence.
F = k(r − x)
u(t) = Kp e(t) = Kp (r(t) − y(t))
so the “force” with which the proportional controller pulls the system’s output toward the setpoint is
proportional to the error, just like a spring.
where Kp is the proportional gain, Kd is the derivative gain, and e(t) is the error at the current
time t.
Kp e(t)
+ e(t) + u(t) y(t)
r(t) Plant
− +
Kd de(t)
dt
A PD controller has a proportional controller for position (Kp ) and a proportional controller for
velocity (Kd ). The velocity setpoint is implicitly provided by how the position setpoint changes over
time. Figure 2.3 shows an example for an elevator.
To prove a PD controller is just two proportional controllers, we will rearrange the equation for a PD
controller.
ek − ek−1
u k = Kp e k + Kd
dt
where uk is the control input at timestep k and ek is the error at timestep k. ek is defined as ek =
rk − xk where rk is the setpoint and xk is the current state at timestep k.
r −r x −x
Notice how k dtk−1 is the velocity of the setpoint. By the same reasoning, k dtk−1 is the system’s
velocity at a given timestep. That means the Kd term of the PD controller is driving the estimated
velocity to the setpoint velocity.
If the setpoint is constant, the implicit velocity setpoint is zero, so the Kd term slows the system
down if it’s moving. This acts like a “software-defined damper”. These are commonly seen on door
closers, and their damping force increases linearly with velocity.
14 Chapter 2. PID controllers
where Kp is the proportional gain, Ki is the integral gain, e(t) is the error at the current time t,
and τ is the integration variable.
The integral integrates from time 0 to the current time t. We use τ for the integration because we
need a variable to take on multiple values throughout the integral, but we can’t use t because we
already defined that as the current time.
Kp e(t)
+ e(t) + u(t) y(t)
r(t) Plant
− +
∫t
Ki 0
e(τ ) dτ
When the system is close to the setpoint in steady-state, the proportional term may be too small to pull
the output all the way to the setpoint, and the derivative term is zero. This can result in steady-state
error, as shown in figure 2.5.
A common way of eliminating steady-state error is to integrate the error and add it to the control
input. This increases the control effort until the system converges. Figure 2.5 shows an example
of steady-state error for a flywheel, and figure 2.6 shows how an integrator added to the flywheel
controller eliminates it. However, too high of an integral gain can lead to overshoot, as shown in
figure 2.7.
R There are better approaches to fix steady-state error like using feedforwards or constraining
when the integral control acts using other knowledge of the system. We will discuss these in
more detail when we get to modern control theory.
2.3 Integral term 15
Figure 2.6: PI controller on a flywheel without Figure 2.7: PI controller on a flywheel with over-
steady-state error shoot from large Ki gain
16 Chapter 2. PID controllers
where Kp is the proportional gain, Ki is the integral gain, Kd is the derivative gain, e(t) is the
error at the current time t, and τ is the integration variable.
Figure 2.8 shows a block diagram for a system controlled by a PID controller.
Kp e(t)
Kd de(t)
dt
R Note: Adding an integral gain to the controller is an incorrect way to eliminate steady-state
error. A better approach would be to tune it with an integrator added to the plant, but this
requires a model. Since we are doing output-based rather than model-based control, our only
option is to add an integrator to the controller.
Beware that if Ki is too large, integral windup can occur. Following a large change in setpoint, the
integral term can accumulate an error larger than the maximal control input. As a result, the system
overshoots and continues to increase until this accumulated error is unwound.
We’ll try to explain this through a bit of math. Let’s say we have a controller u = k(r − x) where
u is the control effort, k is the gain, r is the reference, and x is the current state. Let umax be the
limit of the actuator’s output which is less than the uncapped value of u and kmax be the associated
maximum gain. We will now compare the capped and uncapped controllers for the same reference
and current state.
umax < u
kmax (r − x) < k(r − x)
kmax < k
For the inequality to hold, kmax must be less than the original value for k. This reduced gain is
evident in a system response when there is a linear change in state instead of an exponential one as it
approaches the reference. This is due to the control effort no longer following a decaying exponential
plot. Once the system is closer to the reference, the controller will stop saturating and produce realistic
controller values again.
2.8 Limitations
PID’s heuristic method of tuning is a reasonable choice when there is no a priori knowledge of the
system dynamics. However, controllers with much better response can be developed if a dynamical
model of the system is known. Furthermore, PID only applies to single-input, single-output (SISO)
systems; we’ll cover methods for multiple-input, multiple-output (MIMO) control in part II of this
book.
Treeline by Crown/Merril bus stop at UCSC
3. Transfer functions
This chapter briefly discusses what transfer functions are, how the locations of poles and zeroes affect
system response and stability, and how controllers affect pole locations. As long as you gain under-
standing of those concepts, don’t worry too much about not being able to follow the math presented
here; we won’t use transfer functions in modern control theory (it’s mainly provided for complete-
ness). This chapter is intended to provide a framework within which to understand results from the
mathematical machinery of modern control as well as vocabulary to communicate that understand-
ing.
function are called poles because they make the transfer function approach infinity; on a 3D graph,
these look like the poles of a circus tent (see figure 3.1).
When the factors of the denominator are broken apart using partial fraction expansion into something
like s+a
A B
+ s+b , the constants A and B are called residues, which determine how much each pole
contributes to the system response.
The factors representing poles are each the Laplace transform of a decaying exponential1 . That means
the time domain responses of systems comprise decaying exponentials (e.g., y = e−t ).
R Imaginary poles and zeroes always come in complex conjugate pairs (e.g., −2 + 3i, −2 − 3i).
The locations of the closed-loop poles in the complex plane determine the stability of the system.
Each pole represents a frequency mode of the system, and their location determines how much of
each response is induced for a given input frequency. Figure 3.2 shows the impulse responses in the
time domain for transfer functions with various pole locations. They all have an initial condition of
1.
Location Stability
Left Half-plane (LHP) Stable
Imaginary axis Marginally stable
Right Half-plane (RHP) Unstable
When a system is stable, its output may oscillate but it converges to steady-state. When a system is
1
We are handwaving Laplace transform derivations because they are complicated and neither relevant nor useful.
3.2 Parts of a transfer function 21
Im(jω)
Stable Unstable
t
t
t
LHP RHP
t
t t
Re(σ)
t t t
marginally stable, its output oscillates at a constant amplitude forever. When a system is unstable, its
output grows without bound.
moving to the zero placed next to it. You have the same dynamics as before, but the pole is also stuck
where it is no matter how much feedback gain is applied. For an attempted nonminimum phase zero
cancellation, you have effectively placed an unstable pole that’s unobservable. This means the system
will be going unstable and blowing up, but you won’t be able to detect this and react to it.
Keep in mind when making design decisions that the model likely isn’t perfect. The whole point of
feedback control is to be robust to this kind of uncertainty.
+
X(s) K G Y (s)
−
The transfer function of figure 3.3, a control system diagram with feedback, from input to output is
Y (s) KG
Gcl (s) = = (3.2)
X(s) 1 + KGH
The numerator is the open-loop gain and the denominator is one plus the gain around the feedback
loop, which may include parts of the open-loop gain (see appendix H.1 for a derivation). As another
example, the transfer function from the input to the error is
E(s) 1
Gcl (s) = = (3.3)
X(s) 1 + KGH
The roots of the denominator of Gcl (s) are different from those of the open-loop transfer function
KG(s). These are called the closed-loop poles.
Figure 3.4: Root locus showing pole moving to- Figure 3.5: Root locus showing poles moving to-
ward negative infinity ward asymptotes
If poles are much farther left in the LHP than the typical system dynamics exhibit, they can be con-
sidered negligible. Every system has some form of unmodeled high frequency, nonlinear dynamics,
but they can be safely ignored depending on the operating regime.
To demonstrate this, consider the transfer function for a second-order DC brushed motor (a CIM
motor) from voltage to position
K
G(s) =
s((Js + b)(Ls + R) + K 2 )
Figure 3.7: Root locus of second-order DC Figure 3.8: Root locus of first-order DC brushed
brushed motor plant motor plant
Figure 3.9: Step response of second-order DC Figure 3.10: Step response of first-order DC
brushed motor plant brushed motor plant
frequency stable poles decays rapidly. Unstable poles, on the other hand, represent unstable dynamics
which cause the system output to grow to infinity. Regardless of how slow these unstable dynamics
are, they will eventually dominate the response.
II
Modern control theory
4 Linear algebra . . . . . . . . . . . . . . . . . . . . . . 27
4.1 Vectors
4.2 Linear combinations, span, and basis vectors
4.3 Linear transformations and matrices
4.4 Matrix multiplication as composition
4.5 The determinant
4.6 Inverse matrices, column space, and null space
4.7 Nonsquare matrices as transformations between dimen-
sions
4.8 Eigenvectors and eigenvalues
4.9 Miscellaneous notation
5 State-space controllers . . . . . . . . . . . . . . 29
5.1 From PID control to model-based control
5.2 What is a dynamical system?
5.3 State-space notation
5.4 Controllability
5.5 Observability
5.6 Closed-loop controller
5.7 Pole placement
5.8 Linear-quadratic regulator
5.9 Model augmentation
5.10 Feedforward
5.11 Integral control
6 Digital control . . . . . . . . . . . . . . . . . . . . . . 45
6.1 Phase loss
6.2 s-plane to z-plane
6.3 Discretization methods
6.4 Effects of discretization on controller performance
6.5 Matrix exponential
6.6 Taylor series
6.7 Zero-order hold for state-space
7 Nonlinear control . . . . . . . . . . . . . . . . . . . 57
7.1 Introduction
7.2 Linearization
7.3 Lyapunov stability
7.4 Affine systems
7.5 Further reading
8 State-space applications . . . . . . . . . . . . 61
8.1 Elevator
8.2 Flywheel
8.3 Single-jointed arm
8.4 Pendulum
8.5 Differential drive
8.6 Ramsete unicycle controller
8.7 Linear time-varying unicycle controller (cascaded)
This page intentionally left blank
Grass clearing by Interdisciplinary Sciences building and Thimann Labs at UCSC
4. Linear algebra
Modern control theory borrows concepts from linear algebra. At first, linear algebra may appear very
abstract, but there are simple geometric intuitions underlying it. First, watch 3Blue1Brown’s preview
video for the Essence of linear algebra video series (5 minutes) [4]. The goal here is to provide an
intuitive, geometric understanding of linear algebra as a method of linear transformations.
We would normally include written material here for learning linear algebra, but 3Blue1Brown’s
animated videos are better at conveying the geometric intuition involved than anything we could
write here with static text. Instead of making an inferior copy of his work, we’ll provide bibliography
entries for a selection of his videos.
4.1 Vectors
Watch the “Vectors, what even are they?” video from 3Blue1Brown’s Essence of linear algebra series
(5 minutes) [11].
The matrix denoted by 0m×n is a matrix filled with zeroes with m rows and n columns.
The T in AT denotes transpose, which flips the matrix across its diagonal such that the rows become
columns and vice versa.
The † in B† denotes the Moore-Penrose pseudoinverse given by B† = (BT B)−1 BT . The pseudoin-
verse is used when the matrix is nonsquare and thus not invertible to produce a close approximation
of an inverse in the least squares sense.
Night sky above Dufour Street in Santa Cruz, CA
5. State-space controllers
R Chapters from here on use the frccontrol Python package to demonstrate the concepts dis-
cussed and perform the complex math required. See appendix D for how to install it.
When we want to command a system to a set of states, we design a controller with certain control laws
to do it. PID controllers use the system outputs with proportional, integral, and derivative control
laws. In state-space, we also have knowledge of the system states so we can do better.
Modern control theory uses state-space representation to model and control systems. State-space
representation models systems as a set of state, input, and output variables related by first-order
differential equations that describe how the system’s state changes over time given the current states
and inputs.
where ω is the angular velocity and ωmax is the maximum angular velocity. If DC brushed motors
are said to behave linearly, then why is this?
Linearity refers to a system’s equations of motion, not its time domain response. The equation defin-
ing the motor’s change in angular velocity over time looks like
ω̇ = −aω + bV
where ω̇ is the derivative of ω with respect to time, V is the input voltage, and a and b are constants
specific to the motor. This equation, unlike the one shown before, is actually linear because it only
consists of multiplications and additions relating the input V and current state ω.
Also of note is that the relation between the input voltage and the angular velocity of the output
shaft is a linear regression. You’ll see why if you model a DC brushed motor as a voltage source
and generator producing back-EMF (in the equation above, bV corresponds to the voltage source
and −aω corresponds to the back-EMF). As you increase the input voltage, the back-EMF increases
linearly with the motor’s angular velocity. If there was a friction term that varied with the angular
velocity squared (air resistance is one example), the relation from input to output would be a curve.
Friction that scales with just the angular velocity would result in a lower maximum angular velocity,
but because that term can be lumped into the back-EMF term, the response is still linear.
for linear systems. Including nonzero initial conditions complicates the algebra even more. State-
space representation uses the time domain instead of the Laplace domain, so it can model nonlinear
systems1 and trivially supports nonzero initial conditions.
If modern control theory is so great and classical control theory isn’t needed to use it, why learn
classical control theory at all? We teach classical control theory because it provides a framework
within which to understand results from the mathematical machinery of modern control as well as
vocabulary with which to communicate that understanding. For example, faster poles (poles moved
to the left in the s-plane) mean faster decay, and oscillation means there is at least one pair of complex
conjugate poles. Not only can you describe what happened succinctly, but you know why it happened
from a theoretical perspective.
5.3.3 Definition
Below are the continuous and discrete versions of state-space notation.
ẋ = Ax + Bu (5.1)
y = Cx + Du (5.2)
In the continuous case, the change in state and the output are linear combinations of the state vector
and the input vector. The A and B matrices are used to map the state vector x and the input vector
u to a change in the state vector ẋ. The C and D matrices are used to map the state vector x and
the input vector u to an output vector y.
1
This book focuses on analysis and control of linear systems. See chapter 7 for more on nonlinear control.
32 Chapter 5. State-space controllers
5.4 Controllability
State controllability implies that it is possible – by admissible inputs – to steer the states from any
initial value to any final value within some finite time window.
rank B AB A2 B · · · An−1 B =n (5.5)
where rank is the number of linearly independent rows in a matrix and n is the number of state
variables.
The matrix in equation (5.5) being rank-deficient means the inputs cannot apply transforms along all
axes in the state-space; the transformation the matrix represents is collapsed into a lower dimension.
(C)
The condition number of the controllability matrix C is defined as σσmax
min (C)
where σmax is the max-
imum singular value2 and σmin is the minimum singular value. As this number approaches infinity,
one or more of the states becomes uncontrollable. This number can also be used to tell us which actu-
ators are better than others for the given system; a lower condition number means that the actuators
have more control authority.
5.5 Observability
Observability is a measure for how well internal states of a system can be inferred by knowledge of
its external outputs. The observability and controllability of a system are mathematical duals (i.e., as
controllability proves that an input is available that brings any initial state to any desired final state,
observability proves that knowing enough output values provides enough information to predict the
initial state of the system).
C
CA
rank .. = n (5.6)
.
CAn−1
where rank is the number of linearly independent rows in a matrix and n is the number of state
variables.
The matrix in equation (5.6) being rank-deficient means the outputs do not contain contributions
from every state. That is, not all states are mapped to a linear combination in the output. Therefore,
the outputs alone are insufficient to estimate all the states.
(O)
The condition number of the observability matrix O is defined as σσmax
min (O)
where σmax is the maxi-
2
mum singular value and σmin is the minimum singular value. As this number approaches infinity,
one or more of the states becomes unobservable. This number can also be used to tell us which
2
Singular values are a generalization of eigenvalues for nonsquare matrices.
5.6 Closed-loop controller 33
sensors are better than others for the given system; a lower condition number means the outputs
produced by the sensors are better indicators of the system state.
ẋ = Ax + BK(r − x)
ẋ = Ax + BKr − BKx
ẋ = (A − BK)x + BKr (5.7)
Now for the output equation. Substitute the control law into equation (5.2).
y = Cx + D(K(r − x))
y = Cx + DKr − DKx
y = (C − DK)x + DKr (5.8)
Now, we’ll do the same for the discrete system. We’d like to know whether the system defined by
equation (5.3) operating with the control law uk = K(rk − xk ) converges to the reference rk .
Instead of commanding the system to a state using the vector u directly, we can now specify a vector of
desired states through r and the controller will choose values of u for us over time to make the system
converge to the reference. For equation (5.9) to reach steady-state, the eigenvalues of A − BK must
be in the left-half plane. For equation (5.11) to have a bounded output, the eigenvalues of A − BK
must be within the unit circle.
The eigenvalues of A − BK are the poles of the closed-loop system. Therefore, the rate of conver-
gence and stability of the closed-loop system can be changed by moving the poles via the eigenvalues
of A − BK. A and B are inherent to the system, but K can be chosen arbitrarily by the controller
designer.
where J represents a trade-off between state excursion and control effort with the weighting factors
Q and R. LQR finds a control law u that minimizes the cost function. Q and R slide the cost along
a Pareto boundary between state tracking and control effort (see figure 5.1). Pareto optimality for
this problem means that an improvement in state tracking cannot be obtained without using more
control effort to do so. Also, a reduction in control effort cannot be obtained without sacrificing state
tracking performance. Pole placement, on the other hand, will have a cost anywhere on, above, or to
the right of the Pareto boundary (no cost can be inside the boundary).
The minimum of LQR’s cost function is found by setting the derivative of the cost function to zero
and solving for the control law u. However, matrix calculus is used instead of normal calculus to take
the derivative.
5.8 Linear-quadratic regulator 35
Z∞
min xT Qx + uT Ru dt
u
0
subject to ẋ = Ax + Bu (5.13)
If the system is controllable, the optimal control policy u∗ that drives all the states to zero is −Kx.
To converge to nonzero states, a reference vector r can be added to the state x.
u = K(r − x) (5.14)
This means that optimal control can be achieved with simply a set of proportional gains on all the
states. To use the control law, we need knowledge of the full state of the system. That means we
either have to measure all our states directly or estimate those we do not measure.
See appendix E for how K is calculated in Python. If the result is finite, the controller is guaranteed
to be stable and robust with a phase margin of 60 degrees [25].
R LQR design’s Q and R matrices don’t need discretization, but the K calculated for contin-
uous time and discrete time systems will be different. The discrete time gains approach the
continuous time gains as the sample period tends to zero.
Q and R matrices are chosen based on the maximum acceptable value for each state and actuator.
The nondiagonal elements are zero. The balance between Q and R can be slid along the Pareto
boundary using a weighting factor ρ.
Z∞ " 2 2 # " 2 2 #!
x1 xm u1 un
J= ρ + ... + + + ... + dt
x1,max xm,max u1,max un,max
0
ρ 1
x21,max
0 ... 0 u21,max
0 ... 0
.. ..
0 ρ
. 0 1
.
x22,max u22,max
Q= R =
.. .. .. ..
. . 0 . . 0
ρ 1
0 ... 0 x2m,max
0 ... 0 u2n,max
Small values of ρ penalize control effort while large values of ρ penalize state excursions. Large
values would be chosen in applications like fighter jets where performance is necessary. Spacecrafts
would use small values to conserve their limited fuel supply.
Figure 5.2 shows the response using discrete poles3 placed at (0.1, 0) and (0.9, 0) and LQR with the
following cost matrices.
1
0
Q = 20 2
1 R = 1212
0 402
LQR selected poles at (0.593, 0) and (0.955, 0). Notice with pole placement that as the current pole
moves left, the control effort becomes more aggressive.
3
See section 6.2 for pole mappings of discrete systems (inside unit circle is stable). The pole mappings mentioned so
far (LHP is stable) only apply to continuous systems.
5.9 Model augmentation 37
Figure 5.2: Second-order CIM motor response with pole placement and LQR
current measurements. We just tell the plant what kind of dynamics the term has and the observer
will estimate it for us.
5.9.5 Examples
Snippet 5.1 shows how one packs together the following augmented matrix in Python using concate-
nation.
A B
C D
#!/ usr/bin/env python3
import numpy as np
Snippet 5.2 shows how one packs together the same augmented matrix in Python using array slices.
#!/ usr/bin/env python3
import numpy as np
Section 5.11 demonstrates model augmentation for different types of integral control.
5.10 Feedforward
So far, we’ve used feedback control for reference tracking (making a system’s output follow a desired
reference signal). While this is effective, it’s a reactionary measure; the system won’t start applying
control effort until the system is already behind. If we could tell the controller about the desired
movement and required input beforehand, the system could react quicker and the feedback controller
could do less work. A controller that feeds information forward into the plant like this is called a
feedforward controller.
A feedforward controller injects information about the system’s dynamics (like a model does) or the
desired movement. The feedforward handles parts of the control actions we already know must be
applied to make a system track a reference, then feedback compensates for what we do not or cannot
know about the system’s behavior at runtime.
There are two types of feedforwards: model-based feedforward and feedforward for unmodeled
40 Chapter 5. State-space controllers
dynamics. The first solves a mathematical model of the system for the inputs required to meet desired
velocities and accelerations. The second compensates for unmodeled forces or behaviors directly so
the feedback controller doesn’t have to. Both types can facilitate simpler feedback controllers; we’ll
cover examples of each.
∂xT Ax
Theorem 5.10.1 ∂x = 2Ax where A is symmetric.
Setup
Let’s start with the equation for the reference dynamics
where uk is the feedforward input. Note that this feedforward equation does not and should not take
into account any feedback terms. We want to find the optimal uk such that we minimize the tracking
error between rk+1 and rk .
rk+1 − Ark = Buk
To solve for uk , we need to take the inverse of the nonsquare matrix B. This isn’t possible, but we
can find the pseudoinverse given some constraints on the state tracking error and control effort. To
find the optimal solution for these sorts of trade-offs, one can define a cost function and attempt to
minimize it. To do this, we’ll first solve the expression for 0.
This expression will be the state tracking cost we use in the following cost function as an H2 norm.
Minimization
Given theorem 5.10.1, find the minimum of J by taking the partial derivative with respect to uk and
setting the result to 0.
∂J
= 2BT (Buk − (rk+1 − Ark ))
∂uk
5.10 Feedforward 41
Theorem 5.10.2 — Linear plant inversion. Given the discrete model xk+1 = Axk + Buk ,
the plant inversion feedforward is
where B† is the Moore-Penrose pseudoinverse of B, rk+1 is the reference at the next timestep,
and rk is the reference at the current timestep.
Discussion
Linear plant inversion in theorem 5.10.2 compensates for reference dynamics that don’t follow how
the model inherently behaves. If they do follow the model, the feedforward has nothing to do as
the model already behaves in the desired manner. When this occurs, rk+1 − Ark will return a zero
vector.
For example, a constant reference requires a feedforward that opposes system dynamics that would
change the state over time. If the system has no dynamics, then A = I and thus
uk = B† (rk+1 − Irk )
uk = B† (rk+1 − rk )
uk = B† (rk − rk )
uk = B† (0)
uk = 0
uk = Vapp (5.16)
42 Chapter 5. State-space controllers
where Vapp is a constant. Another feedforward holds a single-jointed arm steady in the presence of
gravity. It has the following form.
uk = Vapp cos θ (5.17)
where Vapp is the voltage required to keep the single-jointed arm level with the ground, and θ is
the angle of the arm relative to the ground. Therefore, the force applied is greatest when the arm is
parallel with the ground and zero when the arm is perpendicular to the ground (at that point, the joint
supports all the weight).
Note that the elevator model could be augmented easily enough to include gravity and still be linear,
but this wouldn’t work for the single-jointed arm since a trigonometric function is required to model
the gravitational force in the arm’s rotating reference frame4 .
ẋI = e = r − Cx
u = K(r − x) − KI xI
r x
u = K KI −
0 xI
ẋ = Ax + B (u + uerror )
ẋ = Ax + Bu + Buerror
ẋ = Ax + Bu + Berror uerror
44 Chapter 5. State-space controllers
where Berror is a column vector that maps uerror to changes in the rest of the state the same way B
does for u. Berror is only a column of B if uerror corresponds to an existing input within u.
Given the above equation, we’ll augment the plant as
x˙ A Berror x B
= + u
uerror 0 0 uerror 0
x
y= C 0 + Du
uerror
Notice how the state is augmented with uerror . With this model, the observer will estimate both the
state and the uerror term. The controller is augmented similarly. r is augmented with a zero for the
goal uerror term.
u = K (r − x) − kerror uerror
r x
u = K kerror −
0 uerror
where kerror is a column vector with a 1 in a given row if uerror should be applied to that input or a
0 otherwise.
This process can be repeated for an arbitrary error which can be corrected via some linear combina-
tion of the inputs.
Chaparral by Merril Apartments at UCSC
6. Digital control
The complex plane discussed so far deals with continuous systems. In decades past, plants and con-
trollers were implemented using analog electronics, which are continuous in nature. Nowadays, mi-
croprocessors can be used to achieve cheaper, less complex controller designs. Discretization con-
verts the continuous model we’ve worked with so far from a differential equation like
ẋ = x − 3 (6.1)
where xk refers to the value of x at the k th timestep. The difference equation is run with some update
period denoted by T , by ∆T , or sometimes sloppily by dt1 .
While higher order terms of a differential equation are derivatives of the state variable (e.g., ẍ in
relation to equation (6.1)), higher order terms of a difference equation are delayed copies of the state
variable (e.g., xk−1 with respect to xk in equation (6.2)).
1
The discretization of equation (6.1) to equation (6.2) uses the forward Euler discretization method.
2
See section C.1 for an explanation of phase margin.
46 Chapter 6. Digital control
at discrete time intervals. As the sample rate of the controller decreases, the phase margin decreases
according to − T2 ω where T is the sample period and ω is the frequency of the system dynamics.
Instability occurs if the phase margin of the system reaches zero. Large amounts of phase loss can
make a stable controller in the continuous domain become unstable in the discrete domain. Here are
a few ways to combat this.
• Run the controller with a high sample rate.
• Designing the controller in the analog domain with enough phase margin to compensate for
any phase loss that occurs as part of discretization.
• Convert the plant to the digital domain and design the controller completely in the digital
domain.
s-plane z-plane
(0, 0) (1, 0)
imaginary axis edge of unit circle
(−∞, 0) (0, 0)
Figure 6.2: Single poles in various locations in Figure 6.3: Complex conjugate poles in various
z-plane locations in z-plane
3 def
The aliases of a frequency f can be expressed as falias (N ) = |f − N fs |. For example, if a 200 Hz sine wave is
sampled at 150 Hz, the observer will see a 50 Hz signal instead of a 200 Hz one.
48 Chapter 6. Digital control
Figure 6.4: The samples of two sine waves can be identical when at least one of them is at a frequency
above half the sample rate. In this case, the 2 Hz sine wave is above the Nyquist frequency 1.5 Hz.
where the function f (tn , yn ) is the slope of y at n and T is the sample period for the discrete system.
Each of these methods is essentially finding the area underneath a curve. The forward and backward
Euler methods use rectangles to approximate that area while the bilinear transform uses trapezoids
(see figures 6.6 and 6.7). Since these are approximations, there is distortion between the real dis-
crete system’s poles and the approximate poles. This is in addition to the phase loss introduced by
discretizing at a given sample rate in the first place. For fast-changing systems, this distortion can
quickly lead to instability.
6.3 Discretization methods 49
well.
Forward Euler z = 1 + Ts z = 1 + Ts
1
Reverse Euler z= 1−T s z = 1 + T s + T 2 s2 + T 3 s3 + . . .
Table 6.2: Taylor series expansions of discretization methods (scalar case). The zero-order hold
discretization method is exact.
∞
X 1 k
eX = X (6.3)
k!
k=0
To understand why the matrix exponential is used in the discretization process, consider the set of
6.6 Taylor series 53
differential equations ẋ = Ax we use to describe systems (systems also have a Bu term, but we’ll
ignore it for clarity). The solution to this type of differential equation uses an exponential. Since we
are using matrices and vectors here, we use the matrix exponential.
x(t) = eAt x0
where x0 contains the initial conditions. If the initial state is the current system state, then we can
describe the system’s state over time as
xk+1 = eAT xk
R Watch the “Taylor series” video from 3Blue1Brown’s Essence of calculus series (22 minutes)
[13] for an explanation of how the Taylor series expansion works.
The definition for the matrix exponential and the approximations below all use the Taylor series
expansion. The Taylor series is a method of approximating a function like et via the summation
of weighted polynomial terms like tk . et has the following Taylor series around t = 0.
∞ n
X t
et =
n!
n=0
where a finite upper bound on the number of terms produces an approximation of et . As n increases,
the polynomial terms increase in power and the weights by which they are multiplied decrease. For
et and some other functions, the Taylor series expansion equals the original function for all values of
t as the number of terms approaches infinity4 . Figure 6.11 shows the Taylor series expansion of et
around t = 0 for a varying number of terms.
We’ll expand the first few terms of the Taylor series expansion in equation (6.3) for X = AT so we
can compare it with other methods.
3
X 1 1 1
(AT )k = I + AT + A2 T 2 + A3 T 3
k! 2 6
k=0
Table 6.3 compares the Taylor series expansions of the discretization methods for the matrix case.
These use a more complex formula which we won’t present here.
Each of them has different stability properties. The bilinear transform preserves the (in)stability of
the continuous time system.
4
Functions for which their Taylor series expansion converges to and also equals it are called analytic functions.
54 Chapter 6. Digital control
Table 6.3: Taylor series expansions of discretization methods (matrix case). The zero-order hold
discretization method is exact.
ẋ = Ac x + Bc u + w
y = Cc x + D c u + v
where w is the process noise, v is the measurement noise, and both are zero-mean white noise sources
with covariances of Qc and Rc respectively. w and v are defined as normally distributed random
variables.
w ∼ N (0, Qc )
v ∼ N (0, Rc )
xk+1 = Ad xk + Bd uk + wk
y k = Cd x k + D d u k + v k
6.7 Zero-order hold for state-space 55
with covariances
wk ∼ N (0, Qd )
vk ∼ N (0, Rd )
A d = e Ac T (6.4)
Z T
Bd = eAc τ dτ Bc = A−1
c (Ad − I)Bc (6.5)
0
Cd = Cc (6.6)
Dd = Dc (6.7)
Z T
T
Qd = eAc τ Qc eAc τ dτ (6.8)
τ =0
1
Rd = Rc (6.9)
T
where a subscript of d denotes discrete, a subscript of c denotes the continuous version of the
corresponding matrix, T is the sample period for the discrete system, and eAc T is the matrix
exponential of Ac .
Qd can be computed as
−Ac Qc
T
ATc −Ad A−1
Φ=e
0
= d Qd
0 ATd
7. Nonlinear control
While many tools exist for designing controllers for linear systems, all systems in reality are inherently
nonlinear. We’ll briefly mention some considerations for nonlinear systems.
7.1 Introduction
Recall from linear system theory that we defined systems as having the following form:
ẋ = Ax + Bu + Γw
y = Cx + Du + v
In this equation, A and B are constant matrices, which means they are both time-invariant and linear
(all transformations on the system state are linear ones, and those transformations remain the same
for all time). In nonlinear and time-variant systems, the state evolution and output are defined by
arbitrary functions of the current states and inputs.
ẋ = f (x, u, w)
y = h(x, u, v)
Nonlinear functions come up regularly when attempting to control the pose of a vehicle in the global
coordinate frame instead of the vehicle’s rotating local coordinate frame. Converting from one to
the other requires applying a rotation matrix, which consists of sine and cosine operations. These
functions are nonlinear.
7.2 Linearization
One way to control nonlinear systems is to linearize the model around a reference point. Then, all
the powerful tools that exist for linear controls can be applied. This is done by taking the Jacobians
of f and h.
∂f (x,u,w) ∂f (x,u,w) ∂h(x,u,v) ∂h(x,u,v)
A= ∂x B= ∂u C= ∂x D= ∂u
58 Chapter 7. Nonlinear control
The subscript denotes a row of f , where each row represents the dynamics of a state.
The Jacobian is the partial derivative of a vector-valued function with respect to one of the vector
arguments. The Jacobian of f has as many rows as f , and the columns are filled with partial deriv-
iatives of f ’s rows with respect to each of the argument’s elements. For example, the Jacobian of f
with respect to x is ∂f
1 ∂f1 ∂f1
∂x1 ∂x2 . . . ∂x m
∂f2 ∂f2 ∂f2
∂f (x, u) ∂x ∂x . . . ∂x
A= = .
1
.
2
. .
m
∂x .. .. .. ..
∂fm ∂fm ∂fm
∂x1 ∂x2 ... ∂xm
∂f1
∂x1 is the partial derivative of the first row of f with respect to the first state, and so on for all rows
of f and states. This has n2 permutations and thus produces a square matrix.
The Jacobian of f with respect to u is
∂f ∂f1 ∂f1
1
...
∂u
∂f2
1 ∂u2
∂f2
∂un
∂f2
∂f (x, u) ∂u ∂u2 ... ∂un
B= = . 1 ..
..
..
∂u .. . . .
∂fm ∂fm ∂fm
∂u1 ∂u2 ... ∂un
∂f1
∂u1 is the partial derivative of the first row of f with respect to the first input, and so on for all rows
of f and inputs. This has m × n permutations and can produce a nonsquare matrix if m ̸= n.
Linearization of a nonlinear equation is a Taylor series expansion to only the first-order terms (that
is, terms whose variables have exponents on the order of x1 ). This is where the small angle approxi-
mations for sin θ and cos θ (θ and 1 respectively) come from.
Higher order partial derivatives can be added to better approximate the nonlinear dynamics. We typ-
ically only linearize around equilibrium points1 because we are interested in how the system behaves
when perturbed from equilibrium. An FAQ on this goes into more detail [17]. To be clear though,
linearizing the system around the current state as the system evolves does give a closer approximation
over time.
Note that linearization with static matrices (that is, with a time-invariant linear system) only works
if the original system in question is feedback linearizable.
1
Equilibrium points are points where ẋ = 0. At these points, the system is in steady-state.
7.4 Affine systems 59
Since the state evolution in nonlinear systems is defined by a function rather than a constant matrix,
the system’s poles as determined by linearization move around. Nonlinear control uses Lyapunov
stability to determine if nonlinear systems are stable. From a linear control theory point of view,
Lyapunov stability says the system is stable if, for a given initial condition, all possible eigenvalues
of A from that point on remain in the left-half plane. However, nonlinear control uses a different
definition.
Lyapunov stability means that the system trajectory can be kept arbitrarily close to the origin by
starting sufficiently close to it. Lyapunov’s direct method uses a function consisting of the energy in
a system or derivatives of the system’s state to prove stability around an equilibrium point. This is
done by showing that the function, and thus its inputs, decay to some ground state.
More than one Lyapunov function can prove stability, and if one function doesn’t prove it, another
candidate should be tried. For this reason, we refer to these functions as Lyapunov candidate func-
tions.
∂f (x, u) ∂f (x, u)
ẋ ≈ f (x0 , u0 ) + δx + δu
∂x x0 ,u0 ∂u x0 ,u0
∂f (x, u) ∂f (x, u)
ẋ = f (x0 , u0 ) + δx + δu
∂x x0 ,u0 ∂u x0 ,u0
An affine system is a linear system with a constant offset in the dynamics. If (x0 , u0 ) is an equilibrium
point, f (x0 , u0 ) = 0, the resulting model is linear, and LQR works as usual. If (x0 , u0 ) is, say, the
current operating point rather than an equilibrium point, the easiest way to correctly apply LQR is
1. Find a control input u0 that makes (x0 , u0 ) an equilibrium point.
2. Obtain an LQR for the linearized system.
3. Add u0 to the LQR’s control input.
For a control-affine system (a nonlinear system with linear control inputs) ẋ = f (x) + Bu, u0 can
be derived via plant inversion as follows.
ẋ = f (x0 ) + Bu0
0 = f (x0 ) + Bu0
Bu0 = −f (x0 )
u0 = −B† f (x0 ) (7.1)
ṙ = f (x) + Bu
Bu = ṙ − f (x)
60 Chapter 7. Nonlinear control
R To use equation (7.2) in a discrete controller, one can approximate ṙ with rk+1 −rk
T where T is
the time period between the two references.
8. State-space applications
Up to now, we’ve just been teaching what tools are available. Now, we’ll go into specifics on how to
apply them and provide advice on certain applications.
The code shown in each example can be obtained from frccontrol’s Git repository at https://github.
com/calcmogul/frccontrol/tree/master/examples. See appendix D for setup instructions.
8.1 Elevator
8.1.1 Continuous state-space model
The position and velocity derivatives of the elevator can be written as
ẋ = v (8.1)
v̇ = a (8.2)
GKt G2 K t
v̇ = V − v
Rrm Rr2 mKv
G2 K t GKt
v̇ = − 2 v+ V (8.3)
Rr mKv Rrm
Augment the matrix equation with the position state x, which has the model equation ẋ = v. The
matrix elements corresponding to v will be 1, and the others will be 0 since they don’t appear, so
ẋ = 0x + 1v + 0V . The existing rows will have zeroes inserted where x is multiplied in.
˙ " #
x 0 1 x 0
= 2K + GKt V
v 0 − RrG2 mK t
v
v Rrm
ẋ = Ax + Bu
y = Cx + Du
x
x= y=x u=V
v
" #
0 1 0
A= G2 Kt B = GKt C= 1 0 D=0 (8.4)
0 − Rr2 mKv Rrm
This will compensate for unmodeled dynamics such as gravity. However, using a constant voltage
feedforward to counteract gravity is preferred over uerror estimation in this case because it results in
a simpler controller with similar performance.
where B is the motor acceleration term from B and uf f is the voltage feedforward.
Buf f = −(−mg)
8.2 Flywheel 63
Buf f = mg
GKt
uf f = mg
Rrm
Rrm2 g
uf f =
GKt
8.1.4 Simulation
Python Control will be used to discretize the model and simulate it. One of the frccontrol examples1
creates and tests a controller for it. Figure 8.1 shows the closed-loop system response.
8.1.5 Implementation
The script linked above also generates two files: ElevatorCoeffs.h and ElevatorCoeffs.cpp. These can
be used with the WPILib StateSpacePlant, StateSpaceController, and StateSpaceObserver classes in
C++ and Java. A C++ implementation of this elevator controller is available online2 .
8.2 Flywheel
8.2.1 Continuous state-space model
By equation (13.18)
G2 Kt GKt
ω̇ = − ω+ V
Kv RJ RJ
Factor out ω and V into column vectors.
˙ h G2 K i GK
ω = − Kv RJt ω + RJt V
1
https://github.com/calcmogul/frccontrol/blob/master/examples/elevator.py
2
https://github.com/calcmogul/allwpilib/tree/state-space/wpilibcExamples/src/main/cpp/examples/
StateSpaceElevator
64 Chapter 8. State-space applications
ẋ = Ax + Bu
y = Cx + Du
This will compensate for unmodeled dynamics such as projectiles slowing down the flywheel.
8.2.3 Simulation
Python Control will be used to discretize the model and simulate it. One of the frccontrol examples3
creates and tests a controller for it. Figure 8.2 shows the closed-loop system response.
Notice how the control effort when the reference is reached is nonzero. This is a plant inversion
feedforward compensating for the system dynamics attempting to slow the flywheel down when no
voltage is applied.
8.2.4 Implementation
The script linked above also generates two files: FlywheelCoeffs.h and FlywheelCoeffs.cpp. These
can be used with the WPILib StateSpacePlant, StateSpaceController, and StateSpaceObserver classes
in C++ and Java. A C++ implementation of this flywheel controller is available online4 .
G2 Kt GKt
ω̇ = − ω+ V
Kv RJ RJ
Therefore,
ẋ = Ax + Bu
y = Cx + Du
Notice that in this model, the output doesn’t provide any direct measurements of the state. To esti-
mate the full state (alsa known as full observability), we only need the outputs to collectively include
linear combinations of every state5 . We’ll revisit this in chapter 9 with an example that uses range
measurements to estimate an object’s orientation.
5
While the flywheel model’s outputs are a linear combination of both the states and inputs, inputs don’t provide new
information about the states. Therefore, they don’t affect whether the system is observable.
66 Chapter 8. State-space applications
The effectiveness of this model’s observer is heavily dependent on the quality of the current sensor
used. If the sensor’s noise isn’t zero-mean, the observer won’t converge to the true state.
where V is the controller’s new input voltage, Vcmd is the old input voltage, Vnominal is the rail voltage
when effects like voltage drop due to current draw are ignored, and Vrail is the real rail voltage.
To drive the model with a more accurate voltage that includes voltage drop, the reciprocal can be
used.
Vrail
V = Vcmd (8.12)
Vnominal
where V is the model’s new input voltage. Note that if both the controller compensation and model
compensation equations are applied, the original voltage is obtained. The model input only drops
from ideal if the compensated controller voltage saturates.
By equation (13.22)
G2 Kt GKt
ω̇arm = − ωarm + V
Kv RJ RJ
Augment the matrix equation with the angle state θarm , which has the model equation θ̇arm = ωarm .
The matrix elements corresponding to ωarm will be 1, and the others will be 0 since they don’t appear,
so θ̇arm = 0θarm +1ωarm +0V . The existing rows will have zeroes inserted where θarm is multiplied
in.
˙ " #
θarm 0 1 θarm 0
= G2 Kt + GKt V
ωarm 0 −K v RJ
ωarm RJ
8.3 Single-jointed arm 67
ẋ = Ax + Bu
y = Cx + Du
θarm
x= y = θarm u = V
ωarm
" #
0 1 0
A= 2
G Kt B = GKt C= 1 0 D=0 (8.15)
0 − Kv RJ RJ
This will compensate for unmodeled dynamics such as gravity or other external loading from lifted ob-
jects. However, if only gravity compensation is desired, a feedforward of the form uf f = Vgravity cos θ
is preferred where Vgravity is the voltage required to hold the arm level with the ground and θ is the
angle of the arm with the ground.
JBuf f = −(r × F)
JBuf f = −(rF cos θ)
Torque is usually written as rF sin θ where θ is the angle between the r and F vectors, but θ in this
case is being measured from the horizontal axis rather than the vertical one, so the force vector is π4
radians out of phase. Thus, an angle of 0 results in the maximum torque from gravity being applied
rather than the minimum.
68 Chapter 8. State-space applications
The force of gravity mg is applied at the center of the arm’s mass. For a uniform beam, this is
halfway down its length, or L2 where L is the length of the arm.
L
JBuf f = − (−mg) cos θ
2
L
JBuf f = mg cos θ
2
RJ , so
GKt
B=
GKt L
J uf f = mg cos θ
RJ 2
RJ L
uf f = mg cos θ
JGKt 2
L Rmg
uf f = cos θ
2 GKt
L
2 can be adjusted according to the location of the arm’s center of mass.
8.3.4 Simulation
Python Control will be used to discretize the model and simulate it. One of the frccontrol examples6
creates and tests a controller for it. Figure 8.3 shows the closed-loop system response.
8.3.5 Implementation
The script linked above also generates two files: SingleJointedArmCoeffs.h and SingleJointedArm-
Coeffs.cpp. These can be used with the WPILib StateSpacePlant, StateSpaceController, and StateS-
6
https://github.com/calcmogul/frccontrol/blob/master/examples/single_jointed_arm.py
8.4 Pendulum 69
paceObserver classes in C++ and Java. A C++ implementation of this single-jointed arm controller
is available online7 .
8.4 Pendulum
8.4.1 State-space model
Below is the model for a pendulum
g
θ̈ = − sin θ
l
where θ is the angle of the pendulum and l is the length of the pendulum.
Since state-space representation requires that only single derivatives be used, they should be broken up
as separate states. We’ll reassign θ̇ to be ω so the derivatives are easier to keep straight for state-space
representation.
g
ω̇ = − sin θ
l
Now separate the states.
θ̇ = ω
g
ω̇ = − sin θ
l
T
This makes our state vector θ ω and our nonlinear model
ω
f (x, u) =
− gl sin θ
Linearization around θ = 0
To apply our tools for linear control theory, the model must be a linear combination of the states and
inputs (addition and multiplication by constants). Since this model is nonlinear on account of the sine
function, we should linearize it.
Linearization finds a tangent line to the nonlinear dynamics at a desired point in the state-space.
The Taylor series is a way to approximate arbitrary functions with polynomials, so we can use it for
linearization.
The taylor series expansion for sin θ around θ = 0 is θ − 61 θ3 + 1 5
120 θ − . . .. We’ll take just the
first-order term θ to obtain a linear function.
θ̇ = ω
g
ω̇ = − θ
l
Now write the model in state-space representation. We’ll write out the system of equations with the
zeroed variables included to assist with this.
θ̇ = 0θ + 1ω
g
ω̇ = − θ + 0ω
l
7
https://github.com/calcmogul/allwpilib/tree/state-space/wpilibcExamples/src/main/cpp/examples/
StateSpaceSingleJointedArm
70 Chapter 8. State-space applications
If we want to linearize around an arbitrary point, we can take the Jacobian with respect to x.
∂f (x, u) 0 1
=
∂x − gl cos θ 0
For full state feedback, knowledge of all states is required. If not all states are measured directly, an
estimator can be used to supplement them.
We may only be measuring θ in the pendulum example, not θ̇, so we’ll need to estimate the latter.
The C matrix the observer would use in this case is
C= 1 0
y = Cx
θ
y= 1 0
ω
y = 1θ + 0ω
y=θ
ẋ = Ax
y = Cx
θ
x= y=θ
ω
0 1
A= g C= 1 0 (8.19)
− l cos θ 0
1 rb2 1 rb2 1 rb2 1 rb2
v̇l = + C1 vl + − C3 vr + + C 2 Vl + − C 4 Vr
m J m J m J m J
1 r2 1 r2 1 r2 1 r2
v̇r = − b C1 vl + + b C3 vr + − b C 2 Vl + + b C 4 Vr
m J m J m J m J
Factor out vl and vr into a column vector and Vl and Vr into a column vector.
˙ 1 rb2 1 rb2 1 rb2 1 rb2
vl m + J C 1 m − J C 3 vl m + J C 2 m − J C 4 Vl
= +
vr r2 r2 r2 rb2
1
− b C1
m J
1
+ b C3 vr
m J
1
− b C2 m
1
+J m J C4
Vr
ẋ = Ax + Bu
y = Cx + Du
vl v V
x= y= l u= l
vr vr Vr
1 rb2 1 r2 1 rb2 1 rb2
+ J C1 − Jb C3 + C 2 m − J C4
A= m m B= m J
1 rb2 1 r2 1 rb2 1 rb2
− m J C1 m + Jb C3 m − J C2 m + J C4 (8.20)
1 0
C= D = 02×2
0 1
G2 K 2
where C1 = − KvlRrt2 , C2 = G l Kt
Rr , C3 = − KGvrRr
Kt
2 , and C4 = Rr .
Gr Kt
Simulation
Python Control will be used to discretize the model and simulate it. One of the frccontrol examples8
creates and tests a controller for it. Figure 8.4 shows the closed-loop system response.
Given the high inertia in drivetrains, it’s better to drive the reference with a motion profile instead of
a step input for reproducibility.
8
https://github.com/calcmogul/frccontrol/blob/master/examples/differential_drive.py
72 Chapter 8. State-space applications
Next, we’ll augment the linear subspace’s state with the global pose x, y, and θ. Here’s the model as
T T
a vector function where x = x y θ vl vr and u = Vl Vr .
vr
2 cos θ + v2l cos θ
vr
2 sin θ + v2l sin θ
vr vl
f (x, u) = 2rb −
2rb (8.21)
1
rb2
1 r2
1 r2
1 rb2
m+ J C1 vl + − Jb C3 vr + m + Jb C2 Vl + m −
J C 4 Vr
m
1 rb2 1 rb2 1 r 2
1 rb2
m − J C1 vl + m + J C3 vr + m − Jb C2 Vl + m + J C 4 Vr
As mentioned in chapter 7, one can approximate a nonlinear system via linearizations around points
of interest in the state-space and design controllers for those linearized subspaces. If we sample
linearization points progressively closer together, we converge on a control policy for the original
nonlinear system. Since the linear plant being controlled varies with time, its controller is called a
linear time-varying (LTV) controller.
If we use LQRs for the linearized subspaces, the nonlinear control policy will also be locally optimal.
We’ll be taking this approach with a differential drive. To create an LQR, we need to linearize
8.5 Differential drive 73
equation (8.21).
0 0 − vl +v
2 sin θ
r 1
2cos θ 1
2 cos θ
0 vl +vr 1 1
0 2 cos θ 2sin θ 2 sin θ
∂f (x, u) 0 0 0 1
− 2rb 1
=
2
2rb
2
∂x 0 0 0 1 r
+ Jb C1 1
−
rb
C 3
m m J
1 r2 1 rb2
0 0 0 m − Jb C1 m + J C3
0 0
0 0
∂f (x, u) 0 0
=
∂u 1 r2 1 rb2
m + Jb C2 m − J C 4
1 rb2 1 rb2
m − J C 2 m + J C 4
Therefore,
ẋ = Ax + Bu
y = Cx + Du
T T T
x = x y θ v l vr y = θ vl vr u = Vl V r
1 1
0 0 −vs 2c 2c 0 0
0 1 1
0 vc 2s 2s 0 0
0 0 0 − 2r1b 1 0 0
A=
r2
2rb
rb2
B=
1 2
rb
2
rb
1 1 1
0 0 0 + Jb C1 − J C3 m+ J C2 − J C4
m m m
1 r2 1 r2 1 rb2 1 r2
0 0 0 m − Jb C1 m + Jb C3 m − J C2 m + Jb C4
0 0 1 0 0
C= 0 0 0 1 0 D = 03×2
0 0 0 0 1
(8.22)
G2l Kt 2
where v = vl +v
2 , c = cos θ, s = sin θ, C1 = −
r
, C2 = GRr l Kt
Kv Rr 2
, C3 = − KGvrRr
Kt
2 , and
The LQR gain should be recomputed around the current operating point regularly due to the high
nonlinearity of this system. A less computationally expensive controller will be developed in later
sections.
We can also use this in an extended Kalman filter as is since the measurement model (y = Cx+Du)
is linear.
With this controller, θ becomes uncontrollable at v = 0 due to the x and y dynamics being equivalent
to a unicycle; it can’t change its heading unless it’s rolling (just like a bicycle). However, a differential
drive can rotate in place. To address this controller’s limitation at v = 0, one can temporarily switch
to an LQR of just θ, vl , and vr ; linearize the controller around a slightly nonzero state; or plan a new
trajectory after the previous one completes that provides nonzero wheel velocities to rotate the robot.
74 Chapter 8. State-space applications
The linearized differential drive model doesn’t track well because the first-order linearization of A
doesn’t capture the full heading dynamics, making the model update inaccurate. This linearization
inaccuracy is evident in the Hessian matrix (second partial derivative with respect to the state vector)
being nonzero.
0 0 − vl +v
2 cos θ 0 0
r
0 0 − vl +vr sin θ 0 0
∂ 2 f (x, u) 0 0
2
2
= 0 0 0
∂x 0 0 0 0 0
0 0 0 0 0
∂f (x, u) 1 ∂ 2 f (x, u)
f (x, u0 ) ≈ f (x0 , u0 ) + (x − x0 ) + (x − x0 )2
∂x 2 ∂x2
To include higher-order dynamics in the linearized differential drive model integration, we recom-
mend using the fourth-order Runge-Kutta (RK4) integration method on equation (8.21). See snippet
8.1 for an implementation of RK4.
def runge_kutta (f, x, u, dt):
""" Fourth order Runge -Kutta integration .
Keyword arguments :
f -- vector function to integrate
x -- vector of states
u -- vector of inputs ( constant for dt)
dt -- time for which to integrate
"""
half_dt = dt * 0.5
k1 = f(x, u)
k2 = f(x + half_dt * k1 , u)
k3 = f(x + half_dt * k2 , u)
k4 = f(x + dt * k3 , u)
return x + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
8.5 Differential drive 75
Figures 8.7 and 8.8 show a simulation using RK4 instead of the first-order model.
The tracking performance of the linearized differential drive controller (figures 8.7 and 8.8) and
Ramsete (figures 8.9 and 8.10) for a given trajectory are similar, but the former’s performance-effort
trade-off can be tuned more intuitively via the Q and R gains. However, if the x and y error cost are
too high, the x and y components of the controller will fight each other, and it will take longer to
converge to the path. This can be fixed by applying a counterclockwise rotation matrix to the global
tracking error to transform it into the robot’s coordinate frame.
R G
ex cos θ sin θ 0 ex
ey = − sin θ cos θ 0 ey
eθ 0 0 1 eθ
76 Chapter 8. State-space applications
where the the superscript R represents the robot’s coordinate frame and the superscript G represents
the global coordinate frame.
With this transformation, the x and y error cost in LQR penalize the error ahead of the robot and cross-
track error respectively instead of global pose error. Since the cross-track error is always measured
from the robot’s coordinate frame, the model used to compute the LQR should be linearized around
θ = 0 at all times.
0 0 − vl +v 2 sin 0
r 1
2 cos 0
1
2 cos 0
0 0 vl +vr cos 0 1 1
2 2 sin 0 2 sin 0
0 0 0 1
− 2rb 1
A= 2r b
1 rb2
1 rb2
0 0 0 + C 1 − C 3
m J
m J
1 rb2 1 rb2
0 0 0 m − J C 1 m + J C 3
1 1
0 0 0 2 2
0 0 v +v
l
2
r
0 0
0 0 0 − 1 1
A= 2r b 2r b
1 r 2
1 r 2
0 0 0 m + J C1
b
m − J C3
b
2
2
1 r 1 r
0 0 0 m − J C1
b
m + J C3
b
This model results in figures 8.11 and 8.12, which show slightly better tracking performance than the
previous formulation.
With the exception of the x gain plot, all functions are a variation of a square root. v = 0 and v = 1
were used to scale each function to produce a close approximation. The sign function9 is used for
symmetry around the origin in the regression for y.
The following A and B matrices of a continuous system are used to compute the LQR.
1 1
0 0 0 2 2 0 0
0 0 v 0 0 0 0
0 0 0 − 2r1b 1 0 0
A=
r2
2rb
r 2
B=
1 rb2
r 2
1 1 1
0 0 0 + Jb C1 − Jb C3 m+ J C2 − Jb C4
m m m
1 r2 1 r2 1 rb2 1 r2
0 0 0 m − Jb C1 m + Jb C3 m − J C2 m + Jb C4
(8.23)
G2 K G r Kt 2
where v = vl +v r l t Gl Kt
2 , C1 = − Kv Rr 2 , C2 = Rr , C3 = − Kv Rr 2 , and C4 = Rr .
G r Kt
The constants
C1 through C4 are from the derivation in section 13.6.
The locally optimal controller for this model is u = K(v)(r − x) where
p
kx k1,2 (v) sgn(v) kθ,1 p|v| k1,4 (v) k2,4 (v)
K(v) = (8.24)
kx −k1,2 (v) sgn(v) −kθ,1 |v| k2,4 (v) k1,4 (v)
p
k1,2 (v) = ky,0 + (ky,1 − ky,0 ) |v| (8.25)
p
k1,4 (v) = kv+ ,0 + (kv+ ,1 − kv+ ,0 ) |v| (8.26)
p
k2,4 (v) = kv− ,0 − (kv+ ,1 − kv+ ,0 ) |v| (8.27)
Figures 8.18 through 8.21 show the responses of the exact and approximate solutions are the same.
9
The sign function is defined as follows:
−1
x<0
sgn(x) = 0 x=0
1 x>0
78 Chapter 8. State-space applications
T
This produces the following linear subspace over x = vl vr xl xr .
1 rb2 1 rb2 1 rb2 1 rb2
+ J C1 − J C3 0 0 + J C2 − J C4
m m m m
1 rb2 1 rb2 1 rb2 1 rb2
A= m − J C1 m + J C3 0 0 B = m − J C2 m + J C4
1 0 0 0 0 0
0 1 0 0 0 0
(8.28)
T
The measurement model for the complete nonlinear model is now y = θ xl xr instead of
T
y = θ vl vr .
U error estimation
As per subsection 5.11.2, we will now augment the model so uerror states are added to the control
inputs.
The plant and observer augmentations should be performed before the model is discretized. After the
controller gain is computed with the unaugmented discrete model, the controller may be augmented.
Therefore, the plant and observer augmentations assume a continuous model and the controller aug-
mentation assumes a discrete controller.
The three uerror states we’ll be adding are uerror,l , uerror,r , and uerror,heading for left voltage error,
right voltage error, and heading error respectively. The left and right wheel positions are filtered
encoder positions and are not adjusted for heading error. The turning angle computed from the
left and right wheel positions is adjusted by the gyroscope heading. The heading uerror state is the
heading error between what the wheel positions imply and the gyroscope measurement.
T
The full state is thus x = x y θ vl vr xl xr uerror,l uerror,r uerror,heading .
The complete nonlinear model is as follows. Let v = vl +v 2 . The three uerror states augment the
r
The left and right voltage error states should be mapped to the corresponding velocity states, so the
system matrix should be augmented with B.
The heading uerror is measuring counterclockwise encoder understeer relative to the gyroscope head-
ing, so it should add to the left position and subtract from the right position for clockwise correction
of encoder positions. That corresponds to the following input mapping vector.
0
0
Bθ = 1
−1
Now we’ll augment the linear system matrix horizontally to accomodate the uerror states.
vl
˙ vr
vl
xl
vr
= A B Bθ xr + Bu (8.30)
xl
uerror,l
xr
uerror,r
uerror,heading
R The process noise for the voltage error states should be how much the voltage can be expected
to drop. The heading error state should be the encoder model uncertainty.
80 Chapter 8. State-space applications
Figure 8.13: Linear time-varying differential Figure 8.14: Linear time-varying differential
drive controller LQR gain regression fit (x) drive controller LQR gain fit regression fit (y)
Figure 8.15: Linear time-varying differential Figure 8.16: Linear time-varying differential
drive controller LQR gain regression fit (θ) drive controller LQR gain regression fit (vl )
Figure 8.18: Linear time-varying differential Figure 8.19: Linear time-varying differential
drive controller x-y plot (approximate) drive controller x-y plot (exact)
Figure 8.20: Linear time-varying differential Figure 8.21: Linear time-varying differential
drive controller response (approximate) drive controller response (exact)
82 Chapter 8. State-space applications
where ex is the tracking error in x, ey is the tracking error in y, and eθ is the tracking error in θ. The
3 × 3 matrix is a rotation matrix that transforms the error in the pose (represented by xd − x, yd − y,
and θd − θ) from the global coordinate frame into the vehicle’s coordinate frame.
We will use the following control laws u1 and u2 for velocity and turning rate respectively.
u1 = −k1 ex
(8.33)
u2 = −k3 eθ − k2 vd sinc(eθ )ey
where k1 , k2 , and k3 are time-varying gains and sinc(eθ ) is defined as sineθeθ . This choice of control
law may seem arbitrary, and that’s because it is. The only requirement on our choice is that there
exist an associated Lyapunov candidate function that proves the control law is globally asymptotically
stable. We’ll provide a sketch of a proof in theorem 8.6.1.
Our velocity and turning rate commands for the vehicle will use the following nonlinear transforma-
tion of the inputs.
v = vd cos eθ − u1
ω = ωd − u2
Theorem 8.6.1 Assuming that vd and ωd are bounded with bounded derivatives, and that
vd (t) → 0 or ωd (t) → 0 when t → ∞, the control laws in equation (8.33) globally asymptotically
stabilize the origin e = 0.
Proof:
To prove convergence, the paper previously mentioned uses the following Lyapunov function.
k2 2 e2
V = (ex + e2y ) + θ
2 2
where k2 is a tuning constant, ex is the tracking error in x, ey is the tracking error in y, and eθ is
the tracking error in θ.
The time derivative along the solutions of the closed-loop system is nonincreasing since
Thus, ∥e(t)∥ is bounded, V̇ (t) is uniformly continuous, and V (t) tends to some limit value. Using
the Barbalat lemma, V̇ (t) tends to zero [24].
v and ω should be the references for a reference tracker for the drivetrain. This can be a linear
controller (see subsection 8.6.3). x, y, and θ are obtained via a pose estimator (see chapter 10 for
84 Chapter 8. State-space applications
how to implement one). The desired velocity, turning rate, and pose can be varied over time according
to a desired trajectory.
With this controller, θ becomes uncontrollable at v = 0. This is fine for controlling unicycles because
that’s already an inherent structural limitation; unicycles can’t change their heading unless they are
rolling (just like bicycles). However, differential drives can rotate in place. To address this controller’s
limitation at v = 0, one can temporarily switch to an LQR of just θ, vl , and vr ; or plan a new trajectory
that provides nonzero desired turning rates to rotate the robot.
G2 K 2
where C1 = − KvlRrt2 , C2 = Gl Kt
Rr , C3 = − KGvrRr
Kt
2 , and C4 = Rr .
G r Kt
To obtain a left/right velocity reference tracker, we can just remove the position states from the
model.
1 rb2 1 rb2 1 rb2 1 rb2
+ J C1 − J C3 + J C2 − J C4
A= m m B= m m
1 rb2 1 rb2 1 rb2 1 rb2
−m J C1 m + J C3 m − J C2 m + J C4 (8.39)
1 0
C= D = 02×2
0 1
G2 K 2
where C1 = − KvlRrt2 , C2 = G l Kt
Rr , C3 = − KGvrRr
Kt
2 , and C4 = Rr .
Gr Kt
See https://github.com/calcmogul/controls-engineering-in-frc/blob/master/modern-control-theory/
ss-applications/ramsete_decoupled.py for an implementation.
8.6 Ramsete unicycle controller 85
Now, we need to solve for v̇c and ω̇. First, we’ll add equation (8.40) to equation (8.41).
2 2
2v̇c = (C1 (vc − ωrb ) + C2 Vl ) + (C3 (vc + ωrb ) + C4 Vr )
m m
1 1
v̇c = (C1 (vc − ωrb ) + C2 Vl ) + (C3 (vc + ωrb ) + C4 Vr )
m m
1 1 1 1
v̇c = (C1 + C3 )vc + (−C1 + C3 )ωrb + C2 Vl + C4 Vr
m m m m
1 rb 1 1
v̇c = (C1 + C3 )vc + (−C1 + C3 )ω + C2 Vl + C4 Vr
m m m m
Next, we’ll subtract equation (8.41) from equation (8.40).
2rb2 2r2
2ω̇rb = (C1 (vc − ωrb ) + C2 Vl ) − b (C3 (vc + ωrb ) + C4 Vr )
J J
rb rb
ω̇ = (C1 (vc − ωrb ) + C2 Vl ) − (C3 (vc + ωrb ) + C4 Vr )
J J
rb rb rb rb
ω̇ = (C1 − C3 )vc + (−C1 − C3 )ωrb + C2 Vl − C4 Vr
J J J J
rb rb2 rb rb
ω̇ = (C1 − C3 )vc + (−C1 − C3 )ω + C2 Vl − C4 Vr
J J J J
Now, just convert the two equations to state-space notation.
" rb
# 1
1 1
m (C1 + C3 ) m (−C1 + C3 ) m C2 m C4
A= rb2 B= rb
rb
J (C1 − C3 ) J (−C1 − C3 ) J C2 − rJb C4
(8.42)
1 −rb
C= D = 02×2
1 rb
G2 K 2
where C1 = − KvlRrt2 , C2 = G l Kt
Rr , C3 = − KGvrRr
Kt
2 , and C4 = Rr .
Gr Kt
86 Chapter 8. State-space applications
See https://github.com/calcmogul/controls-engineering-in-frc/blob/master/modern-control-theory/
ss-applications/ramsete_coupled.py for an implementation.
Reference tracker comparison
Figures 8.22 and 8.23 shows how well each reference tracker performs.
Figure 8.22: Left/right velocity reference tracker Figure 8.23: Linear/angular velocity reference
response to a motion profile tracker response to a motion profile
Figure 8.24: Ramsete controller response with Figure 8.25: Ramsete controller’s left/right ve-
left/right velocity reference tracker (b = 2, ζ = locity reference tracker response
0.7)
Figure 8.26: Ramsete controller response with Figure 8.27: Ramsete controller’s velocity / an-
velocity / angular velocity reference tracker (b = gular velocity reference tracker response
2, ζ = 0.7)
ẋ = v cos θ
ẏ = v sin θ
θ̇ = ω
T T
Here’s the model as a vector function where x = x y θ and u = v ω .
v cos θ
f (x, u) = v sin θ (8.43)
ω
88 Chapter 8. State-space applications
Therefore,
ẋ = Ax + Bu
y = Cx + Du
x
v
x= y y= θ u=
ω
θ
0 0 −v sin θ cos θ 0
A = 0 0 v cos θ B = sin θ 0
(8.44)
0 0 0 0 1
C= 0 0 1 D = 01×2
The LQR gain should be recomputed around the current operating point regularly due to the high
nonlinearity of this system. A less computationally expensive controller will be developed next.
10
The sign function is defined as follows:
−1
x<0
sgn(x) = 0 x=0
1 x>0
8.7 Linear time-varying unicycle controller (cascaded) 89
Figure 8.28: Linear time-varying unicycle con- Figure 8.29: Linear time-varying unicycle con-
troller cascaded LQR gain regression (x) troller cascaded LQR gain regression (y)
Theorem 8.7.2 — Linear time-varying unicycle controller. The following A and B matrices
of a continuous system are used to compute the LQR.
0 0 0 1 0 x
v
A= 0 0 v
B= 0 0
x= y u= (8.45)
ω
0 0 0 0 1 θ
kx 0 0
K(v) = p (8.46)
0 ky (v) sgn(v) kθ |v|
p
ky (v) = ky,0 + (ky,1 − ky,0 ) |v| (8.47)
Stochastic control theory is a subfield of control theory that deals with the existence of uncertainty
either in observations or in noise that drives the evolution of a system. We assign probability distri-
butions to this random noise and aim to achieve a desired control task despite the presence of this
uncertainty.
Stochastic optimal control is concerned with doing this with minimum cost defined by some cost
functional, like we did with LQR earlier. First, we’ll cover the basics of probability and how we rep-
resent linear stochastic systems in state-space representation. Then, we’ll derive an optimal estimator
using this knowledge, the Kalman filter, and demonstrate creative applications of the Kalman filter
theory.
9.1 Terminology
First, we should provide definitions for terms that have specific meanings in this field.
A causal system is one that uses only past information. A noncausal system also uses information
from the future. A filter is a causal system that filters information through a probabilistic model to
produce an estimate of a desired quantity that can’t be measured directly. A smoother is a noncausal
system, so it uses information from before and after the current state to produce a better estimate.
The Luenberger observer and Kalman filter are examples of these, where the latter chooses the pole
locations optimally based on the model and measurement uncertainties.
Computer vision can also be used for localization. By extracting features from an image taken by
the agent’s camera, like a retroreflective target in FRC, and comparing them to known dimensions,
one can determine where the agent’s camera would have to be to see that image. This can be used to
correct our state estimate in the same way we do with an encoder or gyroscope.
Variables denoted with a hat are estimates of the corresponding variable. For example, x̂ is the
estimate of the true state x.
Notice that a Luenberger observer has an extra term in the state evolution equation. This term uses the
difference between the estimated outputs and measured outputs to steer the estimated state toward
the true state. Large values of L trust the measurements more while small values trust the model
more.
9.2 State observers 95
R Using an estimator forfeits the performance guarantees from earlier, but the responses are still
generally very good if the process and measurement noises are small enough. See John Doyle’s
paper Guaranteed Margins for LQG Regulators for a proof.
A Luenberger observer combines the prediction and update steps of an estimator. To run them
separately, use the equations in theorem 9.2.2 instead.
Predict step
x̂− −
k+1 = Ax̂k + Buk (9.5)
Update step
x̂+ − −1
k+1 = x̂k+1 + A L(yk − ŷk ) (9.6)
ŷk = Cx̂−
k (9.7)
For equation (9.8) to have a bounded output, the eigenvalues of A−LC must be within the unit circle.
These eigenvalues represent how fast the estimator converges to the true state of the given model. A
fast estimator converges quickly while a slow estimator avoids amplifying noise in the measurements
used to produce a state estimate.
As stated before, the controller and estimator are dual problems. Controller gains can be found assum-
ing perfect estimator (i.e., perfect knowledge of all states). Estimator gains can be found assuming
an accurate model and a controller with perfect tracking.
The effect of noise can be seen if it is modeled stochastically as
or given that the probability of a given sample is greater than or equal to zero, the sum of probabilities
for all possible input values is equal to one.
9.3 Introduction to probability 97
The mean of a random variable is denoted by an overbar (e.g., x) pronounced x-bar. The expectation
of the difference between a random variable and its mean x − x converges to zero. In other words,
the expectation of a random variable is its mean. Here’s a proof.
Z ∞
E[x − x] = (x − x) p(x) dx
Z−∞∞ Z ∞
E[x − x] = x p(x) dx − x p(x) dx
−∞ −∞
Z ∞ Z ∞
E[x − x] = x p(x) dx − x p(x) dx
−∞ −∞
E[x − x] = x − x · 1
E[x − x] = 0
9.3.3 Variance
Informally, variance is a measure of how far the outcome of a random variable deviates from its mean.
Later, we will use variance to quantify how confident we are in the estimate of a random variable;
98 Chapter 9. Stochastic control theory
we can’t know the true value of that variable without randomness, but we can give a bound on its
randomness. Z ∞
var(x) = σ 2 = E[(x − x)2 ] = (x − x)2 p(x) dx
−∞
The variance of a joint PDF measures how a variable correlates with itself (we’ll cover variances with
respect to other variables shortly).
Z ∞Z ∞
2
var(x) = Σxx = E[(x − x) ] = (x − x)2 p(x, y) dx dy
Z −∞
∞ Z ∞
−∞
9.3.5 Covariance
A covariance is a measurement of how a variable correlates with another. If they vary in the same
direction, the covariance increases. If they vary in opposite directions, the covariance decreases.
Z ∞Z ∞
cov(x, y) = Σxy = E[(x − x)(y − y)] = (x − x)(y − y) p(x, y) dx dy
−∞ −∞
9.3.6 Correlation
Two random variables are correlated if the result of one random variable affects the result of another.
Correlation is defined as
Σxy
ρ(x, y) = p , |ρ(x, y)| ≤ 1
Σxx Σyy
So two variable’s correlation is defined as their covariance over the geometric mean of their variances.
Uncorrelated sources have a covariance of zero.
9.3.7 Independence
Two random variables are independent if the following relation is true.
p(x, y) = p(x) p(y)
This means that the values of x do not correlate with the values of y. That is, the outcome of one
random variable does not affect another’s outcome. If we assume independence,
Z ∞Z ∞
E[xy] = xy p(x, y) dx dy
−∞ −∞
Z ∞Z ∞
E[xy] = xy p(x) p(y) dx dy
−∞ −∞
Z ∞ Z ∞
E[xy] = x p(x) dx y p(y) dy
−∞ −∞
E[xy] = E[x]E[y]
E[xy] = x y
Therefore, the covariance Σxy is zero, as expected. Furthermore, ρ(x, y) = 0, which means they are
uncorrelated.
100 Chapter 9. Stochastic control theory
1
The scale factor C(y ∗ ) is used to scale the area under the PDF to 1.
If x and y are independent, then p(x|y) = p(x), p(y|x) = p(y), and p(x, y) = p(x) p(y).
The elements of x are scalar variables jointly distributed with a joint density p(x1 , x2 , . . . , xn ). The
expectation is
Z ∞
E[x] = x = x p(x) dx
−∞
E[x1 ]
E[x2 ]
E[x] = .
. .
E[xn ]
Z ∞ Z ∞
E[xi ] = ... xi p(x1 , x2 , . . . , xn ) dx1 . . . dxn
−∞ −∞
Z ∞
E[f (x)] = f (x) p(x) dx
−∞
This n × n matrix is symmetric and positive semidefinite. A positive semidefinite matrix satisfies the
relation that for any v ∈ Rn for which v ̸= 0, vT Σv ≥ 0. In other words, the eigenvalues of Σ are
all greater than or equal to zero.
Σxy = cov(x, y)
Σxy = E[(x − x)(y − y)T ]
Σxy = E[xyT ] − E[xyT ] − E[xyT ] + E[xyT ]
Σxy = E[xyT ] − E[x]yT − xE[yT ] + xyT
Σxy = E[xyT ] − xyT − xyT + xyT
Σxy = E[xyT ] − xyT (9.9)
Factor out constants from the inner integral. This includes variables which are held constant for each
inner integral evaluation.
Z Z
T
E[xy ] = p(x) x dx p(y) yT dyT
X Y
Each of these integrals is just the expected value of their respective integration variable.
Σz = cov(z, z)
Σz = E[(z − z)(z − z)T ]
Σz = E[(Ax + By − Ax − By)(Ax + By − Ax − By)T ]
Σz = E[(A(x − x) + B(y − y))(A(x − x) + B(y − y))T ]
Σz = E[(A(x − x) + B(y − y))((x − x)T AT + (y − y)T BT )]
Σz = E[A(x − x)(x − x)T AT + A(x − x)(y − y)T BT +
B(y − y)(x − x)T AT + B(y − y)(y − y)T BT ]
Since x and y are independent, Σxy = 0 and the cross terms cancel out.
Σz = AΣx AT + BΣy BT
E[x] = x
var(x) = σ 2
1 (x−x)2
p(x) = √ e− 2σ 2
2πσ 2
While we could use any random variable to represent a random process, we use the Gaussian random
variable a lot in probability theory due to the central limit theorem.
9.4 Linear stochastic systems 103
Definition 9.3.1 — Central limit theorem. When independent random variables are added,
their properly normalized sum tends toward a normal distribution (a Gaussian distribution or “bell
curve”).
This is the case even if the original variables themselves are not normally distributed. The theorem is a
key concept in probability theory because it implies that probabilistic and statistical methods that work
for normal distributions can be applicable to many problems involving other types of distributions.
For example, suppose that a sample is obtained containing a large number of independent observa-
tions, and that the arithmetic mean of the observed values is computed. The central limit theorem
says that the computed values of the mean will tend toward being distributed according to a normal
distribution.
E[wk ] = 0
E[wk wkT ] = Qk
E[vk ] = 0
E[vk vkT ] = Rk
where Qk is the process noise covariance matrix and Rk is the measurement noise covariance matrix.
We assume the noise samples are independent, so E[wk wjT ] = 0 and E[vk vkT ] = 0 where k ̸= j.
Furthermore, process noise samples are independent from measurement noise samples.
We’ll compute the expectation of these equations and their covariance matrices, which we’ll use later
for deriving the Kalman filter.
Since the error and noise are independent, the cross terms are zero.
The expected value, which is also the maximum likelihood value, is the linear combination of the
prior expected (maximum likelihood) value and the measurement. The expected value is a reasonable
estimator of x.
σ02 σ2
x̂ = E[x|z1 ] = z 1 + x0 (9.13)
σ02 + σ2 2
σ0 + σ 2
x̂ = w1 z1 + w2 x0
Note that the weights w1 and w2 sum to 1. When the prior (i.e., prior knowledge of state) is uninfor-
mative (a large variance)
σ02
w1 = lim 2 =0 (9.14)
σ02 →0 σ0 + σ 2
σ2
w2 = lim 2 =1 (9.15)
σ02 →0 σ0 + σ 2
and x̂ = z1 . That is, the weight is on the observations and the estimate is equal to the measurement.
Let us assume we have a model providing an almost exact prior for x. In that case, σ02 approaches 0
and
σ02
w1 = lim 2 =1 (9.16)
σ02 →0 σ0 + σ 2
σ2
w2 = lim 2 =0 (9.17)
σ02 →0 σ0 + σ 2
The Kalman filter uses this optimal fusion as the basis for its operation.
9.6.1 Derivations
Given the a posteriori update equation x̂+ −
k+1 = x̂k+1 + Kk+1 (yk+1 − ŷk+1 ), we want to find the
value of Kk+1 that minimizes the a posteriori estimate covariance (the error covariance) because this
minimizes the estimation error.
a posteriori estimate covariance update equation
The following is the definition of the a posteriori estimate covariance matrix.
P+ +
k+1 = cov(xk+1 − x̂k+1 )
Substitute in the a posteriori update equation and expand the measurement equations.
P+ −
k+1 = cov(xk+1 − (x̂k+1 + Kk+1 (yk+1 − ŷk+1 )))
106 Chapter 9. Stochastic control theory
P+ −
k+1 = cov(xk+1 − x̂k+1 − Kk+1 (yk+1 − ŷk+1 ))
P+ −
k+1 = cov(xk+1 − x̂k+1 −
Kk+1 ((Ck+1 xk+1 + Dk+1 uk+1 + vk+1 ) − (Ck+1 x̂−
k+1 + Dk+1 uk+1 )))
P+ −
k+1 = cov(xk+1 − x̂k+1 −
Kk+1 (Ck+1 xk+1 + Dk+1 uk+1 + vk+1 − Ck+1 x̂−
k+1 − Dk+1 uk+1 ))
Reorder terms.
P+ −
k+1 = cov(xk+1 − x̂k+1 −
Kk+1 (Ck+1 xk+1 − Ck+1 x̂−
k+1 + Dk+1 uk+1 − Dk+1 uk+1 + vk+1 ))
P+ − −
k+1 = cov(xk+1 − x̂k+1 − Kk+1 (Ck+1 xk+1 − Ck+1 x̂k+1 + vk+1 ))
P+ − −
k+1 = cov(xk+1 − x̂k+1 − Kk+1 (Ck+1 xk+1 − Ck+1 x̂k+1 ) − Kk+1 vk+1 )
P+ − −
k+1 = cov(xk+1 − x̂k+1 − Kk+1 Ck+1 (xk+1 − x̂k+1 ) − Kk+1 vk+1 )
P+ −
k+1 = cov((I − Kk+1 Ck+1 )(xk+1 − x̂k+1 ) − Kk+1 vk+1 )
Covariance is a linear operator, so it can be applied to each term separately. Covariance squares
terms internally, so the negative sign on Kk+1 vk+1 is removed.
P+ −
k+1 = cov((I − Kk+1 Ck+1 )(xk+1 − x̂k+1 )) + cov(Kk+1 vk+1 )
P+ − T T
k+1 = (I − Kk+1 Ck+1 ) cov(xk+1 − x̂k+1 )(I − Kk+1 Ck+1 ) + Kk+1 cov(vk+1 )Kk+1
P+ − T T
k+1 = (I − Kk+1 Ck+1 )Pk+1 (I − Kk+1 Ck+1 ) + Kk+1 Rk+1 Kk+1
P+ − T T
k+1 = (I − Kk+1 Ck+1 )Pk+1 (I − Kk+1 Ck+1 ) + Kk+1 Rk+1 Kk+1
P+ − − T T
k+1 = (Pk+1 − Kk+1 Ck+1 Pk+1 )(I − Kk+1 Ck+1 ) + Kk+1 Rk+1 Kk+1
9.6 Kalman filter 107
P+ − − T T T
k+1 = (Pk+1 − Kk+1 Ck+1 Pk+1 )(I − Ck+1 Kk+1 ) + Kk+1 Rk+1 Kk+1
P+ − T T − T T T
k+1 = Pk+1 (I − Ck+1 Kk+1 ) − Kk+1 Ck+1 Pk+1 (I − Ck+1 Kk+1 ) + Kk+1 Rk+1 Kk+1
Expand terms.
P+ − − T T − − T T
k+1 = Pk+1 − Pk+1 Ck+1 Kk+1 − Kk+1 Ck+1 Pk+1 + Kk+1 Ck+1 Pk+1 Ck+1 Kk+1
Kk+1 Rk+1 KTk+1
P+ − − T T −
k+1 = Pk+1 − Pk+1 Ck+1 Kk+1 − Kk+1 Ck+1 Pk+1 +
Kk+1 (Ck+1 P− T T
k+1 Ck+1 + Rk+1 )Kk+1
P+ − − T T − T
k+1 = Pk+1 − Pk+1 Ck+1 Kk+1 − Kk+1 Ck+1 Pk+1 + Kk+1 Sk+1 Kk+1 (9.18)
tr(P+ − − T T − T
k+1 ) = tr(Pk+1 ) − tr(Pk+1 Ck+1 Kk+1 ) − tr(Kk+1 Ck+1 Pk+1 ) + tr(Kk+1 Sk+1 Kk+1 )
tr(P+ − −T T − T
k+1 ) = tr(Pk+1 ) − tr((Kk+1 Ck+1 Pk+1 ) ) − tr(Kk+1 Ck+1 Pk+1 ) + tr(Kk+1 Sk+1 Kk+1 )
P−
k+1 is symmetric, so we can drop the transpose.
tr(P+ − − T − T
k+1 ) = tr(Pk+1 ) − tr((Kk+1 Ck+1 Pk+1 ) ) − tr(Kk+1 Ck+1 Pk+1 ) + tr(Kk+1 Sk+1 Kk+1 )
The trace of a matrix is equal to the trace of its transpose since the elements used in the trace are on
the diagonal.
tr(P+ − − − T
k+1 ) = tr(Pk+1 ) − tr(Kk+1 Ck+1 Pk+1 ) − tr(Kk+1 Ck+1 Pk+1 ) + tr(Kk+1 Sk+1 Kk+1 )
tr(P+ − − T
k+1 ) = tr(Pk+1 ) − 2 tr(Kk+1 Ck+1 Pk+1 ) + tr(Kk+1 Sk+1 Kk+1 )
Theorem 9.6.1 ∂
∂A tr(ABAT ) = 2AB where B is symmetric.
Theorem 9.6.2 ∂
∂A tr(AC) = CT
108 Chapter 9. Stochastic control theory
Kk+1 = P− T −1
k+1 Ck+1 Sk+1
Kk+1 Sk+1 = P− T
k+1 Ck+1
Kk+1 Sk+1 KTk+1 = P− T T
k+1 Ck+1 Kk+1
P+ − − T T − T
k+1 = Pk+1 − Pk+1 Ck+1 Kk+1 − Kk+1 Ck+1 Pk+1 + Kk+1 Sk+1 Kk+1
P+ − − T T − − T T
k+1 = Pk+1 − Pk+1 Ck+1 Kk+1 − Kk+1 Ck+1 Pk+1 + (Pk+1 Ck+1 Kk+1 )
P+ − −
k+1 = Pk+1 − Kk+1 Ck+1 Pk+1
Factor out P−
k+1 to the right.
P+ −
k+1 = (I − Kk+1 Ck+1 )Pk+1
Predict step
x̂−
k+1 = Ax̂k + Buk (9.19)
P−
k+1 = AP−
kA
T
+ ΓQΓ T
(9.20)
Update step
Kk+1 = P− T − T
k+1 C (CPk+1 C + R)
−1
(9.21)
x̂+
k+1 = x̂−
k+1 + Kk+1 (yk+1 − Cx̂−
k+1 − Duk+1 ) (9.22)
P+k+1 = (I − Kk+1 C)P−
k+1 (9.23)
C, D, Q, and R from the equations derived earlier are made constants here.
R To implement a discrete time Kalman filter from a continuous model, the model and continuous
time Q and R matrices can be discretized using theorem 6.7.1.
Unknown states in a Kalman filter are generally represented by a Wiener (pronounced VEE-ner)
process2 . This process has the property that its variance increases linearly with time t.
2
Explaining why we use the Wiener process would require going much more in depth into stochastic processes and
Itô calculus, which is outside the scope of this book.
110 Chapter 9. Stochastic control theory
9.6.3 Setup
Equations to model
The following example system will be used to describe how to define and initialize the matrices for a
Kalman filter.
A robot is between two parallel walls. It starts driving from one wall to the other at a velocity of
0.8cm/s and uses ultrasonic sensors to provide noisy measurements of the distances to the walls in
front of and behind it. To estimate the distance between the walls, we will define three states: robot
position, robot velocity, and distance between the walls.
xk+1 = xk + vk ∆T (9.24)
vk+1 = vk (9.25)
xw
k+1 = xw
k (9.26)
where the Gaussian random variable wk has a mean of 0 and a variance of 1. The observation model
is
1 0 0
yk = x + θk (9.29)
−1 0 1 k
where the covariance matrix of Gaussian measurement noise θ is a 2 × 2 matrix with both diagonals
10cm2 .
The state vector is usually initialized using the first measurement or two. The covariance matrix
entries are assigned by calculating the covariance of the expressions used when assigning the state
vector. Let k = 2.
Q= 1 (9.30)
10 0
R= (9.31)
0 10
yk,1
x̂ = (yk,1 − yk−1,1 )/dt (9.32)
yk,1 + yk,2
10 10/dt 10
P = 10/dt 20/dt2 10/dt (9.33)
10 10/dt 20
Initial conditions
To fill in the P matrix, we calculate the covariance of each combination of state variables. The
resulting value is a measure of how much those variables are correlated. Due to how the covariance
calculation works out, the covariance between two variables is the sum of the variance of matching
9.6 Kalman filter 111
terms which aren’t constants multiplied by any constants the two have. If no terms match, the variables
are uncorrelated and the covariance is zero.
In P11 , the terms in x1 correlate with itself. Therefore, P11 is x1 ’s variance, or P11 = 10. For
P21 , One term correlates between x1 and x2 , so P21 = 10 dt . The constants from each are sim-
ply multiplied together. For P22 , both measurements are correlated, so the variances add together.
20
Therefore, P22 = dt 2 . It continues in this fashion until the matrix is filled up. Order doesn’t matter
for correlation, so the matrix is symmetric.
Selection of priors
Choosing good priors is important for a well performing filter, even if little information is known.
This applies to both the measurement noise and the noise model. The act of giving a state variable a
large variance means you know something about the system. Namely, you aren’t sure whether your
initial guess is close to the true state. If you make a guess and specify a small variance, you are telling
the filter that you are very confident in your guess. If that guess is incorrect, it will take the filter a
long time to move away from your guess to the true value.
Covariance selection
While one could assume no correlation between the state variables and set the covariance matrix
entries to zero, this may not reflect reality. The Kalman filter is still guarenteed to converge to the
steady-state covariance after an infinite time, but it will take longer than otherwise.
Noise model selection
We typically use a Gaussian distribution for the noise model because the sum of many independent
random variables produces a normal distribution by the central limit theorem. Kalman filters only
require that the noise is zero-mean. If the true value has an equal probability of being anywhere
within a certain range, use a uniform distribution instead. Each of these communicates information
regarding what you know about a system in addition to what you do not.
Process noise and measurement noise covariance selection
Recall that the process noise covariance is Q and the measurement noise covariance is R. To tune
the elements of these, it can be helpful to take a collection of measurements, then run the Kalman
filter on them offline to evaluate its performance.
The diagonal elements of R are the variances of each measurement, which can be easily determined
from the offline measurements. The diagonal elements of Q are the variances of each state. They
represent how much each state is expected to deviate from the model.
Selecting Q is more difficult. If the data is trusted too much over the model, one risks overfitting
the data. One should balance estimating any hidden states sufficiently with actually filtering out the
noise.
Modeling other noise colors
The Kalman filter assumes a model with zero-mean white noise. If the model is incomplete in some
way, whether it’s missing dynamics or assumes an incorrect noise model, the residual y e = y − Cx̂
over time will have probability characteristics not indicative of white noise (e.g., it isn’t zero-mean).
To handle other colors of noise in a Kalman filter, define that color of noise in terms of white noise
and augment the model with it.
9.6.4 Simulation
Figure 9.3 shows the state estimates and measurements of the Kalman filter over time. Figure 9.4
shows the position estimate and variance over time. Figure 9.5 shows the wall position estimate
112 Chapter 9. Stochastic control theory
and variance over time. Notice how the variances decrease over time as the filter gathers more
measurements. This means that the filter becomes more confident in its state estimates.
The final precisions in estimating the position of the robot and the wall are the square roots of the
corresponding elements in the covariance matrix. That is,√0.5188 m and 0.4491 m respectively. They
are smaller than the precision of the raw measurements, 10 = 3.1623 m. As expected, combining
the information from several measurements produces a better estimate than any one measurement
alone.
Figure 9.4: Robot position estimate and variance with Kalman filter
Figure 9.5: Wall position estimate and variance with Kalman filter
114 Chapter 9. Stochastic control theory
import control as ct
import numpy as np
import scipy as sp
Keyword arguments :
sys -- discrete state -space model
Q -- process noise covariance matrix
R -- measurement noise covariance matrix
Returns :
Kalman gain , error covariance matrix .
"""
m = sys.A.shape [0]
return K, P
Snippet 9.1. Steady-state Kalman gain and error covariance matrices calculation in Python
9.7.1 Derivations
9.7.2 State update equation
Let x̂(t) be the state estimate from the forward pass based on samples 0 to t and x̂b (t) be the state
estimate from the backward pass based on samples T to t.
x(t) + A′ x(t) + A′ x
x̂(t|T ) = Ax(t) + Ae eb (t)
9.7 Kalman smoother 115
where x
e(t) represents the error in the forward state estimate and x
eb (t) represents the error in the
backward state estimate.
e(t|T ) = Ae
x x(t) + (I − A)e
xb (t)
Find the minimum of the trace of P(t|T ) by taking the partial derivaitve with respect to A and
setting the result to 0 ((t) has been dropped from covariance matrices for clarity).
I − A = I − Pb (P + Pb )−1
I − A = (P + Pb )(P + Pb )−1 − Pb (P + Pb )−1
I − A = (P + Pb − Pb )(P + Pb )−1
I − A = P(P + Pb )−1 (9.38)
AT = (Pb (P + Pb )−1 )T
AT = ((P + Pb )−1 )T PTb
AT = ((P + Pb )T )−1 PTb
AT = (PT + PTb )−1 PTb
AT = (P + Pb )−1 Pb (9.39)
116 Chapter 9. Stochastic control theory
Now starting from equation (9.36), substitute in equations (9.37) through (9.40).
Apply theorem 9.7.1 to the right sides of each term to combine them.
x̂(t|T ) = (P−1 −1 −1 −1 −1
b + P ) (P x̂(t) + Pb x̂b (t))
9.7.6 Example
We will modify the robot model so that instead of a velocity of 0.8cm/s with random noise, the
velocity is modeled as a random walk from the current velocity.
xk
xk = vk
(9.53)
xw
k
1 1 0 0
xk+1 = 0 1 0 xk + 0.1 wk
(9.54)
0 0 1 0
3
https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/13-Smoothing.ipynb
9.7 Kalman smoother 119
From there, the continuous Kalman filter equations are used like normal to compute the error covari-
ance matrix P and Kalman gain matrix. The state estimate update can still use the function h(x) for
accuracy.
x̂+ − −
k+1 = x̂k+1 + Kk+1 (yk+1 − h(x̂k+1 ))
4
https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/10-Unscented-Kalman-Filter.
ipynb
9.10 Multiple model adaptive estimation 121
The original paper on the unscented Kalman filter is also an option [18]. The equations for van der
Merwe’s sigma point algorithm are in [31]. Here’s a paper on a quaternion-based Unscented Kalman
filter for orientation tracking [21].
5
MMAE section of https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/
14-Adaptive-Filtering.ipynb
This page intentionally left blank
Road next to Stevenson Academic building at UCSC
Pose is a term for the orientation of an agent (a system with a controller). The plant usually includes
the pose in its state vector. We’ll cover several methods for estimating an agent’s pose from local
measurements such as encoders and gyroscope heading.
xk+1 = xk + vk cos θk T
yk+1 = yk + vk sin θk T
θk+1 = θgyro,k+1
where T is the sample period. This odometry approach assumes that the robot follows a straight path
between samples (that is, ω = 0 at all but the sample times).
Associativity means that within an expression containing two or more occurrences in a row of the
same associative operator, the order in which the operations are performed does not matter as long
as the sequence of the operands is not changed. In other words, different groupings of the operations
produces the same result.
Identity, or an identity element, is a special type of element of a set with respect to a binary operation
on that set, which leaves any element of the set unchanged when combined with it. For example,
the additive identity of the set of integers is zero, which means that any integer summed with zero
produces that integer.
Invertibility means there is an element that can “undo” the effect of combination with another given
element. For integers and the addition operator, the inverse element would be the negation.
We use the pose exponential to take encoder measurement deltas and gyro angle deltas (which are in
the tangent space and are thus a twist) and turn them into a change in pose. This gets added to the
pose from the last update.
10.2.4 Derivation
We can obtain a more accurate approximation of the pose than Euler integration by including first
order dynamics for the heading θ.
x
x = y
θ
vx , vy , and ω are the x and y velocities of the robot within its local coordinate frame, which will be
treated as constants.
R There are two coordinate frames used here: robot and global. A superscript on the left side of a
matrix denotes the coordinate frame in which that matrix is represented. The robot’s coordinate
frame is denoted by R and the global coordinate frame is denoted by G.
To transform this into the global frame SE(2), we apply a counterclockwise rotation matrix where θ
changes over time.
G R
dx cos θ(t) − sin θ(t) 0 vx
dy = sin θ(t) cos θ(t) 0 vy dt
dθ 0 0 1 ω
G R
dx cos ωt − sin ωt 0 vx
dy = sin ωt cos ωt 0 vy dt
dθ 0 0 1 ω
Now, integrate the matrix equation (matrices are integrated element-wise). This derivation heavily
utilizes the integration method described in section 11.2.1.
G sin ωt cos ωt
R t
∆x ω ω 0 vx
∆y = − cos ωt sin ωt
0 vy
ω ω
∆θ 0 0 t ω
0
G sin ωt cos ωt−1
R
∆x ω ω 0 vx
∆y = 1−cos ωt sin ωt
0 vy
ω ω
∆θ 0 0 t ω
This equation assumes a starting orientation of θ = 0. For nonzero starting orientations, we can apply
a counterclockwise rotation by θ.
G sin ωt cos ωt−1
R
∆x cos θ − sin θ 0 ω ω 0 vx
∆y = sin θ cos θ 0 1−cos ωt
ω
sin ωt
ω 0 vy (10.1)
∆θ 0 0 1 0 0 t ω
If we factor out a t, we can use change in pose between updates instead of velocities.
G sin ωt cos ωt−1
R
∆x cos θ − sin θ 0 ω ω 0 vx
∆y = sin θ cos θ 0 1−cos ωt sin ωt
0 vy
ω ω
∆θ 0 0 1 0 0 t ω
G sin ωt cos ωt−1
R
∆x cos θ − sin θ 0 ωt ωt 0 vx
∆y = sin θ cos θ 0 1−cos ωt sin ωt
0 vy t
ωt ωt
∆θ 0 0 1 0 0 1 ω
G sin ωt cos ωt−1
R
∆x cos θ − sin θ 0 ωt ωt 0 v x t
∆y = sin θ cos θ 0 1−cos ωt sin ωt
0 vy t
ωt ωt
∆θ 0 0 1 0 0 1 ωt
G R sin ∆θ cos ∆θ−1
R
∆x cos θ − sin θ 0 ∆θ ∆θ 0 ∆x
∆y = sin θ cos θ 0 1−cos ∆θ
∆θ sin ∆θ
∆θ 0 ∆y (10.2)
∆θ 0 0 1 0 0 1 ∆θ
R
The vector ∆x ∆y ∆θ T is a twist because it’s an element of the tangent space (the robot’s
local coordinate frame).
126 Chapter 10. Pose estimation
Equation (10.1) used the current velocity and projected the model forward to the next timestep (into
the future). As such, the prediction must be done after any controller calculations are performed.
With equation (10.2), the locally measured pose deltas can only be measured using past samples, so
the model update must be performed before any controller calculations to advance the model to the
current timestep.
When the robot is traveling on a straight trajectory (ω = 0), some expressions within the equation
above are indeterminate. We can approximate these with Taylor series expansions.
sin ωt t3 ω 2
=t− + ...
ω 6
sin ωt t3 ω 2
∼t−
ω 6
cos ωt − 1 t2 ω t4 ω 3
=− + − ...
ω 2 4
cos ωt − 1 t2 ω
∼−
ω 2
1 − cos ωt t 2 ω t4 ω 3
= − + ...
ω 2 4
1 − cos ωt t2 ω
∼
ω 2
If we let ω = 0, we should get the standard kinematic equations like x = vt with a rotation applied
to them.
sin ωt t 3 · 02
∼t− =t
ω 6
cos ωt − 1 t2 · 0
∼− =0
ω 2
1 − cos ωt t2 · 0
∼ =0
ω 2
As expected, the equations simplify to the first order case with a rotation matrix applied to the veloc-
ities in the robot’s local coordinate frame.
Differential drive robots have vy = 0 since they only move in the direction of the current heading,
which is along the x-axis. Therefore,
G
∆x cos θ − sin θ 0 vx t
∆y = sin θ cos θ 0 0
∆θ 0 0 1 ωt
10.2 Pose exponential 127
G
∆x vx t cos θ
∆y = vy t sin θ
∆θ ωt
or
where G denotes global coordinate frame and R denotes robot’s coordinate frame.
For sufficiently small ω:
Figure 10.1 shows the error in the global pose coordinates over time for the simpler odometry method
compared to the method using twists (uses the Ramsete controller from subsection 8.6).
The highest error is 4cm in x over a 10m path starting with a 2m displacement. The difference
would be even more noticeable for paths with higher curvatures and longer durations. One mitigation
is using a smaller update period. However, there are bigger sources of error like turning scrub on
skid steer robots that should be dealt with before odometry numerical accuracy.
Figure 10.1: Odometry error compared to method using twists (dt = 0.05s)
11 Calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
11.1 Derivatives
11.2 Integrals
11.3 Tables
12 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 135
12.1 Linear motion
12.2 Angular motion
12.3 Vectors
12.4 Curvilinear motion
11. Calculus
This book uses derivatives and integrals occasionally to represent small changes in values over small
changes in time and the infinitesimal sum of values over time respectively. The formulas and tables
presented here are all you’ll need to carry through with the math in later chapters.
If you are interested in more information after reading this chapter, 3Blue1Brown does a fantastic
job of introducing them in his Essence of calculus video series [12]. We recommend watching videos
1 through 3 and 7 through 11 from that playlist for a solid foundation. The Taylor series (presented
in video 11) will be used in chapter 6.
11.1 Derivatives
Derivatives are expressions for the slope of a curve at arbitrary points. Common notations for this
operation on a function like f (x) include
f (x) = xn
f ′ (x) = nxn−1
h(x) = f (x)g(x)
h′ (x) = f ′ (x)g(x) + f (x)g ′ (x)
h(x) = f (g(x))
h′ (x) = f ′ (g(x)) · g ′ (x)
For example
11.2 Integrals
The integral is the inverse operation of the derivative and calculates the area under a curve. Here is
an example of one based on table 11.2.
Z
eat dt
1 at
e +C
a
The arbitrary constant C is needed because when you take a derivative, constants are discarded be-
cause vertical offsets don’t affect the slope. When performing the inverse operation, we don’t have
enough information to determine the constant.
However, we can provide bounds for the integration.
Z t
eat dt
0
t
1 at
e +C
a
0
1 at 1 a·0
e +C − e +C
a a
1 at 1
e +C − +C
a a
1 at 1
e +C − −C
a a
11.3 Tables 133
1 at 1
e −
a a
Let u = ωt.
du = ω dt
1
dt = du
ω
11.3 Tables
11.3.1 Common derivatives and integrals
R
f (x) dx f (x) f ′ (x)
ax a 0
1 2
2 ax ax a
1 a+1
a+1 x xa axa−1
1 ax
ae eax aeax
− cos(x) sin(x) cos(x)
sin(x) cos(x) − sin(x)
cos(x) − sin(x) − cos(x)
− sin(x) − cos(x) sin(x)
12. Dynamics
1
x(t) = x0 + v0 t + at2
2
where x(t) is an object’s position at time t, x0 is the initial position, v0 is the initial velocity, and a is
the acceleration.
where θ(t) is an object’s angle at time t, θ0 is the initial angle, ω0 is the initial angular velocity, and
α is the angular acceleration.
136 Chapter 12. Dynamics
12.3 Vectors
Vectors are quantities with a magnitude and a direction. Vectors in three-dimensional space have a
coordinate for each spatial direction x, y, and z. Let’s take the vector ⃗a = ⟨1, 2, 3⟩. ⃗a is a three-
dimensional vector that describes a movement 1 unit in the x direction, 2 units in the y direction, and
3 units in the z direction.
We define î, ĵ, and k̂ as vectors that represent the fundamental movements one can make the three-
dimensional space: 1 unit of movement in the x direction, 1 unit of movement y direction, and 1 unit
of movement in the z direction respectively. These three vectors form a basis of three-dimensional
space because copies of them can be added together to reach any point in three-dimensional space.
î = ⟨1, 0, 0⟩
ĵ = ⟨0, 1, 0⟩
k̂ = ⟨0, 0, 1⟩
⃗a = ⟨1, 2, 3⟩
î × ĵ = k̂
ĵ × k̂ = î
k̂ × î = ĵ
They proceed in a cyclic fashion through i, j, and k. If a vector is crossed with itself, it produces the
zero vector (a scalar zero for each coordinate). The cross product of the basis vectors in the opposite
order progress backwards and include a negative sign.
î × k̂ = −ĵ
k̂ × ĵ = −î
ĵ × î = −k̂
12.4 Curvilinear motion 137
Given vectors ⃗u = aî + bĵ + ck̂ and ⃗v = dî + eĵ + f k̂, ⃗u × ⃗v is computed using the distributive
property.
+x
+y
Figure 12.1: 2D projection of North-West-Up (NWU) axes convention. The positive z axis is pointed
out of the page toward the reader.
where ⃗vB is the velocity vector at point B, ⃗vA is the velocity vector at point A, ωA is the angular
velocity vector at point A, and ⃗rB|A is the distance vector from point A to point B (also described as
the “distance to B relative to A”).
vl vr
+x
+y
Inverse kinematics
The mapping from v and ω to the left and right wheel velocities vl and vr is derived as follows. Let
⃗vc be the velocity vector of the center of rotation, ⃗vl be the velocity vector of the left wheel, ⃗vr be the
velocity vector of the right wheel, rb is the distance from the center of rotation to each wheel, and ω
is the counterclockwise turning rate around the center of rotation.
Once we have the vector equation representing the wheel’s velocity, we’ll project it onto the wheel
direction vector using the dot product.
First, we’ll derive vl .
⃗vl = vc î + ω k̂ × rb ĵ
⃗vl = vc î − ωrb î
⃗vl = (vc − ωrb )î
Now, project this vector onto the left wheel, which is pointed in the î direction.
î
vl = (vc − ωrb )î ·
∥î∥
⃗vr = vc î + ω k̂ × rb ĵ
⃗vr = vc î + ωrb î
⃗vr = (vc + ωrb )î
Now, project this vector onto the right wheel, which is pointed in the î direction.
î
vr = (vc + ωrb )î ·
∥î∥
vr = vc + ωrb (12.2)
vl = vc − ωrb (12.3)
vr = vc + ωrb (12.4)
Forward kinematics
The mapping from the left and right wheel velocities vl and vr to v and ω is derived as follows.
vr = vc + ωrb
vc = vr − ωrb (12.5)
vl = vc − ωrb
vl = (vr − ωrb ) − ωrb
vl = vr − 2ωrb
2ωrb = vr − vl
vr − vl
ω=
2rb
vc = vr − ωrb
vr − vl
vc = vr − rb
2rb
vr − vl
vc = vr −
2
vr vl
vc = vr − +
2 2
vr + vl
vc =
2
vf l vf r
vrl vrr
+x
+y
Inverse kinematics
First, we’ll derive the front-left wheel kinematics.
î − ĵ
vf l = ((vx − ωrf ly )î + (vy + ωrf lx )ĵ) · √
2
1
vf l = ((vx − ωrf ly ) − (vy + ωrf lx )) √
2
1
vf l = (vx − ωrf ly − vy − ωrf lx ) √
2
1
vf l = (vx − vy − ωrf ly − ωrf lx ) √
2
1
vf l = (vx − vy − ω(rf lx + rf ly )) √
2
1 1 1
vf l = √ vx − √ vy − √ ω(rf lx + rf ly ) (12.8)
2 2 2
î − ĵ
vrr = ((vx − ωrrry )î + (vy + ωrrrx )ĵ) · √
2
1
vrr = ((vx − ωrrry ) − (vy + ωrrrx )) √
2
1
vrr = (vx − ωrrry − vy − ωrrrx ) √
2
1
vrr = (vx − vy − ωrrry − ωrrrx ) √
2
1
vrr = (vx − vy − ω(rrrx + rrry )) √
2
1 1 1
vrr = √ vx − √ vy − √ ω(rrrx + rrry ) (12.11)
2 2 2
142 Chapter 12. Dynamics
Forward kinematics
Let M be the 4×3 inverse kinematics matrix above including the √1 factor. The forward kinematics
2
are
vf l
vx
vy = M† vf r (12.13)
vrl
ω
vrr
v1
v3 v2
+x
v4
+y
Inverse kinematics
⃗vwheel1 = ⃗vrobot + ω
⃗ robot × ⃗rrobot2wheel1
⃗vwheel2 = ⃗vrobot + ω
⃗ robot × ⃗rrobot2wheel2
⃗vwheel3 = ⃗vrobot + ω
⃗ robot × ⃗rrobot2wheel3
⃗vwheel4 = ⃗vrobot + ω
⃗ robot × ⃗rrobot2wheel4
v1x = vx − ωr1y
v2x = vx − ωr2y
v3x = vx − ωr3y
v4x = vx − ωr4y
144 Chapter 12. Dynamics
v1y = vy + ωr1x
v2y = vy + ωr2x
v3y = vy + ωr3x
v4y = vy + ωr4x
To convert from swerve module x and y velocity components to a velocity and heading, use the
Pythagorean theorem and arctangent respectively, as shown below.
q
v = vx2 + vy2 (12.15)
vy
θ = tan−1 (12.16)
vx
Forward kinematics
Let M be the 8 × 3 inverse kinematics matrix above. The forward kinematics are
v1x
v1y
v2x
vx
vy = M† v2y (12.17)
v3x
ω
v3y
v4x
v4y
A model is a set of differential equations describing how the system behaves over time. There are
two common approaches for developing them.
1. Collecting data on the physical system’s behavior and performing system identification with it.
2. Using physics to derive the system’s model from first principles.
We’ll use the second approach in this book.
The models derived here should cover most types of motion seen on an FRC robot. Furthermore,
they can be easily tweaked to describe many types of mechanisms just by pattern-matching. There’s
only so many ways to hook up a mass to a motor in FRC. The flywheel model can be used for spinning
mechanisms, the elevator model can be used for spinning mechanisms transformed to linear motion,
and the single-jointed arm model can be used for rotating servo mechanisms (it’s just the flywheel
model augmented with a position state).
These models assume all motor controllers driving DC brushed motors are set to brake mode instead
of coast mode. Brake mode behaves the same as coast mode except where the applied voltage is zero.
In brake mode, the motor leads are shorted together to prevent movement. In coast mode, the motor
leads are an open circuit.
V
ω
Vemf = Kv
V is the voltage applied to the motor, I is the current through the motor in Amps, R is the resistance
across the motor in Ohms, ω is the angular velocity of the motor in radians per second, and Kv is the
angular velocity constant in radians per second per Volt. This circuit reflects the following relation.
ω
V = IR + (13.1)
Kv
τ = Kt I (13.2)
where τ is the torque produced by the motor in Newton-meters and Kt is the torque constant in
Newton-meters per Amp. Therefore
τ
I=
Kt
τ ω
V = R+ (13.3)
Kt Kv
τ = Kt I
τ
Kt =
I
τstall
Kt = (13.4)
Istall
where τstall is the stall torque and Istall is the stall current of the motor from its datasheet.
13.1 DC brushed motor 147
Resistance R
Recall equation (13.1).
ω
V = IR +
Kv
When the motor is stalled, ω = 0.
V = Istall R
V
R= (13.5)
Istall
where Istall is the stall current of the motor and V is the voltage applied to the motor at stall.
Angular velocity constant Kv
Recall equation (13.1).
ω
V = IR +
Kv
ω
V − IR =
Kv
ω
Kv =
V − IR
where ωf ree is the angular velocity of the motor under no load (also known as the free speed), and
V is the voltage applied to the motor when it’s spinning at ωf ree , and If ree is the current drawn by
the motor under no load.
If several identical motors are being used in one gearbox for a mechanism, multiply the stall torque
by the number of motors.
148 Chapter 13. Model examples
13.2 Elevator
13.2.1 Equations of motion
This elevator consists of a DC brushed motor attached to a pulley that drives a mass up or down.
R ωp
I ωm
Vemf
m v
V G
r
circuit mechanics
Gear ratios are written as output over input, so G is greater than one in figure 13.3.
We want to derive an equation for the carriage acceleration a (derivative of v) given an input voltage
V , which we can integrate to get carriage velocity and position.
First, we’ll find a torque to substitute into the equation for a DC brushed motor. Based on figure 13.3
τm G = τp (13.7)
13.2 Elevator 149
where G is the gear ratio between the motor and the pulley and τp is the torque produced by the
pulley.
rFm = τp (13.8)
where r is the radius of the pulley. Substitute equation (13.7) into τm in the DC brushed motor
equation (13.3).
τp
G ωm
V = R+
Kt Kv
τp ωm
V = R+
GKt Kv
where ωp is the angular velocity of the pulley. The velocity of the mass (the elevator carriage) is
v = rωp
v
ωp = (13.11)
r
Substitute equation (13.11) into equation (13.10).
v
ωm = G (13.12)
r
Solve for Fm .
RrFm G
=V − v
GKt rKv
G GKt
Fm = V − v
rKv Rr
GKt 2
G Kt
Fm = V − v (13.13)
Rr Rr2 Kv
X
F = ma (13.14)
P
where F is the sum of forces applied to the elevator carriage, m is the mass of the elevator carriage
in kilograms, and a is the acceleration of the elevator carriage.
ma = Fm
150 Chapter 13. Model examples
R Gravity is not part of the modeled dynamics because it complicates the state-space model and
the controller will behave well enough without it.
GKt G2 K t
ma = V − v
Rr Rr2 Kv
GKt G2 K t
a= V − v (13.15)
Rrm Rr2 mKv
13.3 Flywheel
13.3.1 Equations of motion
This flywheel consists of a DC brushed motor attached to a spinning mass of non-negligible moment
of inertia.
R ωf
I ωm J
Vemf
V G
circuit mechanics
Gear ratios are written as output over input, so G is greater than one in figure 13.4.
We want to derive an equation for the flywheel angular acceleration ω̇f given an input voltage V ,
which we can integrate to get flywheel angular velocity.
We will start with the equation derived earlier for a DC brushed motor, equation (13.3).
τm ωm
V = R+
Kt Kv
Solve for the angular acceleration. First, we’ll rearrange the terms because from inspection, V is the
model input, ωm is the state, and τm contains the angular acceleration.
R 1
V = τm + ωm
Kt Kv
Solve for τm .
R 1
V = τm + ωm
Kt Kv
R 1
τm =V − ωm
Kt Kv
13.3 Flywheel 151
Kt Kt
τm = V − ωm
R Kv R
τf = J ω̇f (13.17)
where J is the moment of inertia of the flywheel and ω̇f is the angular acceleration. Substitute
equation (13.17) into equation (13.16).
GKt G2 K t
(J ω̇f ) = V − ωf
R Kv R
GKt G2 K t
ω̇f = V − ωf
RJ Kv RJ
G2 K t GKt
ω̇f = − ωf + V
Kv RJ RJ
G2 Kt GKt
ω̇ = − ω+ V (13.18)
Kv RJ RJ
G2 Kt GKt
ω̇ = − ω+ V
Kv RJ RJ
G2 Kt GKt
J ω̇ = − ω+ V
Kv R R
Multiply by Kv R
G2 Kt
on both sides.
JKv R GKt Kv R
ω̇ = −ω + · 2 V
G2 K t R G Kt
JKv R Kv
ω̇ = −ω + V
G2 K t G
152 Chapter 13. Model examples
JKv R Kv
ω=− 2
ω̇ + V (13.19)
G Kt G
m
R ωarm
l
I ωm
Vemf
V G
circuit mechanics
Gear ratios are written as output over input, so G is greater than one in figure 13.5.
We want to derive an equation for the arm angular acceleration ω̇arm given an input voltage V , which
we can integrate to get arm angular velocity and angle.
We will start with the equation derived earlier for a DC brushed motor, equation (13.3).
τm ωm
V = R+
Kt Kv
Solve for the angular acceleration. First, we’ll rearrange the terms because from inspection, V is the
model input, ωm is the state, and τm contains the angular acceleration.
R 1
V = τm + ωm
Kt Kv
Solve for τm .
R 1
V = τm + ωm
Kt Kv
13.4 Single-jointed arm 153
R 1
τm = V − ωm
Kt Kv
Kt Kt
τm = V − ωm
R Kv R
where J is the moment of inertia of the arm and ω̇arm is the angular acceleration. Substitute equation
(13.21) into equation (13.20).
GKt G2 K t
(J ω̇arm ) = V − ωarm
R Kv R
G2 K t GKt
ω̇arm =− ωarm + V
Kv RJ RJ
G2 Kt GKt
ω̇ = − ω+ V (13.22)
Kv RJ RJ
where m is the mass of the arm and l is the length of the arm. Otherwise, a procedure for measuring
it experimentally is presented below.
First, rearrange equation (13.22) into the form y = mx + b such that J is in the numerator of m.
G2 Kt GKt
ω̇ = − ω+ V
Kv RJ RJ
G2 Kt GKt
J ω̇ = − ω+ V
Kv R R
154 Chapter 13. Model examples
Multiply by Kv R
G2 Kt
on both sides.
JKv R GKt Kv R
2
ω̇ = −ω + · 2 V
G Kt R G Kt
JKv R Kv
ω̇ = −ω + V
G2 K t G
JKv R Kv
ω = − 2 ω̇ + V (13.24)
G Kt G
13.5 Pendulum
Kinematics and dynamics are a rather large topics, so for now, we’ll just focus on the basics required
for working with the models in this book. We’ll derive the same model, a pendulum, using three
approaches: sum of forces, sum of torques, and conservation of energy.
θ
θ
y0
y1 θ0
mg sin θ
mg cos θ h
θ
mg
(b) Trigonometry of a pendulum
(a) Force diagram of a pendulum
F = ma
13.5 Pendulum 155
where F is the sum of forces on the object, m is mass, and a is the acceleration. Because we are only
concerned with changes in speed, and because the bob is forced to stay in a circular path, we apply
Newton’s equation to the tangential axis only. The short violet arrow represents the component of the
gravitational force in the tangential axis, and trigonometry can be used to determine its magnitude.
Therefore
−mg sin θ = ma
a = −g sin θ
where g is the acceleration due to gravity near the surface of the earth. The negative sign on the right
hand side implies that θ and a always point in opposite directions. This makes sense because when a
pendulum swings further to the left, we would expect it to accelerate back toward the right.
This linear acceleration a along the red axis can be related to the change in angle θ by the arc length
formulas; s is arc length and l is the length of the pendulum.
s = lθ (13.25)
ds dθ
v= =l
dt dt
d2 s d2 θ
a= 2 =l 2
dt dt
Therefore
d2 θ
l = −g sin θ
dt2
d2 θ g
2
= − sin θ
dt l
g
θ̈ = − sin θ
l
τ =r×F
First start by defining the torque on the pendulum bob using the force due to gravity.
τ = l × Fg
where l is the length vector of the pendulum and Fg is the force due to gravity.
For now just consider the magnitude of the torque on the pendulum.
|τ | = −mgl sin θ
where m is the mass of the pendulum, g is the acceleration due to gravity, l is the length of the
pendulum and θ is the angle between the length vector and the force due to gravity.
Next rewrite the angular momentum.
L = r × p = mr × (ω × r)
156 Chapter 13. Model examples
|L| = mr2 ω
dθ
|L| = ml2
dt
d d 2θ
|L| = ml2 2
dt dt
According to τ = dt ,
dL
we can just compare the magnitudes.
d2 θ
−mgl sin θ = ml2
dt2
g d2 θ
− sin θ = 2
l dt
g
θ̈ = − sin θ
l
Since no energy is lost, the gain in one must be equal to the loss in the other
1
mv 2 = mgh
2
dθ p
v=l = 2gh
dt
dθ 2gh
= (13.26)
dt l
where h is the vertical distance the pendulum fell. Look at figure 13.6b, which presents the trigonom-
etry of a pendulum. If the pendulum starts its swing from some initial angle θ0 , then y0 , the vertical
distance from the pivot point, is given by
y0 = l cos θ0
13.6 Differential drive 157
h = l(cos θ − cos θ0 )
This equation is known as the first integral of motion. It gives the velocity in terms of the location
and includes an integration constant related to the initial displacement (θ0 ). We can differentiate by
applying the chain rule with respect to time. Doing so gives the acceleration.
r
d dθ d 2g
= (cos θ − cos θ0 )
dt dt dt l
d2 θ 1 − 2gl sin θ dθ
2
= q
dt 2 2g (cos θ − cos θ ) dt
l 0
2g
r
d2 θ 1 − l sin θ 2g
= q (cos θ − cos θ0 )
dt2 2 2g (cos θ − cos θ ) l
l 0
d2 θ g
= − sin θ
dt2 l
g
θ̈ = − sin θ
l
ωl ωr
J
r
rb
We want to derive equations for the accelerations of the left and right sides of the robot v̇l and v̇r
given left and right input voltages Vl and Vr .
From equation (13.16) of the flywheel model derivations
GKt G2 Kt
τ= V − ω (13.27)
R Kv R
158 Chapter 13. Model examples
where τ is the torque applied by one wheel of the differential drive, G is the gear ratio of the differ-
ential drive, Kt is the torque constant of the motor, R is the resistance of the motor, and Kv is the
angular velocity constant. Since τ = rF and ω = vr where v is the velocity of a given drive side
along the ground and r is the drive wheel radius
GKt G2 Kt v
(rF ) = V −
R Kv R r
GKt G2 K t
rF = V − v
R Kv Rr
GKt G2 Kt
F = V − v
Rr Kv Rr2
G2 Kt GKt
F =− 2
v+ V
Kv Rr Rr
G2l Kt Gl K t
Fl = − vl + Vl
Kv Rr2 Rr
G2 K t Gr Kt
Fr = − r 2 v r + Vr
Kv Rr Rr
where the l and r subscripts denote the side of the robot to which each variable corresponds.
G2 K 2
Let C1 = − KvlRrt2 , C2 = Gl Kt
Rr , C3 = − KGvrRr
Kt
2 , and C4 = Rr .
G r Kt
Fl = C 1 v l + C 2 V l (13.28)
Fr = C 3 v r + C 4 V r (13.29)
Substitute equation (13.31) back into equation (13.30) to obtain an expression for v̇l .
2 1 rb2
v̇l = (Fl + Fr ) − (Fl + Fr ) + (−Fl + Fr )
m m J
r 2
1
v̇l = (Fl + Fr ) − b (−Fl + Fr )
m J
1 rb2
v̇l = (Fl + Fr ) + (Fl − Fr )
m J
1 1 r2 r2
v̇l = Fl + Fr + b Fl − b Fr
m
m J J
1 rb2 1 r2
v̇l = + Fl + − b Fr (13.33)
m J m J
τ1 = r × F
τ1 = rma
where τ1 is the torque applied by a drive motor during only linear acceleration, r is the wheel radius,
m is the robot mass, and a is the linear acceleration.
τ2 = Iα
where τ2 is the torque applied by a drive motor during only angular acceleration, I is the moment of
inertia (same as J), and α is the angular acceleration. If a constant voltage is applied during both the
linear and angular acceleration tests, τ1 = τ2 . Therefore,
rma = Iα
rmv + C1 = Iω + C2
rmv = Iω + C3
I
v= ω + C3 (13.36)
rm
where v is linear velocity and ω is angular velocity. C1 , C2 , and C3 are arbitrary constants of inte-
gration that won’t be needed. The test procedure is as follows.
1. Run the drivetrain forward at a constant voltage. Record the linear velocity over time using
encoders.
2. Rotate the drivetrain around its center by applying the same voltage as the linear acceleration
test with the motors driving in opposite directions. Record the angular velocity over time using
a gyroscope.
3. Perform a linear regression of linear velocity versus angular velocity. The slope of this line has
the form rm I
as per equation (13.36).
4. Multiply the slope by rm to obtain a least squares estimate of I.
Hills by northbound freeway between Santa Maria and Ventura
For now, this chapter just covers how to obtain state-space models for common mechanisms from a
feedforward model derived from performance data.
162 Chapter 14. System identification
0 = Av + B(Kv v) (14.1)
a = Av + B(Kv v + Ka a) (14.2)
a = B(Ka a)
1 = BKa
1
B=
Ka
ẋ = Ax + Bu
position
x= u = voltage
velocity
0 1 0
ẋ = Kv x + 1 u (14.3)
0 −K a Ka
14.2 Drivetrain velocity model 163
Since the column vectors contain the same element, the elements in the second row can be rearranged.
0 A1 A2 v B1 B2 Kv1 v
= +
0 A1 A2 v B1 B2 Kv1 v
Since the rows are linearly dependent, we can use just one of them.
0 = A1 A2 v + B1 B2 Kv1 v
A1
A2
0 = v v Kv1 v Kv1 v B1
B2
A1
A2
0 = 1 1 Kv1 Kv1 B1
B2
v a
If u = Kv1 v + Ka1 a , then x = and ẋ = .
v a
a A1 A2 v B1 B2 Kv1 v + Ka1 a
= +
a A2 A1 v B2 B1 Kv1 v + Ka1 a
Since the column vectors contain the same element, the elements in the second row can be rearranged.
a A1 A2 v B1 B2 Kv1 v + Ka1 a
= +
a A1 A2 v B1 B2 Kv1 v + Ka1 a
Since the rows are linearly dependent, we can use just one of them.
a = A1 A2 v + B1 B2 Kv1 v + Ka1 a
A1
A2
a = v v Kv1 v + Ka1 a Kv1 + Ka1 a B1
B2
−Kv2 v −v 0
If u = , then x = and ẋ = .
Kv2 v v 0
0 A1 A2 −v B1 B2 −Kv2 v
= +
0 A2 A1 v B2 B1 Kv2 v
164 Chapter 14. System identification
0 −A1 A2 v −B1 B2 Kv2 v
= +
0 −A2 A1 v −B2 B1 Kv2 v
0 −A1 A2 v −B1 B2 Kv2 v
= +
0 A1 −A2 v B1 −B2 Kv2 v
Since the column vectors contain the same element, the elements in the second row can be rearranged.
0 −A1 A2 v −B1 B2 Kv2 v
= +
0 −A1 A2 v −B1 B2 Kv2 v
Since the rows are linearly dependent, we can use just one of them.
0 = −A1 A2 v + −B1 B2 Kv2 v
0 = −vA1 + vA2 − Kv2 vB1 + Kv2 vB2
A1
A2
0 = −v v −Kv2 v Kv2 v
B1
B2
A1
A2
0 = −1 1 −Kv2 Kv2 B1
B2
−Kv2 v − Ka2 a −v −a
If u = , then x = and ẋ = .
Kv2 v + Ka2 a v a
−a A1 A2 −v B1 B2 −Kv2 v − Ka2 a
= +
a A2 A1 v B2 B1 Kv2 v + Ka2 a
−a −A1 A2 v −B1 B2 Kv2 v + Ka2 a
= +
a −A2 A1 v −B2 B1 Kv2 v + Ka2 a
−a −A1 A2 v −B1 B2 Kv2 v + Ka2 a
= +
a A1 −A2 v B1 −B2 Kv2 v + Ka2 a
Since the column vectors contain the same element, the elements in the second row can be rearranged.
−a −A1 A2 v −B1 B2 Kv2 v + Ka2 a
= +
−a −A1 A2 v −B1 B2 Kv2 v + Ka2 a
Since the rows are linearly dependent, we can use just one of them.
−a = −A1 A2 v + −B1 B2 Kv2 v + Ka2 a
−a = −vA1 + vA2 − (Kv2 v + Ka2 a)B1 + Kv2 v + Ka2 a)B2
A1
A2
−a = −v v −(Kv2 v + Ka2 a) Kv2 v + Ka2 a
B1
B2
14.2 Drivetrain velocity model 165
A1
A2
a = v −v Kv2 v + Ka2 a −(Kv2 v + Ka2 a)
B1
B2
Solve for matrix elements with Wolfram Alpha. Let b = Kv1 , c = Ka1 , d = Kv2 , and f = Ka2 .
inverse of {{1, 1, b, b}, {v, v, b v + c a, b v + c a},
{-1, 1, -d, d}, {v, -v, d v + f a, -(d v + f a)}} * {{0}, {a}, {0}, {a}}
A1 −cd − bf
A2
= 1 cd − bf
B1 2cf c + f
B2 f −c
A1 −Ka1 Kv2 − Kv1 Ka2
A2 1 Ka1 Kv2 − Kv1 Ka2
=
B1 2Ka1 Ka2 Ka1 + Ka2
B2 Ka2 − Ka1
To summarize,
ẋ = Ax + Bu
left velocity left voltage
x= u=
right velocity right voltage
A1 A2 B1 B2
ẋ = x+ u
A2 A1 B2 B1
A1 −Ka,linear Kv,angular − Kv,linear Ka,angular
A2 1 Ka,linear Kv,angular − Kv,linear Ka,angular
= (14.4)
B1 2Ka,linear Ka,angular Ka,linear + Ka,angular
B2 Ka,angular − Ka,linear
If Kv and Ka are the same for both the linear and angular cases, it devolves to the one-dimensional
case. This means the left and right sides are decoupled.
A1 −Ka Kv − Kv Ka
A2 1 Ka K v − K v Ka
=
B1 2Ka Ka Ka + Ka
B2 Ka − Ka
166 Chapter 14. System identification
A1 −2Kv Ka
A2 1 0
=
B1 2Ka Ka 2Ka
B2 0
Kv
A1 − Ka
A2 0
= 1
B1
Ka
B2 0
V
Motion planning
If smooth, predictable motion of a system over time is desired, it’s best to only change a system’s
reference as fast as the system is able to physically move. Motion profiles, also known as trajectories,
are used for this purpose. For multi-state systems, each state is given its own trajectory. Since these
states are usually position and velocity, they share different derivatives of the same profile.
These profiles are given their names based on the shape of their velocity trajectory. The trapezoidal
profile has a velocity trajectory shaped like a trapezoid and the S-curve profile has a velocity trajectory
170 Chapter 15. Motion profiles
15.1.1 Jerk
The motion characteristic that defines the change in acceleration, or transitional period, is known as
”jerk”. Jerk is defined as the rate of change of acceleration with time. In a trapezoidal profile, the
jerk (change in acceleration) is infinite at the phase transitions, while in the S-curve profile the jerk
is a constant value, spreading the change in acceleration over a period of time.
From figures 15.1 and 15.2, we can see S-curve profiles are smoother than trapezoidal profiles. Why,
however, do the S-curve profile result in less load oscillation? For a given load, the higher the jerk,
the greater the amount of unwanted vibration energy will be generated, and the broader the frequency
spectrum of the vibration’s energy will be.
This means that more rapid changes in acceleration induce more powerful vibrations, and more vibra-
tional modes will be excited. Because vibrational energy is absorbed in the system mechanics, it may
cause an increase in settling time or reduced accuracy if the vibration frequency matches resonances
in the mechanical system.
the overall transfer time. However, the reduced load oscillation at the end of the move considerably
decreases the total effective transfer time. Trial and error using a motion measurement system is gen-
erally the best way to determine the right amount of “S” because modeling high frequency dynamics
is difficult to do accurately.
Another consideration is whether that “S” segment will actually lead to smoother control of the system.
If the high frequency dynamics at play are negligible, one can use the simpler trapezoidal profile.
where x(t) is the position at time t, x0 is the initial position, v0 is the initial velocity, and a is the
acceleration at time t. The S-curve profile equations also include jerk.
1 1
x(t) = x0 + v0 t + at2 + jt3
2 6
1 2
v(t) = v0 + at + jt
2
a(t) = a0 + jt
where j is the jerk at time t, a(t) is the acceleration at time t, and a0 is the initial acceleration.
More derivations are required to determine when to start and stop the different profile phases. The
derivations for a trapezoid profile are in appendix H.4. We’re omitting an S-curve profile derivation
because the continuous acceleration over time that S-curves provide isn’t required for FRC mecha-
nisms.
and ending accelerations are very high, and there is no “S” phase where the acceleration smoothly
transitions to zero. If load oscillation is a problem, parabolic profiles may not work as well as an
S-curve despite the fact that a standard S-curve profile is not optimized for a stepper motor from the
standpoint of the torque/speed curve.
To contextualize some definitions related to trajectory optimization, let’s use the example of a driv-
etrain in FRC. A path is a set of (x, y) points for the system to follow. A trajectory is a path that
includes both the states (e.g., position and velocity) and control inputs (e.g., voltage) of the system
as functions of time. Trajectory optimization refers to a set of methods for finding the best choice of
trajectory (for some mathematical definition of “best”) by selecting states and inputs of the system
as functions of time.
Matthew Kelly has an approachable introduction to trajectory optimization on his website [20].
This page intentionally left blank
VI
Appendices
A Simplifying block diagrams . . . . . . . . . 177
A.1 Cascaded blocks
A.2 Combining blocks in parallel
A.3 Removing a block from a feedforward loop
A.4 Eliminating a feedback loop
A.5 Removing a block from a feedback loop
H Derivations . . . . . . . . . . . . . . . . . . . . . . . . 215
H.1 Transfer function in feedback
H.2 Zero-order hold for state-space
H.3 Kalman filter as Luenberger observer
H.4 Trapezoidal motion profile
Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . 220
Bibliography . . . . . . . . . . . . . . . . . . . . . . . 221
Online
Miscellaneous
Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223
This page intentionally left blank
Sunset in an airplane over New Mexico
X(s) P1 P2 Y (s)
X(s) P1 P2 Y (s)
+
X(s) P1 Y (s)
±
P2
X(s) P1 ± P2 Y (s)
+
X(s) P1 Y (s)
±
P2
P1 +
X(s) P2 P2 Y (s)
±
+
X(s) P1 Y (s)
∓
P2
P1
X(s) 1±P1 P 2 Y (s)
+
X(s) P1 Y (s)
∓
P2
1 −
X(s) P2 P1 P 2 Y (s)
∓
This appendix uses Laplace transforms and transfer functions to analyze properties of control systems
like steady-state error.
These case studies cover various aspects of PID control using the algebraic approach of transfer
functions. For this, we’ll be using equation (B.1), the transfer function for a PID controller.
Ki
K(s) = Kp + + Kd s (B.1)
s
First, we need to define what Laplace transforms and transfer functions are, which is rooted in the
concept of orthogonal projections.
B.1 Projections
Consider a two-dimensional Euclidean space R2 shown in figure B.1 (each R is a dimension whose
domain is the set of real numbers, so R2 is the standard x-y plane).
Ordinarily, we notate points in this plane by their components in the set of basis vectors {î, ĵ}, where
î (pronounced i-hat) is the unit vector in the positive x direction and ĵ is the unit vector in the positive
y direction. Figure B.2 shows an example vector v in this basis.
How do we find the coordinates of v in this basis mathematically? As long as the basis is orthogonal
(i.e., the basis vectors are at right angles to each other), we simply take the orthogonal projection of
v onto î and ĵ. Intuitively, this means finding “the amount of v that points in the direction of î or ĵ”.
Note that a set of orthogonal vectors have a dot product of zero with respect to each other.
More formally, we can calculate projections with the dot product - the projection of v onto any other
vector w is as follows.
v·w
projw v =
|w|
Since î and ĵ are unit vectors, their magnitudes are 1 so the coordinates of v are v · î and v · ĵ.
182 Appendix B. Laplace domain analysis
2
v = 2î + 1.5ĵ
1.5
ĵ
x
î
We can use this same process to find the coordinates of v in any orthogonal basis. For example,
imagine the basis {î + ĵ, î − ĵ} - the coordinates in this basis are given by v·(√î+2 ĵ) and v·(√î−2 ĵ) . Let’s
“unwrap” the formula for dot product and look a bit more closely.
n
v · (î + ĵ) 1 X
√ =√ vi (î + ĵ)i
2 2 i=0
where the subscript i denotes which component of each vector and n is the total number of compo-
nents. To change coordinates, we expanded both v and î + ĵ in a basis, multiplied their components,
and added them up. Here’s a more concrete example. Let v = 2î + 1.5ĵ from figure B.2. First, we’ll
project v onto the î + ĵ basis vector.
v · (î + ĵ) 1
√ = √ (2î · î + 1.5ĵ · ĵ)
2 2
v · (î + ĵ) 1
√ = √ (2 + 1.5)
2 2
v · (î + ĵ) 3.5
√ =√
2 2
B.1 Projections 183
√
v · (î + ĵ) 3.5 2
√ =
2 2
v · (î + ĵ) √
√ = 1.75 2
2
v · (î − ĵ) 1
√ = √ (2î · î − 1.5ĵ · ĵ)
2 2
v · (î − ĵ) 1
√ = √ (2 − 1.5)
2 2
v · (î − ĵ) 0.5
√ =√
2 2
√
v · (î − ĵ) 0.5 2
√ =
2 2
v · (î − ĵ) √
√ = 0.25 2
2
Figure B.3 shows this result geometrically with respect to the basis {î + ĵ, î − ĵ}.
√
0.25 2
v
î + ĵ √
1.75 2
x
î − ĵ
The previous example was only a change of coordinates in a finite-dimensional vector space. However,
as we will see, the core idea does not change much when we move to more complicated structures.
Observe the formula for the Fourier transform.
Z ∞
ˆ
f (ξ) = f (x)e−2πixξ dx where ξ ∈ R
−∞
This is fundamentally the same formula we had before. f (x) has taken the place of vn , e−2πixξ has
taken the place of (î + ĵ)i , and the sum over i has turned into an integral over dx, but the underlying
concept is the same. To change coordinates in a function space, we simply take the orthogonal pro-
jection onto our new basis functions. In the case of the Fourier transform, the function basis is the
184 Appendix B. Laplace domain analysis
family of functions of the form f (x) = e−2πixξ for ξ ∈ R. Since these functions are oscillatory at a
frequency determined by ξ, we can think of this as a “frequency basis”.
R Watch the “Abstract vector spaces” video from 3Blue1Brown’s Essence of linear algebra series
(17 minutes) [2] for a more geometric introduction to using functions as a basis.
Now, the Laplace transform is somewhat more complicated - as it turns out, the Fourier basis is
orthogonal, so the analogy to the simpler vector space holds almost-precisely. The Laplace transform
is not orthogonal, so we can’t interpret it strictly as a change of coordinates in the traditional sense.
However, the intuition is the same: we are taking the orthogonal projection of our original function
onto the functions of our new basis set.
Z ∞
F (s) = f (t)e−st dt, where s ∈ C
0
Here, it becomes obvious that the Laplace transform is a generalization of the Fourier transform in
that the basis family is strictly larger (we have allowed the “frequency” parameter to take complex
values, as opposed to merely real values). As a result, the Laplace basis contains functions that grow
and decay, while the Fourier basis does not.
1
The Dirac delta function is zero everywhere except at the origin. The nonzero region has an infinitesimal width and
has a height such that the area within that region is 1.
B.2 Fourier transform 185
Im(jω)
Re(σ)
To extend our analogy of each coordinate being represented by some basis, we now have the y coor-
dinate representing the oscillation frequency of the system response (the frequency domain) and also
the x coordinate representing the speed at which that oscillation decays and the system converges to
zero (i.e., a decaying exponential). Figure 3.2 shows this for various points.
If we move the component frequencies in the Fmajor4 chord example parallel to the real axis to
σ = −25, the resulting time domain response attenuates according to the decaying exponential
e−25t (see figure B.7).
Note that this explanation as a basis isn’t exact because the Laplace basis isn’t orthogonal (that is,
the x and y coordinates affect each other and have cross-talk). In the frequency domain, we had a
basis of sine waves that we represented as delta functions in the frequency domain. Each frequency
contribution was independent of the others. In the Laplace domain, this is not the case; a pure
1
exponential is s−a (a rational function where a is a real number) instead of a delta function. This
function is nonzero at points that aren’t actually frequencies present in the time domain. Figure B.8
demonstrates this, which shows the Laplace transform of the Fmajor4 chord plotted in 3D.
Notice how the values of the function around each component frequency decrease according to
√ 1 in the x and y directions (in just the x direction, it would be x1 ).
2
x +y2
B.3 Laplace transform 187
We won’t be computing any Laplace transforms by hand using this formula (everyone in the real world
looks these up in a table anyway). Common Laplace transforms (assuming zero initial conditions)
are shown in table B.1. Of particular note are the Laplace transforms for the derivative, unit step2 ,
and exponential decay. We can see that a derivative is equivalent to multiplying by s, and an integral
is equivalent to multiplying by 1s .
Table B.1: Common Laplace transforms and Laplace transform properties with zero initial conditions
K(s) = Kp
When these are in unity feedback, the transfer function from the input voltage to the error is
E(s) 1
=
V (s) 1 + K(s)G(s)
1
E(s) = V (s)
1 + K(s)G(s)
1
E(s) = V (s)
K
1 + (Kp ) (Js+b)(Ls+R)+K 2
2
The unit step u(t) is defined as 0 for t < 0 and 1 for t ≥ 0.
B.5 Case study: steady-state error 189
1
E(s) = Kp K
V (s)
1+ (Js+b)(Ls+R)+K 2
Notice that the steady-state error is nonzero. To fix this, an integrator must be included in the con-
troller.
Ki
K(s) = Kp +
s
The same steady-state calculations are performed as before with the new controller.
E(s) 1
=
V (s) 1 + K(s)G(s)
1
E(s) = V (s)
1 + K(s)G(s)
1 1
E(s) =
1 + Kp + Ksi K s
(Js+b)(Ls+R)+K 2
1 1
ess = lim s
s→0
1 + Kp + Ksi K s
(Js+b)(Ls+R)+K 2
1
ess = lim
s→0 Ki K
1 + Kp + s (Js+b)(Ls+R)+K 2
1 s
ess = lim
s→0
1 + Kp + Ki K s
s (Js+b)(Ls+R)+K 2
s
ess = lim
s→0 K
s + (Kp s + Ki ) (Js+b)(Ls+R)+K 2
0
ess =
K
0 + (Kp (0) + Ki ) (J(0)+b)(L(0)+R)+K 2
190 Appendix B. Laplace domain analysis
0
ess = K
Ki bR+K 2
For an explanation of where these equations come from, read section 13.1.
First, we’ll solve for dω
dt in terms of V .
Substitute equation (B.6) into equation (B.5).
ω
V = IR +
K
v
τ ω
V = R+
Kt Kv
Solve for dt .
dω
J dω
dt ω
V = R+
Kt Kv
ω J dω
V − = dt R
Kv Kt
dω Kt ω
= V −
dt JR Kv
dω Kt Kt
=− ω+ V
dt JRKv JR
Now take the Laplace transform. Because the Laplace transform is a linear operator, we can take the
Laplace transform of each term individually. Based on table B.1, dω
dt becomes sω and ω(t) and V (t)
become ω(s) and V (s) respectively (the parenthetical notation has been dropped for clarity).
Kt Kt
sω = − ω+ V (B.8)
JRKv JR
Kt Kt
sω = − ω+ V
JRKv JR
Kt Kt
s+ ω= V
JRKv JR
Kt
ω JR
= Kt
V s + JRK v
Kt K t Kp
sω = − ω+ (ωgoal − ω)
JRKv JR
Kt K t Kp Kt Kp
sω = − ω+ ωgoal − ω
JRKv JR JR
Kt Kt Kp Kt Kp
s+ + ω= ωgoal
JRKv JR JR
Kt Kp
ω JR
=
ωgoal s+ Kt
+
Kt K p
JRKv JR
192 Appendix B. Laplace domain analysis
Kt Kp
This has a pole at − Kt
JRKv + JR . Assuming that that quantity is negative (i.e., we are stable),
1
that pole corresponds to a time constant of Kt Kt Kp .
JRKv
+ JR
As can be seen above, a flywheel has a single pole. It therefore only needs a single pole controller to
place all of its poles anywhere.
R This analysis assumes that the motor is well coupled to the mass and that the time constant of
the inductor is small enough that it doesn’t factor into the motor equations. In Austin Schuh’s
experience with 971’s robots, these are pretty good assumptions.
Next, we’ll try a PD loop. (This will use a perfect derivative, but anyone following along closely
already knows that we can’t really take a derivative here, so the math will need to be updated at some
point. We could switch to discrete time and pick a differentiation method, or pick some other way
of modeling the derivative.)
V = Kp (ωgoal − ω) + Kd s(ωgoal − ω)
Kt Kt Kp
K +
So, we added a zero at − Kpd and moved our pole to − JRKvKt KJR
d
. This isn’t progress. We’ve added
1+ JR
more complexity to our system and, practically speaking, gotten nothing good in return. Zeroes
should be avoided if at all possible because they amplify unwanted high frequency modes of the
system and are noisier the faster the system is sampled. At least this is a stable zero, but it’s still
undesirable.
In summary, derivative doesn’t help on an ideal flywheel. Kd may compensate for unmodeled dy-
namics such as accelerating projectiles slowing the flywheel down, but that effect may also increase
B.6 Case study: flywheel PID control 193
recovery time; Kd drives the acceleration to zero in the undesired case of negative acceleration as
well as well as the actually desired case of positive acceleration.
We’ll cover a superior compensation method much later in subsection 5.11.2 that avoids zeroes in the
controller, doesn’t act against the desired control action, and facilitates better tracking.
This page intentionally left blank
Sunset in an airplane over New Mexico
C. Robust control
To install Python packages, run pip3 install --user pkg where pkg should be the name of the
package. Using --user makes installation not require root privileges. Packages can be upgraded
with pip3 install --user --upgrade pkg.
This page intentionally left blank
Sunset in an airplane over New Mexico
E. Linear-quadratic regulator
This appendix will go into more detail on the linear-quadratic regulator’s derivation and interesting
applications.
E.1 Derivation
Let a continuous time linear system be defined as
ẋ = Ax + Bu (E.1)
where J represents a trade-off between state excursion and control effort with the weighting factors
Q, R, and N. Q is the weight matrix for error, R is the weight matrix for control effort, and N is
a cross weight matrix between error and control effort. N is commonly utilized when penalizing the
output in addition to the state and input.
Z∞ T
x Qx + Nu
J= dt
u NT x + Ru
0
Z∞
Qx + Nu
J= xT uT dt
NT x + Ru
0
Z∞
J= xT (Qx + Nu) + uT NT x + Ru dt
0
200 Appendix E. Linear-quadratic regulator
Z∞
J= xT Qx + xT Nu + uT NT x + uT Ru dt
0
Z∞
T T T
T T
J= x Qx + x Nu + x Nu + u Ru dt
0
Z∞
J= xT Qx + 2xT Nu + uT Ru dt
0
Z∞
J= xT Qx + uT Ru + 2xT Nu dt
0
where K is given by
K = R−1 BT S + NT
and S is found by solving the continuous time algebraic Riccati equation defined as
AT S + SA − (SB + N) R−1 BT S + NT + Q = 0
or alternatively
AT S + SA − SBR−1 BT S + Q = 0
with
A = A − BR−1 NT
Q = Q − NR−1 NT
If there is no cross-correlation between error and control effort, N is a zero matrix and the cost
function simplifies to
Z∞
J= xT Qx + uT Ru dt
0
where K is given by
K = R−1 BT S
and S is found by solving the continuous time algebraic Riccati equation defined as
AT S + SA − SBR−1 BT S + Q = 0
The discrete time LQR controller is computed via a slightly different cost function, constraint, and
resulting algebraic Riccati equation. Snippet E.1 computes the optimal infinite horizon, discrete time
LQR controller.
E.1 Derivation 201
import control as ct
import numpy as np
import scipy as sp
.. math :: xdot = A * x + B * u
.. math :: J = \\ int_0 ^\\ infty (x^T Q x + u^T R u + 2 x^T N u) dt
Keyword arguments :
sys -- StateSpace object representing a linear system .
Q -- numpy.array( states x states ), state cost matrix .
R -- numpy.array( inputs x inputs ), control effort cost matrix .
N -- numpy.array( states x inputs ), cross weight matrix .
Returns :
K -- numpy.array( states x inputs ), controller gain matrix .
"""
sys = args [0]
Q = args [1]
R = args [2]
if len(args) == 4:
N = args [3]
else:
N = np.zeros (( sys.A.shape [0], sys.B.shape [1]))
m = sys.A.shape [0]
if sys.dt == None:
P = sp. linalg . solve_continuous_are (a=sys.A, b=sys.B, q=Q, r=R, s=N)
return np. linalg .solve(R, sys.B.T @ P + N.T)
else:
P = sp. linalg . solve_discrete_are (a=sys.A, b=sys.B, q=Q, r=R, s=N)
return np. linalg .solve(R + sys.B.T @ P @ sys.B, sys.B.T @ P @ sys.A +
N.T)
Other formulations of LQR for finite horizon and discrete time can be seen on Wikipedia [14].
MIT OpenCourseWare has a rigorous proof of the results shown above [30].
202 Appendix E. Linear-quadratic regulator
ẋ = Ax + Bu
y = Cx + Du
However, we may not know how to select costs for some of the states, or we don’t care what certain
internal states are doing. Output feedback can accomodate this. Not only can we make our out-
put contain a subset of states, but we can use any other cost metric we can think of as long as it’s
representable as a linear combination of the states and inputs1 .
For output feedback, we want to minimize the following cost functional.
Z ∞
J= yT Qy + uT Ru dt
0
T
x x
Factor out from the left side and from the right side of each term.
u u
Z ∞ T T !
x CT x x 0 0 x
J= Q C D + dt
0 u DT u u 0 R u
Z ∞ T T !
x C 0 0 x
J= T Q C D + dt
0 u D 0 R u
Multiply in Q.
Z ∞ T T !
x C Q 0 0 x
J= C D + dt
0 u DT Q 0 R u
The optimal control policy u∗ is K(r − y) where r is the desired output and y is the measured
output defined as y = Cx + Du. K can be computed via the typical LQR equations based on
the algebraic Ricatti equation.
If the output is only the whole state vector, then C = I, D = 0, and the cost functional simplifies to
that of state feedback LQR.
Z ∞ T
x Q 0 x
J= dt
0 u 0 R u
ẋ = Ax + Bu
y = Cx
ẏ = Cẋ
ẏ = C(Ax + Bu)
ẏ = CAx + CBu
ż = Aref y
ż = Aref Cx
Therefore,
T
x x
Factor out from the left side and from the right side of each term.
u u
Z ∞ T T !
x (CA − Aref C)T x x 0 0 x
J= Q CA − Aref C CB + dt
0 u (CB)T u u 0 R u
Z ∞ T !
x (CA − Aref C)T 0 0 x
J= Q CA − Aref C CB + dt
0 u (CB)T 0 R u
Multiply in Q.
Z ∞ T !
x (CA − Aref C)T Q 0 0 x
J= CA − Aref C CB + dt
0 u (CB)T Q 0 R u
subject to ẋ = Ax + Bu (E.5)
The optimal control policy u∗ is −Kx. K can be computed via the typical LQR equations based
on the algebraic Ricatti equation.
The control law uimf = −Kx makes ẋ = Ax + Buimf match the open-loop response of ż =
Aref z.
If the original and desired system have the same states, then C = I and the cost functional simplifies
to
(A − Aref )T Q(A − Aref ) (A − Aref )T QB
Z ∞ T | {z } | {z }
x Q N x
J= T T dt (E.6)
0 u B Q(A − A ref ) B QB + R u
| {z } | {z }
NT R
This page intentionally left blank
Sunset in an airplane over New Mexico
∂(Ax+b)T C(Dx+e)
Theorem F.1.1 ∂x = AT C(Dx + e) + DT CT (Ax + b)
∂(Ax+b)T C(Ax+b)
Corollary F.1.2 ∂x = 2AT C(Ax + b) where C is symmetric.
Proof:
∂(Ax + b)T C(Ax + b)
= AT C(Ax + b) + AT CT (Ax + b)
∂x
∂(Ax + b)T C(Ax + b)
= (AT C + AT CT )(Ax + b)
∂x
C is symmetric, so
F.2 Setup
Let’s start with the equation for the reference dynamics
rk+1 = Ark + Buk
where uk is the feedforward input. Note that this feedforward equation does not and should not take
into account any feedback terms. We want to find the optimal uk such that we minimize the tracking
error between rk+1 and rk .
rk+1 − Ark = Buk
To solve for uk , we need to take the inverse of the nonsquare matrix B. This isn’t possible, but we
can find the pseudoinverse given some constraints on the state tracking error and control effort. To
find the optimal solution for these sorts of trade-offs, one can define a cost function and attempt to
minimize it. To do this, we’ll first solve the expression for 0.
0 = Buk − (rk+1 − Ark )
This expression will be the state tracking cost we use in our cost function.
Our cost function will use an H2 norm with Q as the state cost matrix with dimensionality states ×
states and R as the control input cost matrix with dimensionality inputs × inputs.
J = (Buk − (rk+1 − Ark ))T Q(Buk − (rk+1 − Ark )) + uTk Ruk
F.3 Minimization
Given theorem 5.10.1 and corollary F.1.2, find the minimum of J by taking the partial derivative
with respect to uk and setting the result to 0.
∂J
= 2BT Q(Buk − (rk+1 − Ark )) + 2Ruk
∂uk
0 = 2BT Q(Buk − (rk+1 − Ark )) + 2Ruk
0 = BT Q(Buk − (rk+1 − Ark )) + Ruk
0 = BT QBuk − BT Q(rk+1 − Ark ) + Ruk
BT QBuk + Ruk = BT Q(rk+1 − Ark )
(BT QB + R)uk = BT Q(rk+1 − Ark )
uk = (BT QB + R)−1 BT Q(rk+1 − Ark )
Theorem F.3.1 — QR-weighted linear plant inversion. Given the discrete model xk+1 =
Axk + Buk , the plant inversion feedforward is
uk = Kf f (rk+1 − Ark )
where Kf f = (BT QB + R)−1 BT Q, rk+1 is the reference at the next timestep, and rk is the
reference at the current timestep.
Figure F.1 shows plant inversion applied to a second-order CIM motor model.
Plant inversion isn’t as effective with both Q and R cost because the R matrix penalized control
effort. The reference tracking with no cost matrices is much better.
F.3 Minimization 209
G. Steady-state feedforward
Steady-state feedforwards apply the control effort required to keep a system at the reference if it is no
longer moving (i.e., the system is at steady-state). The first steady-state feedforward converts desired
outputs to desired states.
xc = Nx yc
Nx converts desired outputs yc to desired states xc (also known as r). For steady-state, that is
The second steady-state feedforward converts the desired outputs y to the control input required at
steady-state.
uc = Nu yc
Nu converts the desired outputs y to the control input u required at steady-state. For steady-state,
that is
uss = Nu yss (G.2)
ẋ = Ax + Bu
y = Cx + Du
0 = Axss + Buss
yss = Cxss + Duss
0 ANx + BNu
= y
yss CNx + DNu ss
0 ANx + BNu
=
1 CNx + DNu
0 A B Nx
=
1 C D Nu
†
Nx A B 0
=
Nu C D 1
0 = (A − I)xss + Buss
yss = Cxss + Duss
0 (A − I)Nx + BNu
= yss
yss CNx + DNu
0 (A − I)Nx + BNu
=
1 CNx + DNu
0 A − I B Nx
=
1 C D Nu
†
Nx A−I B 0
=
Nu C D 1
xss = Nx yss
N†x xss = yss
uss = Nu yss
uss = Nu (N†x xss )
uss = Nu N†x xss
uf f = Nu N†x r
Continuous:
†
Nx A B 0
= (G.3)
Nu C D 1
Discrete:
†
Nx A−I B 0
= (G.4)
Nu C D 1
uf f = Nu N†x r (G.5)
In the augmented matrix, B should contain one column corresponding to an actuator and C should
contain one row whose output will be driven by that actuator. More than one actuator or output can
be included in the computation at once, but the result won’t be the same as if they were computed
independently and summed afterward.
After computing the feedforward for each actuator-output pair, the respective collections of Nx
and Nu matrices can summed to produce the combined feedforward.
If the augmented matrix in theorem G.3.1 is square (number of inputs = number of outputs), the
normal matrix inverse can be used instead.
Figure G.1: Second-order CIM motor response Figure G.2: Second-order CIM motor response
with various feedforwards with plant inversions
Plant inversion isn’t as effective in figure G.1 because the R matrix penalized control effort. If the
R cost matrix is removed from the plant inversion calculation, the reference tracking is much better
(see figure G.2).
Sunset in an airplane over New Mexico
H. Derivations
+ Z(s)
X(s) G(s) Y (s)
−
H(s)
Y (s) = Z(s)G(s)
Z(s) = X(s) − Y (s)H(s)
X(s) = Z(s) + Y (s)H(s)
X(s) = Z(s) + Z(s)G(s)H(s)
Y (s) Z(s)G(s)
=
X(s) Z(s) + Z(s)G(s)H(s)
Y (s) G(s)
= (H.1)
X(s) 1 + G(s)H(s)
where positive feedback uses the top sign and negative feedback uses the bottom sign.
216 Appendix H. Derivations
d At
e = AeAt = eAt A
dt
d −At
e x(t) = e−At Bu(t)
dt
which is an analytical solution to the continuous model. Now we want to discretize it.
def
xk = x(kT )
Z kT
xk = eAkT x(0) + eA(kT −τ ) Bu(τ ) dτ
0
Z (k+1)T
xk+1 = e A(k+1)T
x(0) + eA((k+1)T −τ ) Bu(τ ) dτ
0
Z kT Z (k+1)T
A((k+1)T −τ )
xk+1 = e A(k+1)T
x(0) + e Bu(τ ) dτ + eA((k+1)T −τ ) Bu(τ ) dτ
0 kT
Z kT Z (k+1)T
xk+1 = e A(k+1)T
x(0) + eA((k+1)T −τ ) Bu(τ ) dτ + eA(kT +T −τ ) Bu(τ ) dτ
0 kT
Z kT Z (k+1)T
A(kT −τ )
xk+1 = eAT e AkT
x(0) + e Bu(τ ) dτ + eA(kT +T −τ ) Bu(τ ) dτ
0 kT
| {z }
xk
We assume that u is constant during each timestep, so it can be pulled out of the integral.
Z !
(k+1)T
xk+1 = eAT xk + eA(kT +T −τ ) dτ Buk
kT
H.3 Kalman filter as Luenberger observer 217
The second term can be simplified by substituting it with the function v(τ ) = kT + T − τ . Note
that dτ = −dv.
Z v((k+1)T ) !
xk+1 = eAT xk − eAv dv Buk
v(kT )
Z 0
AT Av
xk+1 = e xk − e dv Buk
T
Z T
xk+1 = eAT xk + Av
e dv Buk
0
AT −1 Av T
xk+1 = e xk + A e |0 Buk
AT −1 AT A0
xk+1 = e xk + A (e −e )Buk
xk+1 = eAT xk + A−1 (eAT − I)Buk
x̂+ −
k+1 = Ax̂k + Buk + L(yk − ŷk ) (H.3)
ŷk = Cx̂−
k (H.4)
where a superscript of minus denotes a priori and plus denotes a posteriori estimate. Combining
equation (H.3) and equation (H.4) gives
x̂+ − −
k+1 = Ax̂k + Buk + L(yk − Cx̂k ) (H.5)
The following is a Kalman filter that considers the current update step and the next predict step
together rather than the current predict step and current update step.
Update step
K k = P− T − T
k C (CPk C + R)
−1
(H.6)
x̂+
k = x̂−
k + Kk (yk − Cx̂−
k) (H.7)
P+k = (I − Kk C)P−
k (H.8)
Predict step
x̂+ +
k+1 = Ax̂k + Buk (H.9)
P−
k+1 = AP+
kA
T
+ ΓQΓ T
(H.10)
x̂+ − −
k+1 = A(x̂k + Kk (yk − Cx̂k )) + Buk
x̂+ − −
k+1 = Axk + AKk (yk − Cx̂k ) + Buk
x̂+ − −
k+1 = Ax̂k + Buk + AKk (yk − Cx̂k )
Let L = AKk .
x̂+ − −
k+1 = Ax̂k + Buk + L(yk − Cx̂k ) (H.11)
218 Appendix H. Derivations
which matches equation (H.5). Therefore, the eigenvalues of the Kalman filter observer can be ob-
tained by
eig(A − LC)
eig(A − (AKk )(C))
eig(A(I − Kk C)) (H.12)
x̂+ − −
k+1 = x̂k+1 + K(yk+1 − Cx̂k+1 )
x̂+ − −1 −
k+1 = x̂k+1 + A L(yk+1 − Cx̂k+1 )
x̂+ − −1
k+1 = x̂k+1 + A L(yk+1 − ŷk+1 )
The predict step is the same as the Kalman filter’s. Therefore, a Luenberger observer run with pre-
diction and update steps is written as follows.
Predict step
x̂− −
k+1 = Ax̂k + Buk (H.13)
Update step
x̂+ − −1
k+1 = x̂k+1 + A L(yk+1 − ŷk+1 ) (H.14)
ŷk+1 = Cx̂−
k+1 (H.15)
Glossary
agent An independent actor being controlled through autonomy or human-in-the-loop (e.g., a robot,
aircraft, etc.).
control effort A term describing how much force, pressure, etc. an actuator is exerting.
control input The input of a plant used for the purpose of controlling it.
control law Also known as control policy, is a mathematical formula used by the controller to deter-
mine the input u that is sent to the plant. This control law is designed to drive the system
from its current state to some other desired state.
control system Monitors and controls the behavior of a system.
controller Applies an input to a plant to bring about a desired system state by driving the difference
between a reference signal and the output to zero.
discretization The process by which a continuous (e.g., analog) system or controller design is con-
verted to discrete (e.g., digital).
disturbance An external force acting on a system that isn’t included in the system’s model.
disturbance rejection The quality of a feedback control system to compensate for external forces
to reach a desired reference.
error Reference minus an output or state.
feedback controller Used in positive or negative feedback with a plant to bring about a desired
system state by driving the difference between a reference signal and the output to zero.
feedback gain The gain from the output to an earlier point in a control system diagram.
feedforward controller A controller that injects information about the system’s dynamics (like a
model does) or the desired movement. The feedforward handles parts of the control actions
we already know must be applied to make a system track a reference, then the feedback
controller compensates for what we do not or cannot know about the system’s behavior at
runtime.
gain A proportional value that shows the relationship between the magnitude of an input signal to
the magnitude of an output signal at steady-state.
220 Glossary
Bibliography
Online
[14] Wikipedia Commons. Linear-quadratic regulator. URL: https : / / en . wikipedia . org / wiki /
Linear%E2%80%93quadratic_regulator (visited on 03/24/2018) (cited on page 201).
[16] Declan Freeman-Gleason. A Dive into WPILib Trajectories. 2020. URL: https://pietroglyph.
github.io/trajectory-presentation/ (cited on page 172).
[17] Sean Humbert. Why do we have to linearize around an equilibrium point? URL: https://www.
cds.caltech.edu/%7Emurray/courses/cds101/fa02/faq/02-10-09_linearization.html (visited
on 07/12/2018) (cited on page 58).
[18] Simon J. Julier and Jeffrey K. Uhlmann. A New Extension of the Kalman Filter to Nonlinear
Systems. URL: https://www.cs.unc.edu/~welch/kalman/media/pdf/Julier1997_SPIE_KF.pdf
(visited on 09/29/2019) (cited on page 121).
[19] Kapitanyuk, Proskurnikov, and Cao. A guiding vector field algorithm for path following control
of nonholonomic mobile robots. URL: https : / / arxiv . org / pdf / 1610 . 04391 . pdf (visited on
08/09/2018) (cited on page 60).
[20] Matthew Kelly. An Introduction to Trajectory Optimization: How to Do Your Own Direct Collo-
cation. URL: https://www.matthewpeterkelly.com/research/MatthewKelly_IntroTrajectoryOptimization_
SIAM_Review_2017.pdf (visited on 08/27/2019) (cited on page 173).
[21] Edgar Kraft. A Quaternion-based Unscented Kalman Filter for Orientation Tracking. URL: https:
//kodlab.seas.upenn.edu/uploads/Arun/UKFpaper.pdf (visited on 07/11/2018) (cited on
page 121).
[22] Charles F. Van Loan. Computing Integrals Involving the Matrix Exponential. URL: https://www.
cs.cornell.edu/cv/ResearchPDF/computing.integrals.involving.Matrix.Exp.pdf (visited on
06/21/2018) (cited on page 55).
[23] Andrew W. Long et al. The Banana Distribution is Gaussian: A Localization Study with Expo-
nential Coordinates. 2008. URL: https://rpk.lcsr.jhu.edu/wp-content/uploads/2014/08/p34_
Long12_The-Banana-Distribution.pdf (cited on page 93).
222 Bibliography
[24] Luca, Oriolo, and Vendittelli. Control of Wheeled Mobile Robots: An Experimental Overview.
URL: https : / / www . dis . uniroma1 . it / ~labrob / pub / papers / Ramsete01 . pdf (visited on
08/09/2018) (cited on pages 82, 83).
[25] MIT OpenCourseWare. Linear Quadratic Regulator. URL: https : / / ocw . mit . edu / courses /
mechanical - engineering / 2 - 154 - maneuvering - and - control - of - surface - and - underwater -
vehicles-13-49-fall-2004/lecture-notes/lec19.pdf (visited on 06/26/2018) (cited on page 35).
[26] Christoph Sprunk. Planning Motion Trajectories for Mobile Robots Using Splines. 2008. URL:
http://www2.informatik.uni-freiburg.de/~lau/students/Sprunk2008.pdf (cited on page 172).
[29] Russ Tedrake. Underactuated Robotics. 2019. URL: http://underactuated.mit.edu/underactuated.
html (cited on page 60).
[30] Russ Tedrake. Chapter 9. Linear Quadratic Regulators. URL: http://underactuated.mit.edu/lqr.
html (visited on 03/22/2020) (cited on page 201).
[31] Eric A. Wan and Rudolph van der Merwe. The Unscented Kalman Filter for Nonlinear Esti-
mation. URL: https://www.seas.harvard.edu/courses/cs281/papers/unscented.pdf (visited on
06/26/2018) (cited on page 121).
Misc
[1] FRC team 254. Motion Planning and Control in FRC. 2015. URL: https://www.youtube.com/
watch?v=8319J1BEHwM (cited on page 172).
[2] 3Blue1Brown. Abstract vector spaces. 2016. URL: https : / / www . youtube . com / watch ? v =
TgKwz5Ikpc8 (cited on page 184).
[3] 3Blue1Brown. Eigenvectors and eigenvalues. 2016. URL: https://www.youtube.com/watch?v=
PFDu9oVAE-g (cited on page 28).
[4] 3Blue1Brown. Essence of linear algebra preview. 2016. URL: https : / / www . youtube . com /
watch?v=kjBOesZCoqc (cited on page 27).
[5] 3Blue1Brown. Inverse matrices, column space, and null space. 2016. URL: https : / / www .
youtube.com/watch?v=uQhTuRlWMxw (cited on page 28).
[6] 3Blue1Brown. Linear combinations, span, and basis vectors. 2016. URL: https://www.youtube.
com/watch?v=k7RM-ot2NWY (cited on page 27).
[7] 3Blue1Brown. Linear transformations and matrices. 2016. URL: https://www.youtube.com/
watch?v=kYB8IZa5AuE (cited on page 27).
[8] 3Blue1Brown. Matrix multiplication as composition. 2016. URL: https://www.youtube.com/
watch?v=XkY2DOUCWMU (cited on page 27).
[9] 3Blue1Brown. Nonsquare matrices as transformations between dimensions. 2016. URL: https:
//www.youtube.com/watch?v=v8VSDg_WQlA (cited on page 28).
[10] 3Blue1Brown. The determinant. 2016. URL: https://www.youtube.com/watch?v=Ip3X9LOh2dk
(cited on page 28).
[11] 3Blue1Brown. Vectors, what even are they? 2016. URL: https://www.youtube.com/watch?v=
fNk_zzaMoSs (cited on page 27).
[12] 3Blue1Brown. Essence of calculus. 2017. URL: https : / / www . youtube . com / playlist ? list =
PLZHQObOWTQDMsr9K-rj53DwVRMYO3t5Yr (cited on page 131).
[13] 3Blue1Brown. Taylor series. 2017. URL: https://www.youtube.com/watch?v=3d6DsjIBzJ4
(cited on page 53).
[15] Brian Douglas. Gain and Phase Margins Explained! 2015. URL: https://www.youtube.com/
watch?v=ThoA4amCAX4 (cited on page 195).
[27] Zach Star. What does the Laplace Transform really tell us? A visual explanation (plus applica-
tions). 2019. URL: https://www.youtube.com/watch?v=n2y7n6jw5d0 (cited on page 19).
[28] Russ Tedrake. 6.832 Underactuated Robotics. 2019. URL: https://www.youtube.com/channel/
UChfUOAhz7ynELF-s_1LPpWg/videos (cited on page 60).
Sunset near Porter Academic building at UCSC
Index