Preface

Does the world really need another introductory textbook on feedback control? For many years we have felt that the many high-quality textbooks that currently exist are more than sufficient. However, several years ago I, the first author (RWB), had an experience that changed my mind and the minds of my colleagues. I had just finished teaching the introductory feedback control course at BYU. I had a need for a new master’s student in my research group and hired one of the top students in the class. We had designed a gimbal system for a small unmanned air vehicle, and we needed a control system implemented on a microcontroller for the gimbal. I tasked this new master’s student to model the gimbal and design the control system. I was surprised by how much the student struggled with this task, especially understanding where to begin and how to model the gimbal. If I had given him a transfer function and asked him to design a PID controller, or if I had given him a state-space model and asked him to design an observer and controller for the system, he would not have had any difficulty. But he did not know how to do an end-to-end design that required developing models for the system, including physical constraints. It was this experience that convinced me that my current approach to teaching feedback control was inadequate.

The next time that I taught the course, rather than focusing on the theory of feedback design, I focused the lectures on the end-to-end design process. Accordingly, I spent significantly more time talking about the modeling process, adding several lectures on the Euler-Lagrange method. I also added several lectures on linearization about set points, and on developing linear models. When I taught feedback design methods, I focused on design specifications and the need to account for saturation, sensor noise, model uncertainty, and external disturbances on the system. Over the next several years, we developed several complete design studies and then used these design studies throughout the course to illustrate the material. That course reorganization and the realization that existing textbooks on feedback control did not fit the pedagogy were the genesis of this textbook.

Therefore, this textbook is unique in several ways. Most importantly, it is focused, in its organization, in its examples, and in its homework problems, on a particular, but fairly general, end-to-end design strategy and methodology. Topics were included in the book only if they aided in the design process. For example, the book does not include a chapter on the root locus (although it is discussed in an appendix). The second unique feature is that the examples and homework problems follow a small set of design studies throughout the book. The text provides complete worked examples for three design studies: a single-link robot arm, a pendulum on a cart, and a satellite attitude control problem. The homework problems parallel the examples with three additional design studies: a mass-spring-damper, a block sliding on a beam, and a planar vertical take-off and landing aircraft, similar to a quadrotor. A third unique feature of the book is that complete solutions to the three example design studies are provided in MATLAB/Simulink and Python at http://controlbook.byu.edu.

The design methodology culminates in computer code that implements the controllers. In the examples provided in the book, the computer code interacts with a simulation model of the plant composed of coupled nonlinear differential equations, rather than the physical plant. However, the example computer code is written in a way that can be directly implemented on a microcontroller, or other computing devices. We have used MATLAB/Simulink because many students are already familiar with these tools. But we also provide Python examples because Python is free of charge and therefore accessible to everyone, and because for many applications, Python can be used to implement the controller in hardware.

One criticism of the book might be that all of the design problems are mechanical systems, and that they are fairly similar to each other. The student may come away with the impression that feedback control and the design methods discussed in the book are only applicable to mechanical systems. Of course this is not true. The tools taught in the book are widely used in many applications including electrical power systems, disk drives, aerospace systems, chemical processing, biological systems, queueing systems, and many more. Although this book is focused on mechanical systems, we believe that the student who masters the material will be able to easily transition to other application domains. Mechanical systems, like those covered in this book, have a compelling pedagogical advantage because the control objectives are intuitive, and the behavior of the system can be easily visualized. Therefore, we feel that they offer the best setting for developing intuition behind feedback design.

We had two other motivations for writing this textbook. The first is the exorbitant prices currently being charged for textbooks. Accordingly, we are providing an electronic version of this book to students free of charge at http://controlbook.byu.edu. The book has required significant time and energy to develop, and so we hope that “free” does not equate to “low quality.” An Amazon version of the book is also available for those desiring a hard copy. Given the amount of high-quality open-source software, students may get the wrong impression about monetary compensation for technical work. Good work should be compensated. However, since our day jobs provide adequate financial remuneration, we ask to be compensated by your feedback about your educational experience with the book, including typos and other errors. This will allow us to update and improve the book. We also hope that you will understand that we will probably not be very responsive to email asking for help with the homework problems. Another motivation for writing the book is our observation that this generation of students seems to prefer learning through electronic resources. The book has therefore been written with the assumption that it will be studied electronically. Therefore, hyperlinks and other electronic aids are embedded throughout the text. We hope that these enhance the educational experience.

1 Introduction

Introduction

Overview

The objective of this book is to prepare the reader to design feedback control systems for dynamic systems and to lay the groundwork for further studies in the area. The design philosophy that we follow throughout the book is illustrated schematically in Fig. 1.1. The objective is to design a control system for the “Physical System” shown in Fig. 1.1. The physical system includes actuators like motors and propellers, and sensors like accelerometers, gyros, GPS, cameras, and pressure sensors. The first step in the design process is to model the physical system using nonlinear differential equations. While approximations and simplifications will be necessary at this step, the hope is to capture, in mathematics, all of the important characteristics of the physical system. While there are many different methods for developing the equations of motion of the physical systems, in this book we will introduce one specific method, the Euler-Lagrange method, that is applicable to a wide variety of mechanical systems. We will, for the most part, abstract the actuators to applied forces and torques on the system. Similarly, the sensors will be abstracted to general measurements of certain physical quantities like position and velocity. The resulting model is called the Simulation Model as shown in Fig. 1.1. In the design process, the simulation model is used for the high fidelity computer simulation of the physical system. Every control design will be tested thoroughly on the simulation model. However, the simulation model is only a mathematical approximation of the physical system, and simply because a design is effective on the simulation model, we should not assume that it will function properly on the physical system. Part I of the book introduces the Euler-Lagrange method and demonstrates how to derive simulation models for mechanical systems using this method.

1.1 Interactive Control Workflow

Click Figure 1.1

The diagram below is the same design-process figure, but the major regions are clickable. Select a region to see how that part of the workflow will be developed later in the book.

Physical system Simulation model Transfer-function design model State-space design model Control design methods Test in simulation Implementation

The simulation model is typically nonlinear and high order and is too mathematically complex to be useful for control design. Therefore, to facilitate design, the simulation model is simplified and usually linearized to create lower-order design models. For any physical system, there may be multiple design models that capture certain aspects of the system that are important for a particular design. In this book we will introduce the two most common design models used for control systems analysis and design: transfer-function models and state-space models. In both cases we will linearize the system about a particular operating condition. The advantage of transfer function models is that they capture the frequency dependence of the input-output relationship in dynamic systems. The advantage of state-space models is that they are more intuitive and are better adapted to numerical implementation. Part II of the book shows how to linearize the equations of motion and how to create transfer function and state-space design models.

The next three parts, Part III–V, introduce three different control design methods. Part III shows how to design Proportional-Integral-Derivative (PID) controllers for second-order systems. The PID method is the most common control technique used in industry. It can be understood from both a time-domain and frequency-domain perspective. The advantage of PID is its simplicity in design and implementation and the fact that it can be used to achieve high performance for a surprising variety of systems. However, the disadvantage of PID control is that stability and performance can typically only be guaranteed for second-order systems. Therefore, Part III is focused exclusively on second-order systems.

In Part IV we introduce the observer-based control method. Observer-based methods can be thought of as a natural extension of PID controllers to higher-order systems. However, observer-based methods are much more general than PID and can be used to obtain higher performance. In addition, most observer-based techniques extend to nonlinear and time-varying systems. The disadvantages of observer-based methods are that they require a higher level of mathematical sophistication to understand and use, and they are primarily time-domain methods.

The final design technique, which is covered in Part V, is the loop shaping control design method. The loop shaping method is a frequency domain technique that explicitly addresses robustness of the system. The loop shaping method is also a natural extension of PID control to higher-order systems. One of the advantages of loop shaping techniques is that they can be used when only a frequency response model of the system is available.

The book is organized around several design studies. For each new concept that is introduced, we will apply that concept to three specific design problems. The design problems are introduced in the following three sections and include a single-link robot arm, a pendulum on a cart, and a simple satellite attitude control problem. The homework problems, which are provided in Part VI, will follow the same format with problems that are identical to those worked in the book, but they will be applied to different physical systems. Our hope is to illustrate the control design process with specific, easy-to-understand problems. We note, however, that the concepts introduced in this book are applicable to a much wider variety of problems than those illustrated in the book. We focus on mechanical systems, but a similar process (with the exception of the modeling section) can be followed for any system that can be modeled using ordinary differential equations. Examples include aerodynamic systems, electrical systems, chemical processes, large structures like buildings and bridges, biological feedback systems, economic systems, queuing systems, and many more.

1.2 Design-Study Preview

1.3 Design Studies

1.3.1 Design Study: Single Link Robot Arm

Figure 1.2 shows a single-link robot arm system. The robot arm has mass \(m\) and length \(\ell\). The angle of the robot is given by \(\theta\) as measured from horizontal to the ground. The angular speed of the arm is \(\dot{\theta}\). There is an applied torque \(\tau\) at the joint. There is also a damping (or viscous friction) torque that opposes the rotation of the joint and which can be modeled as \(-b\dot{\theta}\).

The actual physical parameters of the system are \(m=0.5\) kg, \(\ell=0.3\) m, \(g=9.8\) m/s\(^2\), \(b=0.01\) Nms. The torque is limited by \(\left|\tau\right|\leq 1\) Nm.

A.2 Kinetic energy.

A.3 Equations of motion.

A.4 Linearize equations of motion.

A.5 Transfer function model.

A.6 State-space model.

A.7 Pole placement using PD.

A.8 Second-order design.

A.9 Integrators and system type.

A.P.6 Root locus.

A.10 Digital PID.

A.11 Full state feedback.

A.12 Full state with integrator.

A.13 Observer based control.

A.14 Disturbance observer.

A.15 Frequency response.

A.16 Loop gain.

A.17 Stability margins.

A.18 Loopshaping design.

Note: The notation identifies the design configuration and the chapter in which it is used: Example A.4 refers to the application of approximations in Chapter 4 to linearize design configuration A. Similarly, Example A.P.6 refers to the root locus assignment associated with Appendix A.P.6.

1.3.2 Design Study: Pendulum on a Cart

Figure 1.3 shows the pendulum on a cart system, also called an inverted pendulum. The position of the cart measured from the origin is \(z\) and the linear speed of the cart is \(\dot{z}\). The angle of the pendulum from straight up is given by \(\theta\), and the angular velocity is \(\dot{\theta}\). The rod is of length \(\ell\), and has mass \(m_1\), and is approximated as being infinitely thin. The cart has mass \(m_2\). Gravity acts in the down direction. The only applied force is \(F\), which acts in the direction of \(z\). The cart slides on a frictionless surface, but air friction produces a damping force equal to \(-b\dot{z}\).

The physical constants are \(m_1=0.25\) kg, \(m_2=1.0\) kg, \(\ell=1.0\) m, \(g=9.8\) m/s\(^2\), \(b=0.05\) Ns. The force is limited by \(\left|F\right|\leq 5\) N.

B.2 Kinetic energy.

B.3 Equations of Motion.

B.4 Linearize equations of motion.

B.5 Transfer function model.

B.6 State-space model.

B.8 Successive loop closure.

B.9 Integrators and system type.

B.P.6 Root locus.

B.10 Digital PID.

B.11 Full state feedback.

B.12 Full state with integrator.

B.13 Observer based control.

B.14 Disturbance observer.

B.15 Frequency Response.

B.16 Loop gain.

B.17 Stability margins.

B.18 Loopshaping design.

1.3.3 Design Study: Satellite Attitude Control

Figure 1.4 shows a simplified version of a satellite with flexible solar panels. We will model the flexible panels using a rigid sheet with moment of inertia \(J_p\) connected to the main satellite body by a torsional compliant element with spring constant \(k\) and damping constant \(b\). The moment of inertia of the satellite is \(J_s\). The angle of the satellite body from the inertial reference is \(\theta\) and the angle of the panel from the inertial reference is denoted as \(\phi\). Thrusters are used to command an external torque of \(\tau\) about the axis of the satellite body.

The physical constants are \(J_s=5\) kg m\(^2\), \(J_p=1\) kg m\(^2\), \(k=0.1\) N m, \(b=0.05\) Nms. The torque is limited by \(\left|\tau\right|\leq 5\) Nm.

C.2 Kinetic energy.

C.3 Equations of motion.

C.4 Equilibria.

C.5 Transfer function model.

C.6 State-space model.

C.8 Successive loop closure.

C.9 Integrators and system type.

C.P.6 Root locus.

C.10 Digital PID.

C.11 Full state feedback.

C.12 Full state with integrator.

C.13 Observer based control.

C.14 Disturbance observer.

C.15 Frequency Response.

C.16 Loop gain.

C.17 Stability margins.

C.18 Loopshaping design.

1.4 Knowledge Check

Use these quick checks to test the Chapter 1 workflow vocabulary before moving on. The short-answer checks are intentionally limited to terms that can be graded reliably.

To design efficient feedback control systems, we use mathematical models to understand the system, and then we design controllers that achieve the desired specifications. Developing suitable mathematical models is problem-specific and can be very complicated, requiring specialized knowledge. However, there are some well-known techniques for certain classes of systems. This book covers one general technique that is useful for a large class of rigid-body mechanical systems. The technique is known as the Euler-Lagrange method and requires a mathematical description of the kinetic and potential energy of the system, given its generalized position and velocity. The theory behind the Euler-Lagrange method is deep and well beyond the scope of this book. However, we will give a brief introduction that will enable the reader to apply Euler-Lagrange methods to a surprisingly large class of problems.

As motivation for the Euler-Lagrange equations, consider a particle of mass \(m\) in a gravity field, as shown in Fig. 1.5, where \(g\) is the force of gravity at sea level, \(f_a\) is an externally applied force, and the friction due to air resistance is modeled by a gain \(b\) times the velocity \(\dot{y}\).

Since all of the forces are in the vertical direction, the motion of the vehicle can be completely described by the vertical position of particle \(y(t)\). From Newton’s law, the equations of motion are \[\tag{1.1} m\ddot{y} = f_a-mg-b\dot{y},\] where dots over the variable will be used throughout the book to denote differentiation with respect to time, i.e., \(\dot{y}\stackrel{\triangle}{=}\frac{dy}{dt}\) and \(\ddot{y}\stackrel{\triangle}{=}\frac{d^2y}{dt^2}\). Rearranging 1.1 gives \[\tag{1.2} m\ddot{y} + mg = f_a -b\dot{y}.\] The kinetic energy for the particle is given by \(K(\dot{y})=\frac{1}{2}mv^2=\frac{1}{2}m\dot{y}^2\). The potential energy, which is induced by gravitational pull proportional to the distance from the center of the earth, is given by \(P(y)=P_0+mgy\), where \(P_0\) is the potential energy of the particle when \(y=0\). Note that if the velocity \(\dot{y}\) is thought of as independent of \(y\), then \[\frac{\partial K(\dot{y})}{\partial \dot{y}} = m\dot{y}\] and \[\frac{d}{dt}\left(\frac{\partial K(\dot{y})}{\partial \dot{y}}\right) = m\ddot{y}.\] Also note that \[\frac{\partial P(y)}{\partial y} = mg.\] Therefore, Equation 1.2 can be written as \[\tag{1.3} \frac{d}{dt}\left(\frac{\partial K(\dot{y})}{\partial \dot{y}}\right) + \frac{\partial P(y)}{\partial y} = f_a-b\dot{y}.\] Equation 1.3 is a particular form of the so-called Euler-Lagrange equation. The generalization to a more complicated system takes the form \[\tag{1.4} \frac{d}{dt}\left(\frac{\partial K(y,\dot{y})}{\partial \dot{y}}-\frac{\partial P(y)}{\partial \dot{y}}\right) - \left(\frac{\partial K(y,\dot{y})}{\partial y} - \frac{\partial P(y)}{\partial y}\right) = f_a-b\dot{y},\] which reduces to Equation 1.3 for our particular case since \(K\) is only a function of \(\dot{y}\). To simplify the notation, we define the so-called Lagrangian as \[L(y,\dot{y})\stackrel{\triangle}{=}K(y,\dot{y})-P(y),\] and write Equation 1.4 as \[\tag{1.5} \frac{d}{dt}\left(\frac{\partial L(y,\dot{y})}{\partial \dot{y}}\right) - \left(\frac{\partial L(y,\dot{y})}{\partial y}\right) = f_a-b\dot{y}.\] The position of the system can be generalized to both positions and angles, and the velocity can be generalized to both linear and angular velocities. Throughout the book we will use the notation \(\mathbf{q} = (q_1, q_2, \dots, q_n)^\top\) to denote the general configuration of the mechanical system that we are modeling. The vector \(\mathbf{q}\) is called the generalized-coordinate vector. The kinetic energy of the system is generally a function of both the generalized coordinates \(\mathbf{q}\) and the generalized velocity \(\dot{\mathbf{q}}\). Therefore, we denote the kinetic energy as \(K(\mathbf{q},\dot{\mathbf{q}})\). On the other hand, in this book we will only consider potential fields that are conservative and depend only on the configuration variable \(\mathbf{q}\). Under these assumptions, the general form of the Euler-Lagrange equations is given by \begin{align} \frac{d}{dt}\left(\frac{\partial L(\mathbf{q},\dot{\mathbf{q}})}{\partial \dot{q}_1}\right) - \frac{\partial L(\mathbf{q},\dot{\mathbf{q}})}{\partial q_1} &= \tau_1-b_{11} \dot{q}_1 - \cdots - b_{1n}\dot{q}_n \notag \\ \frac{d}{dt}\left(\frac{\partial L(\mathbf{q},\dot{\mathbf{q}})}{\partial \dot{q}_2}\right) - \frac{\partial L(\mathbf{q},\dot{\mathbf{q}})}{\partial q_2} &= \tau_2-b_{21} \dot{q}_1 - \cdots - b_{2n}\dot{q}_n \tag{1.6} \\ \vdots \\ \frac{d}{dt}\left(\frac{\partial L(\mathbf{q},\dot{\mathbf{q}})}{\partial \dot{q}_n}\right) - \frac{\partial L(\mathbf{q},\dot{\mathbf{q}})}{\partial q_n} &= \tau_n-b_{n1} \dot{q}_1 - \cdots - b_{nn}\dot{q}_n, \notag \end{align} where \(\tau_i\) are the generalized forces and torques that are applied to the system. Or in vector form we can write \[\tag{1.8} \frac{d}{dt}\left(\frac{\partial L(\mathbf{q},\dot{\mathbf{q}})}{\partial \dot{\mathbf{q}}}\right) - \frac{\partial L(\mathbf{q},\dot{\mathbf{q}})}{\partial \mathbf{q}} = \boldsymbol{\tau} -B \dot{\mathbf{q}},\] where \(B\) is a matrix with positive elements and \(\boldsymbol{\tau}=(\tau_1,\dots,\tau_n)^\top\) are the generalized forces. The Euler-Lagrange equation will be discussed in more detail in Chapter 6. However, before showing how to use Equation 1.8 to derive the equations of motion for practical systems, we need to show how to compute the kinetic and potential energy for general (holonomic) rigid body systems.

2 Kinetic Energy of Mechanical Systems

2.1 Theory

Consider a point mass with mass \(m\), traveling at velocity \(\mathbf{v}\), as shown in Fig. 2.1. If the velocity is resolved with respect to a coordinate axis, then \(\mathbf{v}=(v_x, v_y, v_z)^\top\). The kinetic energy of the point mass is given by \[\begin{aligned} K &= \frac{1}{2} m \left|\mathbf{v}\right|^2 \\ &= \frac{1}{2} m \mathbf{v}^\top \mathbf{v} \\ & = \frac{1}{2} m (v_x^2 + v_y^2 + v_z^2). \end{aligned}\]

Now consider two point masses as shown in Fig. 2.2.

The total kinetic energy of the system is given by summing the kinetic energy of each of the point masses as \[\tag{2.1} K = \frac{1}{2} m_1 \mathbf{v}_1^\top \mathbf{v}_1 + \frac{1}{2} m_2 \mathbf{v}_2^\top \mathbf{v}_2.\] If the two masses are connected by a massless rod and the particles are spinning about a point on the rod, as shown in Fig. 2.2, where the angular velocity vector \(\boldsymbol{\omega}\) points out of the page, and the magnitude is given by \(\dot{\theta}\), i.e., \(\boldsymbol{\omega}=(0, 0, \dot{\theta})^\top\), then from elementary physics, the velocity of point \(i\) is given by \[\mathbf{v}_i = \boldsymbol{\omega}\times\mathbf{p}_i,\] where in our case \(\mathbf{p}_1 = (L_1\cos\theta, L_1\sin\theta, 0)^\top\), and \(\mathbf{p}_2 = (-L_2\cos\theta, -L_2\sin\theta, 0)^\top\). Therefore \[\begin{aligned} \mathbf{v}_1 &= \begin{pmatrix} 0 \\ 0 \\ \dot{\theta}\end{pmatrix}\times\begin{pmatrix}L_1\cos\theta \\ L_1\sin\theta \\ 0\end{pmatrix} = \begin{pmatrix} -L_1\dot{\theta}\sin\theta \\ L_1\dot{\theta}\cos\theta \\ 0 \end{pmatrix} \\ \mathbf{v}_2 &= \begin{pmatrix} 0 \\ 0 \\ \dot{\theta} \end{pmatrix}\times\begin{pmatrix}-L_2\cos\theta \\ -L_2\sin\theta \\ 0\end{pmatrix} = \begin{pmatrix} L_2\dot{\theta}\sin\theta \\ -L_2\dot{\theta}\cos\theta \\ 0 \end{pmatrix}. \end{aligned}\] An alternative method for computing the velocities \(\mathbf{v}_1\) and \(\mathbf{v}_2\) is to differentiate \(\mathbf{p}_1\) and \(\mathbf{p}_2\) with respect to time to obtain \[\begin{aligned} \mathbf{v}_1 &\stackrel{\triangle}{=}\dot{\mathbf{p}}_1 = \frac{d\mathbf{p}_1}{dt} = \frac{d}{dt}(L_1\cos\theta, L_1\sin\theta, 0)^\top \\ &= \begin{pmatrix}-L_1\left(\frac{d\theta}{dt}\right)\sin\theta, & L_1\left(\frac{d\theta}{dt}\right)\cos\theta, & 0 \end{pmatrix}^\top \\ &= \begin{pmatrix}-L_1\dot{\theta}\sin\theta, & L_1\dot{\theta}\cos\theta, & 0 \end{pmatrix}^\top \\ \mathbf{v}_2 &\stackrel{\triangle}{=}\dot{\mathbf{p}}_2 = \frac{d\mathbf{p}_2}{dt} = \frac{d}{dt}(-L_2\cos\theta, -L_2\sin\theta, 0)^\top \\ &= \begin{pmatrix}L_2\left(\frac{d\theta}{dt}\right)\sin\theta, & -L_2\left(\frac{d\theta}{dt}\right)\cos\theta, & 0 \end{pmatrix}^\top \\ &= \begin{pmatrix}L_2\dot{\theta}\sin\theta, & -L_2\dot{\theta}\cos\theta, & 0 \end{pmatrix}^\top. \end{aligned}\]

Therefore, using Equation 2.1 the kinetic energy is \[\begin{aligned} K &= \frac{1}{2} m_1 \mathbf{v}_1^\top \mathbf{v}_1 + \frac{1}{2} m_2 \mathbf{v}_2^\top \mathbf{v}_2 \\ &= \frac{1}{2} m_1 \left[ (-L_1\dot{\theta}\sin\theta)^2 + (L_1\dot{\theta}\cos\theta)^2 \right] \\ &\qquad + \frac{1}{2} m_2 \left[ (L_2\dot{\theta}\sin\theta)^2 + (-L_2\dot{\theta}\cos\theta)^2 \right] \\ &= \frac{1}{2} m_1 \left[ L_1^2\dot{\theta}^2\sin^2\theta + L_1^2\dot{\theta}^2\cos^2\theta \right] % + \frac{1}{2} m_2 \left[ L_2^2\dot{\theta}^2\sin^2\theta + L_2^2\dot{\theta}^2\cos^2\theta \right] \\ &= \frac{1}{2} (m_1 L_1^2 + m_2 L_2^2)\dot{\theta}^2. \end{aligned}\]

Now consider the case where there are \(N\) point masses spinning about a general angular velocity vector \(\boldsymbol{\omega}\), where the position of the \(i^{th}\) point mass, \(i=1,\dots,N\), in a specified coordinate system is given by \(\mathbf{p}_i\), as shown in Fig. 2.3.

We will see in the discussion below that it is most convenient to pick the coordinate system to be at the center of mass. The kinetic energy of the \(i^{th}\) particle is \(K_i = \frac{1}{2}m_i\mathbf{v}_i^\top\mathbf{v}_i\). Therefore, the kinetic energy of the total system is given by \[K = \sum_{i=1}^N \frac{1}{2}m_i\mathbf{v}_i^\top\mathbf{v}_i.\] Using the fact that \(\mathbf{v}_i = \boldsymbol{\omega}\times\mathbf{p}_i\) gives \[K = \sum_{i=1}^N \frac{1}{2}m_i\left(\boldsymbol{\omega}\times\mathbf{p}_i\right)^\top\left(\boldsymbol{\omega}\times\mathbf{p}_i\right).\] Recalling the cross product rule \[(\mathbf{a}\times\mathbf{b})^\top(\mathbf{c}\times\mathbf{d}) = (\mathbf{a}^\top\mathbf{c})(\mathbf{b}^\top\mathbf{d})-(\mathbf{a}^\top\mathbf{d})(\mathbf{b}^\top\mathbf{c})\] gives \[\begin{aligned} K &= \frac{1}{2} \sum_{i=1}^N m_i \left[ (\boldsymbol{\omega}^\top\boldsymbol{\omega})(\mathbf{p}_i^\top\mathbf{p}_i)-(\boldsymbol{\omega}^\top\mathbf{p}_i)(\mathbf{p}_i^\top\boldsymbol{\omega})\right] \\ &= \frac{1}{2} \sum_{i=1}^N m_i \boldsymbol{\omega}^\top \left[ (\mathbf{p}_i^\top\mathbf{p}_i) I_3 - \mathbf{p}_i \mathbf{p}_i^\top \right] \boldsymbol{\omega} \\ &= \frac{1}{2} \boldsymbol{\omega}^\top \left[ \sum_{i=1}^N m_i \left((\mathbf{p}_i^\top\mathbf{p}_i) I_3 - \mathbf{p}_i \mathbf{p}_i^\top \right) \right] \boldsymbol{\omega} \\ &\stackrel{\triangle}{=}\frac{1}{2} \boldsymbol{\omega}^\top J \boldsymbol{\omega}. \end{aligned}\] where \(I_3\) is the \(3 \times 3\) identity matrix. The matrix \[J = \sum_{i=1}^N m_i \left((\mathbf{p}_i^\top\mathbf{p}_i) I_3 - \mathbf{p}_i \mathbf{p}_i^\top \right)\] is called the inertia matrix of the mass system. If the mass system is solid then there is an infinite number of particles, and the sum becomes an integral resulting in \[J = \int_{\mathbf{p}\in\text{body}} \left((\mathbf{p}^\top\mathbf{p}) I_3 - \mathbf{p} \mathbf{p}^\top \right)dm(\mathbf{p}),\] where \(dm(\mathbf{p})\) is the elemental mass located at position \(\mathbf{p}\) and the integral is over all elemental masses in the system.

Suppose that \(\mathbf{p}=(x, y, z)^\top\) is expressed with respect to a coordinate system fixed in a rigid body. Then \[\begin{aligned} J &= \int_{(x,y,z)\in\text{body}} \Bigg[ \begin{pmatrix} x & y & z \end{pmatrix}\begin{pmatrix} x \\ y \\ z \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} % - \begin{pmatrix} x \\ y \\ z \end{pmatrix} \begin{pmatrix} x & y & z \end{pmatrix} \Bigg] dm(x,y,z) \\ &= \int_{(x,y,z)\in\text{body}} \Bigg[ (x^2 + y^2 + z^2)\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} % - \begin{pmatrix} x^2 & xy & x z \\ yx & y^2 & yz \\ zx & zy & z^2 \end{pmatrix} \Bigg] dm(x,y,z) \\ &= \int_{(x,y,z)\in\text{body}} \begin{pmatrix} y^2 + z^2 & -xy & -x z \\ -yx & x^2+z^2 & -yz \\ -zx & -zy & x^2+y^2 \end{pmatrix} dm(x,y,z). \end{aligned}\]

As an example, to compute the inertia matrix of the thin rod of mass \(m\) about its center of mass, as shown in Fig. 2.4, the first step is to define a coordinate system at the center of mass, as shown in Fig. 2.4, where \(\hat{i}\) and \(\hat{j}\) are orthogonal unit vectors pointing parallel and perpendicular to the rod, respectively. The unit vector \(\hat{k}\) forms a right-handed coordinate system and points out of the page.

Assuming that the rod is infinitely thin and that the mass is evenly distributed throughout the rod, then the elemental mass is given by \[dm(x) = \frac{m}{\ell}dx.\] Since the rod is infinitely thin in the \(\hat{j}\) and \(\hat{k}\) directions, the inertia matrix is given by \begin{align} J &= \int_{x=-\ell/2}^{\ell/2}\int_{y=0}^0\int_{z=0}^0 \begin{pmatrix} y^2 + z^2 & -xy & -x z \\ -yx & x^2+z^2 & -yz \\ -zx & -zy & x^2+y^2 \end{pmatrix} dz \, dy \, \frac{m}{\ell}dx \notag\\ &= \frac{m}{\ell} \int_{x=-\ell/2}^{\ell/2} \begin{pmatrix} 0 & 0 & 0 \\ 0 & x^2 & 0 \\ 0 & 0 & x^2 \end{pmatrix} dx \notag\\ &= \frac{m}{\ell} \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{x^3}{3}\big|_{x=-\frac{\ell}{2}}^{\frac{\ell}{2}} & 0 \\ 0 & 0 & \frac{x^3}{3}\big|_{x=-\frac{\ell}{2}}^{\frac{\ell}{2}} \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{m\ell^2}{12} & 0 \\ 0 & 0 & \frac{m\ell^2}{12} \end{pmatrix}. \tag{2.2} \end{align} Note that the infinitely thin rod assumption implies that there is only inertia, or rotational mass, about the \(\hat{j}\) and \(\hat{k}\) axes.

As another example, to compute the inertia matrix of the thin plate of mass \(m\) about its center of mass, as shown in Fig. 2.5, the first step is to define a coordinate system at the center of mass, as shown in Fig. 2.5, where \(\hat{i}\) and \(\hat{j}\) are orthogonal unit vectors pointing along the length and width of the plate, respectively, and \(\hat{k}\) forms a right-handed coordinate system and points perpendicular to the plate.

Assuming that the plate is infinitely thin, and that the mass is evenly distributed throughout the plate, then the elemental mass is given by \[dm(x,y) = \frac{m}{\ell w}dx\,dy.\] Since the plate is infinitely thin in the \(\hat{k}\) direction, the inertia matrix is given by \[\begin{aligned} J &= \int_{x=-\ell/2}^{\ell/2}\int_{y=-w/2}^{w/2}\int_{z=0}^0 \begin{pmatrix} y^2 + z^2 & -xy & -x z \\ -yx & x^2+z^2 & -yz \\ -zx & -zy & x^2+y^2 \end{pmatrix} dz \, dy \, \frac{m}{\ell w}\,dx \\ &= \frac{m}{\ell w} \int_{x=-\ell/2}^{\ell/2} \int_{y=-w/2}^{w/2} \begin{pmatrix} y^2 & -xy & 0 \\ -yx & x^2 & 0 \\ 0 & 0 & x^2 +y^2\end{pmatrix} dy\,dx \\ % &= \frac{m}{\ell w} \int_{x=-\ell/2}^{\ell/2} \begin{pmatrix} \frac{w^3}{12} & 0 & 0 \\ 0 & x^2 w & 0 \\ 0 & 0 & x^2 w +\frac{w^3}{12}\end{pmatrix} dx \\ % &= \frac{m}{\ell w} \begin{pmatrix} \frac{w^3}{12}\ell & 0 & 0 \\ 0 & \frac{\ell^3}{12} w & 0 \\ 0 & 0 & \frac{\ell^3}{12}w +\frac{w^3}{12}\ell\end{pmatrix} = \begin{pmatrix} \frac{mw^2}{12} & 0 & 0 \\ 0 & \frac{m\ell^2}{12} & 0 \\ 0 & 0 & \frac{m\ell^2}{12} +\frac{mw^2}{12}\end{pmatrix}. \end{aligned}\]

The diagonal elements of \(J\) are called the moments of inertia and represent the rotational mass about the \(\hat{i}\), \(\hat{j}\), and \(\hat{k}\) axes respectively.

We have shown that if the motion of the system of \(N\) masses is caused purely by spinning about the angular velocity vector \(\boldsymbol{\omega}\), then the kinetic energy is given by \(K=\frac{1}{2}\boldsymbol{\omega}^\top J \boldsymbol{\omega}\). If the system of masses is also translating with velocity \(\mathbf{v}\), then the velocity of the \(i^{th}\) particle is given by \(\mathbf{v}_i = \mathbf{v} + \boldsymbol{\omega}\times\mathbf{p}_i.\) In that case, the kinetic energy is given by \[\begin{aligned} K &= \frac{1}{2}\sum_{i=1}^N m_i \mathbf{v}_i^\top\mathbf{v}_i = \frac{1}{2}\sum_{i=1}^N m_i (\mathbf{v} + \boldsymbol{\omega}\times\mathbf{p}_i)^\top(\mathbf{v} + \boldsymbol{\omega}\times\mathbf{p}_i) \\ &= \frac{1}{2}\sum_{i=1}^N m_i \Big(\mathbf{v}^\top\mathbf{v} + \mathbf{v}^\top(\boldsymbol{\omega}\times\mathbf{p}_i)+(\boldsymbol{\omega}\times\mathbf{p}_i)^\top\mathbf{v} % + (\boldsymbol{\omega}\times\mathbf{p}_i)^\top(\boldsymbol{\omega}\times\mathbf{p}_i) \Big) \\ &= \frac{1}{2}\left(\sum_{i=1}^N m_i\right)\mathbf{v}^\top\mathbf{v} + \frac{1}{2} \sum_{i=1}^N m_i (\boldsymbol{\omega}\times\mathbf{p}_i)^\top(\boldsymbol{\omega}\times\mathbf{p}_i) \\ &\qquad\qquad + \frac{1}{2}\sum_{i=1}^N m_i\left[\mathbf{v}^\top(\boldsymbol{\omega}\times\mathbf{p}_i) + (\boldsymbol{\omega}\times\mathbf{p}_i)^\top\mathbf{v}\right]. \end{aligned}\] Using the cross product property \[\mathbf{a}^\top(\mathbf{b}\times\mathbf{c}) = \mathbf{b}^\top(\mathbf{c}\times\mathbf{a}) = \mathbf{c}^\top(\mathbf{a}\times\mathbf{b})\] we get that \[\begin{gathered} \sum_{i=1}^N m_i\left[\mathbf{v}^\top(\boldsymbol{\omega}\times\mathbf{p}_i) + (\boldsymbol{\omega}\times\mathbf{p}_i)^\top\mathbf{v}\right] = \\ (\sum_{i=1}^N m_i\mathbf{p}_i)^\top(\mathbf{v}\times\boldsymbol{\omega}) + (\sum_{i=1}^N m_i\mathbf{p}_i)^\top(\mathbf{v}\times\boldsymbol{\omega}). \end{gathered}\] If we choose the center of the coordinate system to be the center of mass, then \(\sum_{i=1}^N m_i\mathbf{p}_i=0\). Defining the total mass as \(m=\sum_{i=1}^N m_i\), the kinetic energy is given by \[\begin{aligned} K&= \frac{1}{2}\left(\sum_{i=1}^N m_i\right)\mathbf{v}_{cm}^\top\mathbf{v}_{cm} + \frac{1}{2} \sum_{i=1}^N m_i (\boldsymbol{\omega}\times\mathbf{p}_i)^\top(\boldsymbol{\omega}\times\mathbf{p}_i) \notag \\ &= \frac{1}{2}m\mathbf{v}_{cm}^\top\mathbf{v}_{cm} + \frac{1}{2} \boldsymbol{\omega}^\top J_{cm} \boldsymbol{\omega}. \end{aligned}\]

Therefore, the kinetic energy of a rigid body moving with velocity \(\mathbf{v}_{cm}\) and spinning at angular velocity \(\boldsymbol{\omega}\) is \[\tag{2.3} K=\frac{1}{2}m\mathbf{v}_{cm}^\top\mathbf{v}_{cm} + \frac{1}{2} \boldsymbol{\omega}^\top J_{cm} \boldsymbol{\omega}.\] If the system consists of \(n\) rigid bodies, then the total kinetic energy is computed by summing the kinetic energy for each body as \[K = \sum_{j=1}^n \left[\frac{1}{2}m_j\mathbf{v}_{cm,j}^\top\mathbf{v}_{cm,j} + \frac{1}{2} \boldsymbol{\omega}_j^\top J_{cm,j} \boldsymbol{\omega}_j\right].\]

2.1.1 Example: Mass Spring System

As an example, consider the mass spring system shown in Fig. 2.6 where two masses are moving along a surface and are connected by springs with constants \(k_1\) and \(k_2\) and dampers with constants \(b_1\) and \(b_2\). The kinetic energy of the system is determined by the velocity of the two masses and is therefore given by \[K = \frac{1}{2} m_1 \dot{z}_1^2 + \frac{1}{2} m_2 \dot{z}_2^2.\]

2.1.2 Example: Spinning Dumbbell

As another simple example, consider the spinning dumbbell system shown in Fig. 2.7, where two masses are spinning about a point on the adjoining line between the masses. The angular velocity vector is pointing out of the page with magnitude \(\dot{\theta}\) and the dumbbell system is moving to the right at velocity \(\dot{z}\).

The first step in computing the kinetic energy is to define an inertial coordinate system as shown in Fig. 2.7 where \(\hat{i}\) is a unit vector along the direction of travel, \(\hat{j}\) is perpendicular to the direction of travel, and where \(\hat{k}\) forms a right-handed coordinate system and points out of the page. In this coordinate system, the position of the point of rotation is given by \(\mathbf{p}(t)=(z(t), 0, 0)^\top\). Accordingly, the positions of \(m_1\) and \(m_2\) are given by \[\begin{aligned} \mathbf{p}_1(t) &= (z(t)+\ell_1\cos\theta(t), \ell_1\sin\theta(t), 0)^\top, \\ \mathbf{p}_2(t) &= (z(t)-\ell_2\cos\theta(t), -\ell_2\sin\theta(t), 0)^\top. \end{aligned}\] Since \(z\) and \(\theta\) are functions of time, using the chain rule \[\frac{d}{dt}f(g(t)) = \frac{\partial f}{\partial g}\frac{d g}{dt}\] and differentiating with respect to time, we get the velocities of the masses as \[\begin{aligned} \mathbf{v}_1 &= (\dot{z} - \ell_1\dot{\theta}\sin\theta, \ell_1\dot{\theta}\cos\theta, 0)^\top,\\ \mathbf{v}_2 &= (\dot{z} + \ell_2\dot{\theta}\sin\theta, -\ell_2\dot{\theta}\cos\theta, 0)^\top. \end{aligned}\] The kinetic energy is therefore given by \[\begin{aligned} K &= \frac{1}{2}m_1\mathbf{v}_1^\top \mathbf{v}_1 + \frac{1}{2}m_2\mathbf{v}_2^\top \mathbf{v}_2 \\ &= \frac{1}{2}m_1 \left[(\dot{z} - \ell_1\dot{\theta}\sin\theta)^2 + (\ell_1\dot{\theta}\cos\theta)^2 \right] \\&\qquad + \frac{1}{2}m_2\left[(\dot{z} + \ell_2\dot{\theta}\sin\theta)^2 + (-\ell_2\dot{\theta}\cos\theta)^2 \right] \\ &= \frac{1}{2}m_1 \left[\dot{z}^2 - 2\ell_1\dot{z}\dot{\theta}\sin\theta + \ell_1^2\dot{\theta}^2\sin^2\theta+\ell_1^2\dot{\theta}^2\cos^2\theta \right] \\ &\quad + \frac{1}{2}m_2\left[\dot{z}^2 + 2\ell_2\dot{z}\dot{\theta}\sin\theta + \ell_2^2\dot{\theta}^2\sin^2\theta+\ell_2^2\dot{\theta}^2\cos^2\theta \right] \\ &= \frac{1}{2}(m_1+m_2)\dot{z}^2 + \frac{1}{2} (m_1\ell_1^2+m_2\ell_2^2)\dot{\theta}^2 + (m_2\ell_2-m_1\ell_1)\dot{z}\dot{\theta}\sin\theta. \end{aligned}\] Note that if the pivot point is located at the center of mass, then \(m_1\ell_1=m_2\ell_2\), and the kinetic energy is simply given by \[K = \frac{1}{2}(m_1+m_2)\dot{z}^2 + \frac{1}{2} (m_1\ell_1^2+m_2\ell_2^2)\dot{\theta}^2,\] which is in the form \[K = \frac{1}{2}m\mathbf{v}_{cm}^\top\mathbf{v}_{cm} + \frac{1}{2}\boldsymbol{\omega}^\top J \boldsymbol{\omega},\] where the total mass is \(m=m_1+m_2\), the velocity of the center of mass is \(\mathbf{v}_{cm}=(\dot{z}, 0, 0)^\top\), the angular velocity vector is \(\boldsymbol{\omega} = (0, 0, \dot{\theta})^\top\) and the \((3,3)\) element of \(J\) is given by \(J_z = m_1\ell_1^2 + m_2\ell_2^2\).

2.2 Design Study A: Single Link Robot Arm

A definition of the single-link robot arm is given in Design Study A.

2.2.1 Example Problem 3.1.2

(a) Using the configuration variable \(\theta\), write an expression for the kinetic energy of the system.

(b) Referring to Appendices P.1, P.2, and P.3, write a Python or MATLAB class, or a MATLAB function that creates an animation of the single-link robot arm. Simulate the animation to display a sinusoidal variation on the configuration variable \(q=\theta\).

Design Study A guided exercise

Derive the arm kinetic energy

Work from the center-of-mass motion and the rotational inertia of the slender link.

Use this geometry to identify the center-of-mass position and the angular velocity before deriving kinetic energy.

Step 1. Center-of-mass position

Enter the horizontal and vertical components of \(p_{cm}\).

\(p_x\) \(p_y\)

Check Step 1

Step 2. Link inertia

Enter the moment of inertia about the center of mass for the slender link.

\(J_z\)

Check Step 2

Step 3. Total kinetic energy

Enter the multiplier on \(\dot{\theta}^2\) in the final kinetic energy expression.

Multiplier

Check Step 3

Pseudocode

The complete solution includes an animation class and a simulation script. Before opening it, outline each component separately.

These prompts are ungraded. They are intended to connect the derivation to the implementation code in the complete solution.

Show complete solution

2.3 Design Study B. Pendulum on Cart

A definition of the pendulum on a cart (or inverted pendulum) is given in Design Study 3.2.

2.3.1 Example Problem 3.2.2

(a) Using the configuration variables \(z\) and \(\theta\), write an expression for the kinetic energy of the system.

(b) Referring to Appendices P.1, P.2, and P.3, write a Python or MATLAB class, or a MATLAB function that creates an animation of the inverted pendulum. Simulate the animation to display sinusoidal variations on the configuration variables \(q=(z, \theta)^\top\).

Design Study B guided exercise

Build the pendulum kinetic energy

Use the geometry to identify positions, velocities, and energy terms for the cart and rod.

Use this geometry to identify the cart position, rod center-of-mass position, and pendulum angular velocity.

Step 1. Position variables

Rod center of mass includes Cart position

Check Step 1

Step 2. Velocity variables

Rod velocity depends on Cart velocity

Check Step 2

Step 3. Energy terms

Energy types included

Check Step 3

Pseudocode

The complete solution includes an animation class and a simulation script. Before opening it, outline each component separately.

These prompts are ungraded. They are intended to connect the derivation to the implementation code in the complete solution.

Show complete solution

2.4 Design Study C. Satellite Attitude Control

A definition of the satellite attitude control problem is given in Design Study C.

2.4.1 Example Problem 3.3.2

(a) Using the configuration variables \(\theta\) and \(\phi\), write an expression for the kinetic energy of the system.

(b) Referring to Appendices P.1, P.2, and P.3, write a Python or MATLAB class, or a MATLAB function that creates an animation of the satellite. Simulate the animation to display sinusoidal variations on the configuration variables \(q=(\theta, \phi)^\top\).

Design Study C guided exercise

Build the satellite kinetic energy

The satellite model has two rotational inertias: the body inertia \(J_s\) associated with \(\theta\), and the panel inertia \(J_p\) associated with \(\phi\).

Use inertial angular rates for kinetic energy. The relative deflection \(\phi-\theta\) appears in spring and damping terms later.

Step 1. Match inertias to coordinates

Body inertia coordinate Panel inertia coordinate

Check Step 1

Step 2. Identify the kinetic-energy terms

Number of rotational terms Relative rate used directly?

Check Step 2

Step 3. Enter the total kinetic-energy structure

Body term coefficient Panel term coefficient

Check Step 3

Pseudocode

The complete solution includes satellite animation code. Before opening it, outline how an animation should update the body and panel orientations from \(\theta(t)\) and \(\phi(t)\).

This prompt is ungraded. It connects the energy model to the animation code in the complete solution.

Show complete solution

2.5 Interactive Inertia Check

2.6 Knowledge Check

2.7 Notes and References

The inertia matrix in Equation 2.3 is computed about the center of mass of the object. When an object is composed of several connected rigid bodies and the inertia matrix is known for each body, the parallel axis theorem provides a convenient way to find the total inertia about the center of mass. The theorem states that if \(J_{cm}\) is the inertia matrix about the center of mass, then the inertia matrix about a point located at vector \(\mathbf{r}\) from the center of mass is \[J(\mathbf{r}) = J_{cm} + m\left(\mathbf{r}^\top\mathbf{r} I_3 - \mathbf{r}\mathbf{r}^\top\right).\] For the moment of inertia about a single axis, where the axis is shifted from the center of mass by distance \(d\), the parallel axis theorem implies that \[I(d) = I_{cm} + md^2.\]

3 The Euler-Lagrange Equations

3.1 Theory

In this chapter we show how to define the potential energy for simple mechanical systems and how to define the generalized forces. The Euler-Lagrange equations given in Equation 1.8 can then be used to derive the equations of motion for the system.

3.1.1 Potential Energy

In this book we will be concerned with two sources of potential energy: potential energy due to gravity and potential energy due to springs.

Consider a point mass \(m\) in a gravity field that is at position \(y\) as shown in Fig. 3.1. The potential energy is given by \[P = mgy + P_0,\] where \(P_0\) is the potential energy when \(y=0\). Typically, we desire an expression for potential energy that is zero when the configuration variables are zero.

The second source of potential energy is from translational and rotational springs.

The spring shown in Fig. 3.2 is stretched from rest by a distance of \(z\). The potential energy of a linear stretched/compressed spring is given by \[P = \frac{1}{2}k z^2,\] where \(k\) is the spring constant. Similarly, for a rotational spring, the potential energy is given by \(P=\frac{1}{2}k \theta^2\), where \(\theta\) is the rotational angle and \(\theta=0\) is the angle at which the rotational spring is at rest.

3.1.2 Generalized Coordinates and Generalized Forces

The generalized coordinates are the minimum set of configuration variables. For example, for the mass-spring-damper system shown in Fig. 2.6, the generalized coordinates are \(\mathbf{q}=(z_1, z_2)^\top\). For the single-link robot arm shown in Fig. 2.8, the generalized coordinate is \(\mathbf{q}=\theta\). For the spinning dumbbell system shown in Fig. 2.7, and the pendulum on a cart system shown in Fig. 2.9, the generalized coordinates are \(\mathbf{q}=(z, \theta)^\top\). Finally, for the satellite system shown in Fig. 1.4, the generalized coordinates are \(\mathbf{q}=(\theta,\phi)^\top\).

The generalized forces are the nonconservative forces and torques acting along each generalized coordinate. In other words, they are the applied forces and torques on the system. For example, for the mass-spring-damper system shown in Fig. 2.6, the generalized force acting along \(z_1\) is zero and the generalized force acting along \(z_2\) is the applied force \(F\). Therefore, the generalized force is \(\boldsymbol{\tau} = (0, F)^\top\).

3.1.3 Damping Forces and Torques

Friction and applied dampers also exert forces on the system. Damping forces are externally applied forces that are proportional to the velocity or angular velocity of the system. For example, in the mass-spring-damper system shown in Fig. 2.6, the damping force acting along \(z_1\) is due to the first damper and the damper between the two masses and is given by \(-b_1\dot{z}_1-b_2(\dot{z}_1-\dot{z}_2)\). The damping acting along \(z_2\) is given by \(-b_2(\dot{z}_2-\dot{z}_1)\). Therefore, the damping force on the system is given by \[-B\dot{\mathbf{q}} = -\begin{pmatrix} b_1+b_2 & -b_2 \\ -b_2 & b_2 \end{pmatrix} \begin{pmatrix} \dot{z}_1 \\ \dot{z}_2 \end{pmatrix} = \begin{pmatrix} -b_1\dot{z}_1-b_2(\dot{z}_1-\dot{z}_2) \\ -b_2(\dot{z}_2-\dot{z}_1) \end{pmatrix}.\]

3.1.4 Euler-Lagrange Equations

The equations of motion are given by the Euler-Lagrange equations. Let the Lagrangian be given by \[L(\mathbf{q},\dot{\mathbf{q}}) = K(\mathbf{q},\dot{\mathbf{q}}) - P(\mathbf{q}),\] where \(K(\mathbf{q},\dot{\mathbf{q}})\) is the kinetic energy of the system and is, in general, a function of the generalized coordinates and the generalized velocities, and \(P(\mathbf{q})\) is the potential energy which is, in general, a function of the generalized coordinates. As shown in Equation 1.8 the Euler-Lagrange equations are given in vector form as \[\frac{d}{dt}\left(\frac{\partial L(\mathbf{q},\dot{\mathbf{q}})}{\partial \dot{\mathbf{q}}}\right) - \frac{\partial L(\mathbf{q},\dot{\mathbf{q}})}{\partial \mathbf{q}} = \boldsymbol{\tau}-B \dot{\mathbf{q}}.\]

3.2 Design Study A. Single Link Robot Arm

3.2.1 Example Problem 3.1.6

Find the potential energy for the system.

Define the generalized coordinates.

Find the generalized forces and damping forces.

Derive the equations of motion using the Euler-Lagrange equations.

Referring to Appendices P.1, P.2, and P.3, write a class or s-function that implements the equations of motion. Simulate the system using a variable torque input. The output should connect to the animation function developed in homework A.2.

Design Study A guided exercise

Build the arm equation of motion

Use the kinetic energy from Chapter 2 and gravity potential energy for the center of mass.

Use this geometry to relate the center-of-mass height to \(\theta\) before applying Euler-Lagrange.

Step 1. Potential energy

\(P-P_0\)

Check Step 1

Step 2. Generalized coordinate and force

\(q\) Applied generalized force Damping torque

Check Step 2

Step 3. Equation of motion

Enter the multiplier on \(\ddot{\theta}\) in the final equation.

Inertia coefficient

Check Step 3

Pseudocode

The complete solution includes a dynamics class and a simulation script. Before opening it, outline each component separately.

These prompts are ungraded. They are intended to connect the Euler-Lagrange equation to the implementation code in the complete solution.

Show complete solution

3.3 Design Study B. Inverted Pendulum

3.3.1 Example Problem 3.2.6

Find the potential energy for the system.

Define the generalized coordinates.

Find the generalized forces and damping forces.

Derive the equations of motion using the Euler-Lagrange equations.

Referring to Appendices P.1, P.2, and P.3, write a class for s-function that implements the equations of motion. Simulate the system using a variable force on the cart as an input. The output should connect to the animation function developed in homework B.2.

Design Study B guided exercise

Assemble the pendulum equations of motion

Use energy, generalized coordinates, and generalized forces to form the coupled nonlinear equations.

Use this setup to connect gravity, force, cart motion, and pendulum angle before applying Euler-Lagrange.

Step 1. Potential energy

Potential energy term

Check Step 1

Step 2. Generalized coordinates

Coordinates Generalized force vector

Check Step 2

Step 3. Coupled dynamics structure

Mass matrix size Damping direction

Check Step 3

Pseudocode

The complete solution includes a dynamics class and a simulation script. Before opening it, outline each component separately.

These prompts are ungraded. They are intended to connect the derivation to the implementation code in the complete solution.

Show complete solution

3.4 Design Study C. Satellite Attitude Control

3.4.1 Example Problem 3.3.6

Find the potential energy for the system.

Define the generalized coordinates.

Find the generalized forces and damping forces.

Derive the equations of motion using the Euler-Lagrange equations.

Referring to Appendices P.1, P.2, and P.3, write a class or s-function that implements the equations of motion. Simulate the system using a variable torque on the body as an input. The output should connect to the animation function developed in homework C.2.

Design Study C guided exercise

Assemble the satellite equations of motion

Use the torsional spring, relative damping, body torque, and two angular coordinates to form the coupled satellite-panel equations.

The spring and damper act on relative motion between the body and panel, while the actuator torque is applied only to the body.

Step 1. Potential energy

Spring deflection Energy coefficient

Check Step 1

Step 2. Generalized coordinates and input

Coordinates Generalized force vector

Check Step 2

Step 3. Relative damping

Body damping uses Panel damping uses

Check Step 3

Step 4. Coupled equation structure

Mass matrix size Direct torque equation

Check Step 4

Pseudocode

The complete solution includes dynamics code. Before opening it, outline how the code should compute \(\ddot{\theta}\) and \(\ddot{\phi}\) from state and torque.

This prompt is ungraded. It connects the Euler-Lagrange equations to the implementation code in the complete solution.

Show complete solution

3.5 Interactive Euler-Lagrange Check

3.6 Knowledge Check

3.7 Notes and References

Derivations of the Euler-Lagrange equations can be found in many sources, including [1, 2]. A simple derivation using only calculus can be found in [3].

The equations of motion that we derived in the previous section are often too complicated to facilitate effective engineering design. Accordingly, the standard approach is to simplify the models into so-called design models that capture the essential features of the system. The design techniques studied in this class will require that the design models are linear and time-invariant. In addition, the PID method described in Part III will require that the design models are second-order systems.

Our approach to developing linear time-invariant design models will be to linearize the simulation models developed in Part I and then to convert the linearized models into two standard forms, namely

• Transfer function models and

• State-space models.

Both forms will have advantages and disadvantages that will be explained throughout the remaining parts of the book. In Chapter 7 we will define the important concept of equilibrium and show how the equations of motion can be linearized about the equilibrium. Chapter 8 will show how to transform the linearized equations of motion into transfer function models, and Chapter 9 will show how to put the equations of motion in state-space form.

The remaining parts of the book will use the design models to develop feedback controllers for the system of interest. In Part III we will describe the
Proportional-Integral-Derivative (PID) controllers based on second-order transfer function models. Part IV derives observer-based controllers using state-space models. Finally, Part V describes loopshaping design strategies using general transfer function models.

4 Equilibria and Linearization

4.1 Theory

In this chapter we will describe two methods for linearizing the system, namely Jacobian linearization and feedback linearization. Jacobian linearization is based on the idea of approximating the nonlinear terms in the equations of motion by the first two terms in their Taylor series expansion. Feedback linearization is based on the idea of using the input signal to artificially remove the nonlinear terms.

4.1.1 Jacobian Linearization

If the function \(g:\mathbb{R}\to\mathbb{R}\) is a nonlinear function that can be differentiated an infinite number of times, then its Taylor series expansion about the point \(w_0\) is given by \[g(w) = g(w_0) + \left.\frac{\partial g}{\partial w}\right|_{w_0} (w-w_0) + \frac{1}{2!}\left.\frac{\partial^2 g}{\partial w^2}\right|_{w_0}(w-w_0)^2 + \cdots.\] A pictorial representation of Taylor series expansion is shown in Fig. 4.1, where it can be seen that in a small neighborhood around \(w=w_0\), the function \(g(w)\) is well approximated by the first two terms in the Taylor series expansion, namely \[g(w)\approx g(w_0) + \left.\frac{\partial g}{\partial w}\right|_{w_0} (w-w_0).\]

If \(g\) is a scalar function of several variables, then the Taylor series about \(w_1=w_{10}, w_2=w_{20}, \cdots, w_m=w_{m0}\) is given by \[\begin{gathered} g(w_1,w_2,\dots,w_m) = g(w_{10}, w_{20}, \dots, w_{m0}) + \left.\frac{\partial g}{\partial w_1}\right|_{w_{10},w_{20},\dots,w_{m0}} (w_1-w_{10}) \\ + \left.\frac{\partial g}{\partial w_2}\right|_{w_{10},w_{20},\dots,w_{m0}} (w_2-w_{20}) + \cdots + \left.\frac{\partial g}{\partial w_m}\right|_{w_{10},w_{20},\dots,w_{m0}} (w_m-w_{m0}) + \cdots. \end{gathered}\] For example, given the function \(g(v,w)=v\cos w\), the constant and linear terms in the Taylor series, expanded around \(v_0\) and \(w_0\) are \[\begin{aligned} g(v,w) &\approx g(v_0, w_0) + \left.\frac{\partial g}{\partial v}\right|_{v_{0},w_{0}}(v-v_{0}) + \left.\frac{\partial g}{\partial w}\right|_{v_{0},w_{0}}(w-w_{0}) \\ &= v_0\cos w_0 + \left.\left(\cos w \right)\right|_{v_0,w_0}(v-v_0) % + \left.\left(-v\sin w\right)\right|_{v_0,w_0}(w-w_0) \\ &= v_0\cos w_0 + (\cos w_0)(v-v_0) - (v_0\sin w_0)(w-w_0). \end{aligned}\]

Since approximating a nonlinear function by the first two terms in its Taylor series is only a good approximation in a small neighborhood of \(w_0\), it is important to select \(w_0\) judiciously. In particular, we will be interested in linearizing about equilibrium points of differential equations.

Given the differential equation \(\dot{x} = f(x,u)\), the pair \((x_e,u_e)\) is an equilibrium point if \[f(x_e,u_e)=0.\] In other words, an equilibrium point is a state-input pair where there is not any motion in the system.

For example, suppose that the nonlinear differential equation of a second-order system is given by \[\tag{4.2} \ddot{y} + a\dot{y} + by + g(y,\dot{y}) = u,\] where \(a\) and \(b\) are constants and \(g(y,\dot{y})\) is a known nonlinear function of \(y\) and \(\dot{y}\), and \(u\) is the input signal. Defining \(x=(y,\dot{y})^\top\), we get \[\dot{x} = \begin{pmatrix} \dot{y} \\ \ddot{y} \end{pmatrix} = \begin{pmatrix} \dot{y} \\ u - a\dot{y} - by - g(y,\dot{y}) \end{pmatrix} \stackrel{\triangle}{=}f(x,u).\] The equilibrium is when \(f(x_e,u_e)=0\) or in other words, when \[y_e = \text{~anything~}, \qquad \dot{y}_e=0, \qquad u_e=by_e+g(y_e,0).\] Equivalently, at equilibrium there is no motion in the system, which implies that \(\ddot{y}_e=\dot{y}_e=0\), which from Equation 4.2 implies that \(u_e = by_e+g(y_e,0)\).

Jacobian linearization then proceeds by replacing each term in the nonlinear differential equations describing the system by the first two terms in the Taylor series expansion about the equilibrium point.

For example, since the first term in Equation 4.2 is already linear, it can be expanded exactly as \[\ddot{y} = \ddot{y}_e + \frac{\partial \ddot{y}}{\partial \ddot{y}}\Big|_e (\ddot{y}-\ddot{y}_e) = \ddot{\tilde{y}},\] where we have defined the linearized quantity \(\tilde{y}\stackrel{\triangle}{=}y-y_e\), and we have used the fact that \(\ddot{y}_e=0\). Similarly, the remaining terms in Equation 4.2 are linearized as \[\begin{aligned} a\dot{y} &= a\dot{y}_e + a\frac{\partial \dot{y}}{\partial \dot{y}}\Big|_e (\dot{y}-\dot{y}_e) = a\dot{\tilde{y}} \\ by &= by_e + b\frac{\partial y}{\partial y}\Big|_e (y-y_e) = by_e + b\tilde{y} \\ g(y,\dot{y}) &\approx g(y_e,\dot{y}_e) + \frac{\partial g}{\partial y}\Big|_e(y-y_e) + \frac{\partial g}{\partial \dot{y}}\Big|_e(\dot{y}-\dot{y}_e) \\ &= g(y_e,0) + \frac{\partial g}{\partial y}\Big|_e\tilde{y} + \frac{\partial g}{\partial \dot{y}}\Big|_e\dot{\tilde{y}} \\ u &= u_e + \frac{\partial u}{\partial u}\Big|_e (u-u_e) = u_e + \tilde{u}. \end{aligned}\] Substituting these expressions into Equation 4.2 the resulting linearized equations of motion are \[[\ddot{\tilde{y}}] + [a\dot{\tilde{y}}] + [by_e + b\tilde{y}] + \left[g(y_e,0) + \frac{\partial g}{\partial y}\Big|_e\tilde{y} + \frac{\partial g}{\partial \dot{y}}\Big|_e\dot{\tilde{y}}\right] = [u_e + \tilde{u}]\] Using the equilibrium \(u_e=by_e+g(y_e,0)\) gives \[\ddot{\tilde{y}} + \left(a+\frac{\partial g}{\partial \dot{y}}\Big|_e\right)\dot{\tilde{y}} + \left(b+\frac{\partial g}{\partial y}\Big|_e\right)\tilde{y} = \tilde{u}.\] A block diagram that shows the design model for Equation 4.2 using Jacobian linearization is shown in Fig. 4.2.

4.1.2 Feedback Linearization

It is often the case that the nonlinearities in the system can be directly canceled by a judicious choice for the input variable. Consider again the nonlinear equations of motion in Equation 4.2. If \(u\) is selected as \[u = g(y,\dot{y}) + \tilde{u}\] then the resulting equation of motion is \[\ddot{y} + a\dot{y} + by = \tilde{u},\] which is linear. This technique is called feedback linearization because the feedback signal \(u\) has been used to linearize the system by directly canceling the nonlinearities. The advantage of this method is that the equations of motion are still globally defined and are not restricted to a small neighborhood of an equilibrium point. The disadvantage is that \(g(y,\dot{y})\) may not be precisely known, and therefore may not be completely canceled by \(u\). A block diagram that shows the design model for Equation 4.2 using feedback linearization is shown in Fig. 4.3.

4.2 Design Study A. Single Link Robot Arm

4.2.1 Example Problem 3.1.7

For the single-link robot arm:

(a) Find the equilibria of the system.

(b) Linearize the system about the equilibria using Jacobian linearization.

(c) Linearize the system using feedback linearization.

4.3 Design Study B. Inverted Pendulum

4.3.1 Example Problem 3.2.7

For the inverted pendulum,

(a) Find the equilibria of the system.

(b) Linearize the system about the equilibria using Jacobian linearization.

Note that for the inverted pendulum, since the nonlinearities are not all contained in the same channel as the control force \(F\), the system cannot be feedback linearized.

4.4 Design Study C. Satellite Attitude Control

4.4.1 Example Problem 3.3.7

For the satellite system:

(a) Find the equilibria of the system.

4.5 Interactive Linearization Check

4.6 Knowledge Check

4.7 Notes and References

Jacobian linearization is a standard technique that is well documented in most introductory textbooks on control. For example, see [4–7]. Feedback linearization is a much deeper topic than can be covered in an introductory textbook on control. A discussion that is similar to ours, that assumes that the nonlinearity is in the input channel, is contained in [5]. Feedback linearization is a fairly standard technique for fully actuated mechanical systems modeled by Euler-Lagrange equations, where the model of the system can be reduced to \[M(q)\ddot{q} + c(q,\dot{q})\dot{q}+g(q) = \tau,\] where \(q\in\mathbb{R}^m\) is the generalized coordinate and \(\tau\in\mathbb{R}^m\) is the generalized force, and where \(M(q)\) is an invertible matrix. In that case we can write \[\ddot{q} = M^{-1}(q)\left(\tau - c(q,\dot{q})\dot{q} - g(q)\right).\] Selecting \(\tau\) as \[\tag{4.16} \tau = c(q,\dot{q})\dot{q} + g(q) + M(q)\nu,\] results in the linear system \[\ddot{q} = \nu,\] where the new control variable \(\nu\) can be designed to stabilize the system. The control law in Equation 4.16 is the feedback linearizing control. This is a standard technique for robot manipulators and is often called the method of computed torque [8–10].

When the nonlinearity is not in the input channel, the system can still often be feedback linearized through a suitable change of variables. A good introduction to the topic is contained in [11]. More advanced coverage of the topic can be found in [12–15].

5 Transfer Function Models

5.1 Theory

The (one-sided) Laplace transform of a time domain signal \(y(t)\) is defined as \[\mathcal{L}\{y(t)\} \stackrel{\triangle}{=}Y(s) = \int_0^\infty y(t) e^{-st}\,dt.\] In essence, the Laplace transform indicates the content of the complex signal \(e^{-st}\) that is contained in \(y(t)\). The Laplace transform of the time derivative of \(y(t)\) is given by \[\mathcal{L}\{\dot{y}\} = \int_0^\infty \dot{y} e^{-st}\,dt.\] Using integration by parts \(\int udv = uv-\int vdu\), where \(u=e^{-st}\) and \(dv = \dot{y}dt\), and assuming that \(\lim_{t\to\infty}y(t)e^{-st}=0\), or that \(s\) is in the region of convergence of the Laplace transform, we get \[\begin{aligned} \mathcal{L}\{\dot{y}\} &= \int_0^\infty \dot{y} e^{-st}\,dt \\ &= y(t)e^{-st}\Big|_0^\infty - \int_0^\infty y(t)\left(-se^{-st}dt\right) \\ &= \lim_{t\to\infty}y(t)e^{-st} - y(0) + s\int_0^\infty y(t)e^{-st}\,dt \\ &= sY(s) - y(0). \end{aligned}\] Similarly, the Laplace transform of the second derivative of \(y(t)\) is given by \[\begin{aligned} \mathcal{L}\{\ddot{y}\} &= s\mathcal{L}\{\dot{y}\} - \dot{y}(0) \\ &= s\left(sY(s)-y(0)\right) - \dot{y}(0) \\ &= s^2 Y(s) - sy(0) - \dot{y}(0). \end{aligned}\]

The transfer function is typically defined as the relationship between the input and the output of the system. In that setting the transfer function is found when all initial conditions are set to zero. This is important to remember in the context of modeling. We will see that state-space models explicitly account for initial conditions, but that transfer functions do not.

Given the \(n^{th}\) order differential equation \[\begin{gathered} y^{(n)} + a_{n-1}y^{(n-1)} + \cdots+ a_1 \dot{y} + a_0 y \\ = b_mu^{(m)} + b_{m-1}u^{(m-1)}+\cdots+b_1\dot{u} + b_0 u, \end{gathered}\] the transfer function from \(u\) to \(y\) is found by taking the Laplace transform of both sides of the equation and setting all initial conditions to zero to obtain \[\begin{gathered} s^nY(s) + a_{n-1}s^{n-1}Y(s) + \cdots+ a_1 sY(s) + a_0 Y(s) \\ = b_ms^{m}U(s) + b_{m-1}s^{m-1}U(s)+\cdots+b_1sU(s) + b_0U(s). \end{gathered}\] Factoring \(Y(s)\) and \(U(s)\) we get \[\begin{gathered} \left(s^n + a_{n-1}s^{n-1} + \cdots+ a_1 s + a_0\right) Y(s) \\ = \left(b_ms^{m}+b_{m-1}s^{m-1}+\cdots+b_1s+b_0\right)U(s). \end{gathered}\] Solving for \(Y(s)\) gives \[Y(s) = \left(\frac{b_ms^{m} + b_{m-1}s^{m-1}+\cdots+b_1s + b_0}{s^n + a_{n-1}s^{n-1} + \cdots+ a_1 s + a_0}\right)U(s),\] where \[H(s)= \frac{b_ms^{m} + b_{m-1}s^{m-1}+\cdots+b_1s + b_0}{s^n + a_{n-1}s^{n-1} + \cdots+ a_1 s + a_0}\] is the transfer function from \(u\) to \(y\).

Given the transfer function \(H(s)\), the poles of the transfer function are the roots of the denominator polynomial \[s^n + a_{n-1}s^{n-1} + \cdots+ a_1 s + a_0 = 0.\] Similarly, the zeros of the transfer function are the roots of the numerator polynomial \[b_ms^{m} + b_{m-1}s^{m-1}+\cdots+b_1s + b_0 = 0.\]

As a concrete example, consider the differential equation given by \[y^{(3)} + 2\ddot{y} + 3\dot{y} + 4y = 5\dot{u} + 6u.\] Taking the Laplace transform of both sides, setting the initial conditions to zero, and factoring \(Y(s)\) and \(U(s)\) gives \[(s^3 + 2s^2 + 3s + 4)Y(s) = (5s+6)U(s).\] Solving for \(Y(s)\) gives \[Y(s) = \frac{5s+6}{s^3+2s^2+3s+4} U(s).\] The transfer function from \(u\) to \(y\) is therefore \[H(s) = \frac{5s+6}{s^3+2s^2+3s+4}.\] The poles of the transfer function are the roots of \(s^3+2s^2+3s+4=0\) or \(p_1=-1.6506\), \(p_2=-0.1747+j1.5469\), and \(p_3=-0.1747-j1.5469\). The zeros of the transfer function are the roots of \(5s+6=0\) or \(z_1 = -6/5\).

5.2 Design Study A. Single Link Robot Arm

5.2.1 Example Problem 3.1.8

For the single-link robot arm, find the transfer function of the system from the torque \(\tau\) to the angle \(\theta\).

Design Study A guided exercise

Derive the arm transfer function

Start from the feedback-linearized arm equation and solve for \(\Theta(s)/\tilde{\tau}(s)\).

The completed transfer function should fit this input-output block diagram from \(\tilde{\tau}\) to \(\theta\).

Step 1. Starting differential equation

Enter the coefficient on \(\ddot{\theta}\).

Inertia coefficient

Check Step 1

Step 2. Transfer-function coefficients

Numerator constant Denominator \(s\) coefficient Denominator constant

Check Step 2

Step 3. Pole locations

Integrator pole Damping pole

Check Step 3

Show complete solution

5.3 SIMO Systems and Cascade Approximations

There are many systems that have multiple inputs and/or multiple outputs. In this case, a transfer function exists from each input to each output. For example, suppose that the system is given by the coupled differential equations \[\begin{aligned} \ddot{y}_1 + a_{11}\dot{y}_1 + \ddot{y}_2 + a_{12}\dot{y}_2 &= b_{11}u_1 + b_{12}u_2 + b_{13}\dot{u}_2 \\ \ddot{y}_2 + a_{22}y_2 &= b_{21}u_2. \end{aligned}\] Taking the Laplace transform and arranging in matrix notation gives \[\left(\begin{array}{c|c} s^2+a_{11}s & s^2+a_{12}s \\ \hline 0 & s^2+a_{22}\end{array}\right) \begin{pmatrix}Y_1(s) \\ Y_2(s)\end{pmatrix} \\ = \left(\begin{array}{c|c} b_{11} & (b_{13}s+b_{12}) \\ \hline0 & b_{21}\end{array}\right) \begin{pmatrix} U_1(s) \\ U_2(s)\end{pmatrix}.\] To find the relevant transfer functions, we need to isolate the outputs on the left-hand side by inverting the matrix on the left. Using the inverse formula for a \(2\times 2\) matrix derived in Appendix 31 \[\begin{pmatrix} a & b \\ c & d \end{pmatrix}^{-1} = \frac{\begin{pmatrix} d & -b \\ -c & a\end{pmatrix}}{ad-bc}\] we get \[\begin{aligned} \begin{pmatrix}Y_1 \\ Y_2\end{pmatrix} &= \left(\begin{array}{c|c} s^2+a_{11}s & s^2+a_{12}s \\ \hline 0 & s^2+a_{22}\end{array}\right)^{-1} \left(\begin{array}{c|c} b_{11} & b_{13}s+b_{12} \\ \hline 0 & b_{21}\end{array}\right) \begin{pmatrix} U_1 \\ U_2\end{pmatrix} \\ &= \frac{\left(\begin{array}{c|c} s^2+a_{22} & -s^2-a_{12}s \\ \hline 0 & s^2+a_{11}s\end{array}\right)}{(s^2+a_{11}s)(s^2+a_{22})} \left(\begin{array}{c|c} b_{11} & b_{13}s+b_{12} \\ \hline 0 & b_{21}\end{array}\right) \begin{pmatrix} U_1 \\ U_2\end{pmatrix} \\ &= \frac{\left(\begin{array}{c|c} b_{11}(s^2+a_{22}) & (s^2+a_{22})(b_{13}s+b_{12})-b_{21}(s^2+a_{12}s) \\ \hline 0 & b_{21}(s^2+a_{11}s)\end{array}\right)}{(s^2+a_{11}s)(s^2+a_{22})} \begin{pmatrix} U_1 \\ U_2\end{pmatrix} \\ &= \left(\begin{array}{c|c} \frac{b_{11}}{s^2+a_{11}s} & \frac{b_{13}s^3+(b_{12}-b_{21})s^2+(b_{13}a_{22}-b_{21}a_{12})s+a_{22}b_{12}}{(s^2+a_{11}s)(s^2+a_{22})} \\ \hline 0 & \frac{b_{21}}{s^2+a_{22}}\end{array}\right) \begin{pmatrix} U_1 \\ U_2\end{pmatrix}. \end{aligned}\] Therefore, the transfer function from \(u_1\) to \(y_1\) is \(\frac{b_{11}}{s^2+a_{11}s}\), the transfer function from \(u_2\) to \(y_1\) is \(\frac{b_{13}s^3+(b_{12}-b_{21})s^2+(b_{13}a_{22}-b_{21}a_{12})s+a_{22}b_{12}}{(s^2+a_{11}s)(s^2+a_{22})}\), the transfer function from \(u_1\) to \(y_2\) is \(0\), and the transfer function from \(u_2\) to \(y_2\) is \(\frac{b_{21}}{s^2+a_{22}}\).

For many of the case studies that we will look at in this book, there is a single input, but we are also interested in controlling two or more measurements or outputs. State-space control techniques are able to create control laws to control multiple outputs in a single design step, as we will see in later chapters. However, we also want to be able to create single-input, multiple-output controllers using more classical transfer-function based approaches. Doing so will require more intuition about the systems we are controlling, as well as a clear understanding of our control objectives so that we can make appropriate design decisions and simplifying assumptions.

The majority of our design studies (specifically studies B, C, E, and F) have coupled, higher-order dynamics (greater than order two) with multiple outputs. Instead of coming up with a single feedback loop to control the multiple outputs of interest, we will form a cascade of two lower-order transfer functions and design two independent control loops. Forming the cascade of lower-order transfer functions requires insights into the dynamic behavior of the system, because it may require simplifying assumptions about the dynamic coupling between the two subsystems. Instead of a two-way coupling, the cascade may make the dynamic coupling one-way. To reduce the complexity of the system, we may make simplifying assumptions about the dynamics of one of the subsystems. It is essential for these simplifications to be justifiable given the parameters of the system and our control objectives. For example, in design study C, we form a cascade of two transfer functions: the first is from the torque input to the satellite angular position and the second is from the satellite angular position to the solar panel angular position. Because the satellite is more massive than the solar panels, we will make a valid assumption that while the satellite significantly influences the dynamic behavior of the solar panels, the solar panels do not significantly influence the dynamic behavior of the satellite. As we will see, this will simplify the design of the control system significantly and result in a viable control design that provides good performance.

In forming cascades of transfer functions, other factors may influence the form of the cascade. Our control priorities may play a significant role. For example, in design study B, the balancing of the inverted pendulum is the highest priority and therefore, we will use the force input to control the angle of the pendulum. Of importance, but lower priority, is controlling the cart. These considerations will lead us to form the cascade using the transfer functions from force to pendulum angle and from pendulum angle to cart position. In other situations, the causality of the dynamic behavior is more transparent. In the block-beam system, the force on the beam causes the beam angle to change and the change in beam angle causes the block to slide. We will form the transfer function cascade for the block-beam system accordingly. In the planar VTOL aircraft system, torque from the differential thrust causes the aircraft to pitch, while thrust combined with non-zero pitch causes the aircraft to translate. It is natural, therefore, to think of torque causing pitch, and pitch causing translational motion. We will thus form a cascade of two transfer functions: from torque to pitch angle and from pitch angle to horizontal position.

In forming a cascade of two (or more) transfer functions to represent the dynamics of a single system, it is critically important that we consider the physics of the system and our control objectives as we make simplifying assumptions. Important questions to consider include:

• What is the control input to the system?

• What outputs of the system are to be controlled?

• What outputs can be measured?

• What physically justifiable simplifications can be made about the dynamic coupling between subsystems?

• What physically justifiable simplifying assumptions can be made about the dynamics of the system?

Asking and answering these questions is fundamental to the design of successful control systems. The more insight you have into the physical behavior of the system you are controlling, the better your control system designs will be.

5.4 Design Study B. Inverted Pendulum

5.4.1 Example Problem 3.2.8

For the inverted pendulum:

(a) Start with the linearized equations and use the Laplace transform to convert the equations of motion to the s-domain.

(b) Find the transfer functions from the input \(\tilde{F}(s)\) to the outputs \(\tilde{Z}(s)\) and \(\tilde{\Theta}(s)\). Assume that the damping constant is negligible. How does this simplify the transfer functions? Is this a reasonable assumption?

(c) From the simplified \(\tilde{Z}(s)/\tilde{F}(s)\) and \(\tilde{\Theta}(s)/\tilde{F}(s)\) transfer functions, compute the \(\tilde{Z}(s)/\tilde{\Theta}(s)\) transfer function. Draw a block diagram of the system as a series of two transfer functions from \(\tilde{F}(s)\) to \(\tilde{\Theta}(s)\) and from \(\tilde{\Theta}(s)\) to \(\tilde{Z}(s)\).

(d) Describe how this transfer-function cascade makes sense physically. How does force influence the pendulum angle? How does the pendulum angle influence the cart position?

Design Study B guided exercise

Build the pendulum transfer-function cascade

Start from the linearized equations and identify the force-to-angle and angle-to-position relationship.

Use the cascade structure to separate the fast force-to-angle subsystem from the angle-to-position subsystem.

Step 1. Simplifying assumption

Damping assumption

Check Step 1

Step 2. Force-to-angle subsystem

Plant behavior

Check Step 2

Step 3. Angle-to-position subsystem

Free integrator count Cascade output

Check Step 3

Show complete solution

The PID and loopshaping techniques that we will discuss in Parts III and V will use the approximate transfer function models of the plant derived above. The state-space methods discussed in Part IV will use a state-space model, which does not depend on neglecting the coupling between subsystems and does not depend on a separation into fast and slow subsystems.

5.5 Design Study C. Satellite Attitude Control

5.5.1 Example Problem 3.3.8

For the satellite attitude control problem:

(a) Start with the linearized equations for the satellite attitude problem and use the Laplace transform to convert the equations of motion to the s-domain.

(b) Find the full transfer matrix from the input \(\tau(s)\) to the outputs \(\Phi(s)\) and \(\Theta(s)\).

(c) From the transfer matrix, find the second-order transfer function from \(\Theta(s)\) to \(\Phi(s)\).

(d) Under the assumption that the panel moment of inertia \(J_p\) is significantly smaller than the spacecraft moment of inertia \(J_s\) (specifically, \((J_s+J_p)/J_s \approx 1\)), find the second-order approximation for the transfer function from \(\tau(s)\) to \(\Theta(s)\).

(e) From your results on parts (c) and (d), form the approximate transfer function cascade for the satellite/panel system and justify why it makes sense physically.

Design Study C guided exercise

Build the satellite transfer-function cascade

The transfer-function model separates the body-angle dynamics from the panel-angle dynamics using a physically motivated cascade approximation.

The approximation lets the body angle drive the panel dynamics while neglecting the panel's back-effect on the massive satellite body.

Step 1. Identify the s-domain outputs

Body-angle transform Panel-angle transform

Check Step 1

Step 2. Body angle to panel angle

Transfer-function ratio Model order

Check Step 2

Step 3. Torque to body angle approximation

Approximate inertia Free integrator count

Check Step 3

Step 4. Interpret the approximation

Neglected direction

Check Step 4

Show complete solution

5.6 Interactive Transfer Function Check

5.7 Knowledge Check

5.8 Notes and References

The Laplace transform can be found in introductory textbooks on signals and systems like [16] and in introductory textbooks on differential equations. The transfer function of a linear time-invariant system is also covered in introductory textbooks on signals and systems, and in most introductory textbooks on control like [4].

6 State-Space Models

6.1 Theory

The transfer function models discussed in the previous chapter are a frequency-domain representation of the system because they describe how the system responds to sinusoidal inputs at specific frequencies. In contrast, differential equation models are time-domain representations of the systems that directly model the time-evolution of the system. The state-space model is also a time-domain representation that is essentially a reformatting of the differential equations. In essence, the state-space model represents the system by an input, an output, and a memory element, called the state. In this book, the input will always be denoted as \(u(t)\), the output as \(y(t)\), and the state as \(x(t)\). The state represents all of the memory elements in the system. For example, the state may represent a storage register, the altitude of an aircraft, the velocity of a car, the voltage across a capacitor, or the current through an inductor.

The state-space model is composed of two equations. The first equation is called the state evolution equation and represents how the memory elements, or states, change as a function of the current state and the inputs to the system. The second equation is called the output equation and describes how the current output of the system depends on the current state and the current input. The general state-space equations for a continuous time system are written as \begin{align} \dot{x} &= f(x,u) \tag{6.1} \\ y &= h(x,u), \tag{6.2} \end{align} where Equation 6.1 is the state evolution equation and Equation 6.2 is the output equation. Note that for continuous time systems the state evolution equation is a system of coupled first order nonlinear differential equations. The general state-space equations for a discrete time system are written as \begin{align} x[k+1] &= f(x[k],u[k]) \tag{6.3} \\ y[k] &= h(x[k],u[k]), \tag{6.4} \end{align} where again Equation 6.3 is the state evolution equation, this time given by a system of coupled first order difference equations, and Equation 6.4 is the output equation.

Any nonlinear \(n^{th}\) order differential equation can be transformed into state-space form. For example, consider the nonlinear differential equation \[\dddot{y} + y^2\ddot{y} + \sin(y\dot{y}) + e^y = u,\] where \(u\) is the input and \(y\) is the output. Define the state as \[x=\begin{pmatrix}x_1 \\ x_2 \\ x_3 \end{pmatrix}\stackrel{\triangle}{=}\begin{pmatrix} y \\ \dot{y} \\ \ddot{y} \end{pmatrix}.\] Then the state evolution equation can be derived as \[\begin{aligned} \dot{x} &= \begin{pmatrix} \dot{x}_1 \\ \dot{x}_2 \\ \dot{x}_3 \end{pmatrix} = \begin{pmatrix} \dot{y} \\ \ddot{y} \\ \dddot{y} \end{pmatrix} = \begin{pmatrix} \dot{y} \\ \ddot{y} \\ -y^2\ddot{y} - \sin(y\dot{y}) - e^y + u \end{pmatrix} \\ &= \begin{pmatrix} x_2 \\ x_3 \\ -x_1^2x_3 - \sin(x_1x_2) - e^{x_1} + u \end{pmatrix} \stackrel{\triangle}{=}f(x,u), \end{aligned}\] and the output equation is \[y = x_1 \stackrel{\triangle}{=}h(x,u).\] As another example, consider the coupled differential equations \[\begin{aligned} \ddot{y}_1 + \dot{y}_1y_2 + y_2 &=u \\ \ddot{y}_2 + \dot{y}_2\cos(y_1)+\dot{y_1} &=0. \end{aligned}\] Defining the output as \(y=(y_1, y_2)^{\top}\) and the state as \[x=\begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix} \stackrel{\triangle}{=}\begin{pmatrix} y_1 \\ y_2 \\ \dot{y}_1 \\ \dot{y}_2 \end{pmatrix},\] the state evolution equation can be derived as \[\begin{aligned} \dot{x} &= \begin{pmatrix} \dot{x}_1 \\ \dot{x}_2 \\ \dot{x}_3 \\ \dot{x}_4 \end{pmatrix} = \begin{pmatrix} \dot{y}_1 \\ \dot{y}_2 \\ \ddot{y}_1 \\ \ddot{y}_2 \end{pmatrix} = \begin{pmatrix} \dot{y}_1 \\ \dot{y}_2 \\ -\dot{y}_1y_2 - y_2 +u \\ -\dot{y}_2\cos(y_1)-\dot{y_1} \end{pmatrix} \\ &= \begin{pmatrix} x_3 \\ x_4 \\ -x_3x_2 - x_2 +u \\ -x_4\cos(x_1)-x_3 \end{pmatrix} \stackrel{\triangle}{=}f(x,u), \end{aligned}\] and the output equation is given by \[y = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \stackrel{\triangle}{=}h(x,u).\]

6.1.1 Jacobian Linearization Revisited

Jacobian linearization of the system is more straightforward when the equations are in state-space form. Given the nonlinear state-space equations \begin{align} \dot{x} &= f(x,u) \tag{6.5} \\ y &= h(x,u) \tag{6.6} \end{align} an equilibrium point is defined as any \((x_e,u_e)\) such that \(f(x_e,u_e)=0\). In other words, equilibrium points are combinations of states and inputs where the state stops evolving in time. To linearize about an equilibrium point, Equations 6.5 and 6.6 can be expanded in a Taylor series about \((x_e,u_e)\) as \begin{align} f(x,u) &\approx f(x_e,u_e) + \left.\frac{\partial f}{\partial x}\right|_e (x-x_e) + \left.\frac{\partial f}{\partial u}\right|_e(u-u_e) + H.O.T. \tag{6.7}\\ h(x,u) &\approx h(x_e,u_e) + \left.\frac{\partial h}{\partial x}\right|_e (x-x_e) + \left.\frac{\partial h}{\partial u}\right|_e(u-u_e) + H.O.T., \tag{6.8} \end{align} where the Jacobian matrices are defined as \[\begin{aligned} \frac{\partial f}{\partial x} &\stackrel{\triangle}{=}\begin{pmatrix} \frac{\partial f_1}{\partial x_1} & \dots & \frac{\partial f_1}{\partial x_n} \\ \vdots & & \vdots \\ \frac{\partial f_n}{\partial x_1} & \dots & \frac{\partial f_n}{\partial x_n} \end{pmatrix} \\ \frac{\partial f}{\partial u} &\stackrel{\triangle}{=}\begin{pmatrix} \frac{\partial f_1}{\partial u_1} & \dots & \frac{\partial f_1}{\partial u_m} \\ \vdots & & \vdots \\ \frac{\partial f_n}{\partial u_1} & \dots & \frac{\partial f_n}{\partial u_m} \end{pmatrix} \\ \frac{\partial h}{\partial x} &\stackrel{\triangle}{=}\begin{pmatrix} \frac{\partial h_1}{\partial x_1} & \dots & \frac{\partial h_1}{\partial x_n} \\ \vdots & & \vdots \\ \frac{\partial h_p}{\partial x_1} & \dots & \frac{\partial h_p}{\partial x_n} \end{pmatrix} \\ \frac{\partial h}{\partial u} &\stackrel{\triangle}{=}\begin{pmatrix} \frac{\partial h_1}{\partial u_1} & \dots & \frac{\partial h_1}{\partial u_m} \\ \vdots & & \vdots \\ \frac{\partial h_p}{\partial u_1} & \dots & \frac{\partial h_p}{\partial u_m} \end{pmatrix}. \end{aligned}\] Therefore, \(\frac{\partial f}{\partial x}\) is an \(n\times n\) matrix, \(\frac{\partial f}{\partial u}\) is an \(n\times m\) matrix, \(\frac{\partial h}{\partial x}\) is a \(p\times n\) matrix, and \(\frac{\partial h}{\partial u}\) is a \(p\times m\) matrix.

Defining \(\tilde{x}\stackrel{\triangle}{=}x-x_e\) and \(\tilde{u}=u-u_e\) and letting \(A=\left.\frac{\partial f}{\partial x}\right|_e\) and \(B=\left.\frac{\partial f}{\partial u}\right|_e\), and noting that \(\dot{\tilde{x}}=\dot{x}-\dot{x}_e=\dot{x}\), results in the linearized state evolution equation \[\dot{\tilde{x}}=A\tilde{x} + B\tilde{u}.\] At the equilibria, the output may not necessarily be zero. If we define the equilibrium output to be \(y_e=h(x_e,u_e)\) and the linearized output as \(\tilde{y}=y-y_e\), and define \(C\stackrel{\triangle}{=}\left.\frac{\partial h}{\partial x}\right|_e\) and \(D\stackrel{\triangle}{=}\left.\frac{\partial h}{\partial u}\right|_e\), then using Equation 6.8 in 6.6 gives \[\begin{aligned} \tilde{y} &= y-y_e \\ &\approx h(x_e,u_e)+C\tilde{x} + D\tilde{u} - h(x_e,u_e) \\ &= C\tilde{x} + D\tilde{u}. \end{aligned}\] The linearized state-space equations are therefore \[\begin{aligned} \dot{\tilde{x}} &= A\tilde{x} + B\tilde{u} \\ \tilde{y} &= C\tilde{x} + D\tilde{u}. \end{aligned}\]

6.1.2 Converting State-Space Models to Transfer Function Models

Linear state-space models can be converted to transfer function models using the following technique. Suppose that the linear state-space model is given by \begin{align} \dot{x} &= Ax + Bu \tag{6.9} \\ y &= Cx + Du. \tag{6.10} \end{align} Defining the Laplace transform of a vector to be the Laplace transform of each element, i.e., \[\mathcal{L}\left\{\begin{pmatrix} x_1 \\ \vdots \\ x_n \end{pmatrix}\right\} = \begin{pmatrix} \mathcal{L}\left\{x_1\right\} \\ \vdots \\ \mathcal{L}\left\{x_n\right\}\end{pmatrix},\] and noting that the Laplace transform is a linear operator, which implies that \[\mathcal{L}\{Ax\} = A\mathcal{L}\{x\},\] taking the Laplace transform of Equations 6.9 and 6.10, and setting the initial conditions to zero, gives \begin{align} sX(s) &= AX(s) + BU(s) \tag{6.11} \\ Y(s) &= CX(s) + DU(s). \tag{6.12} \end{align} Solving 6.11 for \(X(s)\) gives \begin{align} & sX(s)-AX(s) = BU(s) \notag \\ \implies & (sI-A)X(s) = BU(s) \notag \\ \implies & X(s) = (sI-A)^{-1}BU(s), \tag{6.13} \end{align} where \(I\) is the \(n\times n\) identity matrix. Substituting 6.13 into 6.12 gives \[\begin{aligned} & Y(s) = C(sI-A)^{-1}BU(s) + DU(s) \\ \implies & Y(s) = \left[ C(sI-A)^{-1}B + D\right]U(s). \end{aligned}\] The transfer function from \(u\) to \(y\) is therefore \[\tag{6.14} H(s) = C(sI-A)^{-1}B + D.\] As an example, suppose that the state-space equations are given by \[\begin{aligned} \dot{x} &= \begin{pmatrix} 0 & 1 \\ -4 & -3 \end{pmatrix} x + \begin{pmatrix} 0 \\ 2 \end{pmatrix} u \\ y &= \begin{pmatrix} 1 & 0 \end{pmatrix} x + (0)u, \end{aligned}\] then the transfer function from \(u\) to \(y\) is given by \[\begin{aligned} H(s) &= C(sI-A)^{-1}B + D \\ &= \begin{pmatrix} 1 & 0\end{pmatrix}\left(s\begin{pmatrix} 1 & 0 \\ 0 & 1\end{pmatrix} - \begin{pmatrix} 0 & 1 \\ -4 & -3 \end{pmatrix}\right)^{-1}\begin{pmatrix} 0 \\ 2\end{pmatrix} + 0 \\ &= \begin{pmatrix} 1 & 0\end{pmatrix}\begin{pmatrix} s & -1 \\ 4 & s+3\end{pmatrix}^{-1}\begin{pmatrix} 0 \\ 2\end{pmatrix}\\ &= \begin{pmatrix} 1 & 0\end{pmatrix}\frac{\begin{pmatrix} s+3 & 1 \\ -4 & s\end{pmatrix}}{s^2+3s+4}\begin{pmatrix} 0 \\ 2\end{pmatrix}\\ &= \frac{\begin{pmatrix} s+3 & 1 \end{pmatrix}}{s^2+3s+4}\begin{pmatrix} 0 \\ 2\end{pmatrix}\\ &= \frac{2}{s^2+3s+4}. \end{aligned}\]

The form of the transfer function in Equation 6.14 allows us to derive a simple relationship between the poles of the transfer function \(H(s)\) and the eigenvalues of \(A\). Recall (see Appendix 31) that for any square invertible matrix \(M\), we have that \[M^{-1} = \frac{\text{adj}(M)}{\text{det}(M)},\] where \(\text{adj}(M)\) is the adjugate of \(M\) defined as the transpose of the cofactors of \(M\), and \(\text{det}(M)\) is the determinant of \(M\). Therefore \[\begin{aligned} H(s) &= C(sI-A)^{-1}B+D \\ &=\frac{C\text{adj}(sI-A)B}{\text{det}(sI-A)}+D \\ &= \frac{C\text{adj}(sI-A)B+D\text{det}(sI-A)}{\text{det}(sI-A)}. \end{aligned}\] This expression implies that the zeros of \(H(s)\) are given by the roots of the polynomial \[C\text{adj}(sI-A)B+D\text{det}(sI-A) = 0,\] and that the poles of \(H(s)\) are given by the roots of the polynomial \[\tag{6.15} \text{det}(sI-A) = 0.\] Recall from linear algebra (see Appendix 31) that Equation 6.15 also defines the eigenvalues of \(A\). Therefore, the poles of \(H(s)\) are equivalent to the eigenvalues of \(A\).

In Chapters 14 and 16 we will show two general techniques for converting SISO transfer function models into state-space models.

6.2 Design Study A. Single Link Robot Arm

6.2.1 Example Problem 3.1.9

For the feedback linearized equations of motion for the single-link robot arm given in Equation 4.7, define the states as \(x=(\theta,\dot{\theta})^{\top}\), the input as \(\tilde{u}=\tilde{\tau}\), and the measured output as \(y=\theta\). Find the linear state-space equations in the form \[\begin{aligned} \dot{x} &= Ax + B\tilde{u} \\ y &= Cx + D\tilde{u}. \end{aligned}\]

6.3 Design Study B. Inverted Pendulum

6.3.1 Example Problem 3.2.9

Defining the states as \(\tilde{x}=(\tilde{z},\tilde{\theta},\dot{\tilde{z}},\dot{\tilde{\theta}})^{\top}\), the input as \(\tilde{u}=\tilde{F}\), and the measured output as \(\tilde{y}=(\tilde{z},\tilde{\theta})^{\top}\), find the linear state-space equations in the form \[\begin{aligned} \dot{\tilde{x}} &= A\tilde{x} + B\tilde{u} \\ \tilde{y} &= C\tilde{x} + D\tilde{u}. \end{aligned}\]

6.4 Design Study C. Satellite Attitude Control

6.4.1 Example Problem 3.3.9

Suppose that a star tracker is used to measure \(\theta\) and a strain gauge is used to approximate \(\phi-\theta\). Defining the states as \(x=(\theta,\phi,\dot{\theta},\dot{\phi})^{\top}\), the input as \(u=\tau\), and the measured output as \(y=(\theta,\phi-\theta)^{\top}\), find the linear state-space equations in the form \[\begin{aligned} \dot{x} &= Ax + Bu \\ y &= Cx + Du. \end{aligned}\]

6.5 Interactive State-Space Check

6.6 Knowledge Check

6.7 Notes and References

There are numerous resources that derive state-space models from differential equations and transfer functions. For example [4–6, 17].

In this part of the book we will introduce feedback control by studying the most commonly used industrial control law, the so-called proportional-integral-derivative control, or PID control. PID controllers are widely used because they are intuitive and easy to understand, they can be designed without any understanding of the underlying physics of the system, and a model of the system is not required, which is especially important in the chemical processing industry. However, the stability and performance of the closed-loop system can only be guaranteed when using PID control on second-order systems. In Part IV on observer-based control and Part V on loopshaping, we will extend the basic ideas covered in this Part to higher-order systems.

Fig. 6.1 shows a basic proportional control strategy. The output \(y\) of the physical system is subtracted from the commanded reference output \(r\) to produce the error \(e=r-y\). Proportional control determines the input to the system \(u\) to be proportional to the error, i.e., \(u=k_P e\). The constant gain \(k_P\) is called the proportional gain.

Proportional control can be augmented with integral control to produce PI control, as shown in Fig. 6.2. The basic idea is to integrate the error so that the system responds to error accrued in the past. The constant gain \(k_I\) is called the integral gain.

Alternatively, proportional control can be augmented with a derivative control as shown in Fig. 6.3. The idea is to respond to how fast the error is changing, or the derivative of the error. The constant gain \(k_D\) is called the derivative gain.

When all three elements are used, as shown in Fig. 6.4, the resulting controller is called a PID control. In Part III of this book, we will explore various aspects of PID control. Many of these concepts can be generalized to other control schemes in later chapters, and so our study of PID control will also help to motivate and understand more advanced concepts.

Intuition for PID control can be gained by considering the step response shown in Fig. 6.5. At the beginning of the response, there is a significant error between the current output \(y(t)\) and the reference command \(r(t)\). The proportional control term will tend to push the system so as to reduce the gap between \(y(t)\) and \(r(t)\). As \(y(t)\) begins to change, a fast response may lead to significant overshoot due to the momentum in the system. Therefore, if it is moving too quickly, we may want to slow down the response, or if it is moving too slowly, we may want to speed up the response. Derivative control is used to affect these changes. After the response has settled into steady state, there may be significant steady-state error, as shown in Fig. 6.5. Integral control adds up this steady error and will eventually act to reduce the steady-state error. A good example is an aircraft in a constant crosswind. An integrator can be used to correct the response of the aircraft so that it crabs into the wind to maintain its desired flight path.

Part III of the book is organized as follows. In Chapter 10 we will show how PD control can be used to specify the desired closed-loop poles of a second-order system. In Chapter 11 we introduce the concepts of rise time, natural frequency, and damping ratio for second-order systems and show how these concepts relate to the location of the closed-loop poles. We also show how these quantities can be used to design PD controllers for second-order systems and cascades of second-order systems. In Chapter 12 we use the final value theorem to show the effect of integrators on the closed-loop system. We introduce the notion of system type in the context of reference tracking and disturbance rejection. Finally, in Chapter 13 we show how to write computer code to implement PID controllers.

7 Pole Placement for Second-Order Systems

7.1 Theory

Suppose that the model for a linearized physical system is given by the second-order transfer function \[\tag{7.1} P(s) = \frac{b_0}{s^2 + a_1s + a_0}.\] The open-loop poles of the plant \(P(s)\) are defined to be the roots of the open-loop characteristic polynomial \[\Delta_{ol}(s) = s^2 + a_1 s + a_0,\] which are given by \[p_{ol} = -\frac{a_1}{2} \pm \sqrt{\left(\frac{a_1}{2}\right)^2-a_0}.\] Note that the open-loop poles are determined by the physical parameters and that the poles may be in the right half-plane. If we use PD control to regulate the output \(y\) to the reference command \(r\), then the block diagram is shown in Fig. 7.1.

The transfer function for the closed-loop system can be derived by noting from the block diagram that \begin{align} & Y(s) = \left(\frac{b_0}{s^2 + a_1s + a_0}\right)\left(k_P + k_D s\right)\left(R(s)-Y(s)\right) \notag \\ \implies & (s^2 + a_1s + a_0)Y(s) = \left(b_0k_P + b_0k_D s\right)\left(R(s)-Y(s)\right) \notag \\ \implies & \left[(s^2 + a_1s + a_0)+\left(b_0k_P + b_0k_D s\right)\right]Y(s) \\ & = \left(b_0k_P + b_0k_D s\right)R(s) \notag \\ \implies & (s^2 + (a_1 + b_0k_D)s + (a_0+b_0k_P))Y(s) \\ & = \left(b_0k_P + b_0k_D s\right)R(s) \notag \\ \implies & Y(s) = \frac{b_0k_D s + b_0k_P}{s^2 + (a_1 + b_0k_D)s + (a_0+b_0k_P)}R(s). \tag{7.4} \end{align} Equation 7.4 constitutes the closed-loop transfer function for the system under the influence of PD control. The closed-loop poles are given by the roots of the closed-loop characteristic polynomial \[\Delta_{cl}(s) = s^2 + (a_1+b_0k_D) s + (a_0+b_0k_P),\] which are \[p_{cl} = -\frac{a_1 + b_0k_D}{2} \pm \sqrt{\left(\frac{(a_1+b_0k_D)}{2}\right)^2-(a_0+b_0k_P)}.\] Note that the closed-loop poles can be specified by the designer by selecting \(k_P\) and \(k_D\). The basic idea is to select desired closed-loop poles \(-p_1^d\) and \(-p_2^d\) and then to form the desired characteristic polynomial \[\Delta_{cl}^d(s) = (s+p_1^d)(s+p_2^d) \stackrel{\triangle}{=}s^2 + \alpha_1 s + \alpha_0.\] By setting \(\Delta_{cl}(s)=\Delta_{cl}^d(s)\) we can equate the leading coefficients of each term of the polynomial and solve for the gains \(k_P\) and \(k_D\) as \[\begin{aligned} k_P &= \frac{\alpha_0-a_0}{b_0} \\ k_D &= \frac{\alpha_1-a_1}{b_0}. \end{aligned}\]

One of the disadvantages to the PD architecture shown in Fig. 7.1 is the introduction of a zero in the closed-loop transfer function. Note that while the open-loop system 7.1 does not have a zero, the closed-loop system 7.4 has a zero at \[z_{cl} = - \frac{k_P}{k_D}.\] The presence of the zero will modify the closed-loop response of the system. A simple trick that removes the zero is to use the control structure shown in Fig. 7.2. In this particular case, the differentiator acts only on the output \(y\) and not on the error \(e\). Notice the change of sign on the signal coming from the differentiator.

In this case the closed-loop transfer function is calculated as \begin{align} & Y(s) = \left(\frac{b_0}{s^2 + a_1s + a_0}\right)\left(k_P (R(s)-Y(s)) - k_D s Y(s) \right) \notag \\ \implies & (s^2 + a_1s + a_0)Y(s) = \left(b_0k_P (R(s)-Y(s)) - b_0k_D s Y(s) \right) \notag \\ \implies & \left[(s^2 + a_1s + a_0)+\left(b_0k_P + b_0k_D s\right)\right]Y(s) = b_0k_P R(s) \notag \\ \implies & (s^2 + (a_1 + b_0k_D)s + (a_0+b_0k_P))Y(s) = b_0k_P R(s) \notag \\ \implies & Y(s) = \frac{b_0k_P}{s^2 + (a_1 + b_0k_D)s + (a_0+b_0k_P)}R(s). \tag{7.5} \end{align} Note that the closed-loop transfer function 7.5 has the same poles as Equation 7.4, but that the zero has been removed. An added benefit to the structure shown in Fig. 7.2 is that when \(r(t)\) is a step input, the derivative can introduce a large signal spike in \(u(t)\). These large spikes are removed by only differentiating \(y(t)\) instead of \(e(t)=r(t)-y(t)\).

7.2 Design Study A. Single Link Robot Arm

7.2.1 Example Problem 3.1.10

(a) Given the open-loop transfer function found in problem 3.1.6, find the open-loop poles of the system, when the equilibrium angle is \(\theta_e=0\).

(b) Using the PD control architecture shown in Fig. 7.2, find the closed-loop transfer function from \(r_{\theta}\) to \(\theta\) and find the closed-loop poles as a function of \(k_P\) and \(k_D\).

(c) Select \(k_P\) and \(k_D\) to place the closed-loop poles at \(p_1=-3\) and \(p_2=-4\).

(d) Using the gains from part (c), implement the PD control for the single-link robot arm in simulation and plot the step response for an amplitude of 30 degrees.

Design Study A guided exercise

Compute the PD gains

Match the closed-loop characteristic polynomial to desired poles at \(-3\) and \(-4\).

Use this PD feedback structure to derive the closed-loop characteristic polynomial.

Step 1. Desired characteristic polynomial

\(s\) coefficient Constant term

Check Step 1

Step 2. PD gains

\(k_P\) \(k_D\)

Check Step 2

Pseudocode

The complete solution includes a PD controller class and a simulation script. Before opening it, outline each component separately.

These prompts are ungraded. They are intended to connect the gain calculation to the implementation code in the complete solution.

Show complete solution

7.3 Interactive PD Pole Placement Check

7.4 Knowledge Check

7.5 Notes and References

The concept of using feedback to place the poles of the system at desired locations is a much deeper topic than is covered in this chapter. While the method discussed in this chapter is only applicable to second-order systems, it is usually possible to arbitrarily assign the poles for an \(n^{th}\) order system. In general, to assign \(n\) poles will require \(n\) gains. We will discuss pole placement in more detail in Chapter 14 of this book.

8 Design Strategies for Second-Order Systems

8.1 Theory

In the previous section we saw that for second-order systems, PD control could be used to exactly specify the pole location of the closed-loop system. This chapter examines the design problem of how to select the pole locations to achieve desired behavior. In Section 11.1.1 we derive the relationship between a single pole location and the step response of the system. In Section 11.1.2 we show that the time-domain response for a second-order system can be understood in terms of natural frequency and damping ratio. Finally, in Section 11.1.3 we show the effect that a zero has on the step response of a second-order system.

8.1.1 First Order Response

Consider a first-order system given by \[Y(s) = \frac{p}{s+p} U(s),\] where \(p>0\). The pole location of the first-order system is at \(-p\). If \(U(s)=\frac{A}{s}\) is a step of magnitude \(A\), then \[Y(s) = \frac{Ap}{s(s+p)}.\] As shown in Appendix 29, the inverse Laplace transform is given by \[\tag{8.1} y(t) = \begin{cases} A(1-e^{-pt}), & t\geq 0 \\ 0, & \text{otherwise} \end{cases}.\] The response is shown in Fig. 8.1, which shows the response to steps of size \(A=1\), \(A=2\), and \(A=3\), respectively.

If we define the rise time \(t_r\) of the system to be the time after the step when the response reaches 90% of its final value, then we note from Fig. 8.1 that the rise time is independent of the step size \(A\). Also, by differentiating Equation 8.1 at time \(t=0\), it is straightforward to show that the slope of the response at \(t=0\) is equal to \(Ap\). This is also shown graphically in Fig. 8.1.

When \(A=1\), the step response for different pole locations is shown in Fig. 8.2. Pole locations corresponding to \(p=1\), \(p=2\), and \(p=4\) are shown. Note that as the pole location moves further into the left half of the complex plane, the rise time decreases.

An important concept for the steady-state response of a system is the DC gain, which is defined as follows.

The DC gain of a transfer function \(H(s)\) is \[\left.H(s)\right|_{\text{DC gain}} = \lim_{s\to 0} H(s).\]

An important fact is that if the poles of \(H(s)\) are in the open left half-plane, then the response of \(H(s)\) to a step of size \(A\) approaches \(A\left.H(s)\right|_{\text{DC gain}}\) as \(t\to\infty\). This is true for first-order systems, as well as general \(H(s)\). For example, the DC gain of the first-order system \(H(s) = \frac{q}{s+p}\) is \(\frac{q}{p}\), as long as \(p>0\).

8.1.2 Second-Order Response

Suppose that the closed-loop system has a second-order transfer function with two real poles and no zeros: \[Y(s) = \frac{K}{(s+p_1)(s+p_2)}R(s).\] To find the output of the system when the input is a step of magnitude \(A\), we find the inverse Laplace transform of \[Y(s) = \frac{K}{(s+p_1)(s+p_2)}\frac{A}{s}.\] Using partial fraction expansion we get \[Y(s) = \frac{\frac{KA}{p_1p_2}}{s} + \frac{\frac{KA}{(-p_1)(p_2-p_1)}}{s+p_1} + \frac{\frac{KA}{(-p_2)(p_1-p_2)}}{s+p_2}.\] Taking the inverse Laplace transform gives \[y(t) = \begin{cases} \frac{KA}{p_1p_2} + \frac{KA}{(-p_1)(p_2-p_1)}e^{-p_1t} + \frac{KA}{(-p_2)(p_1-p_2)} e^{-p_2t}, & t\geq 0 \\ 0, & \text{otherwise} \end{cases}.\] Note that as \(t\to\infty\), \(y(t)\to \frac{KA}{p_1p_2}\), which is equal to the step size \(A\) times the DC gain of the system.

Now suppose that the closed-loop system has two complex poles and no zeros. We will write a generic second-order system with poles in the left half-plane using the notation \[H(s) = \frac{K\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2},\] where \(K\) is the DC gain of the system (typically \(K=1\)), \(\omega_n>0\) is called the natural frequency, and \(\zeta>0\) is called the damping ratio. The poles are given by the roots of the characteristic polynomial \[\Delta(s) = s^2 + 2\zeta\omega_n s + \omega_n^2=0,\] in other words, the poles are \[\begin{aligned} p_{1,2} &= -\zeta\omega_n \pm \sqrt{(\zeta\omega_n)^2-\omega_n^2} \\ &= -\zeta\omega_n \pm \omega_n\sqrt{\zeta^2-1}. \end{aligned}\] If \(0<\zeta<1\), then the roots are complex and given by \[\tag{8.2} p_{1,2} = -\zeta\omega_n \pm j\omega_n\sqrt{1-\zeta^2}.\] Using the Python code

import control as cnt
import matplotlib.pyplot as plt

#start by varying the natural frequency
plt.figure(1), plt.clf()
zeta = 0.707

for wn in [0.5, 1.0, 2.0, 5.0]:
    G = cnt.tf(wn**2, [1, 2*zeta*wn, wn**2])
    t,y = cnt.step_response(G)
    plt.plot(t, y)

plt.plot([0, 20.0] , [1, 1], 'k--')
plt.xlabel('Time (seconds)')
plt.ylabel('Amplitude')
plt.title('Step Response')
plt.legend(('wn=0.5', 'wn=1.0', 'wn=2.0', 'wn=5.0'))

#now we'll see the effect of varying zeta
plt.figure(2), plt.clf()
wn = 1.0

for zeta in [0.2, 0.4, 0.707, 1.0, 2.0]:
    G = cnt.tf(wn**2, [1, 2*zeta*wn, wn**2])
    t,y = cnt.step_response(G)
    plt.plot(t, y)
    
plt.plot([0, 35.0] , [1, 1], 'k--')
plt.xlabel('Time (seconds)')
plt.ylabel('Amplitude')
plt.title('Step Response')
plt.legend(('zeta=0.2', 'zeta=0.4', 'zeta=0.707',
            'zeta=1.0', 'zeta=2.0'))

plt.show() # necessary to make plots appear

we obtain second-order response plots similar to the ones shown in Figures 8.3 and 8.4. Note from Fig. 8.3 that as the natural frequency increases, the rise time decreases. We should note that a similar phenomenon to Fig. 8.1 also occurs, where the rise time is independent of the size of the step. Note from Fig. 8.4 that the damping ratio affects the amount of ringing in the system. For a small damping ratio, e.g., \(\zeta=0.2\), there is a large overshoot and significant ringing in the system. At the other extreme, when \(\zeta\) is larger than one, the poles are real and the response is highly damped. The sweet spot is when \(\zeta=0.707\), where the rise time is small but the overshoot is minimal.

Note from Equation 8.2 that when \(\zeta = 0.707 = \frac{1}{\sqrt{2}}\), the poles are at \[p_{1,2} = -\frac{\omega_n}{\sqrt{2}} \pm j\frac{\omega_n}{\sqrt{2}}.\] In the complex plane, the pole locations are shown in Fig. 8.5, where \(\omega_n\) is the distance from the origin to the poles. When \(\zeta=\frac{1}{\sqrt{2}}\) the angle \(\theta\) equals \(45\) degrees.

In rectangular coordinates, the poles are given by \[\tag{8.3} p_{1,2} = -\sigma \pm j\omega_d,\] as shown in Fig. 8.5. Equations 8.2 and 8.3 give \[\begin{aligned} \sigma &= \zeta \omega_n \\ \omega_d &= \omega_n\sqrt{1-\zeta^2}. \end{aligned}\] From Fig. 8.5 it is also clear that \[\sin\theta = \frac{\sigma}{\omega_n} = \zeta.\] Therefore \[\theta = \sin^{-1} \zeta.\]

Since the pole locations are at \(-\zeta\omega_n \pm \omega_n\sqrt{\zeta^2-1}\), the poles are both real if \(\zeta\geq 1\). When \(\zeta=1\) the poles are both located at \(-\omega_n\). When \(0<\zeta<1\) the poles are imaginary and occur in complex conjugate pairs. As shown in Fig. 8.6, we say that the system is overdamped when \(\zeta>1\), we say that the system is critically damped when \(\zeta=1\), and we say that the system is underdamped when \(0<\zeta<1\).

The response of \(H(s)\) to a step of size \(A\) is the inverse Laplace transform of \[Y(s) = \frac{A\omega_n^2}{s(s^2+2\zeta\omega_n s + \omega_n^2)},\] which, as shown in Appendix 29, when the poles are complex results in \[\tag{8.4} y(t) = A\left[ 1-\frac{1}{\sqrt{1-\zeta^2}}e^{-\sigma t}\cos\left(\omega_d t - \tan^{-1}\left(\frac{\sigma}{\omega_d}\right)\right)\right].\] A plot of this signal is shown in Fig. 8.7.

The rise time is defined to be the time that it takes the output to transition from 10% to 90% of its final value \(A\). A suitable approximation is that the rise time is one half of the peak time, which is the time when \(y(t)\) reaches its first peak. To find the peak time \(t_p\), differentiate \(y(t)\) in Equation 8.4 and then solve for the first instance when the derivative is zero. The result is \[t_p = \frac{\pi}{\omega_d} = \frac{\pi}{\omega_n\sqrt{1-\zeta^2}}.\] Therefore, the rise time is approximately \[t_r \approx \frac{1}{2}t_p = \frac{1}{2}\frac{\pi}{\omega_n\sqrt{1-\zeta^2}}.\] When \(\zeta=0.707\) we get \[\tag{8.5} t_r \approx \frac{2.2}{\omega_n}.\]

Specifications for second-order systems are usually given either in terms of the rise time \(t_r\) and the damping ratio \(\zeta\), or in terms of the natural frequency \(\omega_n\) and the damping ratio \(\zeta\).

8.1.3 Effect of a Zero on the Step Response

In the previous section, we looked at the step response when there were two complex poles and no zeros. In this section we will briefly discuss the effect of zeros on the step response of the system. This is a topic that we could treat at a much deeper level.

A transfer function is called rational if it is the ratio of two polynomials in \(s\). For example, \[H(s) = \frac{s+1}{s^3+3s^2+s+4}\] is a rational transfer function, whereas \[H(s) = \frac{(s+1)e^{-s}}{s^3+3s^2+s+4}\] is not. A rational transfer function can be written as \[H(s) = \frac{N(s)}{D(s)},\] where \(N(s)\) is the numerator polynomial and \(D(s)\) is the denominator polynomial. The zeros of the transfer function are roots of the equation \(N(s)=0\), and the poles of the system are the roots of the equation \(D(s)=0\).

One way to understand the effect of zeros is to consider the effect that zeros have on the partial fraction expansion of a transfer function. For a review of partial fraction expansion, see Appendix 29. First consider a two-pole system with characteristic equation given by \(\Delta(s) = (s+2)(s+3)\). The transfer function with DC gain equal to one and without zeros is given by \[H(s) = \frac{6}{(s+2)(s+3)}.\] The response to a unit step is given by \begin{align} y(t) &= \mathcal{L}^{-1}\left\{\frac{6}{s(s+2)(s+3)}\right\} \notag \\ &= \mathcal{L}^{-1}\left\{ \frac{1}{s} + \frac{-3}{s+2} + \frac{2}{s+3}\right\} \notag \\ &= 1 - 3e^{-2t} + 2e^{-3t}, \quad t\geq 0. \tag{8.6} \end{align} Now consider the step response when a zero has been added at \(z=-2.1\) and the gain adjusted so that the DC gain is still one, i.e., when \[H(s) = \frac{\frac{6}{2.1}(s+2.1)}{(s+2)(s+3)}.\] The response to a unit step is given by \[\begin{aligned} y(t) &= \mathcal{L}^{-1}\left\{\frac{\frac{6}{2.1}(s+2.1)}{s(s+2)(s+3)}\right\} \\ &= \mathcal{L}^{-1}\left\{ \frac{1}{s} + \frac{-0.1429}{s+2} + \frac{-0.8571}{s+3}\right\}\\ &= 1 - 0.1429e^{-2t} - 0.8571e^{-3t}, \quad t\geq 0. \end{aligned}\] Note that \(y(t)\) still settles to one, but the response of the pole at \(-2\) has almost been removed. If, however, the zero is far away from the poles, then the effect is much less. Suppose that the zero is at \(z=-60\) and the gain is adjusted so that the DC gain is still one, i.e., \[H(s) = \frac{\frac{6}{60}(s+60)}{(s+2)(s+3)}.\] The response to a unit step is given by \[\begin{aligned} y(t) &= \mathcal{L}^{-1}\left\{\frac{\frac{6}{60}(s+60)}{s(s+2)(s+3)}\right\} \\ &= \mathcal{L}^{-1}\left\{ \frac{1}{s} + \frac{-2.9}{s+2} + \frac{1.9}{s+3}\right\}\\ &= 1 - 2.9e^{-2t} + 1.9e^{-3t}, \quad t\geq 0. \end{aligned}\] In this case, the response is similar to Equation 8.6.

Another way to understand the effect of zeros on the system is to recall that an \(s\) in the numerator is essentially a differentiator. Therefore \[\begin{aligned} H(s) &= \frac{\omega_n^2(\tau s+1)}{s^2 + 2\zeta\omega_n s + \omega_n^2} \\ &= \tau\frac{\omega_n^2 s}{s^2 + 2\zeta\omega_n s + \omega_n^2} + \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}. \end{aligned}\] Therefore, the step response is given by \[\begin{aligned} y(t) &= \mathcal{L}^{-1}\left\{\tau\frac{\omega_n^2s}{s(s^2 + 2\zeta\omega_n s + \omega_n^2)} + \frac{\omega_n^2}{s(s^2 + 2\zeta\omega_n s + \omega_n^2)}\right\} \\ &= \tau\frac{d}{dt}\mathcal{L}^{-1}\left\{\frac{\omega_n^2}{s(s^2 + 2\zeta\omega_n s + \omega_n^2)}\right\} % + \mathcal{L}^{-1}\left\{\frac{\omega_n^2}{s(s^2 + 2\zeta\omega_n s + \omega_n^2)}\right\}. \end{aligned}\] Therefore, the effect of the zero is to add \(\tau\) times the derivative of the step response to what the step response would have been without the zero. Fig. 8.8 shows the step response for different values of \(\tau\), illustrating that adding a zero has a strong effect on the rise time and damping of the response. Note also that when \(\tau<0\), i.e., when zeros are in the right half-plane, the system response first goes negative before going positive. Right-half-plane zeros are often called nonminimum-phase zeros , and the associated response is often called a nonminimum-phase response. Physical systems that have right-half-plane zeros include the inverted pendulum and aircraft altitude-control systems. For the inverted pendulum, to move the cart to the right from a stopped position, the cart must first move left so that the rod leans right before the cart can move to its desired position. Similarly, an aircraft pitches up to climb. However, when the aircraft pitches up, it loses speed, which causes it to lose lift and altitude. As the aircraft loses altitude, it gains speed and produces the lift necessary to climb.

8.1.4 Successive Loop Closure

The natural frequency \(\omega_n\) and the damping ratio \(\zeta\) are only defined for second-order systems, and therefore, the PD design technique based on \(\omega_n\) and \(\zeta\) are only valid for second-order systems. However, these concepts can be extended to higher-order systems through a technique called successive loop closure, or inner-loop, outer-loop design. When we derived the Laplace transform in Chapter 8, we saw that for many systems, the transfer functions can be arranged in cascade form as shown in Fig. 8.9 where the input to the system is \(u_1\) and the output is \(y_2\), but where \(y_1\) and \(u_2\) are intermediary values that represent the natural physical coupling between internal variables. Many physical systems can be represented by a cascade of systems.

When the system is represented as a cascade, the control strategy can be designed using successive loop closure. The basic idea is to close feedback loops in succession around the open-loop plant dynamics rather than designing a single control system. To illustrate this approach, consider the open-loop system shown in Fig. 8.9. The open-loop dynamics are given by the product of two transfer functions in series: \(P(s)=P_1(s)P_2(s)\). We assume that each of the transfer functions has an output (\(y_1\), \(y_2\)) that can be measured and used for feedback. Typically, each of the transfer functions, \(P_1(s)\) and \(P_2(s)\), is of relatively low order – usually first or second-order. In this case, we are interested in controlling the output \(y_2\). Instead of closing a single feedback loop around \(y_2\), we will close feedback loops around \(y_1\) and \(y_2\) in succession as shown in Fig. 8.10.

We will design the compensators \(C_1(s)\) and \(C_2(s)\) in succession. A necessary condition in the design process is that the inner loop must be much faster than each successive loop. As a rule of thumb, if the rise time of the inner loop is \(t_{r_1}\), the rise time of the next loop should be 5 to 10 times longer, i.e., \(t_{r_2}> W_1 t_{r_1}\), where \(W_1\) is a design parameter, usually on the order of \(5-10\). In part V of this book, we will talk about the frequency response of the system. At that time, we will provide a deeper explanation that can be understood in terms of the bandwidth of each successive loop.

Examining the inner loop shown in Fig. 8.10, the goal is to design a closed-loop system from \(r_1\) to \(y_1\) having a rise time of \(t_{r_1}\). If the rise time of the inner loop is significantly faster than the rise time of the outer loop, and if the DC gain of the inner loop is \(k_{DC_1}\), then relative to the outer loop the inner loop can be effectively modeled as the DC gain. This is depicted schematically in Fig. 8.11. With the inner-loop transfer function modeled as its DC gain, design of the outer loop is simplified because it only includes the plant transfer function \(P_2(s)\) and the compensator \(C_2(s)\) and the DC gain \(k_{DC_1}\).

Because both plant models \(P_1(s)\) and \(P_2(s)\) are first or second-order, conventional PID compensators can be employed. The loopshaping techniques discussed in Part V can also be used for successive loop closure of higher-order systems.

8.1.5 Input Saturation and Limits of Performance

A natural question is to ask if the rise time can be made arbitrarily small. Unfortunately, saturation limits on the actuators place a lower bound on the achievable rise time. All physical systems have saturation limits on their actuators. For example, in the robot arm problem, the torque that can be applied to the arm is limited by the motor that is used to apply that torque. Force and torque will always be limited by physical constraints like current limits. Similarly, for an airplane, the rudder command is limited by how far the rudder can physically move, which is typically on the order of \(\pm 40\) degrees. A natural way to model input constraints is to add a saturation block preceding the physical plant. Fig. 8.12 shows a saturation block in conjunction with a PD control scheme.

Mathematically, the saturation block is given by \[u_{sat} = \begin{cases} u_{\max}, & \text{~if~} u\geq u_{\max} \\ u, & \text{~if~} -u_{\max} \leq u \leq u_{\max} \\ -u_{\max}, & \text{~otherwise} \end{cases},\] where the input to the block is \(u\), the output is \(u_{sat}\), and the saturation limit is \(u_{\max}\). Note that the saturation block is non-linear. Therefore, one design strategy is to select the feedback gains so that the system remains linear by keeping the input out of saturation.

To see how this might be done, consider the second-order system shown in Fig. 8.12 with proportional feedback on the output error and derivative feedback on the output. When the input is not saturated, the closed-loop transfer function is \[Y(s) = \left(\frac{b_0 k_p}{s^2 + (a_1 + b_0 k_D)s + (a_0 + b_0 k_p)}\right) R(s). \tag{8.7}\] We can see that the closed-loop poles of the system are defined by the selection of the control gains \(k_p\) and \(k_D\). Note also that the actuator effort \(u\) can be expressed as \(u=k_p e - k_D\dot{y}\). When \(\dot{y}\) is zero or small, the size of the actuator effort \(u\) is primarily governed by the size of the control error \(e\) and the control gain \(k_p\). If the system is stable, the largest control effort in response to a step input will occur immediately after the step when \(\dot{y}=0\). If the system is to be kept just out of saturation, then immediately after the step we have \(\left|u\right|=u_\text{max} = k_p e_\text{max}\). Rearranging this expression, we find that the proportional control gain can be determined from the maximum anticipated output error and the saturation limits of the actuator as \[k_p = \pm \frac{u_\text{max}}{e_\text{max}}, \tag{8.8}\] where \(u_{\max}\) is the maximum control effort the system can provide, \(e_{\max}\) is the step error that results from a step input of nominal size, and the sign of \(k_p\) is determined by the physics of the system.

Suppose that the desired closed-loop characteristic polynomial is given by \[\Delta_{cl}^d(s) = s^2 + 2\zeta\omega_n s + \omega_n^2,\] then by comparing with the actual closed-loop characteristic polynomial \[\Delta_{cl}(s) = s^2 + (a_1 + b_0 k_D)s + (a_0 + b_0 k_p)\] we have that \[\begin{aligned} \omega_n &= \sqrt{a_0 + b_0 k_p} \leq \sqrt{a_0 + b_0 \frac{u_\text{max}}{e_\text{max}}}, \end{aligned}\] which is an upper limit on the natural frequency of the closed-loop system, ensuring that saturation of the actuator is avoided. The corresponding rise time constraint is \[t_r \geq \frac{2.2}{\omega_n} = \frac{2.2}{\sqrt{a_0 + b_0 \frac{u_\text{max}}{e_\text{max}}}},\] which is the smallest possible rise time for an error step of \(e_{\text{max}}\).

If the input has an equilibrium value, then the saturation constraint used in the above calculation must be modified to account for the equilibrium. For example, suppose that \(u=u_e+\tilde{u}\) where \(u_e\) is the equilibrium value and \(\tilde{u}\) is the output of the PID controller. Also suppose that the input saturation constraint is given by \(u_{\min}\leq u\leq u_{\max}\). The objective is to find a value for \(\tilde{u}_{\max}\) where \(\left|\tilde{u}\right|\leq \tilde{u}_{\max}\) guarantees that \(u_{\min}\leq u\leq u_{\max}\). We have that \[\begin{aligned} & u_{\min}\leq u \leq u_{\max} \\ \implies & u_{\min} \leq u_e + \tilde{u} \leq u_{\max} \\ \implies & u_{\min}-u_e \leq \tilde{u} \leq u_{\max}-u_e. \end{aligned}\] The desired value for \(\tilde{u}_{\max}\) is therefore the minimum (in absolute value) of the right and left-hand sides of this expression. The situation is shown in Fig. 8.13, where it is clear that the value for \(\tilde{u}_{\max}\) can be written as \[\tilde{u}_{\max} = \begin{cases} \left|u_{\max}-u_e\right|, & \text{~if~} u_e \geq \frac{u_{\max}+u_{\min}}{2} \\ \left|u_e - u_{\min}\right|, & \text{~if~} u_e \leq \frac{u_{\max}+u_{\min}}{2}. \end{cases}\]

The takeaway from this discussion is that the rise time is limited by the input saturation constraint and cannot be made arbitrarily small. In practice, the control parameters are tuned so that the closed-loop system has the fastest possible rise time without driving the plant input into saturation.

8.2 Design Study A. Single Link Robot Arm

8.2.1 Example Problem 3.1.11

For the single-link robot arm, do the following:

(a) Suppose that the design requirements are that the rise time is \(t_r \approx 0.8\) seconds, with a damping ratio of \(\zeta=0.707\). Find the desired closed-loop characteristic polynomial \(\Delta_{cl}^d(s)\) and the associated pole locations. Find the proportional and derivative gains \(k_P\) and \(k_D\) to achieve these specifications, and modify the simulation from HW 3.1.8 to verify that the step response satisfies the requirements.

(b) Suppose that the input torque for the robot arm is limited to \(\left|\tau\right|\leq\tau_{\max}=1\) Nm. Modify the simulation to include a saturation on the torque \(\tau\). Using the rise time \(t_r\) and damping ratio \(\zeta\) as tuning parameters, tune the PD control law so that the input just saturates when a step of size \(50\) degrees is placed on \(\tilde{r}_{\theta}\).