Basic Nonlinear Control Methods

Chapter 2 Basic Nonlinear Control Methods Overview This chapter gives a survey of the basic nonlinear control methods. First, the main classes of nonlinear systems are summarized, then some examples are presented how to find the dynamic model of simple systems and an overview of the stability methods is given. It is followed by an introduction to the concepts of some often used nonlinear control methods like input/output linearization, flatness control, backstepping and receding horizon control. 2.1 Nonlinear System Classes At the interface of a nonlinear system Σ , we can observe input and output signals. We can denote the input signal by u(·) and its value at time moment t by u(t). Similarly, y(·) is the output signal and y(t) is its value at time moment t. We can assume that the signals are elements of functions spaces which are usually linear spaces having an appropriate norm and some topological properties (convergence of Cauchy sequences of functions, existence of scalar product etc.). The chosen function space (Banach space, Hilbert space etc.) may depend on the control problem to be solved. The state x of the nonlinear dynamic system is an abstract information which represents the entire past signals of the system observed at the interface, that is, x(τ ) = x represents in abstract form all the signals u(σ ), y(σ ), ∀σ ≤ τ . The system can be characterized by the state transition function x(t) = φ(t, τ, x, u(·)) and the output mapping y(t) = g(t, x(t), u(t)) where t is the actual time, τ and x are the initial time and initial state respectively, and for causal systems u(·) denotes the input signal between τ and t. Specifically, for linear systems the transients can be written as x(t) = Φ(t, τ )x + Θ(t, τ )u(·) and y(t) = C(t)x(t) + D(t)u(t). The state transition for linear systems can be thought as the superposition of the initial state transient and the result of the input excitation. B. Lantos, L. Márton, Nonlinear Control of Vehicles and Robots, Advances in Industrial Control, DOI 10.1007/978-1-84996-122-6_2, © Springer-Verlag London Limited 2011 11 12 2 Basic Nonlinear Control Methods Fig. 2.1 Block scheme of the nonlinear system 2.1.1 State Equation of Nonlinear Systems If the state transition function satisfies some continuity conditions, then the state is the solution of the nonlinear state differential equation with appropriately chosen f (t, x, u) and the system is called continuous time smooth nonlinear system: ẋ = f (t, x, u), (2.1) y = g(t, x, u), where “dot” denote time derivative and we used short forms in notation instead of ẋ(t) = f (t, x(t), u(t)) and y(t) = g(t, x(t), u(t)). We shall consider dominantly smooth systems in the book therefore “smooth” will be omitted and the system is called continuous time nonlinear time variant system (NLTV). If f does not depend explicitly on t, that is, f (x, y), the nonlinear system is time invariant (NLTI). We shell usually assume that x ∈ R n , u ∈ R m , y ∈ R p . The block scheme of the nonlinear system is shown in Fig. 2.1. Linear systems are special cases of nonlinear ones. Continuous time linear time varying (LTV) systems have linear state differential equation: ẋ = A(t)x + B(t)u, y = C(t)x + D(t)u, (2.2) The solution of the state equation for LTV system can be expressed by the fundamental matrix Φ(t, τ ) which is the solution of the matrix differential equation Φ̇(t, τ ) = A(t)Φ(t, τ ), Φ(τ, τ ) = I. (2.3) The fundamental matrix has some important properties: Φ(t, τ )Φ(τ, θ ) = Φ(t, θ ), Φ(t, τ )Φ(τ, t) = I, dΦ(t, τ ) = −Φ(t, τ )A(τ ). dτ (2.4) Based on these properties, the solution of the state differential equation for linear systems can be expressed by the fundamental matrix:  t x(t) = Φ(t, τ )x + Φ(t, θ )B(θ )u(θ ) dθ. (2.5) τ The block scheme of the LTV system is shown in Fig. 2.2. 2.1 Nonlinear System Classes 13 Fig. 2.2 Block scheme of the LTV system Fig. 2.3 Mathematical scheme for solving nonlinear state equation For continuous time linear time invariant (LTI) systems, A, B, C, D are constant matrices and the fundamental matrix is the exponential matrix Φ(t, τ ) = eA(t−τ ) , moreover the initial time can be chosen τ = 0 because the system is time invariant. The LTI system has transfer function G(s). The solution of the state (“differential” will be omitted) equation and the transfer function are  t eA(t−θ) Bu(θ ) dθ, (2.6) x(t) = eAt x + 0 G(s) = C(sI − A)−1 B + D, (2.7) where (sI − A)−1 is the Laplace-transform of eAt . The state equation of nonlinear systems can usually be solved only numerically hence we show some numerical methods to find the approximate solution. Here, we assume that the control u(t) := u0 (t) is known therefore ẋ(t) = f (t, x(t), u0 (t)) =: f ⋆ (t, x(t)), and for simplicity we omit “star” in the notation. Assume x(tn ) is known and we want to determine x(tn + h), see Fig. 2.3. Euler method:   x(tn + h) ≈ x(tn ) + hx ′ (tn ) = x(tn ) + hf tn , x(tn ) . Second order Taylor series: h2 x(tn + h) ≈ x(tn ) + hx ′ (tn ) + x ′′ (tn ) 2   = x(tn ) + hf tn , x(tn ) (2.8) 14 2 + Basic Nonlinear Control Methods      h2  ′  fx tn , x(tn ) · f tn , x(tn ) + ft′ tn , x(tn ) . 2 (2.9) Unfortunately, higher order Taylor series are hard to determine hence Runge and Kutta developed such methods which differ from the n-th order Taylor series only in o(hn ). Now we show two of such methods. Second order Runge–Kutta method:     x(tn + h) ≈ x(tn ) + h a1 f tn , x(tn ) + a2 f tn + αh, x(tn )   + βf tn , x(tn ) , (2.10) where the unknowns are a1 , a2 , α, β. In order to guarantee o(h2 ) (small ordo) precision a1 = 1 − a2 , α = 2a12 , β = 2ah2 should be satisfied. For the remaining parameter, two methods are wide-spread: Heun’s trapezoidal algorithm: a2 = 21 . Modified Euler–Cauchy algorithm: a2 = 1. Fourth order Runge–Kutta method: The essence of the method is that the change of the function x(t) will be estimated from the values of the derivatives determined at four points of the function (once at the left end of the interval, twice in the middle and once at right end), and take the average of them with weights assuring o(h4 ) approximation at the right end point of the interval:   k1 = hf tn , x(tn ) ,   1 h k2 = hf tn + , x(tn ) + k1 , 2 2   1 h (2.11) k3 = hf tn + , x(tn ) + k2 , 2 2   k4 = hf tn + h, x(tn ) + k3 , 1 x(tn + h) ≈ x(tn ) + (k1 + 2k2 + 2k3 + k4 ). 6 The computation is usually performed with variable step size. The goal is to find the solution as quick as possible while guaranteeing the necessary precision for which simple step change strategies are available. Multistep methods: Other type of methods assume equidistant step size and try to find exact solution for polynomial form x(t). The solution is assumed in the form k k xn+1 = i=0 ai xn−i + bi hf (tn−i , xn−i ). i=−1 (2.12) 2.1 Nonlinear System Classes 15 Performing the normalization of the running time . . . , tn−1 , tn , tn+1 to normalized time . . . , −1, 0, 1 the following conditions appear: k k 1= i=0 ai (−i)j + bi j (−i)j −1 , i=−1 j ∈ [1, m], (2.13) k 1= ai . i=0 The accuracy can be influenced by m and k. Wide-spread methods are: • Adams–Bashforth method: k = m − 1, a1 = · · · = ak = 0, b−1 = 0. • Adams–Moulton method: k = m − 2, a1 = · · · = ak = 0, b−1 = 0. • Predictor–corrector method: The first Adams–Bashforth (predictor) step determines the estimation xn+1 because now b−1 = 0 and thus xn+1 does not appear on the right side of the equation system. After this, the last xn+1 can cyclically be applied on the right side of the Adams–Moulton (corrector) steps gradually improving the approximation. 2.1.2 Holonomic and Nonholonomic Systems Nonlinear systems have often to satisfy constraints that restrict the motion of the system by limiting the paths which the system can follow. For simplicity, we assume that the configuration space Q is an open subset of R n with coordinates q = (q1 , . . . , qn )T . More general configuration spaces can be handled by an appropriate choice of local coordinates. A simple example is the case of two particles attached by a massless rod. The configuration of each particle is described by a point ri ∈ R 3 , but all trajectories of the particles have to satisfy r1 − r2 2 = L2 , where L is the length of the rod. The constraints act through the application of constraint forces modifying the motion of the system such that the constraints are always satisfied. The constraint is an example of a holonomic constraint. In the sequel, we shall define holonomic and nonholonomic constraints according to Murray and Li, and Sastry, see [103]. Then we investigate how the motion equation can be simplified by using coordinate transformation. 2.1.2.1 Holonomic Constraints More generally, a constraint is called holonomic if it restricts the motion of the system to a smooth hypersurface in the unconstrained configuration space Q. Holo- 16 2 Basic Nonlinear Control Methods nomic constraints can be represented locally by algebraic constraints on the configuration space: hi (q) = 0, i = 1, . . . , k. (2.14) We shall assume that the constraints are linearly independent, that is, the matrix ∂h ∂q has full row rank. Since holonomic constraints define a smooth hypersurface in the configuration space, it is possible to eliminate the constraints by choosing a set of coordinates for this surface. These new coordinates parametrize all allowable motions of the system and are not subject to any further constraints. The constraint forces are linear combinations of the gradients of the constraint functions hi : Q → R. Let h : Q → R k represent the vector valued constraint function then the constraint forces can be written as Fc = ∂h ∂q T λ, (2.15) where λ ∈ R k is called Lagrange multiplier vector measuring the relative magnitudes of the constraint forces. No work is done if the system is moved along feasible trajectories since FcT q̇ = λT  ∂h d q̇ = λT h(q) = 0. ∂q dt A fundamentally more general constraint occurs in some robotics and mechanical applications (multifingered grasping, rolling disk etc.) where the constraint has the form AT (q)q̇ = 0. (2.16) Such constraint is called Pfaffian constraint. Notice that AT (q) is k × n and q̇ appears linearly in the constraint equations. This linearity gives chance to eliminate the constraints using a nonlinear coordinate transformation based on the fundamental results of differential geometry summarized in Appendix B. On the other hand, we shall consider mechanical systems and will also refer to the results in Appendix A, especially we shall use the Lagrange equation describing the dynamic model. Since Pfaffian constraints restrict the allowable velocity of the system but not necessarily the configurations, we cannot always represent it as an algebraic constraint on the configuration space. The Pfaffian constraint is called integrable if there exists a vector valued function h : Q → R k such that AT (q)q̇ = 0 ⇐⇒ ∂h q̇ = 0. ∂q (2.17) Thus, an integrable Pfaffian constraint is equivalent to a holonomic constraint. No∂h , but only that they define the same subspace tice that we do not require that AT = ∂q of the allowable velocities at every q ∈ Q. 2.1 Nonlinear System Classes 17 2.1.2.2 Nonholonomic Constraints Pfaffian constraint that is not integrable is an example of nonholonomic constraint. For such constraints, the instantaneous velocities of the system are constrained to an n − k dimensional subspace, but the set of reachable configurations is not restricted to some n − k dimensional hypersurface in the configuration space. Since not all Pfaffian constraints are integrable, hence the motion equations have to be extended to this case. It is still possible to speak about constraint forces. Such forces are generated by the Pfaffian constraints to insure that the system does not move in the directions given by the rows of the constraint matrix AT (q). The constraint forces at q ∈ Q are defined by Fc = A(q)λ, where the Lagrange multiplier λ ∈ R k is the vector of relative magnitudes of the constraint forces. If the constraints are integrable, then ∂h T this is identical to the holonomic case since A(q) and [ ∂q ] have the same range space. The constraint forces for Pfaffian constraints prevent motion of the system in directions which would violate the constraints. In order to include these forces in the dynamics, we have to add an additional assumption, namely forces generated by the constraints do no work on the system. This assumption is often referred as d’Alambert principle. By using the Lagrangian L(q, q̇) for the unconstrained system and assuming the constraints in the form AT (q)q̇ = 0 the equation of motion can be written in the form d ∂L ∂L − = A(q)λ + F, (2.18) dt ∂ q̇ ∂q where F represents the nonconservative and externally applied forces. For systems having Lagrangian L(q, q̇) = 12 q̇ T M(q)q̇ − P (q) (kinetic energy minus potential energy), the motion equations can be written as M(q)q̈ + C(q, q̇)q̇ + n(q, q̇) = A(q)λ + F, where n(q, q̇) denotes nonconservative and potential forces and F is the external force. Differentiating the constraint equation yields AT (q)q̈ + ȦT (q)q̇ = 0. (2.19) By solving the motion equation for q̈ and substituting it into the differentiated constraint equation, we obtain: q̈ = M −1 (Aλ + F − C q̇ − n), AT M −1 Aλ = AT M −1 (C q̇ + n − F ) − ȦT q̇, −1  T −1   A M (C q̇ + n − F ) − ȦT q̇ , λ = AT M −1 A where we suppressed the dependence on q and q̇. (2.20) 18 2 Basic Nonlinear Control Methods 2.1.2.3 Lagrange–d’Alambert Equations It is useful to derive the equations without solving for the instantaneous constraint forces present in the system. It means that we project the motion of the system onto the feasible directions by ignoring the forces in the constrained directions. The resulting dynamics is in a form well suited for closed-loop control. At a given configuration q the directions at which the system is allowed to move is given by the null space of AT (q). We adopt the classical notation of virtual displacement δq ∈ R n which satisfies AT (q)δq = 0. If F is a generalized force applied to the system, then we call δW = F · δq the virtual work due to F acting along a virtual displacement δq, where “dot” denotes scalar product. D’Alambert’s principle states that the constraint forces do no virtual work, hence   (2.21) A(q)λ · δq = 0 for AT (q)q̇ = 0. Notice that δq is not the same as q̇. The generalized velocity q̇ satisfies both the velocity constraints and the motion equations, while the virtual displacement only satisfies the constraints. In order to eliminate the constraint forces from the motion equations, we can project the motion onto the linear subspace generated by the null space of AT (q). Since (Aλ) · δq = 0, we obtain the so called Lagrange–d’Alambert equations:   d ∂L ∂L − − F · δq = 0, (2.22) dt ∂ q̇ ∂q AT (q)δq = 0. (2.23) Specifically, if q = (q1T , q2T )T and AT (q) = [AT1 (q) AT2 (q)], where AT2 (q) is square and invertible, then we can use the constraint q̇2 = −(AT2 )−1 AT1 q̇1 to eliminate q̇2 and q̈2 from the motion equations. The method is illustrated in the following example. 2.1.2.4 Dynamics of Rolling Disk Consider a disk rolling on the horizontal plane such that the disk plane remains orthogonal to the horizontal plane during motion. The configuration of the disk is given by the position (x, y) of the contact point, the heading angle θ , and the orientation of the disk with respect to the vertical, φ. We use the notation q = (x, y, θ, φ)T . Let R denote the radius of the disk. We assume as input the driving torque on the wheel, τθ , and the steering torque about the vertical axis, τφ . Let m be the mass of the disk, I1 its inertia about the horizontal (rolling) axis, and I2 its inertia about the vertical axis. We assume that the disk rolls without slipping hence the following velocity constraints have to be satisfied: ẋ − R φ̇Cθ = 0, ẏ − R φ̇Sθ = 0 ⇐⇒ AT (q)q̇ = 1 0 0 1 0 −RCθ q̇ = 0. 0 −RSθ 2.1 Nonlinear System Classes 19 These constraints assume that the disk rolls in the direction in which it is heading. It is evident that the constraints are everywhere linearly independent. The unconstrained Lagrangian is the kinetic energy associated with the disk:  1 1 1  L(q, q̇) = m ẋ 2 + ẏ 2 + I1 θ̇ 2 + I2 φ̇ 2 . 2 2 2 Denote δq = (δx, δy, δθ, δφ)T the virtual displacement of the system. The Lagrange–d’Alambert equations are given by ⎡ ⎤⎞ ⎛⎡ ⎤ m 0 0 0 0 ⎢ 0 ⎥⎟ ⎜⎢ 0 m 0 0 ⎥ 1 0 0 −RCθ ⎢ ⎥⎟ ⎜⎢ ⎥ ⎝⎣ 0 0 I1 0 ⎦ q̈ − ⎣ τθ ⎦⎠ · δq = 0, where 0 1 0 −RSθ δq = 0. 0 0 0 I2 τφ From the constraints, we can solve for δx and δy: δx = RCθ δφ, The scalar product can be written as  0 0 I ẍ + 1 0 mRCθ mRSθ ÿ δy = RSθ δφ. 0 I2 τ θ̈ − θ τφ φ̈  · δθ = 0, δφ and since δθ and δφ are free, we obtain the motion equation 0 mRCθ 0 mRSθ I ẍ + 1 0 ÿ 0 I2 τ θ̈ = θ . τφ φ̈ Differentiating the constraints, we obtain ẍ = RCθ φ̈ − RSθ θ̇ φ̇, ÿ = RSθ φ̈ + RCθ θ̇ φ̇, and substituting the results into the motion equation we obtain the system dynamics as follows: I1 0 0 I2 + mR 2 τ θ̈ = θ . τφ φ̈ Note that in this simple example the system of equations does not depend on x and y, but in general it is not valid. The motion of the x- and y-position of the contact point can be obtained from the first order differential equations ẋ = RCθ φ̇, ẏ = RSθ φ̇, which follow from the constraints. Thus, given the trajectory of θ and φ, we can determine the trajectory of the disk rolling on the horizontal plane. Dividing the motion equations into a set of second-order differential equations in a reduced set of variables plus a set of first-order differential equations representing the constraints is typical for many nonholonomic systems. 20 2 Basic Nonlinear Control Methods 2.1.2.5 Holonomic Property Related to Frobenius Theorem Consider the nonlinear system satisfying constraints in which beside the generalized coordinate q ∈ R n appears also its derivative q̇ ∈ R n : AT (q)q̇ = 0, rank AT (q) = k < n. (2.24) Assume that AT (q) is k × n and denote its null (kernel) space ker AT (q) at q0 . Let us define the function (q0 ) := ker AT (q0 ) which maps q0 into a linear subspace ker AT (q0 ), hence it is a distribution. If this distribution is involutive, that is, it ∂f ∂g is closed for Lie bracket defined by [f, g] = ∂q f − ∂q g for vector fields f (q), g(q) then, based on Frobenius theorem, we may find local coordinates q̄ = (q̄1 , . . . , q̄n )T such that (2.24) can be written as q̄˙ n−k+1 = · · · = q̄˙ n = 0 ⇔ q̄n−k+1 = cn−k+1 , . . . , q̄n = cn (2.25) for appropriately chosen constants cn−k+1 , . . . , cn determined by the initial conditions and thus we can eliminate the coordinates q̄n−k+1 , . . . , q̄n . Therefore, we call the constraints (2.24) holonomic if the distribution  is involutive, otherwise the constraints are nonholonomic. Notice that the application of Frobenius theorem needs to solve a set of partial differential equations which is a hard problem. Hence, we seek for other methods simplifying the motion equation in the presence of constraints. 2.1.2.6 Motion Equation in Hamiltonian Form The equation of motion for the constrained mechanical system can be assumed in Lagrange form     d ∂L ∂L − = A(q)λ, AT (q)q̇ = 0, (2.26) dt ∂ q̇ ∂q where A(q)λ describes the constrained forces, λ ∈ R k are the Lagrange multipliers. For simplicity of notation, the generalized external force was neglected. The Lagrangian L(q, q̇) is defined by L = K − P where K is the kinetic energy and P is the potential energy. The kinetic energy is the quadratic form K = 12 M(q)q̇, q̇. Here, M(q) is the generalized inertia matrix. We can assume that no simplifications were applied hence M(q) is positive definite and hence det ∂ 2L ∂ q̇i ∂ q̇j = 0. (2.27) In order to eliminate the Lagrange multipliers and find the state equation of the mechanical system, the Legendre transformation p := ∂L ∂ q̇ = M(q)q̇ and the Hamil- 2.1 Nonlinear System Classes 21 tonian H = K + P = 2K − L will be applied. Because of K = tonian can be written as n H (q, p) = i=1 pi qi − L(q, q̇), pi = ∂L , ∂ q̇i 1 2 p, q̇, the Hamil- i = 1, . . . , n. (2.28) It follows from the new form of the kinetic energy and the positive definiteness of M(q) that p is an independent system, hence the motion equation of the constrained mechanical system can be written in the Hamiltonian form q̇ = ∂H (q, p) , ∂p ṗ = − ∂H (q, p) + A(q)λ, ∂q ∂H (q, p) A (q) = 0. ∂p (2.29) T The constraint forces A(q)λ can be determined by differentiating the constraint equation by the time and expressing λ, that is, ∂AT (q) ∂H (q, p) ∂q ∂p + AT (q) T ∂H (q, p) ∂p ∂H (q, p) ∂ 2 H (q, p) + A(q)λ = 0. − 2 ∂q ∂p (2.30) Notice that here a matrix AT is differentiated by a vector q therefore it is a problem how to implement it. Fortunately not the derivative has to be determined but the derivative multiplied by a vector ∂H ∂p which is easier to obtain. Let us consider the following mathematical problem. Let A be any matrix of type m × n and the vector x have dimension n. Then A can be assumed as  A = i,j aij Eij where the only nonzero element of the matrix Eij is in the intersection of the i-th row and the j -th column. Hence, differentiating by aij and applying the result for the vector x we obtain Eij x which is a column vector with the only nonzero element in the i-th row which is xj . Repeating the process for all the elements of A we obtain ∂A ∂x x which is an m × (mn) type matrix each column containing the above result for one element of A. T Now we can apply this result for the matrix AT and the vector ∂H ∂p . Since rank A is k, hence it follows from (2.27) and (2.28) det AT (q) ∂ 2 H (q, p) A(q) = 0, ∂p 2 (2.31) thus λ can be expressed from (2.30) and put into (2.29) which results in the motion equation of the constrained (holonomic or nonholonomic) system. 22 2 Basic Nonlinear Control Methods 2.1.2.7 Coordinate Transformation of Constrained Motion Equations Here we try to eliminate the constraint equations without the solution of partial differential equations. The idea is to introduce a coordinate transformation p̃ = T (q)p and partition p̃ as p̃ = (p̃ 1T p̃ 2T )T where p̃ 2 = 0 suits to the constraints. Notice that this is only possible if the constraints are holonomic, otherwise we obtain the motion equation in the variables q and the new p̃ and the constraints are not eliminated. Of course the goal is to eliminate the constraints and assure p̃ 2 = 0. By the assumptions AT (q) has type k × n and rank AT (q) = k thus there exists a matrix S(q) of type n × (n − k) satisfying AT (q)S(q) = 0 and having rank S(q) = n − k. Hence, we are able to define the coordinate transformation according to   AT (q)S(q) = 0 ⇒ S(q) = S 1 (q) . . . S n−k (q) , (2.32) S T (q) , p̃ = T (q)p AT (q)   H̃ (q, p̃) = H q, T −1 (q)p̃ . T (q) = ⇔ p = T −1 (q)p̃, (2.33) (2.34) In order to find the transformed motion equations, we have to determine the derivatives of H̃ and the relations to the derivatives of H . During this process, we can apply the rule about the derivative of any inverse matrix A−1 by a scalar −1 −1 variable a: ∂A∂a = −A−1 ∂A ∂a A . Denote ei the i-th standard unit vector and T T a, b = a b = b a = b, a the scalar product of a and b, respectively. Then   −1 p̃j ej , H̃ (q, p̃) = H q, T (q) j     ∂H (q, p) −1 ∂ H̃ (q, p̃) ∂H (q, p) , T (q)ei = T −T (q) , ei , = ∂ p̃i ∂p ∂p ∂ H̃ (q, p̃) ∂H (q, p) ∂H (q, p) ∂ H̃ (q, p̃) = T −T (q) ⇔ = T T (q) , (2.35) ∂ p̃ ∂p ∂p ∂ p̃   ∂ H̃ (q, p̃) ∂H (q, p) ∂T (q) −1 ∂H (q, p) , −T −1 (q) = + T (q)p̃ , ∂qi ∂qi ∂p ∂qi   ∂ H̃ (q, p̃) ∂H (q, p) ∂H (q, p) ∂T (q) −T , = − T (q) p , ∂qi ∂qi ∂p ∂qi   ∂H (q, p) ∂ H̃ (q, p̃) ∂ H̃ (q, p̃) ∂T (q) , (2.36) = + p . ∂qi ∂qi ∂ p̃ ∂qi p̃ 2 Let us consider now the time derivative of p̃. Since p̃ = T (q)p, p̃ 1 = S T (q)p, = AT (q)p, therefore   dT (q) ∂H (q, p) p̃˙ = p + T (q) − + A(q)λ . (2.37) dt ∂q 2.1 Nonlinear System Classes 23 Since S T (q)A(q) = 0, thus the constraint forces AT (q)A(q)λ appear only in p̃˙ 2 . For simplicity of the next derivations, we suppress writing down the variables q, p. First, using (2.33) and (2.37), it follows dS iT ∂H p̃˙ i1 = p − S iT . dt ∂q (2.38) Second, if the constraints are satisfied then, using (2.35) and (2.38), it yields ∂H ∂ H̃ ∂ H̃ ∂ H̃ = AT [S A] = AT A 2 ⇔ = 0, (2.39) ∂p ∂ p̃ ∂ p̃ ∂ p̃ 2    ∂ H̃ dS iT ∂S i  1 ∂S i ∂H S (q) . . . S n−k (q) p = p, , (2.40) = pT dt ∂q ∂p ∂q ∂ p̃ 1   ∂ H̃ ∂T ∂S ∂A ∂ H̃ ∂S ∂ H̃ , = pT p = pT . (2.41) ∂ p̃ ∂ql ∂ql ∂ql ∂ p̃ ∂ql ∂ p̃ 1 0 = AT q̇ = AT Third, using (2.36) the second term in (2.38) can be determined: S iT S iT ∂H = ∂q ∂H = ∂q Sli l l ∂H = ∂ql ∂ H̃ Sli ∂ql Sli l +p T l   ∂ H̃ ∂S ∂ H̃ , + pT ∂ql ∂ql ∂ p̃ 1   ∂S 1 ∂S n−k i ∂ H̃ Sl ... . ∂ql ∂ql ∂ p̃ 1 (2.42)  ∂S 1 1 i It is easy to show that, for example, l ∂qll Sli = ∂S ∂q S which together with the appropriate term from (2.40) and (2.38) results in the Lie bracket with negative sign   1  i 1 ∂S i ∂S i 1 S − S − S ,S = − ∂q ∂q (2.43) and in the state equations of p̃i1 , i = 1, . . . , n − k,   ∂ H̃ ∂ H̃  T  i 1  + −p S , S · · · − p T S i , S n−k . p̃˙ i1 = −S iT ∂q ∂ p̃ 1 (2.44) Using similar steps the state equations of p̃i2 , i = 1, . . . , k can be obtained:   ∂ H̃ ∂ H̃  T  i 1  + −p A , S · · · − p T Ai , S n−k p̃˙ i2 = −AiT ∂q ∂ p̃ 1 + AiT (q)A(q)λ. (2.45) 24 2 Basic Nonlinear Control Methods Introducing the notation ⎞ ⎛ 0 S ⎛ ⎞ T 1 1 T 1 n−k ⎜ −p [S , S ] ··· −p [S , S ] ⎟ ⎟ ⎜ ⎟⎟ ⎜ −S T ⎜ .. .. .. ⎝ ⎠ ⎟ ⎜ . . . ⎟ ⎜ T n−k 1 T n−k n−k ⎜ , S ] · · · −p [S ,S ] ⎟ −p [S ⎟ ⎜ ˜ J (q, p̃) = ⎜ ⎟ ⎟ ⎜ ⎛ ⎞ T 1 1 T 1 n−k ⎟ ⎜ −p [A , S ] · · · −p [A , S ] ⎟ ⎜ ⎜ ⎜ ⎟ ⎟ .. .. .. ⎠ ⎝−AT ⎝ ⎠ . . . T k 1 T k n−k −p [A , S ] · · · −p [A , S ] the dynamic model will be ⎛ ⎞ ⎛ ⎛ ⎞ ⎞ ∂ H̃ (q,p̃) q̇ 0 ∂q ⎜ ⎟ ⎝ ⎜ ˙ 1⎟ ⎠, 0 ⎠+ ⎝p̃ ⎠ = J˜(q, p̃) ⎝ ∂ H̃ (q,p̃) T (q)A(q)λ 2 A ˙ p̃ ∂ p̃ 1 As can be seen, although ∂ H̃ (q,p̃) ∂ p̃ 2 ∂ H̃ (q, p̃) = 0. ∂ p̃ 2 (2.46) (2.47) = 0, that is, H̃ does not depend on p̃ 2 , there is no guarantee that by using the above transformation the derivative of p̃ 2 is zero and the constraint equation can be omitted. The reason is that, denoting by J˜r (q, p̃ 1 ) the upper part of J˜(q, p̃), ⎞ ⎛ 0 S ⎞ ⎛ −p T [S 1 , S 1 ] ··· −p T [S 1 , S n−k ] ⎟   ⎜ ⎟ ⎜ J˜r q, p̃ 1 := ⎜ T ⎜ ⎟⎟ , (2.48) . . .. .. .. ⎠⎠ ⎝ ⎝−S . −p T [S n−k , S 1 ] · · · −p T [S n−k , S n−k ] a bracket { , }r can be defined for scalar valued functions Fr , Gr mapping (q, p̃ 1 ) in R 1 by  ∂G r    ∂F T ∂F T  ∂q 1 1 r r ˜ {Fr , Gr }r q, p̃ := (q, p̃ ) J . (2.49) r 1 ∂G ∂q ∂ p̃ r ∂ p̃ 1 This bracket is skew-symmetric, satisfies the Leibniz rule, but usually does not satisfy the Jacobi identity. Is also the Jacobi identity satisfied then the bracket { , }r is a so called Poisson bracket, in which case the derivative of p̃ 2 is zero and the constraint equations can be omitted, see [128]. 2.1.3 Differentially Flat Systems In Appendix B, we have shown that under special conditions, perhaps locally, the nonlinear system can be linearized using static or dynamic feedback and then it can 2.1 Nonlinear System Classes 25 be further compensated by using well known linear control methods. However, it remains a serious problem to obtain the symbolic solution of a system of partial differential equations (Frobenius equations) for which no general and efficient methods are available. All these offer reason for seeking special solutions in which range the flatness control is a possible approach. The idea of differentially flat systems arose at the 90th years based on the works of David Hilbert and Elie Cartan at the beginning of the 20th century. We shall discuss differentially flat systems using the results of [39]. Fundamentals of differential geometry and novel results in flatness control can be found in [1] and [80], respectively. First, the idea of dynamic system will be extended to infinite dimensional manifolds which can be obtained from conventional systems through prolongation of vector fields. For infinite dimensional system, the so called Lie–Bäcklund transformation will be defined which is equivalent to an endogenous dynamic feedback and state transformation performed on an appropriate finite dimensional system. Based on this result, the class of differentially flat systems will be defined in such a way that the systems belonging to it can be converted to the trivial series of integrators by the help of the Lie–Bäcklund transformation. First of all, we start with a new interpretation of the vector field. In Appendix B, the idea of the manifold X ∈ R n , its tangent space T X and the Lie derivative Lf h of smooth scalar valued function to a vector field f have already been  h relative ∂h defined. Considering Lf h = ni=1 fi ∂x allow us to identify the vector field f by i n f = (f1 , . . . , fn )T ∼ fi i=1 ∂ ∂xi (2.50) which can be assumed the componentwise or differential operator expression of the vector field f . This interpretation can easily be generalized for infinite dimension. 2.1.3.1 Prolongation of Vector Fields, the Cartan Fields Let us consider the system given by the differential equation ẋ = f (x, u) (2.51) where x is a point of an n dimensional manifold M, u ∈ R m and f is C (∞) vector field on M which can be described by the local coordinates ∂x∂ i (i = 1, 2, . . . , n) in  the form f = ni=1 fi ∂x∂ i . It is assumed u(t) ∈ C (∞) . Let us introduce the infinite dimensional vector  T ξ(t) = [ξ1 ξ2 ξ3 · · · ]T = x(t)T u(t)T u̇(t)T · · · m = M̄, that is, ξ is a point of an infinite for which ξ ∈ M × R m × R m × · · · = M × R∞ m = M̄ derived from f dimensional manifold. Let f¯ be a vector field over M × R∞ 26 2 Basic Nonlinear Control Methods by using ⎡ ⎤ f (ξ1 , ξ2 ) ⎢ ξ3 ⎥ ⎥ ⎢ ξ̇ = f¯(ξ ) = ⎢ ξ4 ⎥ . ⎦ ⎣ .. . (2.52) Notice that, despite of the infinite dimensional vector field, every row depends only on a finite number of elements of ξ . The vector field f¯ can be written in local coordinates n f¯ = fi i=1 ∂ + ∂xi m ∂ (j +1) ui (j ) (2.53) ∂ui j ≥0 i=1 which will be called a Cartan field. It follows from the construction of the vector field f¯ that if a path t → (x(t), u(t)) is a solution of (2.51) then the path t → (x(t), u(t), u̇(t), ü(t), . . .) is a solution of (2.52). In what follows, we shall understand under control system the ordered pair (M̄, f¯) where M̄ is an infinite dimensional manifold and f¯ is a Cartan field of the form (2.53) over M̄. Now consider the linear dynamic system consisting of the series of integrators: (κ ) xr1 1 = u1 (κ ) xr2 2 = u2 .. . r1 = 1; ri = ri−1 + κi−1 ; i κi = n. (κ ) xrmm = um In the above form, ri denotes the dimension, κi is the order (number of integrators). The infinite dimensional vector field in local coordinates belonging to this system can be written in the form m f¯triv = (j +1) xi j ≥0 i=1 ∂ (j ) . (2.54) ∂xi This vector field will be called the trivial vector field. Notice that there is a one-toone relation between the integral curves t → (x(t), u(t)) of the trivial system and the properly regular (differentiable) time functions t → (xr1 (t), . . . , xrm (t)). 2.1.3.2 Lie–Bäcklund Transformation Consider the control systems (M̄, f¯) and (N̄ , ḡ) derived from the differential m equations ẋ = f (x, u) and ẏ = g(y, v) (dim u = dim v) where M̄ = M × R∞ 2.1 Nonlinear System Classes 27 m (dim N = l), that is, the traditional dimension num(dim M = n) and N̄ = N × R∞ ber of the two systems are different, while the number of inputs is the same. Furthermore consider an invertible C (∞) mapping which defines a coordinate transformation between M̄ and N̄ : ξ = Φ(ζ ), ξ ∈ M̄ and ζ ∈ N̄ . Denote Ψ the inverse coordinate transformation which is also C (∞) and let ζ = Ψ (ξ ). Only such coordinate transformations will be considered for which each of their components depends only on a finite number of coordinates, similarly to the prolonged vector fields. Thus ζ1 = ψ1 (ξ̄1 ), ζ2 = ψ2 (ξ̄2 ), . . . , ζi = ψi (ξ̄i ), . . . where ξ̄i = (ξi,1 , . . . , ξi,ri ) and ri are finite. Now consider an integral curve t → ξ(t) of the vector field f¯. Its image over the manifold N̄ is t → ζ (t) = Ψ (ξ(t)). Denote Ψ∗ f¯ that vector field over the manifold N̄ whose integral curve is t → ζ (t) = Ψ (ξ(t)). Differentiating the curve Ψ (ξ(t)) by the time yields  ∂Ψ d  Ψ ξ(t) = · f¯(ξ ) = Lf¯ Ψ (ξ ), dt ∂ξ but since Ψ∗ f¯ : N̄ → T N̄ thus     ∂ + Ψ∗ f¯ = Lf¯ Ψ ◦ Φ = Lf¯ ψ1 ◦ Φ ∂y j ≥0   ∂ Lf¯ ψj +1 ◦ Φ (j ) . ∂v (2.55) Notice that in the above expression all terms are finite since the transformation and each element of the vector field depend on finite number of coordinates. For transformation Ψ of general form, it is not satisfied that the resulting vector field (2.55) is a Cartan field, that is, (2.55) can be derived from a “traditional” form dynamic system over N̄ . To this it has to be required that the prolongation is a series of integrators, that is, dv (j ) = LΨ∗ f¯ v (j ) = v (j +1) = ψj +2 ◦ Φ, dt from which using (2.55) follows   ∂v (j ) LΨ∗ f¯ v (j ) = Lf¯ ψ1 ◦ φ + ∂y k≥0    ∂v (j )  Lf¯ ψk+1 ◦ Φ (k) = Lf¯ ψj +1 ◦ Φ. ∂v Comparing the right sides of the last two equations, we obtain that Ψ∗ f¯ is a Cartan field over N̄ , if and only if, Lf¯ ψj +1 = ψj +2 is guaranteed for every positive j . Coordinate transformations mapping integral curves of Cartan fields into integral curves of Cartan fields will be called Lie–Bäcklund transformations. The control systems (M̄, f¯) and (N̄ , ḡ) are called Lie–Bäcklund equivalent if there exists an invertible and C (∞) coordinate transformation Ψ : M̄ → N̄ , each element of which is function of only a finite number of coordinates and Ψ∗ f¯ = ḡ. The Lie–Bäcklund transformation saves the number of inputs and the equilibrium points. Consequently it follows that any two linear controllable systems of equal number of inputs are Lie–Bäcklund equivalent. 28 2 Basic Nonlinear Control Methods 2.1.3.3 Endogenous State Feedback The state feedback is dynamic if it has own state: ż = β(x, z, v), (2.56) u = α(x, z, v). We shall deal only with such feedback where the number of inputs of the original and fed-back system does not change, that is, dim u = dim v. Using (2.56) for the original system (2.51), the differential equation of the closed loop will be   ẋ = f x, α(x, z, v) , (2.57) ż = β(x, z, v). The dynamic feedback does not necessarily save the controllability properties. An example for this is the common use of state observer and state feedback for linear systems where the closed loop system obtained by using dynamic feedback because the states of the feedback are the state variables x̂ of the observer. The system obtained so is not state controllable since the observer was designed to assure exponential decaying the estimation error according to the differential equation (ẋ − ˙ = F (x − x̂) (where F is stable and quick), hence there is no input making the x̂) behavior of x̂ independent of x. In order to eliminate such phenomenon, we restrict the investigations to so-called endogenous dynamic feedback defined as follows. Let be given a system (M̄, f¯). The feedback given in the form (2.56) is called endogenous dynamic feedback if the prolongation of the closed loop system given by (2.57) is Lie–Bäcklund equivalent to the system (M̄, f¯). Theorem 2.1 Existence condition for endogenous state feedback. Consider the systems (M̄, f¯) and (N̄ , ḡ). If they are Lie–Bäcklund equivalent, then there exist an endogenous state feedback of type (2.56) and a traditional coordinate transformation (diffeomorphism) which transfer the system ẋ = f (x, u) into the system prolonged with finite number of integrators: ẏ = g(y, w), w (r+1) = v. 2.1.3.4 Flat Output Variable The system (M̄, f¯) is differentially flat if it is Lie–Bäcklund equivalent to the trivial system given by the vector field m f¯triv = (j +1) yi j ≥0 i=1 ∂ (j ) ∂yi . (2.58) 2.1 Nonlinear System Classes 29 The variables of the vector y (their number is equal to the number of inputs of (M̄, f¯)) will be called the flat outputs. Corollary Differentially flat systems can be linearized by dynamic state feedback and a coordinate transformation since (2.58) is arisen from a linear system prolonged with a series of integrators. According to this property, y is often called the linearizing output. Theorem 2.2 Differentially flat system can be characterized by the flat output variable. Consider the differentially flat system (M̄, f¯) arisen from the prolongation of the system ẋ = f (x, u). Denote Ψ the Lie–Bäcklund transformation which transfers (M̄, f¯) into a trivial system and denote Φ the inverse transformation. It follows from the definition of differential flatness that there is a finite q such that   x = φ1 y, ẏ, ÿ, . . . , y (q) ,   u = φ2 y, ẏ, ÿ, . . . , y (q+1) , (2.59) and x and u identically satisfy the equation ẋ = f (x, u) describing the system. Furthermore, there is a finite p yielding   y = ψ1 x, u, u̇, . . . , u(p) . (2.60) According to (2.59) and (2.60), the path of x and u can be determined from y and its derivatives and there is a one-to-one correspondence between the integral curves t → (x(t), u(t)) satisfying ẋ = f (x, u) and the properly regular time functions t → y(t). It means that in case of differentially flat systems there exists a set of variables inside of which the differential constraints defining the system disappear. The system (2.51) is differentially flat around a point ⇔ it is possible to find an output variable y, dim y = dim u = m, which smoothly depends around that point on x, u, and a finite number of derivatives of u, such that the resulting square system is left and right input–output invertible with trivial zero dynamics [39]. Corollary In order to admit the differential flatness property, it is “enough” to find the output y and the functions φ1 , φ2 and ψ1 appearing in (2.59) and (2.60). 2.1.3.5 Path Design for Differentially Flat Systems Beside the linearizability (which not always can be exploited because many times not all the components of the state vector can be measured), differential flatness can play important role for path design. Namely, the path design can be performed in the space of the variables of y. Then the path of the state and the inputs can simply be determined using the time derivatives according to (2.59), without the integration 30 2 Basic Nonlinear Control Methods of the differential equations. Hence, during path design it is sufficient to choose a properly regular path satisfying the goals of the control in the space of the variables of y using simple interpolation techniques (Lagrange polynomials etc.). 2.1.3.6 Application Possibilities for Different System Classes Recently, there are strong necessary and sufficient conditions to decide whether the system is differentially flat or not (see, for example, [80, Theorem 6.7]), however the steps to be performed are too complicated for general nonlinear systems. Exceptions are SISO nonlinear systems for which necessary and sufficient conditions are available [56] and the linearizing feedback remains static (if it exists). There are also necessary and sufficient conditions for some kind of driftless systems [90]. For people dealing with robotics, it has been known for a long time that open chain, rigid and in all degrees of freedom actuated mechanical systems can be linearized using static nonlinear feedback (called nonlinear decoupling in robotics) which, as mentioned above, is equivalent to the differential flatness property. On the other hand, there exist also differentially flat mechanical systems that are underactuated or nonholonomic. A counter example is the inverted pendulum which is not differentially flat. 2.2 Dynamic Model of Simple Systems Some simple examples will be shown for finding the dynamic model of simple mechanical systems using physical laws. We try to apply as simple knowledge as possible in the discussions. For example, the kinetic energy and the potential energy of a mass point are K = 21 mv 2 and P = mgh, respectively, where g is the gravity acceleration and h is the hight of the mass point over the horizontal plane. Similarly, the kinetic energy of a mass point rotating around an axis is K = 12 Θω2 where ω is the angular velocity, Θ = ml 2 is the inertia moment to the rotation axis and l is the length of the massless rod containing the mass point at its end. For a bit more complicated case of a single rigid body, we can apply the results in Appendix A. For example, for a rigid body having velocity vc , angular velocity ω and inertia matrix Ic belonging to the center of mass the kinetic energy is K = 21 vc , vc m + 21 Ic ω, ω. However, this result was also derived from the above relations because the rigid body is the limit of infinite many point masses. From the physical laws, we shall apply the Lagrange equations, see also Appendix A. 2.2.1 Dynamic Model of Inverted Pendulum Consider a pendulum mounted on a moving car (see Fig. 2.4). 2.2 Dynamic Model of Simple Systems 31 Fig. 2.4 The model of the inverted pendulum Let L be the length of the rod mounted on the car and M the mass of the rod then the inertial moment of the rod belonging to the center of mass: ! M l 3 !!L M 3 ML2 M Θ= dl = 2L = . = l dm = l 2L 3 !−L 6L 3 −L −L 2L  L 2  L 2 (2.61) The generalized coordinates are x and ϕ. The coordinates of the center of mass of the rod are xM and yM , its velocity is vM . For simplicity, let us introduce the usual notations in robotics: Sϕ = sin(ϕ), Cϕ = cos(ϕ). Then xM = x + LSϕ yM = LCϕ ⇒ ⇒ ẋM = ẋ + ϕ̇LCϕ , ẏM = −ϕ̇LSϕ , 2 2 2 vM = ẋM + ẏM = ẋ 2 + 2LCϕ ẋ ϕ̇ + L2 Cϕ2 ϕ̇ 2 + L2 Sϕ2 ϕ̇ 2 (2.62) = ẋ 2 + 2LCϕ ẋ ϕ̇ + L2 ϕ̇ 2 . The kinetic energy K and the potential energy P are respectively, 1 K = mẋ 2 + 2 1 = mẋ 2 + 2 P = MgLCϕ . 1 1 2 MvM + Θ ϕ̇ 2 2 2   1 1 M ẋ 2 + 2LCϕ ẋ ϕ̇ + L2 ϕ̇ 2 + Θ ϕ̇ 2 , 2 2 (2.63) 32 2 Basic Nonlinear Control Methods The derivatives, needed to Lagrange equations, are as follows: ∂K ∂ ẋ ∂K ∂x ∂P ∂x ∂K ∂ ϕ̇ = mẋ + M ẋ + MLCϕ ϕ̇, = 0, = 0, (2.64) = MLCϕ ẋ + ML2 ϕ̇ + Θ ϕ̇, ∂K = −MLSϕ ẋ ϕ̇, ∂ϕ ∂P = −MgLSϕ , ∂ϕ d dt d dt ∂K = (m + M)ẍ + MLCϕ ϕ̈ − MLSϕ ϕ̇ 2 , ∂ ẋ   ∂K = MLCϕ ẍ − MLSϕ ẋ ϕ̇ + Θ + ML2 ϕ̈. ∂ ϕ̇ (2.65) If the car is pulled with force F , then the Lagrange equations yield: d ∂K ∂K ∂P − + dt ∂ ẋ ∂x ∂x =F (m + M)ẍ + MLCϕ ϕ̈ − MLSϕ ϕ̇ 2 = F, ⇒ (2.66) d ∂K ∂K ∂P − + dt ∂ ϕ̇ ∂ϕ ∂ϕ =0 ⇒   MLCϕ ẍ + Θ + ML2 ϕ̈ − MgLSϕ = 0. (2.67) From here, it can be obtained firstly ẍ = F − MLCϕ ϕ̈ + MLSϕ ϕ̇ 2 , m+M then after some algebraical manipulations F − MLCϕ ϕ̈ + MLSϕ ϕ̇ 2 4 + ML2 ϕ̈ − MgLSϕ = 0, m+M 3   2 MCϕ Cϕ MLSϕ Cϕ ϕ̇ 2 4 ϕ̈ = gSϕ − L − F− , 3 m+M m+M m+M MLCϕ (2.68) 2.2 Dynamic Model of Simple Systems 33 Fig. 2.5 Simplified model of a car active suspension and finally ϕ̈ = gSϕ − Cϕ m+M (F L( 43 − + MLSϕ ϕ̇ 2 ) MCϕ2 m+M ) . (2.69) The state vector is (x, ẋ, ϕ, ϕ̇)T and the input signal is F . The two output signals are x and ϕ. The system is underactuated because there are two output signals and only one input signal. 2.2.2 Car Active Suspension Model Let us consider a vehicle (car, bus etc.) divided into 4 parts in case of 4 wheels, each part is a composite mechanical spring-damper system consisting of the quarter part of the mass of the body (together with the passengers) and the mass of the wheel. The vertical positions are described by upward directed x1 and x2 , see Fig. 2.5. The distance between the road surface and the wheel’s contact point is the disturbance w varying together with the road surface. The connections between the different parts are flexible modeled by viscous friction and spring constant. The suspension is active which means that the actuator produces the force F (control signal). A good suspension system should have satisfactory road holding stability, while providing good traveling comfort when riding over bumps and holes in the road. When the car is experiencing any road disturbance, it should not have large oscillations and the oscillations should dissipate quickly. 34 2 Basic Nonlinear Control Methods Denote m1 and m2 the mass of the (quarter) body and the wheel respectively, while the flexible connections are described by the viscous damping factors (b1 , b2 ) and the spring constants (k1 , k2 ). The motion equations, using also the actionreaction (d’Alambert) principle, are as follows: m1 ẍ1 = F − b1 (ẋ1 − ẋ2 ) − k1 (x1 − x2 ), m2 ẍ2 = −F + b1 (ẋ1 − ẋ2 ) + k1 (x1 − x2 ) − b2 (ẋ2 − ẇ) − k2 (x2 − w). As usual, the second derivatives will be expressed: ẍ1 = − b1 k1 1 (ẋ1 − ẋ2 ) − (x1 − x2 ) + F, m1 m1 m1 b1 k1 b2 k2 1 (ẋ1 − ẋ2 ) + (x1 − x2 ) − (ẋ2 − ẇ) − (x2 − w) − F, m2 m2 m2 m2 m2 y1 := x1 − x2 , ẍ2 = ẍ1 = − ẍ2 = b1 k1 1 ẏ1 − y1 + F, m1 m1 m1 k1 b2 k2 1 b1 ẏ1 + y1 − (ẋ2 − ẇ) − (x2 − w) − F, m2 m2 m2 m2 m2 ÿ1 = ẍ1 − ẍ2     b1 b2 k1 k2 ẏ1 − y1 =− + + m1 m2 m1 m2   k2 1 1 b2 F. (ẋ2 − ẇ) + (x2 − w) + + + m2 m2 m1 m2 In order to eliminate ẇ, the equation of ÿ1 will be integrated:   b1 b2 b2 ẏ1 = − y1 + + (x2 − w) m1 m2 m2   #   "  k1 k2 k2 1 1 − y1 + F dt . + + (x2 − w) + + m1 m2 m2 m1 m2 $ %& ' y2 Denote y2 the integral and let us differentiate it:     k1 k2 k1 1 1 ẏ2 = − y1 + F. + (x2 − w) + + m1 m2 m2 m1 m2 Perform the substitution x2 = x1 − y1 :     k2 k2 k1 k2 1 1 k1 y1 + F− ẏ2 = x1 − + + + w, m2 m1 m2 m2 m1 m2 m2 2.2 Dynamic Model of Simple Systems 35   b1 b2 b2 b2 b2 x1 − + + w, y1 + y2 − m2 m1 m2 m2 m2 # "   k1 b1 b1 b2 b2 1 − y1 + ẍ1 = − + (x1 − y1 − w) + y2 − y1 + F m1 m1 m2 m2 m1 m1  # "  k1 b1 b 2 b1 b1 b2 1 b1 b2 b1 b1 − y1 − = x1 + + + y2 + F+ w. m1 m2 m1 m1 m2 m2 m1 m1 m1 m1 m2 ẏ1 = State equation: x̄˙ = Ax̄ + Bu, y = C x̄, u = (F, w)T , x̄ := (x1 , ẋ1 , y1 , ẏ1 ) , ⎡ 0 1 0 ⎢− b1 b2 0 b1 ( b1 + b1 + b2 ) − ⎢ m1 m2 m1 m1 m2 m2 A=⎢ b1 b1 ⎢ b2 0 −( m1 + m2 + mb22 ) ⎣ m2 k2 m2 T 0 ⎡ −( mk11 + ⎤ k1 m2 + k1 m1 k2 m2 ) 0 0 ⎢ b1 b2 ⎥ 1 ⎢ m1 m1 m2 ⎥ ⎥, B =⎢ ⎢ 0 − mb22 ⎥ ⎣ ⎦ k2 1 1 ( m1 + m2 ) − m2   C= 0 0 1 0 . y = x1 − x2 = y1 , ⎤ 0 − mb11 ⎥ ⎥ ⎥, 1 ⎥ ⎦ 0 Some typical values for a bus in SI units: m1 = 2500, b1 = 350, m2 = 320, b2 = 15020, k1 = 80000, k2 = 500000. Typical requirements in case of bus: In order to assure the traveling comfort, settling time shorter than 10 s and overshoot less that 5% is usually necessary. For example, when the bus runs onto a 10 cm high step the body will oscillate within a range of ±5 mm (i.e., ±0.005 m) and should return to smooth ride within 5 s. 2.2.3 The Model of the 2 DOF Robot Arm Consider the two degree of freedom robot arm moving in the vertical plane, see the block scheme in Fig. 2.6. The generalized coordinates are q = (q1 , q2 )T := (ϕ1 , ϕ2 )T and the generalized forces are the driving torques τ := (τ1 , τ2 )T . The position vectors pointing from 36 2 Basic Nonlinear Control Methods Fig. 2.6 The model of a 2 DOF robot moving in vertical plane the origin of the fixed frame (coordinate system) K0 into the origins of the moving frames K1 and K2 are respectively, p01 = (l1 C1 , l1 S1 , 0)T , p02 = (l1 C1 + l2 C12 , l1 S1 + l2 S12 , 0)T , (2.70) where the notations C1 := cos(q1 ), . . . , S12 = sin(q1 + q2 ) are used. The velocities in the basis of the fixed frame K0 can be obtained by formal differentiation of the position vectors: ⎡ ⎤ −l1 S1 0 v1 = ⎣ l1 C1 ⎦ q̇1 , 0 ⎡ −l1 S1 − l2 S12 0 v2 = ⎣ l1 C1 + l2 C12 0 ⎤ −l2 S12   q̇1 l2 C12 ⎦ . q̇2 0 (2.71) The inertia matrix is constant in the frame of the moving body but is nonconstant in the fixed frame, hence we transform the velocities into the basis of the moving 2.2 Dynamic Model of Simple Systems 37 frame: ⎡ ⎤ ⎤ c0 C1 S 1 0 v1 := 1 v1 = ⎣−S1 C1 0⎦ 0 v1 = ⎣ l1 ⎦ q̇1 , 0 0 0 1 ⎡ ⎡ ⎤ C12 S12 0 l1 S 2 v2 := 2 v2 = ⎣−S12 C12 0⎦ 0 v2 = ⎣l1 C2 + l2 0 0 1 0 ⎡ ⎤ 0   q̇1 l2 ⎦ . q̇2 0 (2.72) If each moving frame will be parallelly transferred to the origin of the center of mass, then the velocity (vc ) and the inertia matrix belonging to the center of mass have to be used in the motion equation. In the block scheme, I1 and I2 are already the right lower elements of the inertia matrices Ic1 and Ic2 , respectively. The velocities vc1 , vc2 can be immediately written down using simple geometrical considerations: ⎡ 0 ⎤ vc1 = ⎣lc1 ⎦ q̇1 , 0 ⎡ ⎤ l1 S 2 0   q̇1 . vc2 = ⎣l1 C2 + lc2 lc2 ⎦ q̇2 0 0 (2.73) From the geometrical image, the angular velocities are simply ω1 = (0, 0, q̇1 )T , ω2 = (0, 0, q̇1 + q̇2 )T , which can also be written in matrix-vector form: ⎡ ⎤ 0 ω1 = ⎣0⎦ q̇1 , 1 ⎤ ⎡ 0 0   q̇1 . ω2 = ⎣0 0⎦ q̇2 1 1 (2.74) Hence, the kinetic and potential energy of the robot arm are respectively, 1 1 vc1 , vc1 m1 + I1 ω1 , ω1  2 2 1 1 + vc2 , vc2 m2 + I2 ω2 , ω2 , 2 2 P = m1 glc1 S1 + m2 g(l1 S1 + lc2 S12 ). K= (2.75) 38 2 Basic Nonlinear Control Methods After simple algebraic manipulations, the kinetic energy can be formulated in the quadratic form K= ) 1 1( H (q)q̇, q̇ = 2 2 Dj k q̇j q̇k , j k where D11 (q) D12 (q) D12 (q) , D22 (q)   2 2 D11 = m1 lc1 + m2 l12 + lc2 + 2l1 lc2 C2 + I1 + I2 , H (q) = (2.76) D12 = m2 (l1 C2 + lc2 )lc2 + I2 , 2 D22 = m2 lc2 + I2 . The dynamic model of the robot arm can be found using the Lagrange equations. Performing the differentiations and introducing the notations D112 (q) = −m2 l1 lc2 S2 , D1 (q) = ∂P = m1 glc1 C1 + m2 g(l1 C1 + lc2 C12 ), ∂q1 D2 (q) = ∂P = m2 glc2 C12 , ∂q2 (2.77) we obtain the following dynamic model of the two degree of freedom robot arm:   D11 (q)q̈1 + D12 (q)q̈2 + D112 (q) 2q̇1 q̇2 + q̇22 + D1 (q) = τ1 , D12 (q)q̈1 + D22 (q)q̈2 − D112 (q)q̇12 + D2 (q) = τ2 . (2.78) The state vector may be x := (q1 , q2 , q̇1 , q̇2 )T , the control signal (the input of the model) is u = (τ1 , τ2 )T . Since, from physical considerations, the matrix H (q) is positive definite, thus invertible, therefore q̈1 and q̈2 can be expressed from (2.77) and (2.78). Hence, the state equation can be written in the input affine form ẋ = f (x) + G(x)u. (2.79) 2.3 Stability of Nonlinear Systems The stability of nonlinear systems will be discussed in state space, in signal spaces and in the class of dissipative systems. The discussion is based on the classical work of Lyapunov [86] and more recent results of LaSalle and Lefschetz [77], Vidyasagar [152], Khalil [62] and van der Schaft [127]. 2.3 Stability of Nonlinear Systems 39 Fig. 2.7 Illustration of Lyapunov stability definition 2.3.1 Stability Definitions Let the state equation of the time varying nonlinear system be ẋ = f (t, x, u), y = g(t, x, u). (2.80) Assuming u(t) = 0 or u(t) = u0 (t) known input signal, the nonlinear system can be brought on the form ẋ = f (t, x). In what follows we shall assume f ∈ (0,1) Ct,x ([a, ∞) × Dx ), that is, f is continuous in t and differentiable in x, furthermore a ∈ R 1 and Dx ⊂ R n is open, for example a = 0 and Dx = R n . Although different stability definitions exist for nonlinear systems, the most widespread is the Lyapunov stability definition. Lyapunov stability Let ξ(t) be the solution of the differential equation ẋ = f (t, x). We say that the solution ξ(t) is stable in Lyapunov sense, if for every t0 ∈ [a, ∞) and for every ε > 0 there exists δ(ε, t0 ) > 0 such that if x(t0 ) − ξ(t0 ) < δ(ε, t0 ) and x(t) is the solution then x(t) − ξ(t) < ε for every t ∈ [t0 , ∞). The definition can also be illustrated so that small perturbing the initial condition everywhere with respect to ξ(t) and starting the system from the perturbed initial condition the new solution remains in the neighborhood of ξ(t). A clear interpretation is shown in Fig. 2.7. The solution ξ(t) may be an equilibrium point of the system, a limit cycle (periodical solution) or any other solution of the state equation (for example during chaos investigations). 40 2 Basic Nonlinear Control Methods The following terminology is usual: (i) The solution ξ(t) of the system is unstable if ξ(t) is not stable. (ii) The solution ξ(t) of the system is uniformly stable if ξ(t) is stable in Lyapunov sense and δ(ε, t0 ) = δ(ε) does not depend on t0 . (iii) The solution ξ(t) of the system is asymptotically stable if ξ(t) is stable in Lyapunov sense and for every t0 ∈ [a, ∞) there exists δ1 (t0 ) > 0 such that if x(t0 ) − ξ(t0 ) < δ1 (t0 ) and x(t) is the solution then limt→∞ x(t) − ξ(t) = 0. The definition of (for the practice particularly important) asymptotic stability means that small perturbing everywhere with respect to ξ(t) and starting the system from the perturbed initial condition then the new solution remains in the neighborhood of ξ(t) and actually approaches ξ(t) if t → ∞. We shall show that the stability investigation can be reduced to the stability investigation of ξ ≡ 0 (equilibrium point) of a modified system satisfying a special condition. Namely, let us introduce the notation x̃(t) := x(t) − ξ(t), that is, x(t) = ξ(t) + x̃(t), then       d x̃(t) dx(t) dξ(t) = − = f t, ξ(t) + x̃(t) − f t, ξ(t) =: f˜ t, x̃(t) , dt dt dt which means that the new form of the stability problem   d x̃(t) = f˜ t, x̃(t) , dt ξ̃ (t) ≡ 0 (equilibrium point) (2.81) ⇔ f˜(t, 0) ≡ 0. However, it should be stressed that the new system has to satisfy the condition f˜(t, 0) ≡ 0. Henceforth, we assume that the transformation has already been performed and omit the notation “ ˜ ”. The stability is a local property in the definition, see δ(ε, t0 ) > 0. It does not say anything directly about the region of attraction for the equilibrium point ξ ≡ 0. Further results are needed to answer what is the region around the equilibrium point from which choosing the initial condition the trajectory (solution) remains in the neighborhood of the equilibrium point or even more it does approach the equilibrium point. We shall speak about global stability if the domain of attraction is the whole space. According to these, we shall use the abbreviations GAS, UGAS etc. meaning that the system is globally asymptotically stable, uniformly globally asymptotically stable etc. 2.3.2 Lyapunov Stability Theorems In the stability investigations, the idea of positive definite function plays an important role, in the applications it serves as Lyapunov function candidate. 2.3 Stability of Nonlinear Systems 41 (i) The scalar valued function V (t, x) is positive definite, if there exists a scalar valued function W (x) ∈ C such that V (t, x) ≥ W (x) > 0 for every x =  0 and V (t, 0) ≡ W (0) = 0. (ii) V (t, x) is negative definite if −V (t, x) is positive definite. 2.3.2.1 Lyapunov’s Direct Method Theorem 2.3 Lyapunov’s first theorem, direct method. Let ξ ≡ 0 an equilibrium point, that is, dx = f (t, x), dt f (t, 0) ≡ 0, (0,1) f (t, x) ∈ Ct,x . (2.82) (1,1) 1. If ∃V (t, x) ∈ Ct,x ([a, ∞) × Dx ) positive definite function (so called Lyapunov function) such that for every x(t) solution along the trajectory yields V̇ (t, x) := dV (t, x(t)) ≤ 0 (V̇ negative semidefinite), dt (2.83) then the equilibrium point ξ ≡ 0 is stable in Lyapunov sense. 2. If additionally V̇ (t, x) := dV (t, x(t)) < 0 (V̇ negative definite), dt (2.84) then the equilibrium point ξ ≡ 0 is asymptotically stable in Lyapunov sense. Notice that dV ∂V dx ∂V ∂V ∂V = + = f (t, x) + . (2.85) dt ∂x dt ∂t ∂x ∂t The following simple proof can seldom be found in control books for this general time variant case. The proof is based on the fact that any continuous function over a compact (bounded and closed) set assumes its extremum value. Proof (Stability part) Let ε > 0 and consider the spherical surface Sε = {x : x = ε} that is evidently compact. Then W (x) ∈ C assumes its minimum over Sε , that is, ∃x ∗ ∈ Sε , such that infx∈Sε W (x) = W (x ∗ ) = α > 0. Let t0 ∈ (a, ∞), then, since V (t, x) is continuous and V (t, 0) = 0, hence ∃0 < δ(t0 , ε) < ε such that ∀{x : x < δ} yields 0 ≤ V (t, x) < α. Let x(t0 ) < δ and x(t) (not identically zero) solution of the state equation. It can be shown that the trajectory x(t) remains in the interior of the spherical surface Sε : x(t) < ε, ∀t ∈ [t0 , ∞). Otherwise, it would exist t1 ∈ [t0 , ∞), such that x(t0 ) < δ < ε, x(t) < ε, ∀t < t1 , but x(t1 ) = ε. This leads to contradiction because dv/dt := V̇ (t, x(t)) ≤ 0, thus v(t) is monotone nonincreasing along the trajectory, hence α ≤ W (x(t1 )) ≤ 42 2 Basic Nonlinear Control Methods V (t1 , x(t1 )) ≤ V (t0 , x(t0 )) < α. Therefore, x(t) < ε, ∀t ∈ [t0 , ∞), and thus the equilibrium point ξ ≡ 0 is stable.  Proof (Asymptotic stability part) We already know that ξ ≡ 0 is stable. Let x(t0 ) < δ1 (t0 ) sufficiently small and x(t) (not identically zero) solution. It has to be proven that limt→∞ x(t) = 0. Now we know that dV /dt < 0 (negative definite) along the trajectory. First, we shall prove that, using the notation limt→∞ V (t, x(t)) = inft V (t, x(t)) =: α, it cannot be satisfied α > 0. If, in contradiction, we would assume α > 0, then it would follow from α > 0 that x(t) ≥ β > 0, ∀t ∈ [t0 , ∞). Otherwise, that is, in case of β = 0, there would exist a sequence {t1 , t2 , . . . , tk , . . .} such that limtk →∞ V (tk , x(tk )) = 0, which would be a contradiction because of limt→∞ V (t, x(t)) = α > 0. Hence, α > 0 ⇒ x(t) ≥ β > 0, ∀t ∈ [t0 , ∞). After this, we shall show that α > 0 leads also to contradiction. Namely, we know that the equilibrium point is stable, thus we can assume there exists 0 < β ≤ x(t) < ε, ∀t ∈ [t0 , ∞). Since V̇ (t, x) < 0 (negative definite), there exists W1 (x) ∈ C positive definite and continuous function such that v̇(t) := V̇ (t, x(t)) ≤ −W1 (x(t)) < 0. Let γ := infβ≤x<ε W1 (x), then γ > 0 and by integration  t  t   v(t) = v(t0 ) + W1 (τ ) dτ ≤ v(t0 ) − γ (t − t0 ). V̇ τ, x(τ ) dτ ≤ v(t0 ) − t0 t0 Hence, ∃t, v(t) < 0 which is a contradiction since V is positive definite, and thus V (t, x(t)) ≥ 0 along the trajectory. Hence, in contradiction to the original assumption, α = 0 and therefore limt→∞ V (t, x(t)) = 0. It has yet to be shown that also the trajectory x(t) approaches zero. But this follows from the fact that V is positive definite, W ∈ C is continuous and positive definite, V (t, x(t)) ≥ W (x(t)) ≥ 0, thus 0 = limt→∞ V (t, x(t)) ≥ limt→∞ W (x(t)) ≥ 0 ⇒ limt→∞ W (x(t)) = 0, and since W ∈ C is positive definite thus limt→∞ x(t) = 0. Hence, the equilibrium point ξ ≡ 0 is asymptotically stable in Lyapunov sense.  In order to use Lyapunov’s direct method, an appropriate Lyapunov function candidate has to be found from which we can deduce the stability of the equilibrium point. Positive functions like norm, energy functions etc. can be suggested as Lyapunov function candidates. In some simple cases other possibilities exist, for example, the question arises whether is it possible to deduce the stability of the nonlinear system from the stability of the linearized system. 2.3.2.2 Lyapunov’s Indirect Method Theorem 2.4 Lyapunov’s second theorem, indirect method. Let ξ ≡ 0 be an equilibrium point, i.e. dx dt = f (t, x), f (t, 0) ≡ 0, f (t, x) ∈ (0,1) Ct,x , and let us introduce the notations A(t) := fx′ (t, x)|x=0 and f1 (t, x) := f (t, x) − A(t)x. 2.3 Stability of Nonlinear Systems 43 If it is satisfied (t,x) (i) limx→0 supt≥0 f1x = 0. (ii) A(·) is bounded. (iii) ẋ = A(t)x linear system is uniformly asymptotically stable, then the equilibrium point ξ ≡ 0 of the original nonlinear system is uniformly asymptotically stable. The sketch of the proof is as follows. Consider the linearized system ẋ = A(t)x and its fundamental matrix Φ(t, τ ). Let  ∞ Φ T (τ, t)Φ(τ, t) dτ, P (t) := t *∞ then P (t)x, x = t Φ(τ, t)x2 dτ > 0. Hence, P (t)x, x is a Lyapunov function candidate for which the following conditions are satisfied: (1) P (t) is positive definite and ∃α, β > 0, such that ∀t, x : αx2 ≤ P (t)x, x ≤ βx2 , furthermore it yields Ṗ (t) + P (t)A(t) + AT (t)P (t) = −I . (2) V (t, x) := P (t)x, x is positive definite function with W (x) := αx2 . (3) V̇ (t, x) = −x2 + 2 P (t)f1 (t, x), x < 0, since 2 P (t)f1 (t, x), x can be made less than x2 in some neighborhood of x = 0 because of assumption (i) and properties (1). Notice that Lyapunov’s indirect method does not give any information about the attraction domain of the equilibrium, hence only the local stability of the equilibrium has been proven. For LTI systems, the classical stability idea is strongly related to the poles of the transfer function which are equal to the eigenvalues of the state matrix. Stability was defined for such systems so that the transients are decaying if the initial conditions are changed from zero to nonzero value. The following theorem shows the connection between classical and Lyapunov stability. Theorem 2.5 Lyapunov equation characterizing stable LTI system. Consider the continuous time linear time invariant system ẋ = Ax for which evidently ξ ≡ 0 equilibrium point. Let P be a constant, symmetrical and positive definite matrix (P > 0) and let V (x) := P x, x be the Lyapunov candidate function. Perform the substitution ẋ = Ax in the derivative of V : V̇ = P ẋ, x + P x, ẋ = (P A + AT P )x, x. Then V̇ < 0 (negative definite), if ∃Q > 0 positive definite matrix such that the Lyapunov equation is satisfied: P A + AT P = −Q, P > 0, Q > 0. (2.86) The following statements are equivalent: (i) For all eigenvalues, si of the matrix A yields Re si < 0. (ii) For every Q > 0, the Lyapunov equation has the positive definite solution P := * ∞ AT t At Qe dt > 0. 0 e The proof is simple see, for example, [152]. 44 2 Basic Nonlinear Control Methods 2.3.2.3 LaSalle’s Method for Time Invariant Systems For the nonlinear time invariant (autonomous) system class ẋ = f (x), f (0) = 0 LaSalle [77] gave a stronger result than Lyapunov’s direct method which is based on the idea of positive limit set and invariant set. Positive limit set. The positive limit set of the trajectory x(t) is Γ + := {y : ∃{tn } sequence, tn → ∞ as n → ∞, such that limn→∞ x(tn ) = y}. The positive limit set of the trajectory x(t) may be for example a limit point limt→∞ x(t) (if it exists), or a limit cycle (isolated, closed and periodical curve) to which the trajectory converges. Invariant set. The set M is invariant if each trajectory (solution of the differential equation ẋ = f (x)) going through any point of M lies entirely in M. If the set E is not invariant then we can speak about its maximal invariant set M contained in E: M = max{H ⊂ E : x(0) ∈ H, x(t) solution ⇒ ∀t : x(t) ∈ H }. The positive limit set Γ + of a bounded trajectory x(t) is always nonempty, compact (bounded and closed) and invariant. Theorem 2.6 LaSalle’s maximal invariant set theorem for autonomous nonlinear system. Assume the system ẋ = f (x) has equilibrium point 0 = f (0) and f (x) ∈ Cx(1) . Assume there exists V (x) Lyapunov function which satisfies the following conditions for some fixed r > 0 over the compact set {x : |V (x)| ≤ r} =: Ωr : (i) V (x) is positive definite, i.e. V (x) ≥ 0 and V (x) = 0 ⇔ x = 0. (ii) V̇ (x) = grad V , f  ≤ 0 (negative semidefinite). (iii) V (x) ∈ Cx(1) (continuously differentiable). Let M be the maximal invariant set contained in E := {x ∈ Ωr : V̇ (x) = 0}. Then for every initial condition x(0) ∈ Ωr the resulting trajectory x(t) asymptotically approaches the maximal invariant set M. More detailed it means that for every ε > 0 there exists T , such that for every t > T the trajectory x(t) remains in the ε-neighborhood of M: There exist y(t, ε) ∈ M + {x : x < ε} such that x(t) − y(t, ε) < ε. Notice that Ωr is an attraction set of the equilibrium point although not necessarily the largest one. Corollary If M = {0} is the maximal invariant set contained in E = {x ∈ Ωr : V̇ (x) = 0}, then for every x(0) ∈ Ωr the corresponding x(t) trajectory asymptotically approaches the equilibrium point 0 = f (0). Corollary Let x = (y T , zT )T and assume that beside the conditions of LaSalle’s theorem it is additionally satisfied: (i) ẏ = ϕ(y, z) and ż = ψ(y, z). 2.3 Stability of Nonlinear Systems 45 (ii) V̇ (x) = V̇ (y, z) = 0 ⇒ y = 0. (iii) ϕ(0, z) = 0 ⇔ z = 0. Then the equilibrium point y = 0, z = 0 is asymptotically stable. It is evident namely, that according to (ii) and the definition of E in each point of M ⊂ E yields y = 0. Since M is an invariant set and y is constant zero over it, therefore ẏ = 0 for each trajectory running in M, and thus by (iii) yields also z = 0. Hence, M consists only of a single point and any trajectory x(t) starting in Ωr asymptotically approaches the equilibrium point. 2.3.2.4 K and KL Comparison Functions for Stability Investigations In order to sharpen the stability theorems in the direction of global asymptotic stability, some further definitions are necessary. (i) The function α : [0, a) → R + is called K-function if α(·) is strictly monotone increasing and α(0) = 0. Especially, if a = ∞ and α(r) → ∞ if r → ∞ then α(r) is a K∞ -function. (ii) The function β : [0, a) × R + → R + is called KL-function if β(r, s) satisfies that β(·, s) is a K-function in the first variable for every fixed s, while β(r, ·) is monotone decreasing in the second variable for every fixed r and β(r, s) → 0 if s → ∞. It is clear that, for example, α(r) = atan(r) is a K-function, but not a K∞ function, however α(r) = r c is a K∞ -function for c > 0, while β(r, s) = r c e−s is a KL-function for c > 0. The following theorems are valid for ẋ = f (t, x), the proofs can be found in [62]. Theorem 2.7 Relations between stability and comparison functions. The equilibrium point ξ ≡ 0 is (i) Uniformly stable ⇔ there exist α(r) K-function and c > 0 independent of t0 such that + + + + +x(t)+ ≤ α +x(t0 )+ , + + ∀t ≥ t0 ≥ 0, ∀+x(t0 )+ < c. (ii) Uniformly asymptotically stable ⇔ there exist β(r, s) KL-function and c > 0 independent of t0 such that + + +  + +x(t)+ ≤ β +x(t0 )+, t − t0 , + + ∀t ≥ t0 ≥ 0, ∀+x(t0 )+ < c. + + +  + +x(t)+ ≤ β +x(t0 )+, t − t0 , ∀t ≥ t0 ≥ 0, ∀x(t0 ) ∈ R n . (iii) Globally uniformly asymptotically stable ⇔ there exist β(r, s) KL-function such that 46 2 Basic Nonlinear Control Methods Theorem 2.8 Relations between Lyapunov function and qualified stability. Let D ⊂ R n be an open set, 0 ∈ D, V : [0, ∞) × D → R continuously differentiable function satisfying W1 (x) ≤ V (t, x) ≤ W2 (x), where W1 (x), W2 (x) are positive definite and continuous, ∂V ∂V f (t, x) + ≤ 0. ∂x ∂t Then the equilibrium point ξ ≡ 0 is (i) Uniformly stable. ∂V (ii) Uniformly asymptotically stable if ∂V ∂x f (t, x) + ∂t ≤ −W3 (x), where W3 (x) is positive definite and continuous. (iii) If beside the previous condition it is also satisfied that r and c are chosen so that Br := {x : x ≤ r} ⊂ D and 0 < c < minx=r W1 (x), then every trajectory starting in {x ∈ Br : W2 (x) ≤ c} satisfies + + + +  +x(t)+ ≤ β +x(t0 )+, t − t0 , ∀t ≥ t0 ≥ 0 for some class KL-function. Specifically, if D = R n and W1 (x) is radially unbounded then the equilibrium is globally uniformly asymptotically stable. Remark The proof of the strong result iii) needs the construction of the class K-functions α1 (x) ≤ W1 (x), W2 (x) ≤ α2 (x), α3 (x) ≤ W3 (x), then α = α3 ◦ α2−1 , y(t) =: σ (y0 , t − t0 ) is the solution of the scalar nonlinear differential equation ẏ = −α(y), y(t0 ) = V (t0 , x(t0 )) and finally β(x(t0 ), t − t0 ) := α1−1 (σ (α2 (x(t0 )), t − t0 )). Till now, the input signal has not appeared in the stability investigations, hence the discussion should be extended to the input-state stability of nonlinear systems. Consider the nonlinear system ẋ = f (t, x, u). We shall assume that the unexcited system has equilibrium point ξ ≡ 0, that is, f (t, 0, 0) ≡ 0, and we shall investigate the influence of the input signal u(t) on the stability. 2.3.2.5 Input-to-State Stability (ISS) The nonlinear system ẋ = f (t, x, u) is called input-to-state stable (ISS), if there exist β(r, s) class KL-function and γ (r) class K-function such that for every initial condition x(t0 ) and every bounded input u(t) there exits the solution x(t), ∀t ≥ t0 satisfying the following inequality:  + + + + + +  +x(t)+ ≤ β +x(t0 )+, t − t0 + γ sup +u(τ )+ . t0 ≤τ ≤t 2.3 Stability of Nonlinear Systems 47 Notice that the above condition assures the boundedness of the trajectory x(t) for every bounded input u(t). If t is increasing the solution x(t) is uniformly bounded by a K-function, in the argument of which stands supt≥t0 u(t), furthermore if t → ∞ yields u(t) → 0 then x(t) → 0 is also valid. Since for u(t) ≡ 0 the inequality reduces to x(t) ≤ β(x(t0 ), t − t0 ), hence the equilibrium point ξ(t) ≡ 0 of the unexcited system is globally uniformly asymptotically stable. The following theorem gives a condition of ISS, the proof can be found in [62]. Theorem 2.9 Sufficient condition for ISS stability. Let V : [0, ∞) × R n → R be a continuously differentiable function satisfying the conditions     α1 x ≤ V (t, x) ≤ α2 x , ∂V ∂V f (t, x, u) + ≤ −W3 (x), ∂x ∂t   ∀x ≥ ρ u > 0 for every (t, x, u) ∈ [0, ∞) × R n × R m , where α1 , α2 are class K∞ -functions, ρ is class K-function and W3 (x) is positive definite over R n . Then ẋ = f (t, x, u) is input-to-state stable with the choice γ := α1−1 ◦ α2 ◦ ρ. In the expressions, −1 denotes the inverse function and circle denotes the composition of functions. 2.3.3 Barbalat Lemmas Functions like sin(t), e−at , log(t) etc. and combinations thereof often appear in dynamic system responses. Unfortunately, it is not an easy thing to deduce system behavior from such signals. Against typical engineering feeling, if f˙ → 0 it does not follow that f (t) has a limit as t → ∞. For example, f (t) = sin(log(t)) has no limit as t → ∞ while f˙(t) = cos(log(t)) → 0 as t → ∞. Similarly, the fact that t limt→∞ f (t) exists, does not imply that limt→∞ f˙(t) = 0. For example, f (t) = e−t sin(e2t ) tends to zero while its derivative f˙(t) = −e−t sin(e2t ) + 2et cos(e2t ) is unbounded. On the other hand, if f (t) is lower bounded and decreasing (f˙(t) ≤ 0), then it converges to a limit. To find reliable results, consider first the definition of the limit value of a function. We know from standard calculus that limt→∞ f (t) = a if for every ε > 0 there exists T (ε) > 0 such that for all t ≥ T follows |f (t) − a| < ε. As a consequence, if limt→∞ f (t) = 0, then there exists a ε0 > 0 such that for all T > 0 we can find a T1 ≥ T satisfying |f (T1 )| ≥ ε0 . Uniformly continuous functions are especially important. We call f (t) uniformly continuous if for every ε > 0 there exists δ(ε) > 0 such that |f (t2 ) − f (t1 )| < ε if |t2 − t1 | < δ. According to the mean value theorem, if the function f : R → R is differentiable, then f (t2 ) − f (t1 ) = f˙(τ )(t2 − t1 ) where τ ∈ (t1 , t2 ). If f˙ ∈ L∞ , that is, f˙ is essentially bounded, then f is uniformly continuous and δ can be chosen ε/f˙∞ . 48 2 Basic Nonlinear Control Methods Lemma 2.1 (Barbalat lemmas) *t 1. Let f : R → R be uniformly continuous on [0, ∞). Suppose that limt→∞ 0 f (τ ) dτ exists and is finite. Then limt→∞ f (t) * t = 0. 2. Suppose that f ∈ L∞ , f˙ ∈ L∞ and limt→∞ 0 f (τ ) dτ exists and is finite. Then limt→∞ f (t) = 0. 3. If f (t) is differentiable and has finite limit limt→∞ f (t) < ∞, and if f˙ is uniformly continuous, then f˙(t) → 0 as t → ∞. Proof 1. If the first statement is not true, then there is a ε0 > 0 such that for all T > 0 we can find T1 ≥ T satisfying |f (T1 )| ≥ ε0 . Since f (t) is uniformly continuous, there is a positive constant δ > 0 such that |f (t + τ ) − f (t)| < ε0 /2 for all t > 0 and 0 < τ < δ. Hence, ! ! ! ! ! ! ! ! !f (t)! = !f (t) − f (T1 ) + f (T1 )! ≥ !f (T1 )! − !f (t) − f (T1 )! > ε0 − ε0 /2 = ε0 /2 for every t ∈ [T1 , T1 + δ]. Therefore, ! T +δ !  ! 1 ! ! != f (t) dt ! ! T1 T1 +δ ! T1 ! !f (t)! dt > ε0 δ/2 where the equality * t holds because f (t) does not change the sign for t ∈ [T1 , T1 + δ]. Thus, 0 f (t) dt cannot converge to a finite limit as t → ∞, which is a contradiction. Hence, the first statement is valid. 2. The second statement follows from the fact that if f˙ ∈ L∞ then f is uniformly continuous and the first statement can be applied. 3. The third statement follows from the fact that because of  ∞ f˙(t) dt f (∞) − f (0) = 0 the right side exists and is finite. Since f˙ is uniformly continuous, we can apply the first statement with f˙ instead of f .  The following lemma gives an immediate corollary of Barbalat’s lemma, which looks like an invariant set theorem in Lyapunov stability analysis of time varying systems. Lemma 2.2 (Lyapunov-like form of Barbalat’s lemma) If a scalar valued function V (x, t) satisfies the following three conditions: (i) V (x, t) is lower bounded. (ii) V̇ (x, t) is negative semidefinite. (iii) V̇ (x, t) is uniformly continuous in t. Then V̇ (x, t) → 0 as t → ∞. 2.3 Stability of Nonlinear Systems 49 Observe that V (x, t) can simply be a lower bounded function instead of being positive definite. Notice that V̇ (x, t) is uniformly continuous if V̈ (x, t) exists and is bounded. 2.3.4 Stability of Interconnected Passive Systems Stability problems can be considered also in function spaces for fixed (typically zero) initial conditions. For p ∈ [1, ∞), * ∞ the Lebesgue measurable function f : [0, ∞) → R 1 belongs to Lp [0, ∞) if 0 |f (t)|p dt < ∞ and then, with the norm *∞ f p = ( 0 |f (t)|p dt)1/p , the space Lp [0, ∞) is a Banach space (linear normed space in which the Cauchy sequences are convergent). If p = ∞, then with the norm f ∞ = ess sup |f (t)| < ∞ the space L∞ [0, ∞) is also a Banach space. The definition can easily be extended to vector valued functions  f = (f1 , . . . , fn )T , fi ∈ Lp [0, ∞), i = 1, . . . , n where, with the norm f p := ( ni=1 fi 2p )1/2 , the space Lnp [0, ∞) is also a Banach space. *∞ For p = 2, a scalar product f, g = 0 f (t), g(t) dt can be explained and the norm can be defined by f  = ( f, f )1/2 . Then Ln2 [0, ∞) is a Hilbert space. Further extensions are possible for functions which are not necessarily bounded but satisfy the above conditions for every finite T > 0 if the function is substituted by zero for ∀t > T . Such extended spaces are denoted by Lnpe . For simplicity, we shall consider only Lnp [0, ∞) signal spaces. 2.3.4.1 Lp -Stable Nonlinear Systems n Let G : Lnp → Lm p a nonlinear system mapping u ∈ Lp [0, ∞) input signals into m y ∈ Lp [0, ∞) output signals. Then the nonlinear system is (i) Lp -stable with finite gain (wfg) if there exist γp , βp ≥ 0 constants such that x ∈ Lnp ⇒ yp ≤ γp xp + βp . (ii) Lp -stable with finite gain without bias (wb) if there exist γp ≥ 0 constant such that x ∈ Lnp ⇒ yp ≤ γp xp . The gain of the nonlinear system is defined by γp (G) = inf{γp } satisfying (i) for wfg and (ii) for wb, respectively. Let us consider now interconnected stable nonlinear systems according to Fig. 2.8 and ask for the stability of the closed loop. The proofs of the following two theorems are simple, see also [152]. 50 2 Basic Nonlinear Control Methods Fig. 2.8 Scheme of nonlinear control system Theorem 2.10 Small gain theorem. Consider the closed loop system according to Fig. 2.8 and let p ∈ [1, ∞]. Assume that G1 , G2 are causal and Lp -stable wb (or wfg) nonlinear systems, γ1p := γ1p (G1 ) and γ2p := γ2p (G2 ). Then the closed loop system is Lp -stable, if γ1p · γ2p < 1. The inputs of nonlinear electrical circuits are voltage and current generators, the outputs are voltages and currents and the instantaneous power is y(t), u(t). The circuit is passive if the input energy (integral of the power) in the system does not exceed the change of the stored * ∞ energy. If the inductive and capacitive elements are initially chargeless, then 0 y(t), u(t) dt ≥ 0 for passive electrical circuits. Passive systems generalize this principle. Notice, scalar product needs that the dimension of input and output signals is the same. 2.3.4.2 Passive Systems Let G : Ln2 → Ln2 a nonlinear system. We say that *∞ (i) G is passive, if ∃β ∈ R 1 such that u, Gu = 0 u(t), (Gu)(t) dt ≥ β, ∀u ∈ L2 . (ii) G is strictly input passive, if ∃α > 0, β ∈ R 1 such that u, Gu ≥ αu22 + β, ∀u ∈ L2 . (iii) G is strictly output passive, if ∃α > 0, β ∈ R 1 such that u, Gu ≥ α(Gu)22 + β, ∀u ∈ L2 . Notice that since u ≡ 0 ∈ L2 thus 0 ≥ β. It can easily be shown that strictly output passive system is L2 -stable wfg. Theorem 2.11 Passivity theorem in Lp space. Consider the closed loop system in Fig. 2.8 with G1 , G2 ∈ Ln2 . (i) If both G1 and G2 are strictly output passive and both u1 and u2 are active then, for the output yy12 , the closed loop is strictly output passive and L2 -stable wfg. (ii) If G1 is strictly output passive, G2 is passive and u2 = 0 (i.e., only u1 is active), then for the output y1 , the closed loop is strictly output passive and L2 -stable wfg. (iii) If G1 is passive, G2 is strictly input passive and u2 = 0 (i.e., only u1 is active), then for the output y1 , the closed loop is strictly output passive and L2 -stable wfg. 2.3 Stability of Nonlinear Systems 51 2.3.4.3 Dissipative Systems Dissipativity can be seen as passivity in state space. The storage function serves as Lyapunov function. The number of inputs and outputs is the same, the input and output signals are from L2 . Consider the nonlinear system Σ : ẋ = f (x, u), y = h(x, u), x ∈ R n , u, y ∈ R m . The real valued function s(u, y) will be called supply rate. The system is called dissipative if there exist S : R n → R + so called storage function such that for every x(t0 ) ∈ R n and for every t1 ≥ t0 the following dissipative inequality yields:     S x(t1 ) ≤ S x(t0 ) +  t1 t0   s u(t), y(t) dt. (2.87) Specifically, if equality holds then the system is called lossless. As an example, consider the supply rate s(u, y) := u, y. Assuming S ≥ 0, then  0 T( )       u(t), y(t) dt ≥ S x(T ) − S x(0) ≥ −S x(0) (2.88) for every x(0), for every T > 0 and for every input u(·) which means that the input– output mapping Gx(0) is passive with the choice β = −S(x(0)). Here, S can be considered to be the stored energy motivating the following definitions. Passive dissipative systems Assume the nonlinear system has the form Σ : ẋ = f (x, u), y = h(x, u), x ∈ R n , u, y ∈ R m , then the system (i) Is passive, if it is dissipative with respect to the supply rate s(u, y) = u, y. (ii) Is strictly input passive if it is dissipative with respect to the supply rate s(u, y) = u, y − δu2 where δ > 0. (iii) Is strictly output passive if it is dissipative with respect to the supply rate s(u, y) = u, y − εy2 where ε > 0. (iv) Is conservative if it is lossless and s(u, y) = u, y. (v) Has L2 -gain ≤ γ if it is dissipative with respect to the supply rate s(u, y) = 1 2 1 2 2 2 γ u − 2 y where γ > 0. For example, in the last case 1 2  0 T +2 + +2  +       γ 2 +u(t)+ − +y(t)+ dt ≥ S x(T ) − S x(0) ≥ −S x(0) , from which after regrouping follows  0 T+ + +y(t)+2 dt ≤ γ 2  0 T+ +   +u(t)+2 dt + 2S x(0) , which is equivalent to Lp -stability wfg ≤ γ . 52 2 Basic Nonlinear Control Methods The dissipativity condition (2.87) can be substituted by its differentiated form which is sometimes easier to check. By dividing the inequality by t1 − t0 and letting t1 → t0 , we obtain the equivalent differential dissipation inequality   Sx′ (x)f (x, u) ≤ s u, h(x, u) , ∀x, u,   (2.89) ∂S ∂S ′ Sx (x) = (x) · · · (x) . ∂x1 ∂xn Specifically, in case of input affine nonlinear system, that is, Σa : ẋ = f (x) + g(x)u, y = h(x), and assuming the system is passive with respect to the supply rate s(u, y) = u, y = u, h(x), then the differential form for Σa amounts to Sx′ (x)[f (x) + g(x)u] ≤ uT h(x), ∀x, u, or equivalently (choosing first u = 0) to the Hill–Moylan conditions Sx′ (x)f (x) ≤ 0, (2.90) Sx′ (x)g(x) = hT (x). Special cases: (i) For LTI system, the conditions are equivalent to the Kalman–Yacubovitch– Popov conditions if the storage function is S(x) = 21 x T P x, where P = P T ≥ 0, namely AT P + P A ≤ 0 and B T P = C. (ii) If the input affine system Σa is strictly output passive with respect to the supply rate s(u, y) = u, y − εy2 where ε > 0, then (choosing first u = 0) yields Sx′ (x)f (x) ≤ −εhT (x)h(x), (2.91) Sx′ (x)g(x) = hT (x). (iii) If the input affine system Σa has finite L2 -gain with respect to the supply rate s(u, y) = 21 γ u2 − 12 y2 , then   1 Sx′ (x) f (x) + g(x)u − γ 2 u2 + 2 + 1+ +h(x)+2 ≤ 0, 2 ∀x, u. (2.92) (iv) The maximum of the left side in the inequality (2.92) is reached at u∗ = ′ 1 T g (x)SxT (x) and substituting it back we obtain the Hamilton–Jacobi– γ2 Bellman inequality Sx′ (x)f (x) + 1 1 ′ ′ Sx (x)g(x)g T (x)SxT (x) + hT (x)h(x) ≤ 0, 2 2 2γ ∀x. (2.93) Hence, Σa has L2 -gain ≤ γ if and only if (2.93) has continuously differentiable solution S ≥ 0. We remark that the equality form is in strong connection with the Hamilton–Jacobi–Bellman equality of dynamic programming. The following theorems were proven in [127]. 2.3 Stability of Nonlinear Systems 53 Theorem 2.12 Extremum property of Sa amongst storage functions. Consider the nonlinear system Σ . The system is dissipative with respect to the supply rate s(u, y) ⇔ for every initial condition x(0) = x yields  T   s u(t), y(t) dt < ∞. (2.94) Sa (x) = sup − u(·),T ≥0 0 Furthermore, every other storage function satisfies Sa (x) ≤ S(x), ∀x. Thus, Sa (x) can be considered the maximal energy which can be extracted from the system Σ if the initial condition is x. If each state x of Σ can be reached from x∗ within finite time, then Σ is dissipative ⇔ Sa (x∗ ) < ∞. Moreover, if Gx∗ is passive wb (β = 0), or has L2 -gain ≤ γ wb, then Sa (x∗ ) = 0. Zero-state detectability. The nonlinear system Σa is called zero-state detectable if u(t) ≡ 0, y(t) ≡ 0, ∀t ≥ 0 ⇒ limt→∞ x(t) = 0. Sa (x∗ ) = 0 is equivalent to the zero-state detectability condition. Theorem 2.13 Connection between passivity and asymptotic stability. Consider the zero-state detectable nonlinear system Σa . Let S ≥ 0 with S(0) = 0 be solution to (2.91) in case of strict output passivity, or (2.92) in case of finite L2 -gain. Then ξ ≡ 0 is asymptotically stable equilibrium point. Moreover if S(x) is proper, that is, {x : S(x) < c} is compact, ∀c > 0, then the equilibrium point is globally asymptotically stable. Consider the closed loop nonlinear system according to Fig. 2.8 but at the place of Gi stands the system Σi : ẋi = fi (xi , ei ), yi = hi (xi , ei ), i = 1, 2. Assume that the strict dissipative conditions  t1 +2  +      ei (t), yi (t) − εi +yi (t)+ dt, i = 1, 2 (2.95) Si xi (t1 ) ≤ Si xi (t0 ) + t0 are satisfied. Theorem 2.14 Stability of interconnected dissipative systems. (i) Suppose both Σ1 and Σ2 are passive (strictly output passive) and (u1 , u2 ) external inputs are active, then the closed loop system with outputs (y1 , y2 ) is output passive (strictly output passive). (ii) Suppose Si satisfies (2.95), is continuously differentiable and has (isolated) local minimum at xi∗ , i = 1, 2, furthermore no external inputs are active, that is, u1 = u2 = 0. Then (x1∗ , x2∗ ) is stable equilibrium point of the closed loop system. (iii) Suppose Σi is strictly output passive and zero-state detectable, Si satisfying (2.95), is continuously differentiable and has (isolated) local minimum at 54 2 Basic Nonlinear Control Methods xi∗ = 0, i = 1, 2, furthermore no external inputs are active, that is, u1 = u2 = 0. Then (0, 0) is asymptotically stable equilibrium point of the closed loop. If additionally Si has global minimum at xi∗ = 0, i = 1, 2, then (0, 0) is globally asymptotically stable equilibrium point of the closed loop system. Theorem 2.15 Small gain theorem for interconnected dissipative systems. (i) Suppose Σi is zero-state detectable, has finite gain γi , i = 1, 2, with γ1 · γ2 < 1. Suppose Si satisfies (2.95), is continuously differentiable and has (isolated) local minimum at xi∗ = 0, i = 1, 2, furthermore no external inputs are active, that is, u1 = u2 = 0. Then (0, 0) is asymptotically stable equilibrium point of the closed loop. If additionally Si has global minimum at xi∗ = 0 and is proper, i = 1, 2, then (0, 0) is globally asymptotically stable equilibrium point of the closed loop system. (ii) Suppose Σai is input affine system, has finite gain γi , i = 1, 2, with γ1 · γ2 < 1. Suppose Si satisfies (2.95), is continuously differentiable and Si (0) = 0 (therefore Si (x) is positive semidefinite), and let xi∗ = 0 be asymptotically stable equilibrium point of the unexcited system ẋi = fi (xi ), i = 1, 2, furthermore no external inputs are active, that is, u1 = u2 = 0. Then (0, 0) is asymptotically stable equilibrium point of the closed loop system. 2.4 Input–Output Linearization It is a meaningful idea to control the nonlinear system in such a way that first the nonlinear system is linearized then the linearized system is compensated by using attractive methods of linear control theory. For linearizing the nonlinear system, state feedback and output feedback can be suggested, called input/state linearization and input/output linearization, respectively. Their theories are well developed, see for example [56, 57, 152], however their attractiveness is limited for low dimensional systems because of computational problems, mainly because of the need for solving partial differential equations, the so called Frobenius equations. The basic results were also surveyed in Appendix B and we assume the reader became acquainted with the main techniques from there. Comparing the two approaches, input/state linearization has the advantage that the number of inputs and outputs may be different, but has the drawback that a large number of sets of Frobenius equations has to be solved. On the other hand, input/output linearization has the disadvantage that the number of inputs and outputs must be equal, but the solution of Frobenius equation can be evaded, which is a benefit, except for SISO system if the input should be eliminated from the zero dynamics. Notice that eliminating the input from the zero dynamics for MIMO nonlinear systems is usually not possible because the arising distribution spanned by the fix part of the nonlinear coordinate transformation needed for linearization is usually not involutive and the Frobenius equation cannot be applied. 2.4 Input–Output Linearization 55 In both cases, the remaining part of the computations can be performed automatically by using symbolic computation tools, like MATLAB Extended Symbolic Toolbox, Mathematica, etc. Now we summarize here the basic results regarding input/output linearization discussed in Appendix B. The system is assumed to be input affine, that is, m ẋ = f (x) + ui gi (x), i=1 yi (x) = hi (x), i = 1, . . . , m, h = (h1 , . . . , hm )T , f, gi ∈ V (X), hi ∈ S(X), i = 1, . . . , m, where f, gi are vector fields and hi are real valued functions, all appropriately smooth. As can be seen, the input does not appear in the output mappings hence the idea is to differentiate the outputs so many times that at the end the input appears in the results. This will define the components of the possible vector degree for the different outputs. For example, differentiating yi by the time yields , m m ẏi = ∇hi ẋ = dhi , f + i=1 gi ui = Lf hi + uj Lgj hi . j =1 If each Lgj hi is zero, then the differentiation is continued and we get , m m ÿi = ∇(Lf hi )ẋ = d(Lf hi ), f + i=1 gi ui = L2f hi + uj Lgj (Lf hi ). j =1 Let ri be the smallest integer for which m (ri ) yi = Lrfi hi + j =1   Lgj Lrfi −1 hi uj . The remaining problem is how to express the inputs from the formulas, which is related to the linear independence of the functions weighting the inputs after differentiation. The MIMO input affine nonlinear system is said to have the relative degree vector r = (r1 , . . . , rm )T at the point x0 if ∃U (x0 ) open set such that: (i) Lgj Lkf hi (x) ≡ 0, ∀x ∈ U, j = 1, . . . , m, k = 0, . . . , ri − 2. r −1 (ii) The matrix S(x) = [sij (x)]m×m is nonsingular at x0 where sij = Lgj Lfi i.e. ⎤ ⎡ Lg1 Lrf1 −1 h1 · · · Lgm Lrf1 −1 h1 ⎥ ⎢ .. .. .. ⎥. S(x) = ⎢ . . . ⎦ ⎣ rm −1 rm −1 Lg1 Lf hm · · · Lgm Lf hm hi , 56 2 Basic Nonlinear Control Methods Fig. 2.9 Block scheme of the MIMO input/output linearized system Since S(x0 ) is nonsingular hence S(x) is also nonsingular in a sufficiently small neighborhood of x0 . The components qi (x) = Lrfi hi (x), i = 1, . . . , m and the column vector q(x) will also be needed to internal compensation, see Fig. 2.9. Theorem 2.16 The MIMO nonlinear system can be input/output linearized in the neighborhood U (x0 ) if it has relative degree vector and S(x0 ) is invertible. There is a smooth nonlinear transformation z = T (x), where z = (z0T , zuT )T , zo is locally observable while zu is locally not observable. By using nonlinear static feedback u(x) = S −1 (v − q(x)), the observable part consists of decoupled integrators in Brunowsky canonical form. The not observable part remains nonlinear and has the form żu = bu (zo , zu ) + Au (zo , zu )u. The observable part can be further compensated by external linear feedback. The full system stability is influenced by the stability of the unobservable part that is called the zero dynamics of the nonlinear system. In Appendix B, formulas can be found for computing the normal form of the input affine system after coordinate transformation T (x). The system components of the observable and not observable parts are given in (B.51a) and (B.51b), respectively. The resulting observable system is given in (B.52a). The resulting zero dynamics is given in (B.52b) containing u and in (B.54) containing v, respectively. The fixed parts of the components of the coordinate transformation T (x) are given in (B.50), these have to be extended with further functions to make the transformation n-dimensional, that is, to a diffeomorphism. This step does not need the solution of any Frobenius equation, see Lemma B.1 in Appendix B, but the price is that the zero dynamics depends also on the input u. An exception may be the SISO case where the distribution  spanned by the single vector g is evidently involutive. Hence, by using Frobenius theorem there are η1 , . . . , ηn−1 ∈ S(X) functions whose differentials are independent at x0 and satisfy dηi , g(x) = 0, ∀x ∈ U . Taking n − r pieces from these functions, for example the first n − r, the coordinate transformation may be T (x) = T T T [dhT , . . . , (dLr−1 f h) , (dη1 ) , . . . , (dηn−r ) ]. By using this coordinate transformation, the input u will not appear in the zero dynamics, i.e. żu = fu (zo , zu ). Further details can also be found in Appendix B regarding the stability test of the zero dynamics and the extension of the minimum phase property to nonlinear systems based on the stability of the zero dynamics. Another possibility is the use of precompensator based on the dynamic extension algorithm if the vector relative 2.5 Flatness Control 57 degree does originally not exist. In this case, the compensator is a dynamic system, see (B.56) and (B.57). 2.5 Flatness Control We illustrate here how can the flatness property be exploited in the path design of underactuated systems, for which not every path is realizable. Assuming underactuated system, it means that for a separately designed path x(t) it may be possible that there is no appropriate control u(t) which can realize the path. On the other hand, for differentially flat systems we know that both the state and the control can be parametrized by the output variable and its derivatives. Hence, choosing a fine sequence of “corner points” for the path in the output variable and interpolating amongst them analytically, that is, determining also the derivatives in the necessary order, we can simultaneously find both the state trajectory and the corresponding nominal control. Flatness properties of cable-driven crane. We shall illustrate the flatness concept and its use in path design for a cable-driven crane. For simplicity, only the two dimensional (2D) planar version is considered. The scheme of the 2D crane is shown in Fig. 2.10. (Since the pulley’s mass in point B has been neglected, hence Ls and its actuator motor is omitted in the discussion.) The planar crane is a 3-DOF system which is underactuated because the only inputs are the two motor torques moving the ropes L1 and L2 + L3 . The model of the planar crane is differentially flat [69]. In simplified formulation, it means that if ẋ = f (x, u) is the dynamic model of the system then there exist a new variable y (the flat output) and finite integers q and p such that x = φ1 (y, ẏ, ÿ, . . . , y (q) ), u = φ2 (y, ẏ, ÿ, . . . , y (q+1) ), y = ψ1 (x, u, u̇, . . . , u(p) ). For the crane, y := (x, z)T where x, z are the coordinates of the load (point C) in the plane of the crane. The coordinates of point B will be denoted by xB , zB . The flatness relations for the 2D crane are summarized as follows:   θ = arctan −ẍ/(z̈ + g) , ẍ(z − (k + l)Cα ) − (z̈ + g)(x − (k + l)Sα ) , (z̈ + g)Sa − ẍCa   β = arcsin (1 − d/ l)Sα+θ ,   γ = π + β − (α + θ ) /2, d= L1 = lSβ /Sγ , L2 = lSγ −β /Sγ , xB = (k + l)Sα − L2 Sα−β , zB = (k + l)Cα − L2 Cα−β , (2.96) 58 2 Basic Nonlinear Control Methods Fig. 2.10 Block scheme of the cable-driven 2D crane . (x − xB )2 + (z − zB )2 ,   T3 = mload (z̈ + g)Cθ − ẍSθ , L3 = J1 L̈1 , ρ1  J2  u2 = T3 ρ2 − L̈2 + L̈3 . ρ2 u1 = T1 ρ1 − If y, ẏ, . . . , y (4) are known, then all system variables can be computed from them except singularities. For prescribed initial y0 and final y1 and transient time TF , five times differentiable y(t) polynomials can be found connecting the points with straight line in the space from which all the system states x(t) and inputs u(t) can be reconstructed. It follows from the equilibrium of the internal forces that the line AB in the direction of the horizontal rope L1 is the bysectrix of the angle 2γ between PB and 2.5 Flatness Control 59 BC. Hence, the magnitudes of the internal forces T3 for L2 and T1 for L1 are related by T1 = 2Cγ T3 . By using the flat output all the variables describing the 2D crane can be determined by (2.96). In order to avoid the first and second order symbolic differentiation of the nonlinear expressions for L1 , L2 , L3 , numerical differentiation can be applied (see spline, unmkpp, mkpp in MATLAB). Finally, the flatness based input torques u1 , u2 for the motors can be determined, where Ji and ρi denote the inertia and the radius of the winch, respectively. Assuming zero derivatives until order five in the initial and final positions, the 11th order polynomial connecting (x0 , z0 ) and (x1 , z1 ) has the form   x(τ ) = x0 + (x1 − x0 )τ 6 × a6 + a7 τ + a8 τ 2 + a9 τ 3 + a10 τ 4 + a11 τ 5 , (2.97) where τ = t/TF is the normalized time and the coefficients for both x(τ ) and z(τ ) are a6 = 0.462, a7 = −1.980, a9 = −3.080, a10 = 1.386, a8 = 3.465, a11 = −0.252. (2.98) Neglecting the feedforward compensated friction and the centripetal and Coriolis effects (because of the slow angular velocities of typical cranes) the basic equations of the underactuated system are J1 L̈1 = T1 ρ1 − u1 , ρ1  J2  L̈2 + L̈3 = T3 ρ2 − u2 . ρ2 (2.99) The system is submitted to constraints with respect to L1 and L2 hence Lagrange multipliers could have been used. Fortunately, instead of this relatively complicated method, the equilibrium between internal forces can be exploited, see Kiss [69], who developed a form of the dynamic model for the choice x = (γ , α − β, L3 , γ̇ , (α̇ − β̇), L̇3 )T where α = const. However, the measured outputs are L1 , L2 + L3 , L̇1 , L̇2 + L̇3 which depend on these state variables in a nonlinear way and would make necessary the use of nonlinear state estimation. Hence, a new model has been developed using x = (L1 , L2 , L3 , L̇1 , L̇2 , L̇3 )T , see [74]. The idea is the use of the cosine theorem by which Cγ = l 2 − L21 − L22 , 2L1 L2 Cβ = l 2 + L22 − L21 . 2lL2 (2.100) For L1 and L2 , the expressions in (2.96) can be used. By elementary geometry and using     ϕ = (π/2) − γ + (α − β) , θ = π − 2γ + (α − β) , it yields x = kSα + L1 Cϕ + L3 Sθ = kSα + L1 Sγ +(α−β) + L3 S2γ +(α−β) , z = kCα + L1 Sϕ − L3 Cθ = kCα + L1 Cγ +(α−β) + L3 C2γ +(α−β) . (2.101) 60 2 Basic Nonlinear Control Methods The main problem is the computation of T3 from the state variables which is needed in (2.99). This can be performed using mẍ = −T3 Sθ , m(z̈ + g) = T3 Cθ for the load, from which it follows   Cθ ẍ + Sθ (z̈ + g) = 0, m −Sθ ẍ + Cθ (z̈ + g) = T3 . (2.102) Notice that the presented flatness based method computes only the nominal path and the nominal (open loop) control but does not solve the feedback control problem. However, the results can be used for feedforward compensation or linearization along the path. 2.6 Backstepping Backstepping control is a popular method for controlling nonlinear systems given in strictly feedback form ẋ1 = G1 (x̄1 )x2 + f1 (x̄1 ), ẋ2 = G2 (x̄2 )x3 + f2 (x̄2 ), .. . (2.103) ẋn−1 = Gn−1 (x̄n−1 )xn + fn−1 (x̄n−1 ), ẋn = Gn (x̄n )u + fn (x̄n ), y = h(x1 ), where xi ∈ R m , i = 1, . . . , n are vector parts of the state vector, y ∈ R m is the output, x̄i = (x1T , . . . , xiT )T is a short notation, h is an invertible function (i.e. its Jacobian hx1 is an invertible matrix), Gi are invertible matrices, i = 1, . . . , n and all functions are appropriately smooth. We assume the task is to realize a prescribed path yd (θ ) ∈ R m for the output variable y. The path is parametrized in θ ∈ R 1 . Let vs (θ, t) ∈ R 1 be the prescribed speed of the path parameter to which a parameter speed error ωs (θ̇ , θ, t) = vs (θ, t) − θ̇ may exist. The result of the path design subtask is the path parameter θ , the path yd (θ ) and the desired path parameter speed vs (θ, t). An extra problem is to make the parameter speed error ωs going to zero using high level control. A special case is θ = t, vs = 1 in which case ωs ≡ 0 is already satisfied and the desired path is yd (t). We shall often need the total derivative by the time of parameters of the form αi−1 (x̄i−1 , θ, t), that is, ∂αi−1 dαi−1 ∂αi−1 ∂αi−1 = θ̇ ẋ1 + · · · + ẋi−1 + dt ∂x1 ∂xi−1 ∂θ θ θ =: σi−1 + αi−1 θ̇ = σi−1 + αi−1 (vs − ωs ). (2.104) 2.6 Backstepping 61 Observe that the derivatives ẋi can be substituted by the right side of the state equations. The first term σi−1 in the total time derivative α̇i−1 is that part which does not contain θ̇ . Notice that any variable in the upper index denote the partial derivative by the variable, see also hx1 . Now we present the backstepping compensation method similarly to Skjetne, Fossen and Kokotovich [130], but we neglect the disturbance effects for simplicity. At every level, we design the virtual control αi−1 , for stabilization purposes we choose Hurwitz matrix Ai , that is, matrix whose eigenvalues are on the open left complex plane, choose positive definite matrix Qi > 0 and solve the Lyapunov equation Pi Ai + ATi Pi = −Qi , Qi > 0, Pi > 0, Ai is Hurwitz, (2.105) to find the symmetric and positive definite matrix Pi > 0 which will be used in the construction of the Lyapunov function Vi . The Lyapunov function for the composite system will be V = Vn . During the process, the original variables xi will be changed to zi and the equilibrium point will be z = (z1T , . . . , znT )T = 0 whose stability has to be guaranteed. The backstepping algorithm consists of the following steps. Step i = 1 In the first step, let z1 = y − yd and z2 = x2 − α1 . Then ż1 = hx1 ẋ1 − ẏd = hx1 (G1 x2 + f1 ) − ẏd = hx1 G1 (z2 + α1 ) + hx1 f1 − ydθ (vs − ωs ). Define the first virtual control by   x −1  1 A1 z1 − hx1 f1 + ydθ vs . α1 := G−1 1 h (2.106) ż1 = hx1 G1 z2 + A1 z1 + ydθ ωs . (2.107) Substituting α1 into ż1 then ż1 = hx1 G1 z2 + hx1 G1 α1 − ydθ (vs − ωs ) and canceling the appropriate terms with positive and negative sign yields Define the first partial Lyapunov function by V1 = z1T P1 z1 and take into consideration that V̇1 = 2z1T P1 ż1 and then  ) (  2z1T P1 A1 z1 = 2 z1 , P1 A1 z1  = z1 , P1 A1 + AT1 P1 z1 = −z1T Q1 z1   V̇1 = 2z1T P1 hx1 G1 z2 + A1 z1 + ydθ ωs = −z1T Q1 z1 + 2z1T P1 hx1 G1 z2 + 2z1T P1 ydθ ωs , $ %& ' τ1 (2.108) 62 2 Basic Nonlinear Control Methods τ1 = 2z1T P1 ydθ . Step i = 2 Since the desired path appears only in the definition of z1 , therefore some differences will be in ż2 . In the second step, let z3 = x3 − α2 . Then ż2 = ẋ2 − α̇1 and thus ż2 = G2 (z3 + α2 ) + f2 − σ1 − α1θ (vs − ωs ) = G2 z3 + G2 α2 + f2 − σ1 − α1θ (vs − ωs ). Define the second virtual control by   −1 T  x1 T α2 := G−1 P1 z1 − f2 + σ1 + α1θ vs . 2 A2 z2 − P2 G1 h (2.109) Substituting α2 into ż2 and canceling the appropriate terms with positive and negative sign yields  T ż2 = G2 z3 + A2 z2 − P2−1 GT1 hx1 P1 z1 + α1θ ωs . (2.110) Define the second partial Lyapunov function by V2 = V1 + z2T P2 z2 and take into consideration that V̇2 = V̇1 + 2z2T P2 ż2 and 2z2T P2 A2 z2 = −z2T Q2 z2 then V̇2 = −z1T Q1 z1 + 2z1T P1 hx1 G1 z2 + τ1   T  + 2z2T P2 G2 z3 + A2 z2 − P2−1 GT1 hx1 P1 z1 + α1θ ωs   = −z1T Q1 z1 − z2T Q2 z2 + 2z2T P2 G2 z3 + τ1 + 2z2T P2 α1θ ωs , (2.111) %& ' $ τ2 τ2 = τ1 + 2z2T P2 α1θ . Steps i = 3, . . . , n − 1 Notice that α1 , ż1 , V̇1 , α2 , ż2 contain hx1 but V̇2 does not. Hence, in the next steps some modifications are necessary: zi = xi − αi−1 ,   θ αi = Gi−1 Ai zi − Pi−1 GTi−1 Pi−1 zi−1 − fi + σi−1 + αi−1 vs , θ żi = Gi zi+1 + Ai zi − Pi−1 GTi−1 Pi−1 zi−1 + αi−1 ωs , i Vi = zjT Pj zj , j =1 i V̇i = − j =1 zjT Qj zj + 2ziT Pi Gi zi+1 + τi ωs , θ τi = τi−1 + 2ziT Pi αi−1 . (2.112) 2.6 Backstepping 63 Step i = n θ , from which the In the last step, let zn = xn − αn−1 , żn = Gn u + fn − σn−1 − αn−1 control law can be determined:   θ −1 T u = G−1 n An zn − Pn Gn−1 Pn−1 zn−1 − fn + σn−1 + αn−1 vs , θ ωs , żn = An zn − Pn−1 GTn−1 Pn−1 zn−1 + αn−1 n Vn = zjT Pj zj , (2.113) j =1 n V̇n = − j =1 zjT Qj zj + τn ωs , θ τn = τn−1 + 2znT Pn αn−1 . The closed loop system can be brought on the form ż = A(x̄n , θ, t)z + b(x̄n−1 , θ, t)ωs , (2.114) where A(x̄n , θ, t) ⎡ A1 ⎢−P −1 GT (hx1 )T P1 1 ⎢ 2 ⎢ 0 =⎢ ⎢ .. ⎣ . 0 ⎛ ⎞ ydθ ⎜ θ ⎟ ⎜ α ⎟ ⎜ 1 ⎟ ⎜ ⎟ b(x̄n−1 , θ, t) = ⎜ α θ ⎟ . ⎜ 2 ⎟ ⎜ .. ⎟ ⎝ . ⎠ θ αn−1 hx1 G1 A2 −P3−1 GT2 P2 .. . 0 G2 A3 .. . 0 0 G3 .. . ··· ··· ··· .. . 0 0 0 .. . 0 ··· 0 −Pn−1 GTn−1 Pn−1 An ⎤ ⎥ ⎥ ⎥ ⎥, ⎥ ⎦ (2.115) Remark • The system stability depends on τn ωs where ωs is the speed error depending on the path parameter. However, if a high level controller makes ωs convergent to zero then the derivative of the Lyapunov function V̇ = V̇n becomes negative definite. • If the path is designed as function of the time, that is, θ = t and the path is yd (t), then ωs ≡ 0, and the closed loop is globally asymptotically stable. 64 2 Basic Nonlinear Control Methods θ θ̇ hides the complexity of the computation. In • The notation α̇i−1 = σi−1 + αi−1 reality, the formulas are early exploding with increasing n. The use of symbolic computation tools is suggested to find the derivatives for n > 2. 2.7 Sliding Control Sliding control is an important method in the class of variable structure control approaches [8, 150]. Its importance is justified by the fact that sliding controller does not need the precise model of the process to be controlled, it is enough to know a simplified model and some upper bounds of the modeling errors. On the other hand, important other methods exist which ever are not sliding controller but integrate in their algorithms the sliding control principle too. Such a method is for example self-tuning adaptive control. We try to introduce the sliding control concept as simple as possible. Assume the state variable is x = (y, y ′ , . . . , y (n−1) )T and the SISO system has input affine model and error signal y (n) = f (x) + g(x)u, (2.116) ỹ(t) = yd (t) − y(t), (2.117) respectively. Suppose xd (0) = x(0) otherwise the error will be composed by integration. The switching surface s(x, t) = 0, called sliding surface, is defined by the rule  n−1 d s(x, t) := ỹ(t), (2.118) +λ dt where λ > 0. It is evident that if the trajectory remains on the sliding surface (sliding mode) then the error goes exponentially to zero. Notice that 21 |s|2 can be considered as Lyapunov function hence it is useful to require 1 d 2 |s| = ss ′ ≤ −η|s|, 2 dt (2.119) where η > 0. If the above conditions are satisfied but in the beginning the trajectory is not on the sliding surface, that is, s(0) = 0, then (since, for example, in case of s(0) > 0 yields s ′ ≤ η) the trajectory surely reaches the sliding surface s = 0 within the time ts ≤ |s(0| . η (2.120) The time for approaching the sliding surface can be influenced by the parameter η while the time of decaying by the parameter λ, which simplify the design. In consequence, the main task of every control system, that is, eliminating the error, is conceptually fulfilled. 2.7 Sliding Control 65 2.7.1 Sliding Control of Second Order Systems The remaining question is how to choose the control law in order to satisfy the above conditions. We illustrate it through examples for second order systems of increasing complexity. 2.7.1.1 Proportional Control Let the system to be controlled ÿ = f (x) + u, the error signal is ỹ(t), n = 2, and by definition s = ỹ˙ + λỹ, (2.121) ˙ s ′ = ÿd − ÿ + λỹ˙ = ÿd − f − u + λỹ. (2.122) Assume that instead of f (x) only its estimation fˆ(x) is available. Let the form of the control u = û + k(x) sign(s). (2.123) It is clear that on the sliding surface s = 0, sign(0) = 0 by definition, hence u = û. Reaching the sliding surface the best concept is to remain on it since then the error goes exponentially to zero. But if s = 0 (constant), then s ′ = 0, too. The only information is fˆ(x) therefore the strategy may be s = 0, fˆ(x) ⇒ ˙ û = ÿd − fˆ(x) + λỹ. s ′ = 0, u = û, (2.124) (2.125) However, if the sliding surface has not been reached yet, that is, s = 0, then in order to reach the sliding surface within finite time, ss ′ ≤ −η|s| has to be satisfied, hence s ′ = fˆ − f − k(x) sign(s), ! !   ss ′ = f − fˆ s − k(x)|s| ≤ !f − fˆ!|s| − k(x)|s| ≤ −η|s|, ! ! k(x) ≥ !f (x) − fˆ(x)! + η. (2.126) (2.127) (2.128) Let the upper bound of the deviation between model and estimation be given by the known function F (x), ! ! !f (x) − fˆ(x)! ≤ F (x) (2.129) then k(x) can be chosen k(x) = F (x) + η. (2.130) 66 2 Basic Nonlinear Control Methods 2.7.1.2 Integral Control the system be again ÿ = f (x) + u, but now the error is given by the integral *Let t ỹ(τ ) dτ , for which formally n = 3 can be assumed. Then by definition 0 s=  d +λ dt 2  0 t ỹ(τ ) dτ = ỹ˙ + 2λỹ + λ2  t ỹ(τ ) dτ, (2.131) 0 s ′ = ÿd − ÿ + 2λỹ˙ + λ2 ỹ = ÿd − f − u + 2λỹ˙ + λ2 ỹ. (2.132) Hence, similarly to the previous case, the strategy may be s = 0, fˆ(x) ⇒ s ′ = 0, u = û yielding u = û + k(x) sign(s), û = ÿd − fˆ(x) + 2λỹ˙ + λ2 ỹ, k(x) = F (x) + η. (2.133) (2.134) (2.135) It is an advantage of integral control that s(0) = 0 can also be assured, that is, the control can already start on the sliding surface eliminating ts , if we choose  t 2 ˙ ˙ ỹ(τ ) dτ − ỹ(0) − 2λỹ(0). (2.136) s(x, t) := ỹ + 2λỹ + λ 0 Remark The difference between proportional and integral control is in the computation of s and û, however k(x) is the same. Henceforth, we shall deal only with proportional control, and in case of integral control we modify the final s and û according to the above results. 2.7.1.3 Proportional Control in Case of Variable Gain In the above examples, we had unit gain. Now, we consider again the second order system ÿ = f (x) + g(x)u but with state dependent gain having strictly positive sign. Suppose that the estimation ĝ(x) of the state dependent gain and a parameter β > 0 (may also be state dependent) is known, furthermore the relation between real and estimated gain satisfies the special condition β −1 ≤ ĝ ≤ β. g (2.137) By definition s=   d + λ ỹ = ỹ˙ + λỹ, dt ˙ s ′ = ÿd − ÿ + λỹ˙ = ÿd − f − gu + λỹ. (2.138) (2.139) 2.7 Sliding Control 67 Let the form of the controller   u = b̂−1 û + k(x) sign(s) , (2.140) û = ÿ − fˆ(x) + λỹ. (2.141) and derive û from the strategy s = 0, fˆ(x), b̂(x) ⇒ s ′ = 0, u = û, then Considering the general case s = 0 and supposing the sliding surface has not been reached yet, then   s ′ = fˆ − f + 1 − g ĝ −1 û − g ĝ −1 k(x) sign(s), (2.142) ! ! ! ! ss ′ ≤ !f − fˆ!|s| + !1 − g ĝ −1 !|û||s| − g ĝ −1 k(x)|s| ≤ −η|s|, (2.143) ! ! ! ! (2.144) k(x) ≥ ĝg −1 !f − fˆ! + !ĝg −1 − 1!|û| + ĝg −1 η. According to the bounds of ĝg −1 the nonlinear gain k(x) can be chosen k(x) = β(F + η) + (β − 1)|û|. (2.145) Now let us consider the more practical case if minimal and maximal bounds of the gain are known: gmin ≤ g ≤ gmax . (2.146) This case can be transformed to the previous one by using the geometric average of the minimal and maximal gains: √ ĝ := gmin · gmax , (2.147) / gmax β := , (2.148) gmin β −1 ≤ ĝ ≤ β. g (2.149) Notice that bmin , bmax may be state dependent. 2.7.2 Control Chattering Unfortunately, the realization of sliding control has an unpleasant problem. Since k(x) is usually nonzero on the sliding surface and ⎧ ⎪ ⎨+1, if s > 0, sign(x) = 0, if s = 0, ⎪ ⎩ −1, if s < 0 68 2 Basic Nonlinear Control Methods Fig. 2.11 Boundary layer for second order system hence the control is discontinuous, the jump is ±k(x). When the trajectory reaches the sliding surface s(x, t) = 0, the control should perform ±k(x) jump within zero time, which is impossible. The trajectory steps over to the other side of the sliding surface, then returns back to the sliding surface, but cannot stop within zero time there, hence steps over to the first side, and repeats the process. Meanwhile between the two sides of the sliding surface the change of the control is ±k(x) hence the control is chattering with high frequency (control chattering). This high frequency excitation may be harmfully, especially if the system has nonmodeled structural frequency ωstruct of resonance character near to the frequency of chattering, because the magnitude may be amplified which can damage the system. 2.7.2.1 Compensation of Control Chattering Using Boundary Layer For eliminating high frequency control chattering, first the sign function will be changed to the sat function: ⎧ ⎪ ⎨+1, if s > φ, sat(s/φ) = s/φ, if − φ ≤ s ≤ φ, ⎪ ⎩ −1, if s < −φ. The price is that a boundary layer neighboring should be formed around the sliding surface: !   ! B(x, t) = x : !s(x, t)! ≤ φ . Outside the boundary layer we save the original control while inside it u will be interpolated by sat(s/φ). For example, in case of known gain, proportional control and yd = const the boundary layer around the sliding surface s = −ẏ + λ(yd − y) can be characterized by its thickness φ and precision ε, see Fig. 2.11. 2.7 Sliding Control 69 From theoretical point of view, it is an important consequence that both ỹ and ỹ (i) are bounded within the boundary layer. Denoting the differential operator (Laplace d then for x(0) = 0 initial condition transform variable) by p = dt ỹ = ỹ (i) 1 s, (p + λ)n−1 i  λ pi 1 1− = s. s= (p + λ) (p + λ)n−1 (p + λ)n−1−i By repeated application of the convolution theorem, it follows within the boundary layer ! t ! ! ! ! e−λτ s(t − τ ) dτ ! ≤ φ/λ ⇒ |ỹ| ≤ φ =: ε, ! ! λn−1 0 ! !  t ! ! !1 − λe−λτ dτ !! ≤ 2, ! 0 ! (i) ! !ỹ ! ≤ φ 1 2i φ = (2λ)i n−1 = (2λ)i ε. λn−1−i λ If the initial conditions differ from zero, then the error remains still bounded because of the stability of the filter: x = xd + O(ε). Now consider the design of the controller if a boundary layer is present. The thickness φ of the boundary layer will be chosen time varying in order to decrease the tracking error ỹ. Hence, the condition for 12 s 2 will be modified to   1 d 2 |s| = ss ′ ≤ φ̇ − η |s|, 2 dt (2.150) which means that we prescribe stronger condition for φ̇ < 0 than for φ̇ ≥ 0. 2.7.2.2 Boundary Layer and Time Varying Filter for Second Order System Consider first proportional control with constant gain g = 1. Evidently φ̇ and η have opposite effects, hence, assuming k(x) belongs to η without boundary layer, then in case of boundary layer it needs k̃(x) = k(x) − φ̇, u = û + k̃(x) sat(s/φ). (2.151) (2.152) Therefore, within the boundary layer yields s ′ = fˆ − f − k̃(x)s/φ, (2.153) 70 2 Basic Nonlinear Control Methods from which it follows s ′ = −k̃(xd )s/φ + f (xd ) + O(ε). (2.154) Notice that s ′ behaves as a first order linear system along the trajectory xd (t) under the effect of bounded disturbance. Let choose the time constant equal to λ then we obtain the balance condition k̃(xd ) := λ, φ (2.155) from which it follows for the varying boundary layer φ̇ + λφ = k(xd ). (2.156) Consider now the case of state dependent gain. Then using (2.144) and η := η − φ̇, it follows ĝ k̃(x) ≥ k(x) − φ̇, g φ̇ > 0 : (2.157) k̃(x) ≥ k(x) −   ĝ φ̇ = k(x) − β −1 φ̇, g min k̃(xd ) = k(xd ) − β −1 φ̇,     ĝ φ̇ < 0 : k̃ ≥ k(x) + −φ̇ = k(x) − β φ̇, g max (2.158) k̃(xd ) = k(xd ) − β φ̇. Prescribe for variable gain the balance condition k̃(xd ) β := λ, φ (2.159) then we obtain for the thickness of the boundary layer the law φ̇ > 0 : φλ = k(xd ) − β −1 φ̇ β φ̇ < 0 : φλ = k(xd ) − β φ̇ β ⇒ ⇒ φ̇ + λφ = βk(xd ), λ 1 φ̇ + 2 φ = k(xd ). β β (2.160) 2.7.3 Sliding Control of Robot Sliding control elaborated for SISO system can be applied for robots considering each degree of freedom as a separate SISO system and modeling the coupling effects 2.8 Receding Horizon Control 71 as disturbances [8]. Robot: H q̈ + h = τ. Controller: τ = Ĥ u + ĥ, u = ? sliding control.   Closed loop: q̈ = H −1 ĥ − h + H −1 Ĥ u. All the functions are state dependent, x = (q T , q̇ T )T is the state, q is the joint variable. Suppose the estimation Ĥ , ĥ is good enough, especially the elements in the main diagonal of H −1 Ĥ = I + H −1 (Ĥ −1 − H ) are strictly positive (i.e., the gains are positive). Denote (H −1 )i the i-th row of H −1 , let h = ĥ − h and x + = (|x1 |, . . . , |xn |)T for any vector x. Develop for each subsystem a sliding controller:   −1     H Ĥ ij uj , q̈i = H −1 i h + H −1 Ĥ ii ui +  −1  ⎫ ⎪ ⎪ ⎪ ⎬ i=j gi,min ≤ H Ĥ ii ≤ gi,max . functions of q. βi = gi,max /gi,min ⎪ ⎪ ⎪ √ ⎭ Gi = ĝi−1 = 1/ gi,max · gi,min Then the control law for the integrating sliding controller yields   ui = Gi (q) ûi + k̃i (q, q̇) sat(si /φi ) , ûi = q̈di + 2λq̃˙i + λ2 q̃i ,  t si = q̃˙i + 2λq̃i + λ2 q̃i (θ ) dθ − q̃˙i (0) − 2λq̃i (0), 0 ki (q, q̇) ≥ βi "   +  1 − βi−1 |ûi | + H −1 i h+ + j =i # ! −1  ! ! ! Gj |q̈dj | H Ĥ ij + ηi . Notice that during the computation of ki the approximation u ≈ Gj q̈dj was applied in order to avoid the solution of an equation system. For computing k̃i and φi , the formulas elaborated above can be used. 2.8 Receding Horizon Control Model Predictive Control (MPC) is a popular method especially in the process industry where relatively slow process models allow online optimization. Linear predictive control is well elaborated both in frequency (operator) domain [25] and state space [87]. Depending on the type of constraints, optimal prediction leads to 72 2 Basic Nonlinear Control Methods Quadratic (QP) or Nonlinear Programming (NP) which are well supported by existing softwares (e.g., MATLAB Optimization Toolbox). 2.8.1 Nonlinear Receding Horizon Control For nonlinear systems, recent approaches are usually based on new optimum seeking methods suited for the predictive control problem or traditional analytical optimum conditions and gradient based optimization techniques [4, 65]. Basis for the later ones is a general form of the Lagrange multiplicator rule which is also valid in Banach spaces [73]. Here a version of the Receding Horizon Control (RHC) approach of model predictive control will be presented. Using RHC, the optimization is performed in the actual horizon in open loop, then the first element of the optimal control sequence will be performed in closed loop, the horizon is shifted to the right and the process is repeated for the new horizon. For simplification of the notation, the start of the moving horizon will be normalized to zero in every step. It follows from the RHC principle that stability problems can arise if the objective function does not contain any information about the behavior after the horizon. To improve the stability properties, the techniques of Frozen Riccati Equation (FRE) and Control Lyapunov Function (CLF) can be suggested [161]. From engineering point of view, the increase of the horizon length N can improve the stability chance however at the same time the real-time optimization problem will be more complex. Critical may be the presence of constraints in form of limited control set or even more in form of state variable constraints. Although efficient numerical optimization techniques exist for constrained optimization like for example, active set method or interior point method, however in real-time environment simpler gradient based techniques are more realistic. They apply penalty function and control projection instead of the constraints and use conjugate gradient or Davidon–Fletcher–Power (DFP) method for unconstrained optimization. They are especially important for quick systems like robots and vehicles where the necessary sampling time is typically in the order of 10–20 ms for the (high level) control system. Hence, we will concentrate here only on quadratic objective function however the system may be nonlinear. Typical finite horizon nonlinear predictive control problems in discrete time lead to optimization in finite dimensional space where the variables are x = {xi }N i=0 and N −1 u = {ui }i=0 , Qi ≥ 0 and Ri > 0 punish large states and control signals respectively, the optimality criterion to be minimized is 1 F0 (x, u) = 2 F0 (x, u) =: N−1 i=0 N −1 i=0   1 Qi xi , xi  + Ri ui , ui  + QN xN , xN , 2 Li (xi , ui ) + Φ(xN ), (2.161) 2.8 Receding Horizon Control 73 the constraints are the discrete time state equation ϕ(xi , ui ) − xi+1 = 0, (2.162) the control set ui ∈ M and the initial condition a − x0 = 0. If (x ∗ , u∗ ) is the optimal solution then Jx′ (x ∗ , u∗ )x + Ju′ (x ∗ , u∗ )u is the derivative of ( ) J (x, u) = F0 (x, u) + λ0 , a − x0  + λ1 , ϕ(x0 , u0 ) − x1 + · · · ) ( + λN , ϕ(xN−1 , uN−1 ) − xN . (2.163) By introducing the Hamiltonian functions ( ) Hi = λi+1 , ϕ(xi , ui ) + Li (xi , ui ), and using the results of [73] the necessary condition of the optimality are obtained as λN = QN xN , (2.164) λi = ∂Hi /∂xi = Qi xi + (∂ϕ/∂xi )T λi+1 , T ∂Hi /∂ui = Ri ui + (∂ϕ/∂ui ) λi+1 , dJ = N −1 i=0 ( ) ∂Hi /∂ui , u∗i − ui ≤ 0, (2.165) (2.166) ∀ui ∈ M. (2.167) For the control design within the actual horizon, first the initial condition x0 and the approximation of u are needed. The latter may be the solution in the previous horizon shifted to the left and supplying the missing term uN (for example, repeating the last term etc.). Notice that for the very first horizon no initial control sequence is present yet thus special method is necessary to supply it. The optimization repeats the following steps: • • • • Solution of the state equations in x = {xi }N i=0 . Computation of the Lagrange multipliers λi . Computation of the derivatives ∂Hi /∂ui . Numerical optimization based on gradient type methods (gradient, conjugate gradient, DFP etc.) to find u = {ui }N−1 i=0 . Nonpredictive design method should be used to find the initial approximation for the first horizon. If the original system is a continuous time one, then first it can be approximated by a discrete time one, for example, ẋ = fc (x, u) ⇒ xi+1 = xi + Tfc (xi , ui ) =: ϕ(xi , ui ), (2.168) 74 2 Basic Nonlinear Control Methods where T is the sampling time. If the full state can not be measured then x0 can be approximated by using extended Kalman filter or other state estimator/observer [30]. If y = Cx is the system output and ỹ = yd − y is the error, then the cost function can be modified as ỹ = yd − Cx, ( ) 2Li (xi , ui ) = Q̃i ỹi , ỹi + Si xi , xi  + Ri ui , ui , ) ( 2Φ(xN ) = Q̃N ỹN , ỹN , (2.169) (2.170) (2.171) and the derivatives can be computed by ∂Li /∂xi = −C T Q̃i ỹi + Si xi , (2.172) ∂Φ/∂xN = −C Q̃N ỹN . (2.173) T Input constraints are enforced by projecting ui into the constraints set. State constraints can be taken into consideration as an additional penalty added to L(xi , ui ) in the cost function. The weighting in the term Φ(xN ) has great influence on the stability and dynamic behavior of the system under predictive control. The appropriate weighting term can be chosen in experiments. The control design strategy can be summarized in the following steps: • Development of the nonlinear dynamic model of the system. • Optimal (suboptimal, flatness-based etc.) open loop control signal design used for initial approximation of the control sequence in the horizon. • Elaboration of the disturbance model reduced on the system input. • Development of the model based nonlinear predictive controller and its use in closed loop control. The first element of the control sequence in the actual horizon is completed by the feedforward compensations of the disturbance. Remark As can be seen from (2.165) and (2.166), ∂ϕ/∂xi and ∂ϕ/∂ui are needed during numerical optimization therefore their computation time has large influence to the real-time properties. In mechanical systems often appear terms of the form f (x) = A−1 (x)b(x) at the right side of the state equation. See, for example, robots where q̈ = −H −1 (q)h(q, q̇) + H −1 (q)τ . The following formula can speed up the computation: ∂f ∂A −1 ∂b = −A−1 A b + A−1 . ∂xi ∂xi ∂xi (2.174) 2.8.2 Nonlinear RHC Control of 2D Crane We can assume that the path and the nominal control for the differentially flat 2D crane have already been designed, see Sect. 2.1.3. Now we shall concentrate only 2.8 Receding Horizon Control 75 on the steps which can accelerate the RHC computation based on the above remark. Since ϕ(xi , ui ) = xi + Tfc (xi , ui ), hence it is enough to find fc (x, u) and then solve its differentiation. First, the Jacobian (dx) and Hessian (D2x) have to be found for which carefully performed symbolic differentiation can be suggested. By using them the derivatives, ẋ, ẍ, ż, z̈ can be determined. For example, in the case of x(t) yields  T ẋ = dx · L̇1 L̇2 L̇3 , T   T   ẍ = dx · L̈1 L̈2 L̈3 + L̇1 L̇2 L̇3 · D2x · L̇1 L̇2 L̇3 . (2.175) The resulting continuous time dynamic model in MATLAB convention for indexing is the following: A(1, :) = Cθ dx + Sθ dz,   T  b(1) = L̇1 L̇2 L̇3 (−Cθ D2x − Sθ D2z) L̇1 L̇2 L̇3 − Sθ g, A(2, :) = 2Cγ m(Sθ dx − Cθ dz)ρ1 , A(2, 1) = A(2, 1) + J1 /ρ1 ,   T  b(2) = L̇1 L̇2 L̇3 2Cγ m(−Sθ D2x + Cθ D2z)ρ1 L̇1 L̇2 L̇3 + 2Cγ mCθ gρ1 , A(3, :) = m(Sθ dx − Cθ dz)ρ2 , A(3, 2 : 3) = A(3, 2 : 3) + [J2 /ρ2 J2 /ρ2 ],   T  b(3) = L̇1 L̇2 L̇3 m(−Sθ D2x + Cθ D2z)ρ2 L̇1 L̇2 L̇3 + mCθ gρ2 , T  L̈1 L̈2 L̈3 = A−1 b − A−1 (0u1 u2 )T . For RHC A and b have to be differentiated once more, see the above remark. However, this can also be performed by using symbolic computation tools. Notice that symbolic computations can be performed offline, only the final formulas have to be implemented in real-time. An important question remains the appropriate choice of the sampling time T so that the approximation of the continuous time system by the discrete time one using Euler formula is accurate enough. Do not forget that the sampling time has influence on the horizon length because the effective horizon length in time units is N T which has to be fitted to the speed of the system. For large N , the use of basis functions (splines etc.) can be suggested in order to decrease the number of variables to be optimized. Numerical results with RHC control of a small size model of a real cable-driven crane of US Navy can be found in [74]. 76 2 Basic Nonlinear Control Methods 2.8.3 RHC Based on Linearization at Each Horizon Till now we discussed RHC based on the nonlinear dynamic model. Of course there exist the possibility to linearize the system at the beginning of each horizon if the linear approximation is accurate enough in the entire horizon. Applying RHC for the LTV or LTI model, valid in the horizon, simple control implementations can be found. Especially, if the objective function is quadratic and the constraints are linear equalities and/or linear inequalities then in each horizon a Quadratic Programming (QP) problem has to be solved for which efficient active set or interior point methods are available. Especially, if there are no constraints then analytical solution is possible which enormously simplifies the real-time conditions. 2.9 Closing Remarks The main goal of the chapter was to give a survey of the basic nonlinear control methods. Central question was the introduction of the main control ideas and the presentation of the results supporting the solutions of nonlinear control problems. In the presentation, a compromise was sought between theoretical precision and engineering view. Nonlinear control systems have to be tested before application for which simulation methods can be applied. During simulation, the real system is substituted by its model. For modeling of continuous time dynamic systems, nonlinear and linear models can be used. Important problem is the computation of the transients for which in case of nonlinear systems numerical methods should be used. For this purpose, simple (Euler) and more accurate (Runge–Kutta and predictor–corrector) methods were shown. For linear systems, the computation can be simplified by using the fundamental matrix (LTV systems) or the exponential matrix or transfer function (LTI systems). Nonlinear systems can be divided into different classes. Nonlinear systems often have to satisfy kinematic constraints in the space of generalized coordinates and their derivatives. The constraints can be divided into holonomic and nonholonomic ones. The kinematic constraints are usually given as a set of equations linear in the derivatives of the generalized coordinates but nonlinear in the generalized coordinates (Pfaffian constraints). If there exist a nonlinear coordinate transformation (diffeomorphism) after which the constraints can be eliminated and the number of variables can be reduced, then the system is holonomic, otherwise it is nonholonomic. In the dynamic model of nonholonomic system appear constraint forces, too, described by parts of the kinematic constraints and Lagrange multipliers. A method was shown how to determine the Lagrange multipliers from the dynamic model. Based on energy concept, the Hamiltonian form was introduced and an experiment was performed to simplify the Hamiltonian using a coordinate transformation which is related to the kernel space of the kinematic constraints. Conditions for omitting the kinematic constraints were given based on Poisson bracket. 2.9 Closing Remarks 77 An other nonlinear system class is the class of differentially flat systems which are Lie–Bäcklund equivalent to a trivial system. Differentially flat systems can be linearized by dynamic state feedback and a coordinate transformation. For such system, there is a flat output variable, having the same dimension as the system input, such that the state and the input can be expressed by the flat output variable and its derivatives of finite order. The resulting square system is left and right input–output invertible with trivial zero dynamics. Path design for differentially flat systems is a simple problem because it can be reduced to smooth interpolation between points in the space of the flat output variable. The result is not only the path but also the nominal control along the path. The nominal (open loop) control can be used for feedforward compensation or may be the basis for linearization along the path. The differential flatness property is especially useful for underactuated systems where not every path can be realized by control. Modeling nonlinear systems is an often arising problem. Three simple examples were considered to illustrate how to find the dynamic model of physical systems. The first example was the dynamic modeling of an inverted pendulum mounted on a moving car. To find the model the kinetic and potential energy was found and the Lagrange equations were applied. The inverted pendulum is an underactuated system. The second example was the dynamic model of a quarter car using active suspension. The flexible connections between body and wheel were described by viscous damping factors and spring constants. The resulting model is an LTI system with active force input. The third model was a 2 DOF robot having revolute joints and moving in vertical plane. First the position, velocity and angular velocity of the links were determined then the dynamic model was derived from the Lagrange equations. The coupling and centripetal effects were determined from the kinetic energy by using differentiation. The dynamic model is nonlinear. Stability is a central problem of every control system, especially for nonlinear ones. First the stability, uniform stability and asymptotic stability of any solution of the time varying nonlinear state equation were defined. Then the solution was transformed to the equilibrium point at zero satisfying some conditions. The idea of positive definite Lyapunov function was introduced and based on it Lyapunov’s first theorem (direct method) was formulated for stability and asymptotic stability. For time varying nonlinear systems, Lyapunov’s second theorem (indirect method) was formulated allowing to conclude from the asymptotic stability of the linearized system to the asymptotic stability of the nonlinear one. For LTI system, the stability is equivalent to the solvability of the Lyapunov equation. For time invariant nonlinear system, LaSalle’s theorem was formulated which gives also information about the attraction domain of the stability and brings into relation the asymptotic stability with the maximal invariant set. The classes of K and KL functions were introduced and two theorems were given for showing the relation to the uniformly asymptotically and the globally uniformly asymptotically stable systems. The input-to-state stability (ISS) was defined and a sufficient condition was formulated for ISS. Three forms of the Barbalat’s lemma were formulated which give sufficient conditions for zero limit value of scalar-valued time functions as the time goes to infinity. As an immediate corollary of Barbalat’s lemma, a Lyapunov-like condition 78 2 Basic Nonlinear Control Methods was given which looks like an invariant set theorem in Lyapunov stability analysis of time varying systems. The system behavior in Lp function spaces was characterized with the nonlinear gain. A small gain theorem was formulated for closed loop systems in such function spaces. The discussions implicitly assumed zero initial conditions. Passive systems, amongst them passive, strictly input passive and strictly output passive systems were introduced and a passivity theorem in Lp spaces was formulated which allows to conclude from the stability of the open loop subsystems to the stability of the closed loop if only one external input is active. Dissipative systems allow to integrate the Lyapunov theory and the stability in L2 spaces in some sense. Special is that the number of inputs and outputs must be equal because the scalar product plays a central role. Dissipativity is defined by using the idea of storage function and supply rate which have to satisfy the dissipativity inequality. Special classes of dissipative systems are the passive, the strictly input passive, the strictly output passive systems and the system having finite L2 gain. The dissipative property can be expressed also in differential form which can be checked easier. Three theorems were formulated giving connection between passivity and asymptotic stability, stability conditions of interconnected systems without external input, and a small gain theorem for interconnected dissipative systems without external inputs. These theorems will play important role in the stability investigations of multiagent systems moving in formation. The control of nonlinear dynamic systems is a complex problem. In order to simplify the problem, it is a useful idea to linearize the nonlinear system using appropriate nonlinear coordinate transformation and then compensate it further using linear control methods. Sufficient condition for the input–output linearization of the MIMO nonlinear input affine system is the existence of vector relative degree and the nonsingularity of a special matrix composed from repeated Lie derivatives. If the conditions are satisfied, then a nonlinear coordinate transformation has to be determined which can be formed by completing a set of functions appearing during the process. Then an internal static nonlinear state feedback can be applied which brings the system onto Brunovsky canonical form equivalent to decoupled integrators. However, beside of this linear system appears also a nonlinear one called zero dynamics which is not observable (has no influence on the output). The zero dynamics has to be stable in order to assure closed loop stability. For SISO system, the zero dynamics is independent from the system input, but in this case the Frobenius equation has to be solved in order to construct the nonlinear coordinate transformation. For MIMO systems, this process cannot be applied and the zero dynamics depends also on the system input. The number of inputs and outputs must be the same in the construction, however the output may differ from the physical system output. The weak point of the method is the large number of symbolic computations, the solution of the Frobenius (partial differential) equation in SISO case or the computation of the nonlinear coordinate transformation in the MIMO case, and the computation of the inverse coordinate transformation which is also needed to find the zero dynamics. Some other nonlinear control methods were also discussed in detail, especially because of the complexity of input–output linearization. In flatness control, only the 2.9 Closing Remarks 79 path planning problem for differentially flat underactuated systems was considered. It was demonstrated for the planar (2D) cable-driven crane. The flat output variable was determined and it was shown that the state and the control can be computed from the flat output and its derivatives. Then the path planning was reduced to polynomial interpolation. Backstepping control is a popular method for controlling nonlinear systems given in strictly feedback form. The dimension of the portions of the state vector, the input vector and the physical output must be the same and some invertability conditions have to be satisfied. It was assumed that the path was designed in the space and parametrized in a scalar path parameter which has to satisfy a parameter speed constraint, too. The design is decomposed in steps, in every step a nonlinear virtual control is developed, the state variable portion is changed to a new variable and a partial Lyapunov function is designed. At the end, the last Lyapunov function assures the closed loop stability, but high level control has to eliminate the parameter speed error. Especially, if the path is given as a function of time then the closed loop is already stable without high level control. Backstepping plays also important role in formation stabilization. Sliding control is a kind of variable structure control. The sliding variable is the output of a stable filter with equal time constants driven by the output error. The sliding surface belongs to the zero value of the sliding variable. If the system remains on the sliding surface, then the error goes exponentially to zero. If the sliding variable is initially not on the sliding surface, then the sliding surface can be reached within limited time. A method was presented for the design of sliding controller satisfying the above conditions for SISO systems including proportional and integral control with constant or state dependent gain. About the functions in the system model only upper bounds were assumed. Since the method can cause high frequency control chattering in the neighborhood of the sliding surface therefore a boundary layer was introduced whose thickness determines the precision. Hence, the thickness was chosen state dependent in order to decrease the error. It was shown that with some modifications the sliding control developed for SISO systems can also be applied for robots considering coupling effects amongst the joints as disturbances. Receding Horizon Control (RHC) is a kind of optimal control considered in discrete time in order to make the problem finite dimensional. Using RHC, the optimization is performed in the actual horizon in open loop, then the first element of the optimal control sequence is performed in closed loop, the horizon is shifted to the right and the process is repeated for the new horizon. It follows from the RHC principle that stability problems can arise if the objective function does not contain any information about the behavior after the horizon. It is clear that the increase of the horizon length can improve the stability chance however at the same time the real-time optimization problem becomes more complex. Critical is the presence of constraints for the control and/or state variables. In order to avoid time consuming constrained optimization in real-time, the method of penalty function, control projection and unconstrained optimization using gradient based techniques (conjugate gradient, DFP etc.) were suggested. Under these assumptions, the RHC algorithm for nonlinear discrete time system was elaborated. The critical steps of the compu- 80 2 Basic Nonlinear Control Methods tations were illustrated in the example of the RHC control of the planar (2D) crane which is differentially flat but underactuated. Another RHC strategy may be the linearization at the beginning of each horizon. The linear approximation should be accurate enough in the entire horizon. Considering the actual LTV or LTI model and assuming the objective function is quadratic and the constraints are linear equalities and/or linear inequalities, then in each horizon a Quadratic Programming (QP) problem has to be solved for which efficient numerical methods are available. This approach can also be applied in safety critical systems where high level RHC computes optimal reference signals for the low level control system. http://www.springer.com/978-1-84996-121-9