Spacecraft and Rockets, Notas de estudo de Engenharia Aeroespacial

Baixe Spacecraft and Rockets e outras Notas de estudo em PDF para Engenharia Aeroespacial, somente na Docsity! ADVANCED CONTROL OF AIRCRAFT, SPACECRAFT AND ROCKETS Aerospace Series List Basic Helicopter Aerodynamics: Third Edition Seddon and Newman July 2011 Cooperative Path Planning of Unmanned Aerial Vehicles Tsourdos et al November 2010 Principles of Flight for Pilots Swatton October 2010 Air Travel and Health: A Systems Seabridge et al September 2010 Perspective Design and Analysis of Composite Kassapoglou September 2010 Structures: With applications to aerospace Structures Unmanned Aircraft Systems: UAVS Design, Development and Deployment Austin April 2010 Introduction to Antenna Placement and Installations Macnamara April 2010 Principles of Flight Simulation Allerton October 2009 Aircraft Fuel Systems Langton et al May 2009 The Global Airline Industry Belobaba April 2009 Computational Modelling and Simulation of Aircraft and the Environment: Volume 1 - Platform Kinematics and Synthetic Environment Diston April 2009 Handbook of Space Technology Ley, Wittmann April 2009 and Hallmann Aircraft Performance Theory and Swatton August 2008 Practice for Pilots Surrogate Modelling in Engineering Forrester, Sobester August 2008 Design: A Practical Guide and Keane Aircraft Systems, Third Edition Moir and Seabridge March 2008 Introduction to Aircraft Aeroelasticity And Loads Wright and Cooper December 2007 Stability and Control of Aircraft Systems Langton September 2006 Military Avionics Systems Moir and Seabridge February 2006 Design and Development of Aircraft Systems Moir and Seabridge June 2004 Aircraft Loading and Structural Layout Howe May 2004 Aircraft Display Systems Jukes December 2003 Civil Avionics Systems Moir and Seabridge December 2002 The intellectual or logical man, rather than the understanding or observant man, set himself to imagine designs – to dictate purposes to God. – Edgar Allan Poe Contents Series Preface xiii Preface xv 1 Introduction 1 1.1 Notation and Basic Definitions 1 1.2 Control Systems 3 1.2.1 Linear Tracking Systems 7 1.2.2 Linear Time-Invariant Tracking Systems 9 1.3 Guidance and Control of Flight Vehicles 10 1.4 Special Tracking Laws 13 1.4.1 Proportional Navigation Guidance 13 1.4.2 Cross-Product Steering 16 1.4.3 Proportional-Integral-Derivative Control 19 1.5 Digital Tracking System 24 1.6 Summary 25 Exercises 26 References 28 2 Optimal Control Techniques 29 2.1 Introduction 29 2.2 Multi-variable Optimization 31 2.3 Constrained Minimization 33 2.3.1 Equality Constraints 34 2.3.2 Inequality Constraints 38 2.4 Optimal Control of Dynamic Systems 41 2.4.1 Optimality Conditions 43 2.5 The Hamiltonian and the Minimum Principle 44 2.5.1 Hamilton–Jacobi–Bellman Equation 45 2.5.2 Linear Time-Varying System with Quadratic Performance Index 47 2.6 Optimal Control with End-Point State Equality Constraints 48 2.6.1 Euler–Lagrange Equations 50 2.6.2 Special Cases 50 2.7 Numerical Solution of Two-Point Boundary Value Problems 52 2.7.1 Shooting Method 54 2.7.2 Collocation Method 57 2.8 Optimal Terminal Control with Interior Time Constraints 61 2.8.1 Optimal Singular Control 62 2.9 Tracking Control 63 viii Contents 2.9.1 Neighboring Extremal Method and Linear Quadratic Control 64 2.10 Stochastic Processes 69 2.10.1 Stationary Random Processes 75 2.10.2 Filtering of Random Noise 77 2.11 Kalman Filter 77 2.12 Robust Linear Time-Invariant Control 81 2.12.1 LQG/LTR Method 82 2.12.2 H2/H∞ Design Methods 89 2.13 Summary 96 Exercises 98 References 101 3 Optimal Navigation and Control of Aircraft 103 3.1 Aircraft Navigation Plant 104 3.1.1 Wind Speed and Direction 110 3.1.2 Navigational Subsystems 112 3.2 Optimal Aircraft Navigation 115 3.2.1 Optimal Navigation Formulation 116 3.2.2 Extremal Solution of the Boundary-Value Problem: Long-Range Flight Example 119 3.2.3 Great Circle Navigation 121 3.3 Aircraft Attitude Dynamics 128 3.3.1 Translational and Rotational Kinetics 132 3.3.2 Attitude Relative to the Velocity Vector 135 3.4 Aerodynamic Forces and Moments 136 3.5 Longitudinal Dynamics 139 3.5.1 Longitudinal Dynamics Plant 142 3.6 Optimal Multi-variable Longitudinal Control 145 3.7 Multi-input Optimal Longitudinal Control 147 3.8 Optimal Airspeed Control 148 3.8.1 LQG/LTR Design Example 149 3.8.2 H∞ Design Example 160 3.8.3 Altitude and Mach Control 166 3.9 Lateral-Directional Control Systems 173 3.9.1 Lateral-Directional Plant 173 3.9.2 Optimal Roll Control 177 3.9.3 Multi-variable Lateral-Directional Control: Heading-Hold Autopilot 180 3.10 Optimal Control of Inertia-Coupled Aircraft Rotation 183 3.11 Summary 189 Exercises 192 References 194 4 Optimal Guidance of Rockets 195 4.1 Introduction 195 4.2 Optimal Terminal Guidance of Interceptors 195 4.3 Non-planar Optimal Tracking System for Interceptors: 3DPN 199 4.4 Flight in a Vertical Plane 208 4.5 Optimal Terminal Guidance 211 Contents xi B.3.2 Global Asymptotic Stability 406 B.3.3 Lyapunov’s Theorem 407 B.3.4 Krasovski’s Theorem 408 B.3.5 Lyapunov Stability of Linear Systems 408 References 408 Appendix C: Control of Underactuated Flight Systems 409 C.1 Adaptive Rocket Guidance with Forward Acceleration Input 409 C.2 Thrust Saturation and Rate Limits (Increased Underactuation) 415 C.3 Single- and Bi-output Observers with Forward Acceleration Input 417 References 432 Index 433 Series Preface The field of aerospace is wide ranging and multi-disciplinary, covering a large variety of products, disciplines and domains, not merely in engineering but in many related supporting activities. These combine to enable the aerospace industry to produce exciting and technologically advanced products. The wealth of knowledge and experience that has been gained by expert practitioners in the various aerospace fields needs to be passed onto others working in the industry, including those just entering from University. The Aerospace Series aims to be practical and topical series of books aimed at engineering profes- sionals, operators, users and allied professions such as commercial and legal executives in the aerospace industry. The range of topics is intended to be wide ranging, covering design and development, manufac- ture, operation and support of aircraft as well as topics such as infrastructure operations and developments in research and technology. The intention is to provide a source of relevant information that will be of interest and benefit to all those people working in aerospace. Flight Control Systems are an essential element of all modern aerospace vehicles, for instance: civil aircraft use a flight management system to track optimal flight trajectories, autopilots are used to follow a set course or to land the aircraft, gust and manoeuvre load alleviation systems are often employed, military aircraft are often designed for enhanced manoeuvrability with carefree handling in flight regimes that have statically unstable open loop behaviour, Unmanned Autonomous Vehicles (UAVs) are able to follow a pre-defined task without any human intervention, rockets require guidance and control strategies that are able to deal with rapidly time varying parameters whilst ensuring that rocket’s orientation follows the desired flight path curvature, and spacecraft need to be guided accurately on missions that may last many years whilst also maintaining the desired attitude throughout. This book, Advanced Control of Aircraft, Spacecraft and Rockets, considers the application of optimal control theory for the guidance, navigation and control of aircraft, spacecraft and rockets. A uniform approach is taken in which a range of modern control techniques are described and then applied to a number of non-trivial multi-variable problems relevant to each type of vehicle. The book is seen as a valuable addition to the Aerospace Series. Peter Belobaba, Jonathan Cooper, Roy Langton and Allan Seabridge Preface There are many good textbooks on design of flight control systems for aircraft, spacecraft, and rockets, but there is none that covers all of them in a similar manner. The objective of this book is to show that modern control techniques are applicable to all flight vehicles without distinction. The main focus of the presentation is on applications of optimal control theory for guidance, navigation, and control of aircraft, spacecraft, and rockets. Emphasis is placed upon non-trivial, multi-variable applications, rather than those that are handled by single-variable methods. While the treatment can be considered advanced and requires a basic course in control systems, the topics are covered in a self-contained manner, without the need frequently to refer to control theory textbooks. The spectrum of flight control topics covered in the book is rather broad, ranging from the optimal roll control of aircraft and missiles to multi-variable, adaptive guidance of gravity-turn rockets, and two-point boundary value, optimal terminal control of aircraft, spacecraft, and rockets. However, rotorcraft and flexible dynamics have been excluded for want of space. The task of designing an optimal flight controller is best explained by the parable of a pilgrim (cf. The Pilgrim’s Progress by John Bunyan), who, in order to attain a certain religious state, xf , beginning from an initial state, xi, must not stray too far from the righteous (optimal) path. In other words, the pilgrim must be aware of the desired trajectory to be followed, xd(t), ti ≤ t ≤ tf , and must also have an in-built tracking system that detects deviations therefrom and applies corrective control inputs, u(t), in order to return to the optimal path. A successful tracking system is the one which keeps deviations from the desired trajectory, || x(t) − xd(t) ||, ti ≤ t ≤ tf , to small values. A linear tracking system can be designed for maintaining small deviations by the use of linear optimal feedback control, and is guaranteed to be successful provided the initial deviation is not too large. The tracking system must only require control inputs, u(t), that are within the capability of a pilgrim to generate. Furthermore, the tracking system must be robust with respect to small (but random) external disturbances (like the various enticements and drawbacks experienced on a pilgrimage), which would otherwise cause instability, that is, a total departure from the righteous path. The job of generating a desired trajectory, xd(t), ti ≤ t ≤ tf , which will be followed by the tracking system from a given initial state, xi, to the desired terminal state, xf , is by no means an easy one. It requires an understanding of the pilgrim’s behavior with time, x(t), which is modeled for example by a nonlinear, ordinary differential equation, ẋ(t) = f[x(t), u(t), t], and should be solved subject to the boundary conditions at the two points, x(ti) = xi, x(tf ) = xf . Now, it is quite possible that there exist several solutions to such a boundary value problem, each depending upon the applied control inputs, u(t). One must, however, select the solution which requires 1 Introduction 1.1 Notation and Basic Definitions Throughout this book, we shall represent scalar quantities by italic letters and symbols (a, α), vectors in boldface (a, α), and matrices in boldface capitals (A). Unless otherwise mentioned, the axes of a frame are specified by unit vectors in the same case as that of respective axis labels, for example, the axis ox would be represented by the unit vector, i, while OX is given by I. The axes are labeled in order to constitute a right-handed triad (e.g., i × j = k). Components of a vector have the same subscripts as the axis labels, for example, a vector a resolved in the frame oxyz is written as a = axi + ayj + azk, (1.1) or alternatively as a = ⎧⎨ ⎩ ax ay az ⎫⎬ ⎭ . (1.2) An overdot represents a vector (or matrix) derived by taking the time derivative of the components in a frame of reference, for example, ȧ .= ⎧⎪⎨ ⎪⎩ dax dt day dt daz dt ⎫⎪⎬ ⎪⎭ = ⎧⎨ ⎩ ȧx ȧy ȧz ⎫⎬ ⎭ . (1.3) The vector product of two vectors a, b is often expressed as the matrix product a × b = S(a)b, where S(a) is the following skew-symmetric matrix of the components of a: S(a) = ⎛ ⎝ 0 −az ay az 0 −ax −ay ax 0 ⎞ ⎠ . (1.4) Advanced Control of Aircraft, Spacecraft and Rockets, Ashish Tewari. © 2011 John Wiley & Sons, Ltd. Published 2011 by John Wiley & Sons, Ltd. Appendix A Linear Systems The theory of linear systems refers to the mathematical framework of a discretized (finite-dimensional) system linearized about a particular solution. For details we refer the reader to Kailath (1980). A.1 Definition A system with state vector, ξ(t), and input vector, η(t), governed by vector state equation, ξ̇ = f(ξ, η, t), (A.1) initial condition, ξ(0) = ξ0, (A.2) and output equation, y = g(ξ, η, t), (A.3) is said to be linear if its output vector resulting from the applied input vector, η(t) = c1η1(t) + c2η2(t), (A.4) is given by y(t) = c1y1(t) + c2y2(t), (A.5) where y1(t) and y2(t) are the output vectors of the system to the inputs η1(t) and η2(t), respectively, and c1, c2 are arbitrary scalar constants. By inspecting the governing equations of a system, equations (A.1) and (A.3), it is possible to determine whether it is linear. If the functionals f(.), g(.) are continuous and do not contain nonlinear functions of the state and input variables, then the system is linear. Advanced Control of Aircraft, Spacecraft and Rockets, Ashish Tewari. © 2011 John Wiley & Sons, Ltd. Published 2011 by John Wiley & Sons, Ltd. 392 Appendix A: Linear Systems A.2 Linearization Let a reference state vector, ξr(t), and a corresponding reference input vector, ηr(t), satisfy the system’s governing vector state equation (A.1), ξ̇r = f(ξr, ηr, t), (A.6) subject to initial condition, ξr(0) = ξr0. (A.7) Then, let x(t) and u(t) be state and control deviations, respectively, from the reference solution, (ξr, ηr), such that the perturbed solution is given by ξ(t) = ξr(t) + x(t) η(t) = ηr(t) + u(t), (A.8) subject to initial conditions (A.2) and (A.7), implying x(0) = ξ(0) − ξr(0) = ξ0 − ξr0 = x0. (A.9) Assuming that the vector functional, f(.), possesses continuous derivatives with respect to state and control variables up to an infinite order at the reference solution, (ξr, ηr), we expand the state equation about the reference solution by neglecting the quadratic and higher-order terms as follows: ξ̇ − ξ̇r = ẋ = f(ξr + x, ηr + u, t) − f(ξr, ηr, t), (A.10) where f(ξr + x, ηr + u, t)  f(ξr, ηr, t) + ∂f ∂ξ (ξr, ηr, t)x + ∂f ∂η (ξr, ηr, t)u. Substitution of equation (A.11) into equation (A.10) yields the following linearized state equation about the reference solution: ẋ = A(t)x + B(t)u, (A.11) subject to initial condition (A.9), where A(t) and B(t) are the following Jacobian matrices: A(t) = ∂f ∂ξ (ξr, ηr, t) B(t) = ∂f ∂η (ξr, ηr, t). (A.12) A.3 Solution to Linear State Equations The solution to the linear state equation, ẋ = A(t)x + B(t)u, (A.13) subject to initial condition, x(0) = x0, (A.14) is expressed as the sum of homogeneous and particular solutions. Appendix A: Linear Systems 395 A.5 Linear Time-Invariant Stability Criteria Equation (A.28) for U(s) = 0 represents an eigenvalue problem, whose solution yields the eigenvalues, s, and eigenvectors, X(s). The eigenvalues of the linear system are obtained by the following characteristic equation: | sI − A |= 0. (A.32) The n generally complex roots of the characteristic equation (eigenvalues of A) signify an important system property, called stability (Appendix B). Considering that an eigenvalue is generally complex, its imaginary part denotes the frequency of oscillation of the characteristic vector about the equilibrium point, and the real part signifies the growth (or decay) of its amplitude with time. We can now state the following criteria for the stability of an LTI system: (a) If all eigenvalues have negative real parts, the system is asymptotically stable and regains its equi- librium in the steady state. (b) A system having complex eigenvalues with zero real parts (and all other eigenvalues with negative real parts) displays oscillatory behavior of a constant amplitude, and is said to be stable (but not asymptotically stable). (c) If at least one eigenvalue has a positive real part, its contribution to the system’s state is an exponen- tially growing amplitude, and the system is said to be unstable. (d) If a multiple eigenvalue of multiplicity p is at the origin (i.e., has both real and imaginary parts zero), its contribution to the system’s state has terms containing the factors ti , i = 0, . . . , p − 1, which signify an unbounded behavior with time. Hence, the system is unstable. A.6 Controllability of Linear Time-Invariant Systems Controllability is the property of a system where it is possible to move it from any initial state, x(ti), to any final state, x(tf ), solely by the application of the input vector, u(t), acting in a finite control interval, ti ≤ t ≤ tf . The words “any” and “finite” are emphasized because it may be possible to move an uncontrollable system only between some specific states by applying the control input, or to take an infinite time to change an uncontrallable system between arbitrary states. For a system to be controllable, all its state variables must be influenced, either directly or indirectly, by control inputs. If there is a subsystem that is unaffected by the control inputs, then the entire system is uncontrollable. The necessary and sufficient condition for controllability of an LTI system (A, B) is that the following test matrix must have rank n, the order of the system: P = (B, AB, A2B, . . . , An−1B) . (A.33) If a system is unstable but controllable, it can be stabilized by a feedback control system. It is often possible to decompose an uncontrollable LTI system into controllable and uncontrollable subsystems. A system that is both unstable and uncontrollable could also be stabilized, provided its uncontrollable subsystem is stable. In such a case, the system is said to be stabilizable. A.7 Observability of Linear Time-Invariant Systems Observability is the property of an unforced (homogeneous) system where it is possible to estimate any initial state, x(ti), of the system solely by a finite record, t ≥ ti, of the output vector, y(t). For a system to be observable, all of its state variables must contribute, either directly or indirectly, to the output vector. If there is a subsystem that leaves the output vector unaffected, then the entire system is unobservable. 396 Appendix A: Linear Systems An unforced LTI system, ẋ = Ax, whose output is related to the state vector by y = Cx, (A.34) is observable1 if and only if the following test matrix has rank n, the order of the system: N = (CT , AT CT , (AT)2CT , . . . , (AT )n−1CT ) . (A.35) It is often possible to decompose an unobservable LTI system into observable and unobservable subsys- tems. A system whose unobservability is caused by a stable subsystem is said to be detectable. A.8 Transfer Matrix An LTI system’s response to a specified input is most efficiently analyzed and designed with the the concept of the transfer matrix. The transfer matrix, G(s), of such a system is defined as the following relationship between the output’s Laplace transform, Y(s), and that of the input vector, U(s), subject to zero initial conditions, y(0) = ẏ(0) = ÿ(0) = . . . = 0: Y(s) = G(s)U(s). (A.36) Clearly, the transfer matrix is the Laplace transform of the impulse response matrix subject to zero initial conditions. The roots of the common denominator polynomial of the transfer matrix are called the poles of the system and are same as the eigenvalues of the system’s state dynamics matrix, A. The roots of the numerator polynomials of the transfer matrix are called the zeros of the system. The transfer matrix of an LTI system can be expressed in terms of its state-space coefficients as G(s) = C (sI − A)−1 B + D. (A.37) When s = iω the transfer matrix becomes the frequency response matrix, G(iω), whose elements denote the steady-state response of an output variable to a simple harmonic input variable, subject to zero initial conditions and all other inputs being zero. A.9 Singular Value Decomposition Singular values of a transfer matrix, G(iω), at a given frequency, ω, are the positive square roots of the eigenvalues of the square matrix GT (−iω)G(iω). The computation of singular values is carried out by the singular value decomposition G(iω) = U(iω)VT (−iω), (A.38) where U, V are unitary complex matrices satisfying U(iω)UT (−iω) = I, V(iω)VT (−iω) = I, (A.39) and  contains the singular values of G(iω), denoted by σk {G(iω)} , k = 1, 2, . . . , n, 1 This definition of observability can be extended to a forced linear system by requiring in addition that the applied inputs, u(t) should be known in the period of observation, t ≥ ti. Appendix A: Linear Systems 397 as its diagonal elements. The largest of all the singular values is denoted by σ̄ and the smallest by σ. The frequency spectrum of the singular values indicates the variation of the “magnitude” (or gain) of the frequency response matrix, which is supposed to lie between σ(ω) and σ̄(ω). The range of frequencies (0 ≤ ω ≤ ωb) over which σ̄(ω) stays above 0.707σ̄(0) is called the system’s bandwidth and is denoted by the frequency, ωb. The bandwidth indicates the highest-frequency signal to which the system has an appreciable response. For higher frequencies (ω > ωb), the singular values of a strictly proper system generally decline rapidly with increasing frequency, called roll-off. A system with a steep roll-off of the largest singular value, σ̄, has a reduced sensitivity to high-frequency noise, which is a desirable property. The highest magnitude achieved by the largest singular value over the system bandwidth is the H∞ norm of the transfer matrix, given by || G ||∞= supω [σ̄ {G(iω)}] , (A.40) where supω(.) is the supremum (the maximum value) with respect to the frequency. A.10 Linear Time-Invariant Control Design Consider the plant dynamics expressed in a linear state-space form as follows: ẋ = Ax + Bu, (A.41) where x(t) is the state vector, u(t) the control input vector, and A(t), B(t) are the state space coefficient matrices, which could be time-varying. However, for the time being, we will confine our attention to a linear plant with constant coefficient matrices, A, B (i.e., an LTI system). A.10.1 Regulator Design by Eigenstructure Assignment It is possible to design a linear state feedback regulator for the LTI plant of equation (A.41) with the control law u = −Kx (A.42) by assigning a structure for the eigenvalues and eigenvectors of the closed-loop dynamics matrix, A − BK. For single-input plants, this merely involves selecting the locations for the closed-loop poles (pole placement) by the following Ackermann formula that yields the desired closed-loop characteristics (Tewari 2002): K = (ad − a)(PP′)−1. (A.43) Here a is the row vector formed by the coefficients, ai, of the plant’s characteristic polynomial in de- scending order, |sI − A| = sn + ansn−1 + an−1sn−2 + . . . + a2s + a1; (A.44) ad is the row vector formed by the characteristic coefficients of the closed-loop system in descending order, |sI − A + BK| = sn + adnsn−1 + ad(n−1)sn−2 + . . . + ad2s + ad1; (A.45) 400 Appendix A: Linear Systems The unmeasurable part, x2, needs estimation by a reduced-order observer, and can be expressed as follows: xo2 = Ly + z, (A.58) where z is the state vector of the reduced-order observer with state equation ż = Fz + Hu + Gy, (A.59) whose coefficient matrices are determined from the requirement that the estimation error, eo = x2 − xo2, should go to zero in the steady state, irrespective of the control input and the output (Tewari 2002): F = A22 − LCA12, G = FL + (A21 − LCA11) C−1, (A.60) H = B2 − LCB1, with the observer gain L selected by eigenstructure assignment or Kalman filter (Chapter 2) approach, such that all the eigenvalues of the oberver dynamics matrix, F, are in the left-half s-plane. References Kailath, T. (1980) Linear Systems. Prentice Hall, Englewood Cliffs, NJ. Kautsky, J., Nichols, N.K., and van Dooren, P. (1985) Robust pole assignment in linear state feedback. International Journal of Control 41, 1129–1155. Tewari, A. (2002) Modern Control Design with MATLAB and Simulink. John Wiley & Sons, Ltd, Chichester. Appendix B Stability B.1 Preliminaries Stability is broadly defined as the tendency of an unforced system to remain close to a given reference state despite small initial disturbances. Consider an unforced system with state vector, ξ(t), governed by the state equation ξ̇ = f(ξ, t), (B.1) with initial condition ξ(0) = ξ0. (B.2) Let the system have a reference solution, ξr(t), t ≥ 0, to the state equation (B.1), subject to the given initial condition (B.2): ξ̇r = f(ξr, t) (B.3) and ξr(0) = ξ0. (B.4) If there is a change in the initial condition such that ξ(0) = ξ0 + δ, (B.5) then the perturbed solution, ξ(t), t ≥ 0, to the state equation (B.1) may depart from the reference solution such that ξ(t) = ξr(t) + (t), t ≥ 0. (B.6) Then stability is concerned with the boundedness of (t), t ≥ 0, if | δ | is small. The unforced system may have a constant reference solution, ξe, called an equilibrium point, or a state of rest, that is, ξ̇e = f(ξe, t) = 0. (B.7) For an equilibrium solution it must be true that ξ(0) = ξe. (B.8) Advanced Control of Aircraft, Spacecraft and Rockets, Ashish Tewari. © 2011 John Wiley & Sons, Ltd. Published 2011 by John Wiley & Sons, Ltd. 402 Appendix B: Stability Equations (B.7) and (B.8) indicate that by redefining the system’s state as the deviation from an equilib- rium solution, x(t) = ξ(t) − ξe, t ≥ 0, (B.9) we have the transformed state equation ẋ = f(x, t), (B.10) with the initial condition x(0) = 0. (B.11) For a mechanical system, stability is usually defined with reference to an equilibrium point. Further- more, the functional, f(x, t), of a mechanical system is continuous and possesses continuous derivatives with respect to the state and time. There two broad mathematical definitions of stability applied to mechanical systems, namely stability in the sense of Lagrange (de la Grange 1788) and stability in the sense of Lyapunov (Lyapunov 1992). Both address the issue of the size of the solution, that is, the Euclidean norm of the state deviation from the origin, | x |= √√√√ n∑ i x2i , (B.12) where xi, i = 1, . . . , n are the state variables of the nth-order system. B.2 Stability in the Sense of Lagrange A system described by equation (B.10) is said to be stable about the equilibrium point, xe = 0 in the sense of Lagrange if, for each real and positive number, δ, there exists another real and positive number, , such that | x(0) |≤ δ (B.13) implies that | x(t) |≤ , t ≥ 0. (B.14) In other words, stability in the sense of Lagrange guarantees a bounded solution about the equilibrium point at the origin, resulting from every bounded perturbation from the original initial condition. Example B.1 Consider a simple pendulum consisting of a bob suspended by an inextensible string of length, . The unforced equation of motion describing angular displacement, θ(t), from the vertically downward position of the bob is θ̈ + g  sin θ = 0, (B.15) where g is the acceleration due to gravity (constant). The state equation with the choice of state variables, x1(t) = θ(t) and x2(t) = θ̇(t), is expressed as ẋ = ( x2 − g  sin x1 ) . (B.16) Appendix B: Stability 405 0 10 20 30 40 50 −4 −2 0 2 4 x 1 0 10 20 30 40 50 −10 −5 0 5 10 t x 2 Figure B.3 Initial response of a van der Pol oscillator with a = 2 and initial condition x1(0) = 0, x2(0) = 0.01 Example B.2 Consider a van der Pol oscillator with state equations ẋ1 = x2, ẋ2 = −x1 + ( a − x21 ) x2, (B.20) where a is a real, positive constant. It is clear that the origin, (0, 0), is an equilibrium point of the system. With any small, initial perturbation from the equilibrium point the system departs from the equilibrium, but exhibits a bounded oscillation. An example of this behavior is seen in Figures B.3 and B.4 for a = 2 and initial condition x1(0) = 0, x2(0) = 0.01, computed with the following MATLAB statements: >> [t,x]=ode45(@vanderpol,[0 50],[0 0.01]); >> subplot(211),plot(t,x(:,1)),ylabel(’x_1’),hold on, >> subplot(212),plot(t,x(:,2)),xlabel(’t’),ylabel(’x_2’) >> figure >> plot(x(:,1),x(:,2)),xlabel(’x_1’),ylabel(’x_2’) where vanderpol.m is the following file specifying the state equation of the unforced pendulum: function xdot=vanderpol(t,x) a=2; xdot=[x(2,1); -x(1,1)+(a-x(1,1)ˆ2)*x(2,1)]; Figure B.3 shows that the van der Pol oscillator eventually reaches a constant (but large) amplitude oscillation called a limit cycle. The plot of x2 = ẋ1 vs. x1 in Figure B.4 (called the phase portrait) illustrates the limit cycle behavior at the outer boundary of the closed curve (called the attracting boundary). Since even a small perturbation makes the system depart from the small neighborhood of the equilib- rium point, the origin of the van der Pol oscillator is unstable in the sense of Lyapunov. However, since the initial response to a bounded perturbation is always bounded, the origin of the van der Pol oscillator is said to be stable in the sense of Lagrange. 406 Appendix B: Stability −3 −2 −1 0 1 2 3 −6 −4 −2 0 2 4 6 x 1 x 2 Attracting Boundary Figure B.4 Phase portrait of a van der Pol oscillator with a = 2 and initial condition x1(0) = 0, x2(0) = 0.01 In all our applications in this book, we will regard a stable equilibrium point to be the one that satisfies the stability criterion in the sense of Lyapunov. B.3.1 Asymptotic Stability A system described by equation (B.10) is said to be asymptotically stable about the origin, xe = 0, if, it is stable in the sense of Lyapunov and if for each real and positive number, , however small, there exist real and positive numbers, δ and τ, such that | x(0) |< δ (B.21) implies that | x(t) |< , t > τ. (B.22) Asymptotic stability is thus possessed by a special class of Lyapunov stable systems whose slightly perturbed solutions approach the origin asymptotically for large times. Thus, the equilibrium is eventually restored. Clearly, asymptotic stability guarantees stability in the sense of Lyapunov (but not vice versa). B.3.2 Global Asymptotic Stability A system described by equation (B.10) is said to be globally asymptotically stable about the origin, xe = 0, if it is stable in the sense of Lyapunov and if, for each real and positive pair, (δ, ), there exists a real and positive number, τ, such that | x(0) |< δ (B.23) Appendix B: Stability 407 implies that | x(t) |< , t > τ. (B.24) Global asymptotic stability thus refers to asymptotic stability with respect to all possible initial con- ditions, x(0), such that | x(0) |< δ, and not just a few specific ones. There are further classifications of asymptotic stability, such as uniform asymptotic stability, and asymptotic stability in the large and whole (Loría and Panteley 2006), which we shall not consider here. B.3.3 Lyapunov’s Theorem Let V (x, t) be a continuously differentiable scalar function of the time as well as of the state variables of a system described by equation (B.10), whose equilibrium point is xe = 0. If the following conditions are satisfied: V (0, t) = 0, V (x, t) > 0, dV dt (x) < 0, for all x /= 0, (B.25) || x ||→ ∞ implies V (x, t) → ∞, (B.26) then the origin, xe = 0, is globally asymptotically stable. Proof of Lyapunov’s theorem (Slotine and Li 1991) is obtained from the unbounded, positive definite nature of V (x, t), and negative definite nature of V̇ (x, t), implying that for any initial perturbation from the origin, x(0) /= 0, the resulting solution satisfies V (x(t), t) ≤ V (x(0), t), t > 0 (i.e., remains in a bounded neighborhood of the origin). Furthermore, the same also implies that V (x(t2), t2) ≤ V (x(t1), t1), t2 > t1, which means a convergence of every solution to the origin. Lyapunov’s theorem gives a test of global asymptotic stability without the necessity of solving the system’s state equations for every possible initial condition, and is thus a powerful tool in nonlinear control design. It merely requires finding a suitable Lyapunov function of the state variables. Example B.3 Consider a system described by ẋ1 = −x1 − x2e−t , ẋ2 = x1 − x2. (B.27) Clearly, the system has an equilibrium point at the origin, x1 = x2 = 0. Let us select the Lyapunov function V (x, t) = x21 + x22(1 + e−t). (B.28) We have V (x, t) > 0, for all x /= 0, (B.29) and dV dt = V̇ (x, t) = ∂V ∂t + ∂V ∂x f(x, t) = −x22e−t + 2x1ẋ1 + 2x2ẋ2(1 + e−t) = −2(x21 + x22 − x1x2) − 3x22e−t (B.30) ≤ −2(x21 + x22 − x1x2) = −x21 − x22 − (x1 − x2)2 < 0, for all x /= 0. (B.31) 410 Appendix C: Control of Underactuated Flight Systems one can use the LQR design approach with variable gains, as demonstrated in Chapter 4. However, an adaptive scheduling of the gains is now necessary in order to have a stable closed-loop system for the underactuated plant. Furthermore, the linear feedback law requires a precise thrust modulation that may only be possible in the future. In the following example, we assume that a required thrust magnitude is instantly available to guide the rocket. Example C.1 Repeat the gravity-turn trajectory guidance of Example 4.9 with a closed-loop tracking system for regulating the altitude, speed, and flight-path angle using full-state feedback and forward acceleration control input. A regulator with the time-varying, state feedback gain matrix K(t) is designed by the infinite-time LQR method and the closed-loop eigenvalues computed with constant cost parameter matrices, Q = I, R = 1, using the following MATLABr statements: n=size(t,2); mu=398600.4; r0=6378.14; g0=mu/r0ˆ2; Eg=[];K=[]; for i=1:n; h=x(1,i)-r0; r=r0+h; v=x(2,i); phi=x(3,i); A=[0 sin(phi) v*cos(phi); 2*g0*r0ˆ2*sin(phi)/rˆ3 0 -g0*r0ˆ2*cos(phi)/rˆ2; (-v/rˆ2+2*g0*r0ˆ2/(v*rˆ3))*cos(phi) ... (1/r+g0*r0ˆ2/(vˆ2*rˆ2))*cos(phi) ... -(v/r-g0*r0ˆ2/(v*rˆ2))*sin(phi)]; B=[0 1 0]’; [k,S,E]=lqr(A,B,eye(3),1); E=sort(E); K(i,:)=k; Eg(:,i)=E; end The resulting closed-loop eigenvalues are plotted in Figure C.1, indicating an asymptotically stable closed-loop system with coefficients frozen at all times along the nominal trajectory. The variation of the LQR regulator gains, kh(t), kv(t), kφ(t), is plotted in Figure C.2. While the assumption kh = 1 s−2 can be readily made as a possible adaptation law, the other two gains require a scheduling with the increasing flight speed according to the following simple adaptation laws (Figure C.2) in the given control interval: kv(t) = 1.5v−1/6, kφ(t) = 0.5v3.01. In order to verify how accurately the adaptation laws capture the LQR designed dynamics, we compare the closed-loop eigenvalues of the two systems in Figure C.3. Clearly, the adapted system poles have marginally decreased real-part magnitudes during the initial phase of flight, when compared to the system with LQR gains. However, as the time progresses, the difference between the two systems decreases, and near the end of the trajectory, the adapted system eigenvalues have a slightly larger real-part magnitude. In order to check the closed-loop system’s response with adapted gains, we consider an initial response to a speed disturbance, v = 0.1 km/s, applied at t = 54.5 s (Figure C.4), with the coefficients frozen at their respective values also at t = 54.5 s. Both the control systems are seen to mitigate the disturbance Appendix C: Control of Underactuated Flight Systems 411 Figure C.1 Closed-loop guidance system eigenvalues along the nominal, gravity-turn trajectory of a launch vehicle with forward acceleration input completely in about 10 s. While Figure C.4 is not an accurate picture of how the time-varying system would behave, it indicates a close agreement between the adapted and the LQR design systems, thereby validating the choice of the adaptation laws. To simulate the actual nonlinear time-varying system with adaptive feedback gains we construct a Simulink model as shown in Figure C.5, with the previously programmed nonlinear subsystem block, Rocket Translational Dynamics (Chapter 4). Note the adaptation loops in Figure C.5 specified with Figure C.2 Feedback gains for LQR and adaptive guidance systems along the nominal, gravity-turn trajectory of a launch vehicle with forward acceleration input 412 Appendix C: Control of Underactuated Flight Systems Figure C.3 Closed-loop eigenvalues comparison of LQR and adapted guidance systems along the nominal, gravity- turn trajectory of a launch vehicle with forward acceleration input user-defined mathematical function blocks. The results of the simulation – carried out with a fourth-order Runge–Kutta scheme – are plotted in Figures C.6–C.9. The initial error of 0.1 rad (5.7◦) in the flight-path angle – and the attendant undershoot in flight speed – is seen to be attenuated in about 400 s (Figure C.6), compared to about 50 s for the normal acceleration guidance case of Chapter 4, without affecting the Figure C.4 Closed-loop initial response comparison of LQR and adapted guidance systems with frozen parameters and forward acceleration input Appendix C: Control of Underactuated Flight Systems 415 Figure C.9 Closed-loop simulated speed with initial flight-path perturbation with nonlinear adaptive rocket guidance with forward acceleration input C.2 Thrust Saturation and Rate Limits (Increased Underactuation) Modern liquid propellant engines can be throttled up to some extent by regulating the supply of propellants through valves. However, a propellant flow rate below a certain limit cannot sustain combustion, resulting in flame-out (zero thrust). Even when throttling is possible, the engine response is highly nonlinear due to limits on the rate of change of thrust. Thus the rocket engine can be regarded as a nonlinear actuator with output saturation and rate limits, and care must be taken to account for the engine behavior when simulating a closed-loop guidance and control system involving thrust modulation. Due to difficulty in throttling a large rocket engine, most launch vehicles employ small liquid propellant rockets called vernier engines dedicated to the task of guidance and attitude control. The vernier rockets – having a faster throttle response – produce axial and normal acceleration inputs with some alacrity required for controlling the flight path and attitude, while the main engines operate at nearly full power until burn-out. However, this does not allow precise control of the flight speed. Having magnitude and rate limits on the thrust results in a thrust-controlled vehicle becoming further underactuated due to the actuator restrictions not allowing the full level of control to be exercised, as demanded by a state feedback regulator. Furthermore, now the system’s underactuated behavior is non- linear due to the presence of actuation limits. While there are certain procedures for designing nonlinear, underactuated control systems (such as the describing function method (Gibson 1963), which applies a frequency domain approach through elaborate quasi-linear gains for each class of nonlinearity), we will refrain from such a design. Instead, we will test our time-varying adaptive control design for robustness in the presence of magnitude and rate limits. Example C.2 Taking into account the following saturation and rate limits of the forward acceleration provided by the rocket engine, simulate the closed-loop response of the gravity-turn guidance system designed in Example C.1: 6 ≤ u1 ≤ 40 m/s2, ∣∣∣∣du1dt ∣∣∣∣ ≤ 30 m/s3. 416 Appendix C: Control of Underactuated Flight Systems Figure C.10 Simulink block diagram for simulating the nonlinear, closed-loop, adaptive rocket guidance system with thrust saturation and rate limits in a gravity-turn launch with initial perturbation We insert the nonlinear saturation and rate limit blocks as constituting the engine actuator into the Simulink block diagram of Figure C.5, as shown in Figure C.10, and repeat the simulation of Example C.1. The results of the simulation are plotted in Figures C.11 and C.12. As a consequence of the full level of control not being exercised, there is a terminal error (excess speed of 0.01km/s and excess altitude of 0.057km) due to the inability of the engine to throttle below u1 = 6 m/s2 as required by the guidance system designed in Example C.1. As a result, the vehicle is sent into a slightly higher orbit than necessary. However, the design can be considered conservative as the mission requirements are slightly exceeded, Figure C.11 Closed-loop trajectory error with an adaptive rocket guidance system simulated with thrust saturation and rate limits Appendix C: Control of Underactuated Flight Systems 417 Figure C.12 Closed-loop forward acceleration input with an adaptive rocket guidance system simulated with thrust saturation and rate limits and the desired orbit can be subsequently attained by minor orbital corrections (Chapter 6). When the atmospheric drag is taken into account, the final error can be seen to be much smaller than that observed in Figure C.11. It is interesting to compare the closed-loop forward acceleration inputs plotted in Figures C.12 and C.8. Clearly, the curtailment of both the magnitude and the rate of actuation results in the engine switching between its upper and lower limits in order to drive the system toward a smaller error. Such a behavior of nonlinear underactuated, asymptotically stable closed-loop systems is called switching. It is also as if the control system were sliding along the nominal trajectory by alternately increasing and decreasing the input magnitude between the actuator limits. Such a nonlinear control strategy is called sliding mode control (Slotine and Li 1991). Although we have not designed a sliding mode controller deliberately, it arises automatically as a consequence of putting the thrust limitations in the closed-loop regulator. The ability of the nonlinear, time-varying regulator to exercise control even in the presence of substantial magnitude and rate limitations attests to its robustness. C.3 Single- and Bi-output Observers with Forward Acceleration Input We now extend the state feedback controller to an observer-based guidance of the underactuated rocket. Let us first consider the single output to be the speed error, y = v, with the rearranged state vector, x = (x1, x2), where x1 = v, x2 = (h, φ), and the following state-space coefficient matrices: A = ⎛ ⎝ 0 | a31 a34 a13 | 0 a14 a43 | a41 a44 ⎞ ⎠ , B = ⎛ ⎝ 1 0 0 ⎞ ⎠ . (C.3) The output equation, y = x1, and the partitioning of A, B produce the following sub-matrices required for the reduced-order observer: C = 1, A11 = 0, A12 = (a31 a34), B1 = 1, (C.4) 420 Appendix C: Control of Underactuated Flight Systems Table C.1 % Program ’rocket_obs2.m’ for generating the closed-loop % state equations for a launch vehicle to circular orbit % while guided by a forward acceleration controlled, % speed output based reduced-order compensator. % (c) 2010 Ashish Tewari % To be used by Runge-Kutta solver ’ode45.m’ % t: current time (s); x(1,1): current radius; % x(2,1): current speed % x(3,1): current flight-path angle; % x(4:5,1): reduced-order observer state vector. function xdot=rocket_obs2(t,x) global tn; %Nominal time listed in the calling program. global rn; %Nominal radius listed in the calling program. global vn; %Nominal speed listed in the calling program. global phin;%Nominal flt-path angle listed in the calling program. global un; %Nominal control input listed in the calling program. mu=398600.4;%Earth’s gravitational constant (kmˆ3/sˆ2) zeta=1/sqrt(2);w=.001/zeta; %Observer damp.-ratio & nat. frequency. r=interp1(tn,rn,t); %Interpolation for current nominal radius. v=interp1(tn,vn,t); %Interpolation for current nominal speed. phi=interp1(tn,phin,t); %Interpolation current nominal f.p.a. u=interp1(tn,un,t); %Interpolation for current nominal control input. % Linearized plant’s state coefficient matrix: A=[0 sin(phi) v*cos(phi); 2*mu*sin(phi)/rˆ3 0 -mu*cos(phi)/rˆ2; (-v/rˆ2+2*mu/(v*rˆ3))*cos(phi) ... (1/r+mu/(vˆ2*rˆ2))*cos(phi) -(v/r-mu/(v*rˆ2))*sin(phi)]; % Regulator gains by adaptation law: k=[1.39959e-2 1.627*vˆ(-1/6) .65367*vˆ3.01]; % Observer coefficients follow: a1=A(1,2);a2=A(1,3);a6=A(2,1);a7=A(2,3);a8=A(3,1); a9=A(3,2);a10=A(3,3); L1=(a2*a8+wˆ2-a2*a6*(a10+2*zeta*w)/a7)/(a7*a8-a6*a10-a2*a6ˆ2/a7); L2=(2*zeta*w+a10-L1*a6)/a7; L=[L1;L2]; F=[0 a2;a8 a10]-L*[a6 a7]; H=-L; G=F*L+[a1;a9]; % Trajectory error variables: dr=x(1,1)-r; dv=x(2,1)-v; dphi=x(3,1)-phi; % Estimated trajectory error variables: z=x(4:5,1); xo2=L*dv+z; % Feedback control input: du=-k(1,2)*dv-k(1,1)*xo2(1,1)-k(1,3)*xo2(2,1); if t>600 u=0; Appendix C: Control of Underactuated Flight Systems 421 Table C.1 (Continued) else u=u+du; end % Trajectory state equations: xdot(1,1)=x(2,1)*sin(x(3,1)); xdot(2,1)=u-mu*sin(x(3,1))/x(1,1)ˆ2; xdot(3,1)=(x(2,1)/x(1,1)-mu/(x(1,1)ˆ2*x(2,1)))*cos(x(3,1)); % Observer’s state equation: xdot(4:5,1)=F*z+H*du+G*dv; Table C.2 % Program ’rocket_obs_u2.m’ for generating the closed-loop % control input for a launch vehicle to circular orbit guided by % forward acceleration controlled, speed output based reduced-order % compensator, using post-processing of simulated trajectory (T,x). % T: (nx1) vector containing simulation ’n’ time instants. % x: (nx5) vector containing trajectory and observer state. function u=rocket_obs_u2(T,x) global tn; global rn; global vn; global phin; global un; mu=398600.4; zeta=1/sqrt(2);w=.001/zeta; n=size(T,1); for i=1:n t=T(i,1); r=interp1(tn,rn,t); v=interp1(tn,vn,t); phi=interp1(tn,phin,t); U=interp1(tn,un,t); A=[0 sin(phi) v*cos(phi); 2*mu*sin(phi)/rˆ3 0 -mu*cos(phi)/rˆ2; (-v/rˆ2+2*mu/(v*rˆ3))*cos(phi) ... (1/r+mu/(vˆ2*rˆ2))*cos(phi) -(v/r-mu/(v*rˆ2))*sin(phi)]; a1=A(1,2);a2=A(1,3);a6=A(2,1);a7=A(2,3);a8=A(3,1); a9=A(3,2);a10=A(3,3); L1=(a2*a8+wˆ2-a2*a6*(a10+2*zeta*w)/a7)/(a7*a8-a6*a10-a2*a6ˆ2/a7); L2=(2*zeta*w+a10-L1*a6)/a7; L=[L1;L2]; F=[0 a2;a8 a10]-L*[a6 a7]; H=-L; G=F*L+[a1;a9]; dr=x(i,1)-r; dv=x(i,2)-v; dphi=x(i,3)-phi; z=x(i,4:5); xo2=L’*dv+z; (Continued) 422 Appendix C: Control of Underactuated Flight Systems Table C.2 (Continued) if t>600 u(i,1)=0; else k=[1.39959e-2 1.627*vˆ(-1/6) .65367*vˆ3.01]; du=-k(1,2)*dv-k(1,1)*xo2(1,1)-k(1,3)*xo2(1,2); u(i,1)=U+du; end end Figure C.13 Closed-loop compensated trajectory during launch phase of a guided launch vehicle with a speed output based reduced-order observer and adaptive regulator y = x1 = (h, v)T , x2 = φ, and the following partitioned coefficient matrices: A = ⎛ ⎝ 0 a13 | a14 a31 0 | a34 a41 a43 | a44 ⎞ ⎠ , B = ⎛ ⎝ 0 1 0 ⎞ ⎠ . (C.12) The sub-matrices required for the reduced-order observer are thus the following: C = I, A11 = ( 0 a13 a31 0 ) , A12 = ( a14 a34 ) , B1 = ( 0 1 ) , (C.13) A21 = (a41 a43), A22 = a44, B2 = 0. (C.14) The state equation of the single-order observer with estimated state xo2 = Ly + z and gain matrix L = (L1, L2) is ż = Fz + Hu + Gy, (C.15) where F = A22 − LCA12 = a44 − L1a14 − L2a34, H = B2 − LCB1 = −L2, (C.16) Appendix C: Control of Underactuated Flight Systems 425 While we can directly utilize the adaptive regulator gains of Example C.1, it is better to fine-tune the regulator in order to avoid the usual degradation of performance with observer-based feedback. Hence, a better choice of the regulator gains for use in the reduced-order, bi-output compensator is kh(t) = 1, kv(t) = v−1/2, kφ(t) = 1.0 × 10−5. For closing the loop of the compensated system, we have u1 = khh + kvv + kh + kφxo2. The simulation of the gravity-turn guidance system with the bi-output observer is carried out by a fourth-order Runge–Kutta method using a modification of the MATLAB code, rocket obs2.m (Table C.1), listed in Table C.3. A separate code (not listed here for brevity) similar to rocket obs u2.m (Table C.2) computes the closed-loop control input for the bi-output case. The results of the simulation are plotted in Figures C.17 and C.18, showing a successful launch even in the presence of a large initial angular error while requiring a control input slightly smaller in magnitude compared to the state feedback regulator of Example C.1. The best feature of the bi-output compensator is its greater robustness with respect to variation in controller gains, kh, kv, kφ, L2, enabling it to have a practical application. Table C.3 % Program ’rocket_obs.m’ for generating the closed-loop % state equations for a launch vehicle to circular orbit while % guided by a forward acceleration controlled, bi-output % reduced-order compensator based upon measurement and feedback of % speed and altitude % (c) 2010 Ashish Tewari % To be used by Runge-Kutta solver ’ode45.m’ % t: current time (s); x(1,1): current radius; % x(2,1): current speed % x(3,1): current flight-path angle; % x(4,1): reduced-order observer state vector. % function xdot=rocket_obs(t,x) global tn; %Nominal time listed in the calling program. global rn; %Nominal radius listed in the calling program. global vn; %Nominal speed listed in the calling program. global phin; %Nominal f.p.a. listed in the calling program. global un; %Nominal control listed in the calling program. mu=398600.4; %Earth’s gravitational constant (kmˆ3/sˆ2) r=interp1(tn,rn,t); %Interpolation for current nominal radius. v=interp1(tn,vn,t); %Interpolation for current nominal speed. phi=interp1(tn,phin,t); %Interpolation for current nominal f.p.a. u=interp1(tn,un,t); %Interpolation for current nominal control. % Linearized plant’s state coefficient matrix: A=[0 sin(phi) v*cos(phi); 2*mu*sin(phi)/rˆ3 0 -mu*cos(phi)/rˆ2; (-v/rˆ2+2*mu/(v*rˆ3))*cos(phi) ... (1/r+mu/(vˆ2*rˆ2))*cos(phi) -(v/r-mu/(v*rˆ2))*sin(phi)]; % Regulator gains by adaptation law: k=[1 vˆ(-1/2) 1e-5]; 426 Appendix C: Control of Underactuated Flight Systems Table C.3 (Continued) % Observer coefficients follow: a1=A(1,2);a2=A(1,3);a6=A(2,1);a7=A(2,3);a8=A(3,1); a9=A(3,2);a10=A(3,3); F=-0.01; L2=(a10-F)/a7; G=F*[0 L2]+[a8 a9]-[L2*a6 0]; H=-L2; % Trajectory error variables: dr=x(1,1)-r; dv=x(2,1)-v; dphi=x(3,1)-phi; % Estimated trajectory error variables: z=x(4,1); xo2=L2*dv+z; % Feedback control input: du=-k(1,2)*dv-k(1,1)*dr-k(1,3)*xo2; if t>600 u=0; else u=u+du; end % Trajectory state equations: xdot(1,1)=x(2,1)*sin(x(3,1)); xdot(2,1)=u-mu*sin(x(3,1))/x(1,1)ˆ2; xdot(3,1)=(x(2,1)/x(1,1)-mu/(x(1,1)ˆ2*x(2,1)))*cos(x(3,1)); % Observer’s state equation: xdot(4,1)=F*z+H*du+G*[dr dv]’; Figure C.17 Closed-loop compensated trajectory of a guided launch vehicle with a bi-output, reduced-order, observer-based adaptive compensator Appendix C: Control of Underactuated Flight Systems 427 Figure C.18 Closed-loop forward acceleration input of a guided launch vehicle with a bi-output, reduced-order, observer-based adaptive compensator Example C.5 Compute the drag, mass, and thrust required for the observer-based gravity-turn trajectory guidance system with forward acceleration input of Example C.4 using the two-stage rocket’s drag and mass data given in Chapter 4. The results of the code rocket drag mass.m run for the closed-loop trajectory simulated in Example C.4 are plotted in Figures C.19 and C.20. Note that the thrust required is always several times greater than the Figure C.19 Mass, atmospheric drag, and required thrust of a guided launch vehicle with a bi-output, reduced-order, observer-based adaptive compensator 430 Appendix C: Control of Underactuated Flight Systems Figure C.23 Simulink Observer Coefficients subsystem block for computing adaptive observer gains required by the block diagram of Figure C.21 Figure C.24 Simulink Adaptive Observer subsystem block for adaptive observer dynamics required by the block diagram of Figure C.21 Appendix C: Control of Underactuated Flight Systems 431 Figure C.25 Closed-loop simulated trajectory of nonlinear adaptive rocket guidance system with adaptive observer and magnitude and rate limits on thrust acceleration input the nominal trajectory is defined in the MATLAB workspace by row vectors tn, rn, vn, phin, and the nominal control input by the matrix un1 = [tn; u1]. Furthermore, the mass and drag computed by the code rocket drag mass.m (Chapter 4) along the nominal trajectory are stored as column vectors m and D, respectively. Note the presence of table look-up blocks for nominal mass and drag as well as thrust acceleration saturation and rate limit blocks in Figure C.21. The results of the simulation carried out by the fourth-order Runge–Kutta method with a relative tolerance of 10−6 are plotted in Figures C.25 and C.26. 0 100 200 300 400 500 600 0 10 20 30 40 50 u 1 (m /s 2 ) 0 100 200 300 400 500 600 0 5 10 15 Time (s) f T (× 10 5 N ) Perturbed Nominal Figure C.26 Closed-loop simulated forward acceleration input and required thrust of nonlinear adaptive rocket guidance system with adaptive observer and magnitude and rate limits on thrust acceleration input 432 Appendix C: Control of Underactuated Flight Systems There is no difference in the closed-loop trajectory compared to that observed in Figure C.6 for the case of state feedback regulator (Example C.1), even though an observer and rate and magnitude control limits are being implemented. Figure C.26 shows that thrust saturation occurs only during the initial phase of the flight due to the reduced feedback gains compared to those in Example C.2. Furthermore, the saturation limits cause a massive reduction (almost 50%) in the required peak thrust from 2.3 × 106 N (Figure C.19) to only 1.25 × 106 N (Figure C.26). Clearly, the bi-output, observer-based adaptive compensator is much more efficient than the adaptive full-state feedback regulator for the present example. References Gibson, J.E. (1963) Nonlinear Automatic Control. McGraw-Hill, New York. Slotine, J.-J.E. and Li, W. (1991) Applied Nonlinear Control. Prentice Hall, Englewood Cliffs, NJ. Tewari, A. (2002) Modern Control Design with MATLAB and Simulink. John Wiley & Sons, Ltd, Chichester. Index 435 loop-transfer recovery, 85 low-pass filter, 77 Lyapunov function, 93 Lyapunov’s theorem, 407 Mach number, 137, 168 Markov process, 73 matrix exponential, 9, 394 matrix Riccati equation, 48, 66, 80 maximum overshoot, 6 mean anomaly, 302 mean motion, 301 mean position, 311 mean value, 70 mean value theorem, 311 measurement noise, 3 minimum principle, 45 miss distance, 16, 206 modified Rodrigues parameters, 185 navigation, 11 necessary condition, 31 neighboring optimal solution, 46 Newton’s method, 315 nominal trajectory, 3 normal distribution, 71 Nyquist frequency, 25 objective function, 31 observability, 395 observable, 4 observer, 4, 78 optimal point, 32 optimal return function, 45 orbit equation, 299 orbital mechanics, 297 orbital perturbation, 334 orbital plane, 298 orbits, 298 order, 2, 3 output, 2 output coefficient matrix, 9, 393 output equation, 3 outside air temperature, 168 parabola, 299 periapsis, 299 perifocal frame, 299 phase portrait, 405 PID tuning, 21 plant, 2 polar coordinates, 298, 299 position error, 13 positive definiteness, 32 power spectral density matrix, 75 probability, 69 probability density function, 69 probability theory, 69 process noise, 2, 3 projected interception point, 13 proportional navigation guidance, 5 proportional-derivative, 23 proportional-integral-derivative control, 21 pure rolling mode, 177 quaternion, 185, 188, 366 random, 69 rectilinear motion, 298, 299 regulator, 4 retrograde, 315 return ratio at input, 85 return ratio at output, 85 Reynolds number, 137 right ascension of the ascending node, 300 Robust Control Toolbox, 86 robust design, 9 robustness, 69 rocket, 10 roll-off, 77, 87 saddle point, 32 sample size, 69 sampling, 24 sampling interval, 73 semi-major axis, 299 semi-perimeter, 308 sensitivity matrix, 84 sensor, 2 separation principle, 9 servo, 145 sideslip angle, 175 singular control, 62 singular point, 32 singular value, 84, 396 singular value decomposition, 396 six-degrees-of-freedom motion, 185 sliding mode control, 417 small perturbation, 8 spacecraft, 10 436 Index specific angular momentum, 298 specific impulse, 207 spiral mode, 177 stability, 7, 395 stability augmentation system, 177 stability in the sense of Lyapunov, 7 stabilizability, 82, 395 standard deviation, 70 state, 2 state equation, 3 state space, 2 state transition matrix, 8, 305, 338 stationary point, 32 stationary process, 72 steady wind, 110 steady-state response, 6 step response, 394 stochastic, 2, 69 strictly proper, 9, 83, 393 Strouhal numbers, 137 Stumpff function, 313, 314 sufficient condition, 32 Sun direction, 382 switching, 417 switching condition, 63 switching function, 63 symplectic matrix, 67, 305, 338, 355 system, 2 Taylor series, 302 terminal control, 4, 29 terminal guidance, 16 terminal state, 4, 29 three-dimensional proportional navigation, 199 thrust, 10 time of periapsis, 301 time step, 24 time to go, 14 time-shift, 74 tolerance, 302, 315 trace, 91 tracking control, 4, 29 transition matrix, 65 transitional probability, 73 transmission matrix, 84 true anomaly, 299 two-point boundary value problem, 4, 29, 313 underactuated, 409 universal gravitational constant, 297 universal variable, 313 unstable, 395 van der Pol oscillator, 404 variance, 70 velocity error, 13 velocity to be gained, 13 vernal equinox, 300 white noise, 9, 73 winds, 110 winds-aloft, 110 z-transform, 25 zero-mean Gaussian white noise, 76 2 Advanced Control of Aircraft, Spacecraft and Rockets Table 1.1 Control system variables Symbol Variable Dimension u(t) control input vector m × 1 û(t) optimal control input vector m × 1 w(t) measurement noise vector l × 1 x(t) state vector n × 1 x̂(t) optimal state vector n × 1 y(t) output vector l × 1 z(t) state vector for augmentation q × 1 ν(t) process noise vector p × 1 The generic term planet is used for any celestial body about which the flight is referenced (Earth, Moon, Sun, Jupiter, etc.). The Euclidean (or L2) norm of a vector, a = (ax, ay, az)T , is written as || a ||= √ a2x + a2y + a2z . (1.5) Standard aerospace symbols define relevant flight parameters and variables as and when used. The nomenclature for control system variables is given in Table 1.1. Any departure from this labeling scheme, if necessary, will be noted. Control is the name given to the general task of achieving a desired result by appropriate adjustments (or manipulations). The object to be controlled (a flight vehicle) is referred to as the plant, while the process that exercises the control is called the controller. Both the plant and the controller are systems, defined as self-contained sets of physical processes under study. A system has variables applied to it externally, called the input vector, and produces certain variables internally, called the output vector, which can be measured. In modeling a system, one must account for the relationship between the input and output vectors. This relationship generally takes the form of a set of differential and algebraic equations, if the system is governed by known physical laws. A system having known physical laws is said to be deterministic, whereas a system with unknown (or partially known) physical laws is called non-deterministic or stochastic. Every system has certain unwanted external input variables – called disturbance inputs – that cannot be modeled physically and are thus treated as stochastic disturbances. The disturbances are generally of two types: (i) process noise that can arise either due to unwanted external inputs, or internally due to uncertainty in modeling the system, and (ii) measurement noise that results from uncertainty in measuring the output vector. The presence of these external and internal imperfections renders all practical systems stochastic. The condition, or state, of a system at a given time is specified by a set of scalar variables, called state variables, or, in vector form, the state vector. The vector space spanned by the state vector is called a state space. The state of a system is defined as a collection of the smallest number of variables necessary to completely specify the system’s evolution in time, in absence of external inputs. The number of state variables required to represent a system is called order of the system, because it is equal to the net order of differential equations governing the system. While the size of the state space (i.e., the order of the system) is unique, any given system can be described by infinitely many alternative state-space representations. For example, a flight vehicle’s state can be described by the position, r(t), velocity, v(t), angular velocity, ω(t), and orientation, ξ(t), relative to a frame of reference. Thus, the state vector of a flight vehicle’s motion is x(t) = {r(t), v(t), ω(t), ξ(t)}T . However, x(t) can be transformed into any number of different state vectors depending upon the choice of the reference frame. A system consisting of the plant and the controller is called a control system. The controller manipulates the plant through a control input vector, which is actually an input vector to the plant but an output of the controller. In physical terms, this output can take the form of either a force or a torque (or both) applied Introduction 5 − Actuatorsd )(tx o(t)x Plant Sensors Flight Dynamics Observer y(t)u(t)e(t) Feedforward Controller Regulator Figure 1.1 Tracking control system with an observer RT(t), respectively, with respective velocities V(t) and VT(t) relative to a stationary frame of reference, (SXYZ). The instantaneous position of the target relative to the missile is given by (the vector triangle Soo′ in Figure 1.2) r(t) = RT(t) − R(t), (1.11) while the target’s relative velocity is obtained by differentiation as follows: v(t) = dr dt = dRT dt − dR dt = VT(t) − V(t). (1.12) Without considering the equations of motion of the missile and the target (to be derived in Chapter 4), we propose the following control law for missile guidance: V(t) = K(t) [RT(t) − R(t)] = K(t)r(t), (1.13) where K(t) is a time-varying gain matrix. A linear feedback control law of the form given by equation (1.13) is called a proportional navigation guidance law (PNG), whose time derivative gives the required Y X Z R(t) RT(t) VT(t) V(t) r(t) Target Missile S o o' Figure 1.2 Missile guidance for interception of an aerial target 6 Advanced Control of Aircraft, Spacecraft and Rockets acceleration control input, u(t), to be applied to the missile: u(t) = dV dt = K(t) [VT(t) − V(t)] + dK dt r (1.14) = Kv(t) + K̇r(t). For a successful interception of the target, the relative separation, r, must vanish at some time, T = tf − t, without any regard to the relative velocity, v, prevailing at the time. A likely choice of state vector for the interception problem is x(t) = [r(t), v(t)]T , (1.15) which yields the following linear feedback control law of the tracking system: u(t) = [K̇(t), K(t)] x(t). (1.16) The main advantage of the control law given by equation (1.14) is the linear relationship it provides between the required input, u, and the measured outputs, (r, v), even though the actual plant may have a nonlinear character. Thus the PNG control law is quite simple to implement, and nearly all practical air-to-air missiles are guided by PNG control laws. As the missile is usually rocket powered, its thrust during the engagement is nearly constant. In such a case, PNG largely involves a rotation of the missile’s velocity vector through linear feedback relationships between the required normal acceleration com- ponents (that are generated by moving aerodynamic fins and/or by thrust vectoring) and the measured relative coordinates and velocity components of the target (obtained by a radar or an infrared sensor mounted on the missile). The proportional navigation gain matrix, K(t), must be chosen such that the interception [r(tf ) → 0] for the largest likely initial relative distance, || r(0) ||, takes place within the al- lowable maximum acceleration, || u ||≤ um as well as the maximum time of operation, tf , of the engines powering the missile. We shall have occasion to discuss the PNG law a little later. A tracking system with a time-varying reference state, xd(t), can be termed successful only if it can maintain the plant’s state, x(t), within a specified percentage error, || xd(t) − x(t) ||≤ δ, of the desired reference state at all times. The achieved error tolerance (or corridor about the reference state), δ, thus gives a measure of the control system’s performance. The control system’s performance is additionally judged by the time taken by the plant’s state to reach the desired error tolerance about the reference state, as well as the magnitude of the control inputs required in the process. The behavior of the closed-loop system is divided into the response at large times, t → ∞, called the steady-state response, and that at small values of time when large deviations (called overshoots) from the desired state could occur. A successful control system is one in which the maximum overshoot is small, and the time taken to reach within a small percentage of the desired state is also reasonably small. Example 1.2 Consider a third-order tracking system with the state vector x = (x1, x2, x3)T . A plot of the nominal trajectory, xd(t) is shown in Figure 1.3. The tracking error corridor is defined by the Euclidean norm of the off-nominal state deviation as follows: || xd(t) − x(t) ||= √ (x1 − x1d)2 + (x2 − x2d)2 + (x3 − x3d)2 ≤ δ, (1.17) where δ is the allowable error tolerance. The actual trajectory is depicted in Figure 1.3 and the maximum overshoot from the nominal is also indicated. Introduction 7 Nominal Trajectory x2 x3 2δ x1 Error Tolerance Corridor Actual Trajectory Maximum Overshoot xd(t)x(t) Figure 1.3 Nominal trajectory and tracking error corridor with tolerance, δ 1.2.1 Linear Tracking Systems With the definition for a successful control system given above, one can usually approximate a nonlinear tracking system by linear differential equations resulting from the assumption of small deviations from the reference state. A first-order Taylor series expansion of the control system’s governing (nonlinear) differential equations about the reference state thus yields a linear tracking system (Appendix A), and the given reference solution, xd(t), is regarded as the nominal state of the resulting linear system. A great simplification occurs by making such an assumption, because we can apply the principle of linear superposition to a linearized system, in order to yield the total output due to a linear combination of several input vectors. Linear superposition also enables us to utilize operational calculus (such as Laplace and Fourier transforms) and linear algebraic methods for design and analysis of control systems. Appendix A briefly presents the linear systems theory, which can be found in detail in any textbook on linear systems, such as Kailath (1980). Let the control system without disturbance variables be described by the state equation ξ̇ = f[ξ(t), η(t), t], (1.18) where ξ is the state vector, and η, the input vector. The nonlinear vector functional, f(.), is assumed to possess partial derivatives with respect to state and input variables in the neighborhood of the reference, nominal trajectory, ξ0(t), which is a solution to equation (1.18) and thus satisfies ξ̇0(t) = f[ξ0(t), η0(t), t], ti ≤ t ≤ tf , (1.19) where η0(t) is the known input (called the nominal input) applied to the system in the interval (ti ≤ t ≤ tf ). In order to maintain the system’s state close to a given reference trajectory, the tracking system must possess a special property, namely stability about the nominal reference trajectory. While stability can be defined in various ways, for our purposes we will consider stability in the sense of Lyapunov (Appendix B), which essentially implies that a small control perturbation from the nominal input results in only a small deviation from the nominal trajectory. 10 Advanced Control of Aircraft, Spacecraft and Rockets K LTI Plant LTI Regulator z• = Fz + Gy + Hu xo = z + Ly LTI Observer xo(t) x• = Ax + Bu y = Cx + Du − y(t)u(t) Figure 1.4 Linear time-invariant tracking system The regulator gain matrix, K, the observer gain matrix, L, and the constant coefficient matrices of the observer, F, G, H, must satisfy the asymptotic stability criteria for the overall closed-loop system, which translates into separately guaranteeing asymptotic stability of the regulator and observer by the separation principle (Kailath 1980). Furthermore, in order to be practical, the design of both the regulator and the observer must be robust with respect to modeling uncertainties (process noise) and random external disturbances (measurement noise), both of which are often small in magnitude (thereby not causing a departure from the linear assumption), but occur in a wide spectrum of frequencies. A block diagram of an LTI tracking system with observer is shown in Figure 1.4. The design of such systems for achieving stability and robustness with respect to uncertainties and random disturbances is the topic of several control systems textbooks such as Kwakernaak and Sivan (1972), Maciejowski (1989), Skogestad and Postlethwaite (1996), and Tewari (2002). Examples of the LTI design methods for achieving robustness are linear quadratic Gaussian (LQG) and H∞ techniques that are now standard in modern control design and involve the assumption of a zero-mean Gaussian white noise1 process and measurement noise. A linear observer ensuring the maximum robustness in the state estimation process with respect to a zero-mean Gaussian white noise and measurement noise is called a Kalman filter, or linear quadratic estimator, and is a part of the LQG controller. Since LQG control, H∞ control, and the Kalman filter naturally arise out of the optimal control theory, we shall have occasion to consider them in Chapter 2. 1.3 Guidance and Control of Flight Vehicles Vehicles capable of sustained motion through air or space are termed flight vehicles, and are broadly classified as aircraft, spacecraft, and rockets. Aircraft flight is restricted to the atmosphere, and spacecraft flight to exo-atmospheric space, whereas rockets are equally capable of operating both inside and outside the sensible atmosphere. Of these, aircraft have the lowest operating speeds due to the restriction imposed by aerodynamic drag (an atmospheric force opposing the motion and proportional to the square of the speed). Furthermore, the requirement of deriving lift (a force normal to the flight direction necessary for balancing the weight) and thrust (the force along flight direction to counter drag) from the denser regions of the atmosphere restricts aircraft flight to the smallest altitudes of all flight vehicles. In contrast, spacecraft have the highest speeds and altitudes due to the requirement of generating centripetal acceleration entirely 1 White noise is a statistical process with a flat power spectrum, that is, the signal has the same power at any given frequency. A Gaussian process has a normal (bell-shaped) probability density distribution. Introduction 11 by gravity, while the operating speeds and altitudes of rockets (mostly employed for launching the spacecraft) fall somewhere between those of the aircraft and spacecraft. All flight vehicles require manipulation (i.e., adjustment or control) of position, velocity, and attitude (or orientation) for successful and efficient flight. A transport aircraft navigating between points A and B on the Earth’s surface must follow a flight path that ensures the smallest possible fuel consumption in the presence of winds. A fighter aircraft has to maneuver in a way such that the normal acceleration is maximized while maintaining the total energy and without exceeding the structural load limits. A spacecraft launch rocket must achieve the necessary orbital velocity while maintaining a particular plane of flight. A missile rocket has to track a maneuvering target such that an intercept is achieved before running out of propellant. An atmospheric entry vehicle must land at a particular point with a specific terminal energy without exceeding the aero-thermal load limits. In all of these cases, precise control of the vehicle’s attitude is required at all times since the aerodynamic forces governing an atmospheric trajectory are very sensitive to the body’s orientation relative to the flight direction. Furthermore, in some cases, attitude control alone is crucial for the mission’s success. For example, a tumbling (or oscillating) satellite is useless as an observation or communications platform, even though it may be in the desired orbit. Similarly, a fighter (or bomber) aircraft requires a stable attitude for accurate weapons delivery. Flight vehicle control can be achieved by either a pilot, an automatic controller (or autopilot), or both acting in tandem. The process of controlling a flight vehicle is called flying. A manned vehicle is flown either manually by the pilot, or automatically by the autopilot that is programmed by the pilot to carry out a required task. Unmanned vehicles can be flown either remotely or by onboard autopilots that receive occasional instructions from the ground. It is thus usual to have some form of human intervention in flying, and rarely do we have a fully automated (or autonomous) flight vehicle that selects the task to be performed for a given mission and also performs it. Therefore, the job of designing an automatic control system (or autopilot) is quite simple and consists of: (i) generating the required flight tasks for a given mission that are then stored in the autopilot memory as computer programs or specific data points; and (ii) putting in place a mechanism that closely performs the required (or reference) flight tasks at a given time, despite external disturbances and internal imperfections. In (i), the reference tasks are generally a set of positions, velocities, and attitudes to be followed as functions of time, and can be updated (or modified) by giving appropriate signals by a human being (pilot or ground controller). The result is an automatic control system that continually compares the actual position, velocity, and attitude of the vehicle with the corresponding reference values (i) and makes the necessary corrections (ii) in order that the vehicle moves in the desired manner. Any flight vehicle must have two separate classes of control systems: first, control of position and linear velocity relative to a planet fixed frame, called trajectory control (or more specifically guidance that results in the vehicle’s navigation2 from one position to another; and second, control of vehicle’s orientation (attitude control) with respect to a frame of reference. The desired position and velocity – usually derived from the solution of a trajectory optimization problem – could be stored onboard at discrete times, serving as nominal (reference) values against which the actual position and velocity can be compared. A guidance system continually compares the vehicles’s actual position and velocity with the nominal ones, and produces linear acceleration commands in order to correct the errors. Most flight vehicles require reorientation (rotation) is realized in practice. Vehicle rotation is performed by applying an angular acceleration external/internal torques. Thus, in such a case, attitude control system becomes subservient (actuator) to the guidance system. In layman terms, the guidance system can be said to “drive the vehicle on an invisible highway in the sky” by using the attitude control system to twist and turn the vehicle. In many flight conditions, a natural stability is inherent in the attitude dynamics so that 2 The term navigation is sometimes applied specifically to the determination of a vehicle’s current position. However, in this book, we shall employ the term in its broader context, namely the process of changing the vehicle’s position between two given points. 12 Advanced Control of Aircraft, Spacecraft and Rockets ⎧ ⎨ ⎧ ⎨⎧ R (t ) − G ui da nc e Fe ed ba ck C on tr ol le r d d (t ) )t ⎫ ⎧R ⎬ ⎨ ⎭ ⎩V ( G ui da nc e Fe ed fo rw ar d C on tr ol le r d d ) ( ) (t t ⎫ ⎬ ⎭ ⎩ω σ A tti tu de Fe ed ba ck C on tr ol le r A tti tu de Fe ed fo rw ar d C on tr ol le r Fl ig ht D yn am ic s ) ( ) (t t ⎫ ⎬ ⎭ ⎩ω σ V (t )⎫ ⎨ ⎬ ⎭ ⎩ 1 ) (t w 2 ) (t w υ (t) − ) t u( F ig ur e 1. 5 B lo ck di ag ra m of a ty pi ca lfl ig ht co nt ro ls ys te m Introduction 15 x y z o Figure 1.7 Geometry of interception of a ball by a dragonfly vector due to gravity, ag, is always vertical. The time of flight is small enough that a frame of refer- ence, (oxyz), fixed to the planetary surface can be considered inertial by neglecting planetary rotation and revolution. A dragonfly initially at the origin of the frame (oxyz), is trying to alight on a cricket ball that has been hit at an initial velocity, 10i + 20k m/s, by a batsman at initial coordinates, (−10, −5, 0.5) m (Figure 1.7). Devise a proportional navigation law for the dragonfly such that it makes an in-flight rendezvous with the ball. Neglect atmospheric drag on the ball. The acceleration of the ball (the target) can be expressed as R̈T = ag = −gk, (1.40) with initial position and velocity vectors, RT(0) and VT(0), respectively. Since the acceleration is constant, the instantaneous position of the ball is obtained by integration to be RT(t) = RT(0) + VT(0)t − g 2 t2k. (1.41) The total time of flight, tf , is easily calculated as follows: zT(tf ) = 0 = zT(0) + żT(0)tf − g 2 t2f , (1.42) or tf = żT(0) g + √ ż2T(0) g2 + 2zT(0) g = 4.1023 s. Thus, it is necessary to have interception at a time T < tf . Let us select T = 30/g, which yields the following initial velocity for the dragonfly, by virtue of equation (1.15) and the dragonfly’s initial position, R(0) = 0: V(0) = VT(0) + 1 T RT(0) = ⎛ ⎝ 6.7300 −1.6350 20.1635 ⎞ ⎠ m/s. The dragonfly must be traveling with this much velocity initially in order to successfully intercept the ball by the PNG approach. 16 Advanced Control of Aircraft, Spacecraft and Rockets For the simulation of the flight of the dragonfly by a fourth-order Runge–Kutta method, we wrote a MATLAB program called run dragonfly.m that has the following statements: global g; g=9.81;% m/sˆ2 global T; T=30/g;% projected intercept time (s) global RT0; RT0=[-10 -5 0.5]’; % target’s initial position (m) global VT0; VT0=[10 0 20]’; % target’s initial velocity (m/s) V0=VT0+RT0/T % interceptor’s initial velocity (m/s) %Runge-Kutta integration of interceptor’s equations of motion: [t,x]=ode45(@dragonfly,[0 3],[0 0 0 V0’]); % Calculation of instantaneous position error: RT=RT0*ones(size(t’))+VT0*t’+.5*[0 0 -g]’*(t.ˆ2)’; error=sqrt((x(:,1)-RT(1,:)’).ˆ2+(x(:,2)-RT(2,:)’).ˆ2+(x(:,3)-RT(3,:)’).ˆ2); The equations of motion of the dragonfly are specified in the following file named dragonfly.m (which is called by run dragonfly.m): function xdot=dragonfly(t,x) global g; global T; global RT0; global VT0; rT=RT0+VT0*t+0.5*[0 0 -g]’*tˆ2; r=rT-x(1:3,1); xdot(1:3,1)=x(4:6,1); eps=1e-8; % Avoiding singularity at t=T: if abs(T-t)>eps xdot(4:6,1)=[0 0 -g]’+r/((T-t)*(T-t+1)); else xdot(4:6,1)=[0 0 -g]’; end In order to avoid the singularity in equation (1.39) at precisely t = T , we specify a zero corrective acceleration at that time, so that the dragonfly is freely falling in the small interval, T −  < t < T + . The results of the simulation are plotted in Figures 1.8–1.10. Note that due to the non-planar nature of its flight path, the dragonfly is never able to actually sit on the ball, with the smallest position error, || r || (called miss distance) being 0.0783 m at t = 2.295 s. The miss distance is a characteristic of the PNG approach, which almost never results in a zero error intercept. A relatively small miss distance (such as in the present example) is generally considered to be a successful interception. When applied to a missile, the PNG law is often combined with another approach called terminal guidance that takes over near the point of interception in order to further reduce the position error. The maximum speed reached by the dragonfly is its initial speed of 21.3198 m/s, which appears to be within its physical capability of about 80 km/hr. 1.4.2 Cross-Product Steering It was seen above in the dragonfly example that non-coplanar target and interceptor trajectories caused a near miss instead of a successful rendezvous. A suitable navigational strategy for a rendezvous would then appear to be the one that simultaneously achieves a zero miss distance and a zero relative speed. Intuition suggests that a simple approach for doing so is to turn the velocity vector, V(t), such that the Introduction 17 −10 0 10 20 −10 −5 0 0 5 10 15 20 25 x (m)y (m) z (m ) Ball Dragonfly Figure 1.8 Three-dimensional plot of the flight paths taken by the ball and the dragonfly for PNG guidance instantaneous velocity error (or velocity to be gained), v(t) = VT(t) − V(t), becomes aligned with the acceleration error, a = dv dt = dVT dt − dV dt = V̇T − V̇. (1.43) This implies that the cross-product of velocity and acceleration errors must vanish: v × a = 0. (1.44) Such a navigational law is termed cross-product steering. 0 0.5 1 1.5 2 2.5 3 0 2 4 6 8 10 12 Time (s) || r || (m ) Figure 1.9 Separation error vs. time between the ball and the dragonfly with PNG guidance 20 Advanced Control of Aircraft, Spacecraft and Rockets 0 0.5 1 1.5 2 2.5 3 3.5 0 5 10 15 || r || (m ) 0 0.5 1 1.5 2 2.5 3 3.5 0 10 20 30 40 Time (s) || v || (m /s ) Figure 1.12 Separation and relative speed error vs. time between the ball and the dragonfly with cross-product steering 0 0.5 1 1.5 2 2.5 3 3.5 5 10 15 20 25 30 35 Time (s) V ( m /s ) Figure 1.13 Time history of the required speed of the dragonfly for rendezvous by cross-product steering Introduction 21 (generalized velocity), q̇(t), both of which constitute the nonlinear system’s state vector, ξ(t) = { q(t) q̇(t) } . (1.48) For ensuring that the error from a desired (or nominal) generalized trajectory, qd(t), defined by e(t) = qd(t) − q(t), (1.49) is always kept small, a likely feedback control law is the one that applies a corrective control input, u(t), proportional to the error, the time integral of the error, and its rate as follows: u(t) = Kpe(t) + Ki ∫ t 0 e(τ)dτ + Kd ė(t), t ≥ 0, (1.50) where (Kp, Ki, Kd) are (usually constant) gain matrices. In this way, not only can a correction be applied based upon the instantaneous deviation from the trajectory, but also the historically accumulated error as well as its tendency for the future can be corrected through the integral and derivative terms, respectively. Such a control law is termed proportional-integral-derivative (PID) control and requires a measurement and feedback of not only the error vector, but also its time integral and time derivative, as shown in the schematic block diagram of Figure 1.14. The process of determination of the gain matrices (often from specific performance objectives) is called PID tuning. While the PID control was originally intended for single-variable systems, its application can be extended to multi-variable plants as follows. − ( )teqd(t) Flight Dynamics Plant u(t)Kp q(t) Kd Ki 0 ( )d t τ τ∫ e { }d ( ) d t t e PID Controller Figure 1.14 Schematic block diagram of PID control system Since a successful application of PID control would result in a small error vector, the plant’s governing equations can be generally linearized about the nominal trajectory resulting in the following governing error equation: ë = A1e + A2ė + Bu, (1.51) with the small initial conditions, e(0) = e0, ė(0) = ė0. (1.52) 22 Advanced Control of Aircraft, Spacecraft and Rockets Here, (A1, A2, B) are Jacobian matrices of the linearized plant (Appendix A) and could be varying with time along the nominal trajectory. Substitution of equation (1.50) into equation (1.51) with the error state vector x(t) = ⎧⎪⎨ ⎪⎩ ∫ t 0 e(τ)dτ e(t) ė(t) ⎫⎪⎬ ⎪⎭ (1.53) results in the following dynamic state equation of the closed-loop system: ẋ = Ax, (1.54) where the closed-loop dynamics matrix is given by A = ⎡ ⎣ 0 I 0 0 0 I BKi A1 + BKp A2 + BKd ⎤ ⎦ . (1.55) Clearly, for the error to remain small at all times t ≥ 0 beginning from a small initial state of x(0) = ⎧⎨ ⎩ 0 e0 ė0 ⎫⎬ ⎭ , (1.56) the closed-loop error dynamics (equation (1.54)) must be stable in the sense of Lyapunov (Appendix B). The PID gain matrices must therefore be selected in a way that ensures stability of the closed-loop system. However, such a derivation of the gain matrices can be termed practical only for a time-invariant system (i.e., when (A1, A2, B) are constant matrices). Example 1.5 Consider an axisymmetric, rigid spacecraft of principal moments of inertia, Jxx = Jyy and Jzz, equipped with a rotor that can apply a gyroscopic, internal torque input, u = (ux, uy)T , normal to the axis of symmetry. The spacecraft has the following equations of motion: u = Jω̇ + ω × Jω̇, (1.57) where J = ⎛ ⎝ Jxx 0 0 0 Jxx 0 0 0 Jzz ⎞ ⎠ (1.58) and ω = ⎧⎨ ⎩ ωx ωy ωz ⎫⎬ ⎭ . (1.59) When the torque input is zero, the spacecraft has an equilibrium state with pure spin about the axis of symmetry, ωT = (0, 0, n). If an initial disturbance with small, off-axis angular velocity components, ωx(0), ωy(0), is applied at t = 0, the spacecraft has the following LTI dynamics linearized about the equilibrium state: { ω̇x ω̇y } = ( 0 −k k 0 ){ ωx ωy } + 1 Jxx ( 1 0 0 1 ){ ux uy } , (1.60) Introduction 25 − xo(kΔt) Flight Dynamics Plant u(t) Digital Observer y(t) Digital Feedforward Controller Digital Regulator A/D Converter D/A Converter u(kΔt)A/D Converter d(t)x Clock Figure 1.16 Block diagram of a digital flight control system interval. At the arrival of the clock signal, the A/D block releases its output signal that was held at either a constant value (zero-order hold), or an interpolated value of the continuous time input signals received during the sampling interval. In contrast, the D/A block continuously produces its output signal as a weighted average of discrete input signals received (and stored) over the past few sampling intervals. The digital sampling and holding process involved in A/D conversion mathematically transforms the governing linear differential equations in time into linear algebraic, difference equations that describe the evolution of the system’s state over a sampling interval. In the frequency domain (Appendix A), the transform method for a discrete time, linear time-invariant system (equivalent to continuous time Laplace transform) is the z-transform. For more details about digital control systems, see Chapter 8 of Tewari (2002), or Phillips and Nagle (1995). A digital control system is unable to respond to signals with frequency greater than the Nyquist frequency, defined as half of the sampling frequency and given by π/ t. Hence, a digital control system has an inherent robustness with respect to high-frequency noise that is always present in any physical system. This primary advantage of a digital control system over an equivalent analog control system is responsible for the replacement of electronic hardware in almost all control applications in favor of digital electronics – ranging from the common digital video recorders to aircraft and spacecraft autopilots. Since a digital control system is designed and analyzed using techniques equivalent to a continuous time control system, one can easily convert one into the other for a sufficiently small sampling interval, t. While many important statistical processes and noise filters are traditionally described as discrete rather than continuous time signals, their continuous time analogs can be easily derived. Therefore, it is unnecessary to discuss both continuous time and discrete time flight control design, and we shall largely restrict ourselves to the former in the remainder of this book. 1.6 Summary A control system consists of a plant and a controller, along with actuators and sensors that are usu- ally clubbed with the plant. Closed-loop (automatic) control is the only practical alternative in reach- ing and maintaining a desired state in the presence of disturbances. There are two distinct types of automatic controllers: (a) terminal controllers that take a nonlinear system to a final state in a given time, and (b) tracking controllers which maintain a system close to a nominal, reference trajectory. A tracking controller generally consists of an observer, a regulator, and a feedforward controller. While most plants are governed by nonlinear state equations, they can be linearized about a nominal, refer- ence trajectory, which is a particular solution of the plant’s state equations. Linearized state equations 26 Advanced Control of Aircraft, Spacecraft and Rockets are invaluable in designing practical automatic controllers as they allow us to make use of linear sys- tems theory, which easily lends itself to solution by transition matrix as well as providing systematic stability criteria. All flight vehicles require control of position, velocity, and attitude for a successful mission. Separate control subsystems are required for guiding the vehicle’s center of mass along a specific trajectory, and control of vehicle’s attitude by rotating it about the center of mass. Due to their greatly different time scales, the guidance system and the attitude control system are usually designed separately and then combined into an overall flight control system. While optimal control techniques are generally required for designing the guidance and attitude control systems, certain intuitive techniques have been successfully applied for practical guidance and control systems in the past, and continue to be in use today due to their simplicity. Examples of intuitive methods include proportional navigation, cross-product steering, and proportional-integral-derivative control. Many practical control systems are implemented as discrete time (or digital) rather than continuous time (or analog) systems. Due to a bandwidth naturally limited by the sampling rate, a digital control system has an inbuilt robustness with respect to high-frequency noise signals that is not present in analog systems, which has led to the replacement of analog systems by digital systems in all modern control hardware. Exercises (1) A system has the following state equations: ẋ1 = −2x2 + u1, ẋ2 = x3 + u2, ẋ3 = −3x1. (a) If the initial condition at t = 0 is x1(0) = 10, x2(0) = 5, x3(0) = 2 and u(t) = (t, 1)T , solve the state equations for the first 10 s. (b) Is the system controllable? (See Appendix A.) (c) Is the system observable with x1 and x2 as outputs? (See Appendix A.) (2) For a system with the following state-space coefficient matrices: A = ⎛ ⎝ −2 −3 0 0 0 1 0 1 0 ⎞ ⎠ , B = ⎛ ⎝ 1 0 0 0 0 −1 ⎞ ⎠ (a) Determine the response to a zero initial condition and a unit impulse function applied at t = 0 as the first input. (b) What are the eigenvalues of the matrix A? (Ans. 1, −1, −2.) (c) Analyze the controllability of the system when only the first input is applied. (d) Analyze the controllability of the system when only the second input is applied. (3) Can the following be the state transition matrix of a homogeneous linear system with state x(t) = (−1, 0, −2)T at time t = 1? (t, 0) = ⎛ ⎝ 1 − sin t 0 −te−2t 0 e−t cos t 0 −2t 0 1 ⎞ ⎠ Why? (Ans. No.) Introduction 27 (4) A homogeneous linear system has the following state transition matrix: (t, 0) = ( cos t − sin t sin t cos t ) . If the state at time t = 1 is x(t) = (−1.3818, −0.3012)T , what was the state at t = 0? (5) Compute and plot the internal control torque components for the rotor-stabilized spacecraft of Example 1.5. What is the maximum torque magnitude required? (6) Repeat the simulation of Example 1.3 by adding the following drag deceleration on the ball: V̇T = −0.0016(VTVT) (m/s2) Is there a change in the miss distance and the maximum speed of the dragonfly? (7) Repeat the simulation of Example 1.4 with the drag term given in Exercise 4. What (if any) changes are observed? (8) Consider a controller for the plant given in Exercise 2 based upon proportional-integral-derivative feedback from the first state variable to the second control input: u2(t) = kpx1(t) + ki ∫ t 0 x1(τ)dτ + kdẋ1(t). (a) Draw a block diagram of the resulting closed-loop system. (b) Is it possible to select the PID gains, kp, ki, kd, such that closed-loop system has all eigenvalues in the left-half s-plane? (9) Replace the PID controller of Exercise 8 by a full-state feedback to the second control input given by u2(t) = −k1x1(t) − k2x2(t) − k3x1(t). (a) Draw a block diagram of the resulting closed-loop system. (b) Select the regulator gains, k1, k2, k3, such that closed-loop system has all eigenvalues at s = −1. (Ans. k1 = −1/3, k2 = −2, k3 = −1.) (c) Determine the resulting control system’s response to a zero initial condition and a unit impulse function applied at t = 0 as the first input, u1(t) = δ(t). (10) Rather than using the full-state feedback of Exercise 9, it is decided to base the controller on the feedback of the first state variable, x1(t). In order to do so, a full-order observer (Appendix A) must be separately designed for estimating the state vector from the measurement of x1(t) as well as the applied input, u2(t), which is then supplied to the regulator designed in Exercise 9. The state equation of the observer is given by ẋo = (A − LC)xo + B2u2 + Lx1, where C = (1, 0, 0), B2 = ⎛ ⎝ 0 0 −1 ⎞ ⎠ , L = ⎛ ⎝ 1 2 3 ⎞ ⎠ . 30 Advanced Control of Aircraft, Spacecraft and Rockets where ud(t) is the nominal control input. If the state vector is available for measurement and feedback, a state feedback control law called the regulator, of the form u = u − ud = h(e(t), t), (2.5) is derived by ensuring that the tracking error, e(t) = x(t) − xd(t), (2.6) is either driven to zero, or in any case kept small, in the presence of modeling errors and external disturbances. However, it is rarely possible to measure the entire state vector; instead only a few output variables, y(t), can be sensed by available sensors. The plant’s output is mathematically given by the output equation, y = g[x(t), u(t), t]. (2.7) In order to be successful, a tracking controller must have a separate subsystem called the observer that generates an estimate of the plant’s state vector, xo(t), and supplies it to the regulator according to the following observer state equation: ẋo = fo[xo(t), y(t), u(t), t], (2.8) which is designed to minimize the estimation error, eo(t) = xo(t) − x(t). (2.9) Thus, any flight control task can be practically handled by first deriving a reference trajectory and the corresponding nominal input history, ud(t), in the open loop by the solution of a nonlinear, two-point boundary value problem, and then designing a feedback controller for tracking the reference trajectory despite disturbances. The resulting control system therefore has both feedforward and feedback paths described by u = ud + h[e(t) + eo(t), t] (2.10) and schematically shown in Figure 2.1. Note the presence of process noise, ν(t), as a disturbance input, which arises due to errors in the plant model. In attempting to measure and feed back the outputs, an additional disturbance called the measurement noise, w(t), is invariably introduced. Both process and measurement noise inputs are assumed to be statistically known, and hence excluded from the nominal plant model used to derive the terminal controller. However, the tracking controller usually incorporates statistical models for filtering out the noise signals through an appropriately designed observer. Owing to the fact that a successful tracking control keeps the plant’s state rather close to the reference trajectory, the assumption of small deviations from the nominal can usually be employed, leading to a linearized plant and great simplification in the controller design. This is the most common approach in deriving flight controllers and takes full advantage of the linear systems theory (Appendix A). Typically, one writes the linearized state equation of the plant as ė = Ae + Bu + Fν, (2.11) where ν(t) is the process noise vector, u = u − ud, (2.12) Optimal Control Techniques 31 o ),+ teh(e d )(tx Terminal Controller )(x,u,tx=f Plant ),(x,u, ),( t t x=f x,u,y=g ν ν )tw( − ooo ),u,y,(x tf=x d(t)u )ty( )tu()tu(Δ o )(tx (t)ν (Tracking Control) Observer Regulator ii ff ) = x ] = 0),s[ t tt x( x( Interior Constraints )]),u([ tt ≤ 0x(α o )() tt + ee( Figure 2.1 A general control system the corrective tracking control, and A(t), B(t), and F(t) are the Jacobian matrices (Appendix A) derived from a first-order Taylor series expansion of equation (2.1) about the reference trajectory, xd(t). Taking advantage of the linear systems theory, we can easily derive linear regulator laws of the u = −Ke (2.13) kind (in the case of state feedback) for driving the error vector to zero in the steady state, provided the plant satisfies certain conditions such as controllability (Appendix A). A suitable modification of equation (2.13) for the case of observer-based output feedback tracking control is easily obtained through a linear observer, provided the plant is observable with the given outputs. Optimal control theory (Athans and Falb 2007; Bryson and Ho 1975) offers a powerful method of deriving both terminal and tracking controllers, and provides a systematic framework for modern control design methods. Optimal control refers to the time history of the control inputs that take the system from an initial state to a desired final state while minimizing (or maximizing) an objective function. The evolution of the system from the initial to final states under the application of optimal control inputs is called an optimal trajectory. Finding an optimal trajectory between two system states is the primary task of terminal control that is generally carried out in open loop. A small variation in initial condition – or the presence of small external disturbances – causes a perturbation of the system from its optimal trajectory, tracking which requires a modification of the optimal control for a neighboring optimal trajectory, generally in a closed loop. Thus open-loop optimal control can yield a terminal control input, superimposed on which is the tracking control input calculated by an optimal feedback control law. Therefore, optimization of control inputs – either in an open loop or a closed loop – is integral to modern control design. 2.2 Multi-variable Optimization We begin with a few preliminaries necessary for studying multi-variable optimization problems. For a scalar function, L(x, u), called the objective function, to have a minimum with respect to a set of m 32 Advanced Control of Aircraft, Spacecraft and Rockets mutually independent, unconstrained scalar variables, u = (u1, u2, . . . , um)T , called control variables, the following necessary condition must be satisfied: ( ∂L ∂u )T = ⎧⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎩ ∂L ∂u1 ∂L ∂u2 ... ∂L ∂um ⎫⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎭ = 0. (2.14) A point u∗ that satisfies the necessary condition for a minimum is called an extremal point (or stationary point) of L. Note that an extremal point is not guaranteed to be a minimum point, as it might also be a maximum point, or neither. In order to ensure that an extremal point is also a minimum point, we require a separate condition to be satisfied by the objective function L at the extremal point. Such a condition would require an increase of L while moving away from the extremal point in any direction, and is therefore related to the sign of the second-order partial derivatives of L(x, u) with respect to u. If the square matrix ∂2L ∂u2 = ⎛ ⎜⎜⎜⎜⎜⎜⎝ ∂2L ∂u1 2 ∂2L ∂u1∂u2 · · · ∂2L ∂u1∂um ∂2L ∂u2∂u1 ∂2L ∂u2 2 · · · ∂2L∂u2∂um ... ... ... ... ∂2L ∂um∂u1 ∂2L ∂um∂u2 · · · ∂2L ∂um∂um ⎞ ⎟⎟⎟⎟⎟⎟⎠ , (2.15) called the Hessian, is positive definite (i.e., has all real, positive eigenvalues) at an extremal point, then the point is guaranteed to be a minimum of L(x, u) with respect to u. In other words, the positive definiteness of the Hessian at an extremal point is a sufficient condition for the existence of a minimum of L with respect to u, and is expressed as ∂2L ∂u2 > 0, (2.16) Similarly, if the extremal point satisfies ∂2L ∂u2 < 0, (2.17) then it is a maximum point L(x, u) with respect to u. But there can be extremal points that are neither a minimum nor a maximum point. If the Hessian has both positive and negative eigenvalues at an extremal point, then the point concerned is said to be a saddle point. Furthermore, if any of the eigenvalues of the Hessian vanishes at an extremal point then we have a singular point. A point û that satisfies the necessary and sufficient conditions for a minimum is called an optimal point (or minimum point) of L. It is possible that we may only be able to find an extremal point, but not a minimum point for a given function, L(x, u). If L(x, û) ≤ L(x, u) for all u /= û, then the optimal point û is said to be a global (or absolute) min- imum of L(x, u) with respect to u. If a minimum point is not a global minimum then it is said to be a local minimum. Optimal Control Techniques 35 which requires a non-singular Jacobian, ∂f/∂x, at an extremal point (that is also necessary for determining x∗ (and thus L∗) from u∗). Once the extremal Lagrange multipliers are determined, the second part of equation (2.20) yields ∂L ∂u + (λ∗)T ∂f ∂u = 0, (2.22) a solution of which yields the extremal point, u∗. Generally, equations (2.18), (2.21), and (2.22) require an iterative solution for the extremal values, (x∗, λ∗, u∗). Example 2.2 For the function minimization of Example 2.1, determine the minimum points, subject to the following equality constraints: x1 − 2x2 + 5x3 = 2u1 − cos u2, −2x1 + x3 = 2 sin u2, −x1 + 3x2 + 2x3 = 0. The constraint equations expressed in the form of equation (2.18) are f = ⎧⎨ ⎩ f1 f2 f3 ⎫⎬ ⎭ = ⎧⎨ ⎩ x1 − 2x2 + 5x3 − 2u1 + cos u2 −2x1 + x3 − 2 sin u2 −x1 + 3x2 + 2x3 ⎫⎬ ⎭ . Now we write ( ∂L ∂x )T = ⎧⎨ ⎩ 2x1 −2x2 3 ⎫⎬ ⎭ , ∂f ∂x = ⎛ ⎜⎝ ∂f1 ∂x1 ∂f1 ∂x2 ∂f1 ∂x3 ∂f2 ∂x1 ∂f2 ∂x2 ∂f2 ∂x3 ∂f3 ∂x1 ∂f3 ∂x2 ∂f3 ∂x3 ⎞ ⎟⎠ = ⎛ ⎜⎝ 1 −2 5 −2 0 1 −1 3 2 ⎞ ⎟⎠ . Since the determinant | ∂f/∂x |= −39 /= 0, we can invert the Jacobian matrix ∂f/∂x to derive the extremal Lagrange multipliers as follows: (λ∗)T = −∂L ∂x ( ∂f ∂x )−1 or λ∗ = − 1 39 ⎧⎨ ⎩ 6x∗1 + 6x∗2 + 18 −38x∗1 + 14x∗2 + 3 4x∗1 − 22x∗2 + 12 ⎫⎬ ⎭ . 36 Advanced Control of Aircraft, Spacecraft and Rockets Next, we derive the extremal control variables, u∗, from equation (2.22) as follows:( ∂L ∂u∗ )T = { 2u∗1 −2 cos u∗2 } , ∂f ∂u = ⎛ ⎜⎝ ∂f1 ∂u1 ∂f1 ∂u2 ∂f2 ∂u1 ∂f2 ∂u2 ∂f3 ∂u1 ∂f3 ∂u2 ⎞ ⎟⎠ = ⎛ ⎜⎝ −2 − sin u2 0 −2 cos u2 0 0 ⎞ ⎟⎠ . Thus, equation (2.22) yields the following extremal control variables: u∗1 = − 1 39 (6x∗1 + 6x∗2 + 18), u∗2 = tan−1 { 72 + 38x∗1 − 14x∗2 18 + 6x∗1 + 6x∗2 } . Finally, the extremal values, x∗ and u∗, are related by the equality constraints as follows: x∗ = ( ∂f ∂x )−1 ⎧⎨ ⎩ 2u∗1 − cos u∗2 2 sin u∗2 0 ⎫⎬ ⎭ or ⎧⎨ ⎩ x∗1 x∗2 x∗3 ⎫⎬ ⎭ = 139 ⎛ ⎝ 3 −19 2−3 −7 11 6 1 4 ⎞ ⎠ ⎧⎨ ⎩ 2u∗1 − cos u∗2 2 sin u∗2 0 ⎫⎬ ⎭ . Clearly, we have to solve for (x∗, λ∗, u∗) by an iterative scheme. Fortunately, the constraints are linear in x, allowing us to uniquely determine x∗ from u∗ by matrix inversion. If it were not so, a more sophisticated, nonlinear programming algorithm would be necessary (such as the Newton–Raphson method). A simple iterative solution procedure for problems with explicitly solvable constraints is depicted in Figure 2.3, where the convergence criterion is that the Euclidean norm of the variation of the control vector between two successive steps is less than or equal to a specified tolerance, : || u∗ − u∗0 || ≤ . The iterative solution depends upon the starting point of the iteration. We demonstrate the solution procedure for this problem, beginning with an initial guess of u∗1 = 0, u∗2 = π/2 which is a minimum point for the unconstrained minimization (Example 2.1). The computational steps with  = 10−10 are programmed in the following MATLAB statements: >> eps=1e-10; >> i=0; >> u10=0; u20=pi/2; % initial guess >> A=[1 -2 5;-2 0 1;-1 3 2]; % Jacobian, df/dx >> Ain=inv(A); Optimal Control Techniques 37 * 0Initial guess : u ** * ) = 0,uSolve : ⇒ f(x x )( *Compute : T L −1 −= f xx λ ∂∂ ∂∂ )( * * Solve : TL ⎞⎞ + =⎟⎟ ⎠⎠ ⎛ ⎜⎝ ⎛ ⎜⎝ ⎞ ⎟⎠ ⎛ ⎜⎝ ⎞ ⎟⎠ ⎛ ⎜⎝ ⇒ 0 uu u λ ∂f∂ ∂ ε ∂ ** 0 ?≤− uu Stop Yes * * 0 = uu No Figure 2.3 A simple iteration procedure for solving the static, constrained minimization problem with explicit constraints, f(x, u) = 0 >> x=Ain*[2*u10-cos(u20);2*sin(u20);0]; %solving the constraints for x >> u1=-(6*x(1,1)+6*x(2,1)+18)/39; >> u2=atan((72+38*x(1,1)-14*x(2,1))/(18+6*x(1,1)+6*x(2,1))); >> du=[u1-u10;u2-u20]; >> while norm(du)>eps >> u10=u1;u20=u2; >> x=Ain*[2*u10-cos(u20);2*sin(u20);0]; >> u1=-(6*x(1,1)+6*x(2,1)+18)/39; >> u2=atan((72+38*x(1,1)-14*x(2,1))/(18+6*x(1,1)+6*x(2,1))); >> du=[u1-u10;u2-u20] >> i=i+1 >> end The resulting solution, ⎧⎨ ⎩ x∗1 x∗2 x∗3 ⎫⎬ ⎭ = ⎧⎨ ⎩ −1.00114096910 −0.28584773541 −0.07179888144 ⎫⎬ ⎭ ,