CHAPTER 9
NONLINEAR AND MULTIVARIABLE CONTROL

The Effect of Joint Flexibility

In Chapter 8 we considered the effect of joint flexibility and showed for a lumped model of a single-link robot that a PD control could be designed for set-point tracking. In this section we will discuss the analogous result in the general case of an n-link manipulator.

We first derive a model similar to Equation (9.6) to represent the dynamics of an n-link robot with joint flexibility. For simplicity, assume that the joints are revolute and are actuated by permanent magnet DC motors. We model the flexibility of the i^th joint as a linear torsional spring with stiffness constant k_i, for i = 1, …, n.

Referring to Figure 9.1, we let q₁ = [θ₁, ⋅⋅⋅, θ_n]^T be the vector of DH-joint variables and q₂ = be the vector of motor shaft angles (reflected to the link side of the gears). Then q₁ − q₂ is the vector of elastic joint deflections. Neglecting the effect of link motion on the kinetic energy of the actuators, as we did in Equation (9.6), the kinetic and potential energies of the manipulator are given by

(9.15)

(9.16)

where J and K are diagonal matrices of inertia and stiffness constants, respectively.

(9.17)

Note that we have simply augmented the kinetic energy of the rigid joint robot with the kinetic energy of the actuators and we have likewise augmented the potential energy due to gravity of the rigid joint model with the elastic potential energy of the linear springs at the joints. It is now straightforward to compute the Euler–Lagrange equations for this system as (Problem 9–2)

(9.18)

For the problem of set-point tracking with PD control, consider a control law of the form

(9.19)

where and q^d is a vector of constant set points, and suppose again that the gravity vector g(q₁) = 0. To show asymptotic tracking for the closed-loop system consider the Lyapunov function candidate

(9.20)

We leave it as an exercise (Problem 9–3) to show that an application of LaSalle’s theorem proves global asymptotic stability of the system. Consequently, the motor and link angles are equal in the steady state. Thus, one may choose the set point q^d₂ as a vector of desired DH variables.

If gravity is present, then it is not apparent how one can achieve gravity compensation in a manner similar to the rigid joint case. We will address this question later in the context of feedback linearization in Chapter 12.

9.3.1 Joint Space Inverse Dynamics

Consider again the dynamic equations of an n-link rigid robot in matrix form

(9.21)

The idea of inverse dynamics is to seek a nonlinear feedback control law

(9.22)

which, when substituted into Equation (9.21), results in a linear closed-loop system. For general nonlinear systems, such a control law may be quite difficult or impossible to find. In the case of the manipulator dynamics given by Equations (9.21), however, the problem is actually easy. By inspecting Equation (9.21) we see that if we choose the control u according to the equation

(9.23)

then, since the inertia matrix M is invertible, the combined system given by Equations (9.21)–(9.23) reduces to

(9.24)

The term a_q represents a new input that is yet to be chosen. Equation (9.24) is known as the double integrator system as it represents n uncoupled double integrators. Equation (9.23) is called the inverse dynamics control and achieves a rather remarkable result, namely that the system given by Equation (9.24) is linear and decoupled. This means that each input can be designed to control a SISO linear system. Moreover, assuming that is a function only of q_k and , then the closed-loop system will be decoupled.

Since a_q can now be designed to control a linear second-order system, an obvious choice is to set

(9.25)

where , , K₀, K₁ are diagonal matrices with diagonal elements consisting of position and velocity gains, respectively and the reference trajectory

(9.26)

defines the desired time history of joint positions, velocities, and accelerations. Note that Equation (9.25) is nothing more than a PD control with feedforward acceleration as defined in Chapter 8.

Substituting Equation (9.25) into Equation (9.24), results in

(9.27)

A simple choice for the gain matrices K₀ and K₁ is

(9.28)

which results in a decoupled closed-loop system with each joint response equal to the response of a critically damped linear second-order system with natural frequency ω_i. As before, the natural frequency ω_i determines the speed of response of the joint, or equivalently, the rate of decay of the tracking error.

The inverse dynamics approach is important as a basis for control and it is worthwhile to examine it from alternative viewpoints. We can give a second interpretation of the control law (9.23) as follows. Consider again the manipulator dynamics (9.21). Since M(q) is invertible for all , we may solve for the acceleration of the manipulator as

(9.29)

Suppose we were able to specify the acceleration as the input to the system. That is, suppose we had actuators capable of producing directly a commanded acceleration (rather than indirectly by producing a force or torque). Then the dynamics of the manipulator, which is after all a position control device, would be given as

(9.30)

where a_q(t) is the input acceleration vector. This is again the familiar double integrator system. Note that Equation (9.30) is not an approximation in any sense; rather it represents the actual open-loop dynamics of the system provided that the acceleration is chosen as the input. The control problem for the system (9.30) is now easy and the acceleration input a_q can be chosen as before according to Equation (9.25).

In reality, however, such “acceleration actuators” are not available to us and we must be content with the ability to produce a generalized force (torque) u_i at each joint i. Comparing Equations (9.29) and (9.30) we see that the torque input u and the acceleration input a_q of the manipulator are related by

(9.31)

Solving Equation (9.31) for the input torque u(t) yields

(9.32)

which is the same as the previously derived expression (9.23). Thus, the inverse dynamics can be viewed as an input transformation that transforms the problem from one of choosing torque input commands to one of choosing acceleration input commands.

Note that the implementation of this control scheme requires the real-time computation of the inertia matrix and the vectors of Coriolis, centrifugal, and gravitational generalized forces. An important issue therefore in the control system implementation is the design of the computer architecture for the above computations.

Figure 9.2 illustrates the so-called inner-loop/outer-loop control architecture. By this we mean that the computation of the nonlinear control law (9.23) is performed in an inner loop with the vectors , and a_q as its inputs and u as output. The outer loop in the system is then the computation of the additional input term a_q. Note that the outer-loop control a_q is more in line with the notion of a feedback control in the usual sense of being error driven. The design of the outer-loop feedback control is in theory greatly simplified since it is designed for the plant represented by the dotted lines in Figure 9.2, which is now a linear system.

9.3.2 Task Space Inverse Dynamics

As an illustration of the importance of the inner-loop/outer-loop control architecture, we will show that tracking in task space can be achieved by modifying the outer-loop control a_q in Equation (9.24) while leaving the inner-loop control unchanged. Let represent the end-effector pose using any minimal representation of SO(3). Since X is a function of the joint variables , we have

(9.33)

(9.34)

where J = J_a is the analytical Jacobian of Equation (4.88). Given the double integrator system (9.24) in joint space we see that if a_q is chosen as

(9.35)

the result is a double integrator system in task space coordinates

(9.36)

Given a task space trajectory X^d(t) satisfying the same smoothness and boundedness assumptions as the joint space trajectory q^d(t), we may choose a_X as

(9.37)

so that the task space tracking error satisfies

(9.38)

Therefore, a modification of the outer-loop control term a_q achieves a linear and decoupled system directly in the task space coordinates without the need to compute a joint space trajectory and without the need to modify the nonlinear inner-loop control.

Note that we have used a minimal representation for the orientation of the end effector in order to specify a trajectory . In general, if the end-effector coordinates are given in SE(3), then the Jacobian J in the above formulation will be the geometric Jacobian of Equation (4.41). In this case

(9.39)

and the outer-loop control

(9.40)

applied to Equation (9.24) results in the system

(9.41)

(9.42)

(9.43)

Although, in this latter case, the dynamics have not been linearized to a double integrator, the outer-loop terms a_v and a_ω may still be used to define control laws to track end-effector trajectories in SE(3).

In both cases we see that nonsingularity of the Jacobian is necessary to implement the outer-loop control. Avoiding singularities, therefore, is an important aspect of motion planning and trajectory planning that we have considered in Chapter 7. Also, if the robot has more or fewer than six joints, then the Jacobians are not square. In this case, other schemes have been developed using, for example, the Jacobian pseudoinverse. See [32] for details.

9.3.3 Robust Inverse Dynamics

The inverse dynamics approach relies on exact cancellation of nonlinearities in the robot equations of motion. The practical implementation of inverse dynamics control requires consideration of various sources of uncertainties such as modeling errors, unknown loads, and computation errors. Let us return to the Euler–Lagrange equations of motion

(9.44)

and write the inverse dynamics control input u as

(9.45)

where the notation represents the computed or nominal value of ( · ) and indicates that the theoretically exact inverse dynamics control is not achievable in practice due to the uncertainties in the system. The error or mismatch is a measure of one’s knowledge of the system parameters.

If we substitute Equation (9.45) into Equation (9.44) we obtain, after some algebra (Problem 9–8),

(9.46)

where the nonlinear function is given by

(9.47)

and is called the uncertainty. Henceforth, we will drop the arguments of for simplicity. However, it is important to keep in mind that η is, in general, a function of both the state and the reference trajectory.

We define E = E(q) as

(9.48)

which allows us to express the uncertainty η as

(9.49)

The system (9.46) is still nonlinear and coupled due to the uncertainty and, therefore, we have no guarantee that the outer-loop control given by Equation (9.25) will satisfy desired tracking performance specifications.

In this section we show how to modify the outer-loop control (9.25) to guarantee global convergence of the tracking error for the system (9.46). The method we outline below is often called Lyapunov redesign. In this approach we set the outer-loop control a_q as

(9.50)

where δa is an additional term to be designed. In terms of the tracking error

(9.51)

we may write Equations (9.46) and (9.50) as

(9.52)

where

(9.53)

Thus, the double integrator is first stabilized by the linear feedback term , and the additional control term δa should be designed to overcome the potentially destabilizing effect of the uncertainty η. The basic idea is to assume that we are able to compute a bound ρ(e, t) ⩾ 0 on the uncertainty η as

(9.54)

and design the additional input term δa to guarantee ultimate boundedness of the error trajectory e(t) in Equation (9.52). Note that the bound ρ is in general a function of the tracking error e and time reflecting the fact that η likewise depends on e as well as on t through a time-varying reference trajectory.

Returning to our expression for the uncertainty η and substituting for a_q from Equation (9.50) we have

(9.55)

Let us assume that we can find constants α < 1, γ₁, and γ₂, together with possibly time-varying γ₃ such that

(9.56)

We note that the condition determines how close our estimate must be to the true inertia matrix. Suppose that we have constant bounds on M^{− 1} as

(9.57)

Remark 9.1.

If we choose the estimated inertia matrix as

(9.58)

then it can be shown that

(9.59)

The point is that there is always a choice for that satisfies the condition ‖E‖ < 1 provided upper and lower bounds and can be determined.

Next, assume, for the moment, that ‖δa‖ ⩽ ρ(e, t) which must then be checked a posteriori. It follows that

(9.60)

which, since α < 1, defines ρ as

(9.61)

Since K₀ and K₁ are chosen so that the matrix A in Equation (9.52) is Hurwitz, we may choose Q > 0 and let P > 0 be the unique symmetric positive definite matrix satisfying the Lyapunov equation

(9.62)

Defining the control δa according to

(9.63)

it follows that the Lyapunov function

(9.64)

satisfies along solution trajectories of Equation (9.52). To show this result, we compute

(9.65)

For simplicity set w = B^TPe and consider the second term w^T{δa + η} in the above expression. If w = 0 this term vanishes, and hence,

(9.66)

independent of the choice of δa. For w ≠ 0 we have

(9.67)

and hence, using the Cauchy–Schwartz inequality we have

(9.68)

since ‖η‖ ⩽ ρ. Therefore

(9.69)

and the result follows. Finally, note from Equation (9.67) that ‖δa‖ ⩽ ρ as required.

Since the above control term δa is discontinuous on the subspace defined by B^TPe = 0, solution trajectories on this subspace are not well defined in the usual sense. One may define solutions in a generalized sense, the so-called Filippov solutions [46]. In practice, the discontinuity in the control results in the phenomenon of chattering, where the control switches rapidly between the control values in (9.63).

To avoid control chattering, we may implement a continuous approximation to the above discontinuous control as

(9.70)

In this case, since the control signal given by Equation (9.70) is continuous and differentiable except at ‖B^TPe‖, a solution to the system (9.52) exists for any initial condition and we can prove the following result.

Theroem 9.1.

Proof: As before, choose V(e) = e^TPe and compute

(9.71)

(9.72)

with ‖w‖ = ‖B^TPe‖ as above. For ‖w‖ ⩾ ε the argument proceeds as above and . For ‖w‖ ⩽ ε the second term in (9.72) becomes

This expression attains a maximum value of when . Thus, we have

(9.73)

provided

(9.74)

Using the relationship (see B.5)

(9.75)

from Equation (B.5), where λ_min(Q), λ_max(Q) denote the minimum and maximum eigenvalues, respectively, of the matrix Q, we have that if

(9.76)

or, equivalently

(9.77)

Let S_δ denote the smallest level set of V containing B(δ), the ball of radius δ and let B_r denote the smallest ball containing S_δ. Then all solutions of the closed-loop system are u.u.b. with respect to B_r. The situation is shown in Figure 9.3. All trajectories will eventually reach the boundary of S_δ since is negative definite outside of S_δ and hence remain inside the ball B_r.

Note that the radius of the ultimate boundedness set, and hence, the magnitude of the steady-state tracking error, is proportional to the product of the uncertainty bound ρ and the constant ε.

9.3.4 Adaptive Inverse Dynamics

This section is concerned with the adaptive control problem in which an estimation scheme is used to generate estimated values of parameters that are then used in lieu of the true parameters.

Once the linear parametrization property for manipulators became widely known in the mid-1980’s, the first globally convergent adaptive control results began to appear. These first results were based on the inverse dynamics approach discussed above. Consider the plant given by Equation (9.44) and the control given by Equation (9.45) as above, but now suppose that the parameters appearing in Equation (9.45) are not fixed as in the robust control approach, but are time-varying estimates of the true parameters. Substituting Equation (9.45) into Equation (9.44) and setting

(9.78)

it can be shown (Problem 9–11), using the linear parametrization property, that

(9.79)

where is the regressor function defined in Equation (6.122) and , where is the estimate of the parameter vector θ. In state space we write the system (9.79) as

(9.80)

where

(9.81)

with K₀ and K₁ chosen as before as diagonal matrices of positive gains so that A is a Hurwitz matrix. Let P be the unique symmetric, positive definite matrix P satisfying the matrix Lyapunov equation

(9.82)

and choose the parameter update law as

(9.83)

where Γ is a constant, symmetric, positive definite matrix. Then, global convergence to zero of the tracking error with all internal signals remaining bounded can be shown using the Lyapunov function

(9.84)

To see this we calculate as (Problem 9–12)

(9.85)

the latter term following since θ is constant, that is, . Using the parameter update law (9.83) we have

(9.86)

Remark 9.2.

From Equation (9.86) it follows that the position tracking errors converge to zero asymptotically and the parameter estimation errors remain bounded. We will not go through the details of the proof in this section. The argument is similar to that of the passivity-based, adaptive control approach of the next section so we will defer the details until later. In order to implement this adaptive inverse dynamics scheme, however, one notes that the acceleration is needed in the parameter update law and that must be invertible. The need for the joint acceleration in the parameter update law presents a serious challenge to its implementation. Acceleration sensors are noisy and introduce additional cost whereas calculating the acceleration by numerical differentiation of position or velocity signals is not feasible in most cases. The invertibility of can be enforced in the algorithm by resetting the parameter estimate whenever would otherwise result in becoming singular. The passivity-based approaches that we treat next remove both of these impediments.

9.4.1 Passivity-Based Robust Control

In this section we use the passivity-based approach above to derive an alternative robust control algorithm that exploits both the skew symmetry property and linearity in the parameters and leads to a much easier design in terms of computation of uncertainty bounds. We modify the control input (9.88) as

(9.95)

where K, Λ, v, a, and r are given by (9.89). In terms of the linear parametrization of the robot dynamics, the control (9.95) becomes

(9.96)

and the combination of Equation (9.95) with Equation (9.87) yields

(9.97)

We now choose the term in Equation (9.96) as

(9.98)

where θ₀ is a fixed nominal parameter vector and δθ is an additional control term. The system (9.97) then becomes

(9.99)

where is a constant vector and represents the parametric uncertainty in the system. If the uncertainty can be bounded by finding a nonnegative constant ρ ⩾ 0 such that

(9.100)

then the additional term δθ can be designed according to

(9.101)

Using the same Lyapunov function candidate from Equation (9.91) above, we can show uniform ultimate boundedness of the tracking error. Carrying out the details of the calculation of yields

(9.102)

Uniform ultimate boundedness of the tracking error follows with the control δθ from Equation (9.101) exactly as in the proof of Theorem 9.5. The details are left as an exercise (Problem 9–14).

Comparing this approach with the approach in Section 9.3.3 we see that finding a constant bound ρ for the constant vector is much simpler than finding a time-varying and state-dependent bound for η in Equation (9.47). The bound ρ in this case depends only on the inertia parameters of the manipulator, while ρ(x, t) in Equation (9.54) depends on the manipulator state vector, the reference trajectory and, in addition, requires some assumptions on the estimated inertia matrix .

9.4.2 Passivity-Based Adaptive Control

In the adaptive approach the vector in Equation (9.96) is taken to be a time-varying estimate of the true parameter vector θ. Combining the control law, Equation (9.95), with Equation (9.87) yields

(9.103)

The parameter estimate may be computed using standard methods of adaptive control such as gradient or least squares. For example, using the gradient update law

(9.104)

together with the Lyapunov function

(9.105)

results in global convergence of the tracking errors to zero and boundedness of the parameter estimates. Computing along trajectories of the system (9.103), we obtain

(9.106)

Substituting the expression for from the gradient update law (9.104) into Equation (9.106) yields

(9.107)

where e and Q are defined as before, showing that the closed-loop system is stable in the sense of Lyapunov. It follows that the tracking error e(t) converges asymptotically to zero and the estimation error is bounded. In order to show this we note that since is nonincreasing from Equation (9.107), the value of V(t) can be no greater than its value at t = 0. Since V consists of a sum of nonnegative terms, this means that each of the terms r, , and are bounded as functions of time.

With regard to the tracking error, , , we also note that is quadratic in the error vector e(t). Integrating both sides of Equation (9.107) gives

(9.108)

As a consequence, the tracking error vector e(t) is a square integrable function. Such functions, under some mild additional restrictions, must tend to zero as t → ∞. Specifically, we may appeal Barbalat’s Lemma from Appendix C. Since both and have already been shown to be bounded, it follows that is also bounded. Therefore, we have that is square integrable and its derivative is bounded. Hence, the tracking error as t → ∞.

To show that the velocity tracking error also converges to zero, one must appeal to the equations of motion (9.103), from which one may argue that the acceleration is bounded. It follows that the velocity error asymptotically converges to zero provided that the reference acceleration is bounded.