Nonlinear Receding-Horizon Control of Rigid Link Robot Manipulators

Abstract

The approximate nonlinear receding-horizon control law is used to treat the trajectory tracking control problem of rigid link robot manipulators. The derived nonlinear predictive law uses a quadratic performance index of the predicted tracking error and the predicted control effort. A key feature of this control law is that, for their implementation, there is no need to perform an online optimization, and asymptotic tracking of smooth reference trajectories is guaranteed. It is shown that this controller achieves the positions tracking objectives via link position measurements. The stability convergence of the output tracking error to the origin is proved. To enhance the robustness of the closed loop system with respect to payload uncertainties and viscous friction, an integral action is introduced in the loop. A nonlinear observer is used to estimate velocity. Simulation results for a two-link rigid robot are performed to validate the performance of the proposed controller.

1. Introduction

During recent years much emphasis has been placed on flexible manufacturing processes where the most important factors are quality, costs and time. Both fast motion in unconstrained space and mechanical interaction with the environment are required in most manufacturing systems. Industrial robots are often used to meet this demand and to perform various tasks such as material assembling, painting or welding. To accomplish these tasks efficiently and accurately, several control approaches have been proposed in the literature. Among these, a simple PD-control scheme that achieves satisfactory performance (Spong M. W. & Vidyasagar M. 1989) in the absence of gravity. However, robot manipulator is highly nonlinear system with coupling between joints and the gravity effects. The computed torque control or feedback linearization control has been also used to achieve best tracking performance. The implantation of the computed torque controller requires exact knowledge of the robot dynamics. Unfortunately, model uncertainties are frequently encountered in robotics due to unknown or changing payload and friction. These model uncertainties may decrease significantly the performance of this method in terms of tracking accuracy. Therefore, to achieve acceptable performance, even when all kinds of uncertainties are encountered, numerous robust control algorithms have been used like the variable structure approach (Slotine J. J. E. & Sastry S. S. 1983), robust adaptive approach (Ortega R. & Spong M. W. 1989), (Lee K. W. & Khalil H. K. 1997), (Canudas C. W. & Fixot N. 1992), (Spong M. W. 1992) and nonlinear H_∞ approach (Chen B. S. et al 1994). A first survey of early results in robust control has been compiled in (Abdullah C. et al; 1991) and the second survey of recent results has been given in (Sage H. G. et al, 1999). Finally the robotic applications require effective control laws that achieve accurate tracking of fast motion despite the variations of inertia and gravitational load of the manipulator during operation.

Model predictive control of linear systems has received considerable attention in the last decade due to its robustness with respect to model uncertainties. However, many systems are inherently nonlinear. Since linear models are often inadequate to describe accurately the process dynamics, then nonlinear models should be used. Much effort has been made to extend linear predictive control to nonlinear systems (Michalska H & Mayne D. Q. 1993). The disadvantages of the proposed approach is the heavy online computation burden that causes two important problems in implementation of the nonlinear predictive control. One is the computation delay that cannot be ignored and the other is the global solution that cannot be guaranteed in each optimization problem. The application of these nonlinear control laws to nonlinear systems characterized by fast dynamics (such robotics) sound like unusual proposal. To overcome the computation burden, several nonlinear predictive laws have been developed in (Ping L. 1995), (Singh S. M. 1995), (Souroukh M. & Kravaris C. 1996), (Chen W. H. et al, 2003), where the one step ahead predictive output error is obtained by expanding the output signal and the is used to derive the offline control laws. In this paper, the nonlinear receding-horizon controller proposed in (Ping L. 1998) is applied to robot manipulator to achieve position angular tracking objectives. To derive the control law, the predictive tracking error and the predicted control effort are minimized over a fixed time horizon. This approximate nonlinear controller is given in a closed form and thus no online optimization is required. Moreover, to increase the robustness of the control algorithm with regard to model uncertainties, we propose to introduce an integral action in the loop. The well-known Lyapunov based theory is used to show the asymptotic stability of the closed loop system in matched or mismatched case.

The major drawback of the proposed schemes is the requirement of measurement of motor speed. Speed measurements increase cost and impose constraints on the achievable bandwidth. Thus, to overcome this problem a nonlinear observer is used to estimate position and velocity angular of robot manipulator.

The outline of this paper proceeds as follows. In the next section, a dynamic model of robot manipulator is presented. In section 3, the approximate receding-horizon control scheme is developed to allow position angular tracking of a desired references trajectory. Stability analysis and robustness are treated in section 4. The high gain observer used to estimate unmeasured output elements (velocity angular) is presented in section 5. Simulation results are given in section 6. In the last, we conclude with some remarks.

2. Dynamic model of rigid link robot manipulators

The Euler–Lagrange equations are a tool from analytical mechanics that can be used to derive the equations of motion for a mechanical system. In this approach the joint q(t) are considered as generalized coordinates. The kinetic energy of a robot manipulator with n degrees of freedom can be calculated as:

$\Gamma (\mathbf{q},\stackrel{.}{\mathbf{q}})=\frac{1}{2}\stackrel{.}{\mathbf{q}}{(t)}^T\mathbf{D}(\mathbf{q})\stackrel{.}{\mathbf{q}}(t),$

where D(q) is the inertia matrix. Let U(q) : ℜⁿ → ℜ be a continuously differentiable function, called the potential energy. The Lagrangian function is defined (Spong M. W. & Vidyasagar M. 1989) by:

$L(\mathbf{q},\stackrel{.}{\mathbf{q}})=\Gamma (\mathbf{q},\stackrel{.}{q})-U(\mathbf{q}).$

The dynamics of the manipulator are described by Lagrange's equations:

$\frac{d}{dt}\frac{\partial L(\mathbf{q},\stackrel{.}{\mathbf{q}})}{\partial {\stackrel{.}{\mathbf{q}}}_k}-\frac{\partial L(\mathbf{q},\stackrel{.}{\mathbf{q}})}{\partial {\mathbf{q}}_k}={\mathbf{u}}_k,\mathrm{k}=1,\dots \dots ,n,$

reference signal in a r_i^th order Taylor series, ri is the relative degree of the i^th element of the output. Then, the continuous minimization of the predictive tracking errors

where u₁, u₂,……u_n represent generalized input torques. Inserting the kinetic energy and the potential energy for the Lagrange L(q, q̇) above leads to the matrix description:

$\mathbf{D}(\mathbf{q})\stackrel{..}{\mathbf{q}}+\mathbf{C}(\mathbf{q},\stackrel{.}{\mathbf{q}})\stackrel{.}{\mathbf{q}}+\mathbf{G}(\mathbf{q})+{\mathbf{f}}_r={\mathbf{u}}_r,$

where q(t) ε ℜⁿ is the vector of the generalized coordinates representing the angular joint positions and controlled with the driving torques u _r ε ℜⁿ, D(q) ε ℜ^nxn, D(q)=D(q)^T>0, is the link inertia matrix, C(q, q̇) q̇ ε ℜⁿ is the vector of the coriolis and centripetal torques, G(q) ε ℜⁿ is the vector of gravitational torques and f _r represents friction torques acting on the joints. This is described in [9], when only the mechanical parts of actuators dynamics are included. The dynamic model of a rigid robot manipulator becomes:

$\mathbf{M}(\mathbf{q})\stackrel{..}{\mathbf{q}}+\mathbf{C}(\mathbf{q},\stackrel{.}{\mathbf{q}})\stackrel{.}{\mathbf{q}}+\mathbf{G}(\mathbf{q})+\mathbf{f}=\mathbf{u},$

(1)

with:

$\begin{array}{ll}\mathbf{u}=N{\mathbf{u}}_m,\hfill & \mathbf{M}(\mathbf{q})=\mathbf{D}(\mathbf{q})+N^2{\mathbf{J}}_m=\mathbf{D}(\mathbf{q})+\mathbf{J}\hfill \end{array}$

and f = f _r + N f _m ,

where:

N is the diagonal matrix of the gear ratios.

u _m is the vector of torque supplied by the actuators.

f _m is the vector friction torque acting on the motors.

J _m is the diagonal matrix containing the effective motors' inertia.

It is assumed that the position q(t) is available for measurement.

Control Objective: The desired reference trajectory for the control object to follow is assumed to be available as bounded functions of time in terms of generalized position q_ref (t). That is, there exist three positive constants r_i, i=0,1,2 such that the following inequalities hold:

$\Vert {\mathbf{q}}_{ref}(t)\Vert \le r_{\mathit{0}},\Vert {\stackrel{.}{\mathbf{q}}}_{ref}(t)\Vert \le r_{\mathit{1}}\hspace{0.17em}\text{and}\Vert {\stackrel{..}{\mathbf{q}}}_{ref}(t)\Vert \le r_{\mathit{2}}$

(2)

State space representation: The dynamic equation of n link robot manipulator (1) can be written in the state space representation as:

$\{\begin{array}{l}{\stackrel{.}{\mathbf{x}}}_{\mathit{1}}(t)={\mathbf{x}}_{\mathit{2}},\hfill \\ {\stackrel{.}{\mathbf{x}}}_{\mathit{2}}(t)=\mathbf{f}({\mathbf{x}}_{\mathit{1}},{\mathbf{x}}_{\mathit{2}})+\mathbf{P}({\mathbf{x}}_{\mathit{1}})\mathbf{u}(t),\hfill \\ \mathbf{y}(t)={\mathbf{x}}_{\mathit{1}},\hfill \end{array}$

(3)

where x(t) = [x ₁ x ₂ ] ^T = [q q̇]^T ε ℜ²ⁿ is the state vector. u(t) ε ℜⁿ represents the control torque vector and y(t) ε ℜⁿ is the output vector (position angular). f(x ₁ , x₂) = –M(q)⁻¹ (C(q, q̇) q̇ + G(q))) ε ℜⁿ and P(x₁) = M(x₁)⁻¹ ε ℜ^nxn.

Properties (Spong M. W. & Vidyasagar M. 1989):

P₁. The matrix M(x₁) is symmetric definite positive, then there exist two positive constants: M and M̄ such that: M ≤ ||M(x₁)|| ≤ M̄.

P₂. ∃ μ > 0 such that ||C(q, x)||≤ μ ||x||, ∀ x ε ℜ ⁿ .

P₃ The vector function f(x₁,x₂) is Lipschitz with respect to x₂. Thus there exists κ > 0 such that:

$\begin{array}{l}\Vert \mathbf{f}({\mathbf{x}}_{\mathit{1}},{\mathbf{x}}_{\mathit{2}})-\mathbf{f}({\mathbf{x}}_{\mathit{1}},{\stackrel{.}{\mathbf{q}}}_{ref})\Vert \le \kappa \Vert {\mathbf{x}}_{\mathit{2}}-{\stackrel{.}{\mathbf{q}}}_{ref}\Vert =\kappa \Vert {\mathbf{e}}_2\Vert ,\hfill \\ \forall (x_{\mathit{1}},x_{\mathit{2}})\in {\Re }^n\times {\Re }^n.\hfill \end{array}$

3. Receding-horizon control law

In the receding-horizon control strategy, the following control problem is solved at each t>0 and x(t):

$\begin{array}{l}\underset{u(t)}{Min}J(x(t),t,\mathbf{u}(t))=\hfill \\ \underset{u(t)}{Min}\frac{1}{2}\underset{t}{\overset{t+h}{\int }}\left[\mathbf{x}{(\tau )}^T\mathbf{Q}\mathbf{x}(\tau )+\mathbf{u}{(\tau )}^T\mathbf{R}\hspace{0.17em}\mathbf{u}(\tau )\right]d\tau \hfill \end{array}$

(4)

subject to the equation (3) and x(t + h) = 0 for some h>0, where Q is positive definite and R positive semi-definite. Denote the optimal control to the above problem by u*(τ), τ ε [t, t + h]. The currently applied control is u(t) set equal to u*(t). This process is repeated for every next t for stabilization of the system at the origin. However, to solve a nonlinear dynamic optimization problem with equality constraints is highly computationally intensive, and in many cases it is impossible to be performed within a reasonable time limit. Furthermore, the global optimization solution cannot be guaranteed in each optimization procedure since, in general, it is a non-convex, constrained nonlinear optimization problem.

In order to find the current control that improves tracking error along a fixed interval, the output tracking error e(τ) = q(τ) – q_ref(τ) is used instead of the state vector x(τ) in the above receding control problem:

$J(\mathbf{e},\mathbf{u},t)=\frac{1}{2}{\int }_t^{t+T}\left({\mathbf{e}}^T(\tau )\mathbf{Q}\mathbf{e}(\tau )+\mathbf{u}{(\tau )}^T\mathbf{R}\mathbf{u}(\tau )\right)d\tau$

(5)

where Qεℜ^nxn is positive definite, Rεℜ^nxn positive semi-definite, T is the predicted tracking horizon.

To avoid the computational burden, we shall approximate the above receding- horizon control problem by Simpson's rule (Atkinson K. E. 1978):

$\begin{array}{l}J=\frac{1}{2}\underset{t}{\overset{t+T}{\int }}L(\tau )d\tau \hfill \\ =\frac{T}{6}\left[L(t)+4L\left(t+\frac{T}{2}\right)+L(t+T)\right]\hfill \\ =\frac{h}{3}\left[L(t)+4L(t+h)+L(t+2h)\right]\hfill \end{array}$

with T=2h is the prediction horizon and

$L(\tau )={\mathbf{e}}^T(\tau )\mathbf{Q}\hspace{0.17em}\mathbf{e}(\tau )+{\mathbf{u}}^T(\tau )\mathbf{R}\hspace{0.17em}\mathbf{u}(\tau ).$

A simple and effective way of predicting the cost function L(.) is to expand the predicted tracking error in a first order Taylor series, in the following way: q(t + h) = q(t) + h q̇(t) and the reference trajectory is predicted as follow: q _ref (t + h) = q _ref (t) + h q̇ _ref .

The predicted tracking error is then given by:

$\mathbf{e}(t+h)=\mathbf{e}(t)+h\stackrel{.}{\mathbf{e}}(t).$

Predict e (t+2h) by another first-order Taylor series expansion at e (t+h) to have:

$\begin{array}{l}\mathbf{e}(t+2h)=\mathbf{q}(t+2h)-{\mathbf{q}}_{ref}(t+2h)\hfill \\ =\mathbf{e}(t)+2h\stackrel{.}{\mathbf{e}}(t)+h^2(\mathbf{f}-{\stackrel{..}{\mathbf{q}}}_{ref})+h^2\mathbf{P}\mathbf{u}(t)\hfill \end{array}$

(6)

where $\mathbf{q}(t+2h)=\mathbf{q}(t)+2h\stackrel{.}{\mathbf{q}}(t)+h^2\mathbf{f}+h^2\mathbf{P}\hspace{0.17em}\mathbf{u}(t)$ and ${\mathbf{q}}_{ref}(t+2h)={\mathbf{q}}_{ref}(t)+2h{\stackrel{.}{\mathbf{q}}}_{ref}(t)+h^2{\stackrel{..}{\mathbf{q}}}_{ref}(t).$ .

Thus, the performance index (5) can be approximated as:

$\begin{array}{l}J=\frac{h}{3}[{\mathbf{e}}^T(t)\mathbf{Q}\hspace{0.17em}\mathbf{e}(t)+{\mathbf{u}}^T(t)\mathbf{R}\hspace{0.17em}\mathbf{u}(t)\hfill \\ +4{\mathbf{e}}^T(t+h)\mathbf{Q}\hspace{0.17em}\mathbf{e}(t+h)+4{\mathbf{u}}^T(t+h)\mathbf{R}\hspace{0.17em}\mathbf{u}(t+h)\hfill \\ +{\mathbf{e}}^T(t+2h)\mathbf{Q}\hspace{0.17em}\mathbf{e}(t+2h)+{\mathbf{u}}^T(t+2h)\mathbf{R}\hspace{0.17em}\mathbf{u}(t+2h)]\hfill \end{array}$

(7)

We can rewrite the performance index (7) in the conventional quadratic form by using the predicted tracking error given above, as:

$\bar{J}=\frac{3}{2h}J=\frac{1}{2}{\mathbf{U}}^T\mathbf{\theta }(x)\mathbf{U}+{\mathbf{G}}^T(\mathbf{x})\mathbf{U}+m(\mathbf{e},\stackrel{.}{\mathbf{e}})$

where $\mathbf{\theta }(x)=\left|\begin{array}{lll}\hfill \mathbf{R}+h^4{\mathbf{P}}^T({\mathbf{x}}_1)\mathbf{Q}\hspace{0.17em}\mathbf{P}({\mathbf{x}}_1)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 4\mathbf{R}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \mathbf{R}\hfill \end{array}\right|$ is positive definite matrix, m(e, ė): terms that are independent of U(t) where U(t)^T = [u(t)^T u(t+h)^T u(t+2h)^T],

${\mathbf{G}}^T(\mathbf{x})=\left|h^2{\left(\mathbf{e}+2h\stackrel{.}{\mathbf{e}}+h^2(\mathbf{f}-{\stackrel{..}{\mathbf{q}}}_{ref}\right)}^T\mathbf{Q}\hspace{0.17em}\mathbf{P}(\mathbf{x})\hspace{0.17em}0\hspace{0.17em}0\right|.$

The receding-horizon control problem that minimizes the cost function J̄ is: U(t) = –θ(x)⁻¹ G(x). The applied control signal to nonlinear system at time t is given by:

$\begin{array}{l}\mathbf{u}(t)=-h^2\mathbf{M}({\mathbf{x}}_1){\left(h^4\mathbf{Q}+\mathbf{M}({\mathbf{x}}_1)\mathbf{R}\hspace{0.17em}\mathbf{M}({\mathbf{x}}_1)\right)}^{-1}\times \hfill \\ \mathbf{Q}\left(\mathbf{e}+2h\stackrel{.}{\mathbf{e}}+h^2\left(\mathbf{f}({\mathbf{x}}_1,{\mathbf{x}}_2)-{\stackrel{..}{\mathbf{q}}}_{ref}\right)\right)\hfill \end{array}$

(9)

Note that with R=0, the above nonlinear predictive control law leads to the well known computed torque controller.

4- Stability analysis and robustness issues

In this section, we will investigate the stability and the robustness of the closed loop system with respect to model uncertainties.

1- Stability Analysis

Let Q = qI _n and R = rI _n (we give the same penalty to all joints), the tracking error of the nonlinear system (3) closed by the nonlinear feedback (9) is given by:

$\{\begin{array}{l}{\stackrel{.}{\mathbf{e}}}_1={\mathbf{e}}_2\hfill \\ {\stackrel{.}{\mathbf{e}}}_2=-q{h}^2{\bar{\mathbf{P}}}^{-1}{\mathbf{e}}_1-2q{h}^3{\bar{\mathbf{P}}}^{-1}{\mathbf{e}}_2+r{\bar{\mathbf{P}}}^{-1}\mathbf{M}{({\mathbf{x}}_1)}^2\times \hfill \\ \left(\mathbf{f}({\mathbf{x}}_1,{\mathbf{x}}_2)-{\stackrel{..}{\mathbf{q}}}_{ref}\right)\hfill \end{array}$

(10)

Where ${\mathbf{e}}_1=\mathbf{q}(t)-{\mathbf{q}}_{ref}(t);\hspace{0.17em}\hspace{0.17em}{\mathbf{e}}_2=\stackrel{.}{\mathbf{q}}(t)-{\stackrel{.}{\mathbf{q}}}_{ref}(t)$ and $\bar{\mathbf{P}}=q{h}^4{\mathbf{I}}_n+r\mathbf{M}{({\mathbf{x}}_1)}^2.$ .

This equation can be written in compact form as:

$\stackrel{.}{\mathbf{e}}=\mathbf{A}(h,{\mathbf{x}}_1)\mathbf{e}+\mathbf{B}\mathbf{\eta }$

(11)

Where $\mathbf{A}(h,{\mathbf{x}}_1)=\left|\begin{array}{ll}0\hfill & I_n\hfill \\ -q{h}^2{\bar{\mathbf{P}}}^{-1}\hfill & -2q{h}^3{\bar{\mathbf{P}}}^{-1}\hfill \end{array}\right|$ ; $\mathbf{B}=\left|\begin{array}{l}0\hfill \\ {\mathbf{I}}_n\hfill \end{array}\right|$ ; $\mathbf{\eta }=r{\bar{\mathbf{P}}}^{-1}\mathbf{M}{({\mathbf{x}}_1)}^2\left(\mathbf{f}({\mathbf{x}}_1,{\mathbf{x}}_2)-{\stackrel{..}{\mathbf{q}}}_{ref}\right).$ .

Lemma-1: The matrix A(h,x₁) is Hurwitz.

Proof: Both matrix P̄ and its inverse are symmetric positives definite. Let x̄ ε ℜ ⁿ and $\bar{\lambda }\in \Re$ are the eigenvector and the correspondent eigenvalue of the inverse of the above matrix. Thus, we have the equality:

$\mathbf{A}(h,{\mathbf{x}}_1)\left|\begin{array}{l}\hfill \bar{\mathbf{x}}\hfill \\ \hfill \lambda \bar{\mathbf{x}}\hfill \end{array}\right|=\left|\begin{array}{l}\hfill \lambda \bar{\mathbf{x}}\hfill \\ \hfill -q{h}^2\bar{\lambda }\bar{\mathbf{x}}-2q{h}^3\lambda \bar{\lambda }\bar{\mathbf{x}}\hfill \end{array}\right|=\left|\begin{array}{l}\hfill \lambda \bar{\mathbf{x}}\hfill \\ \hfill {\lambda }^2\bar{\mathbf{x}}\hfill \end{array}\right|=\lambda \left|\begin{array}{l}\hfill \bar{\mathbf{x}}\hfill \\ \hfill \lambda \bar{\mathbf{x}}\hfill \end{array}\right|\mathrm{N}$

ote that λ is the solution of equation:

${\lambda }^2+q{h}^2\bar{\lambda }+2q{h}^3\bar{\lambda }\lambda =0$

(12)

Therefore, λ is the eigenvalue of the matrix A(h, x₁) and $\left|\begin{array}{l}\hfill \bar{\mathbf{x}}\hfill \\ \hfill \lambda \bar{\mathbf{x}}\hfill \end{array}\right|$ the correspondent eigenvector. Set λ₁ and λ₂ the solution of the equation (12), we have the relations:

$\begin{array}{l}{\lambda }_1+{\lambda }_2=-2q{h}^3\bar{\lambda }\hfill \\ {\lambda }_1{\lambda }_2=q{h}^2\bar{\lambda }\hfill \end{array}$

Since the eigenvalue $\bar{\lambda }$ is positive, then λ₁ and λ₂ have a negative real part (end of the proof).

Since the matrix A(h,x₁) is a Hurwitz matrix, then for any symmetric positive definite matrix Q_A(h, x₁), there exists a symmetric positive definite matrix P_A(h, x₁) solution of the lyapunov equation:

$\begin{array}{l}{\stackrel{.}{\mathbf{P}}}_{\mathbf{A}}(h,{\mathbf{x}}_1)+\mathbf{A}{(h,{\mathbf{x}}_1)}^T{\mathbf{P}}_{\mathbf{A}}(h,{\mathbf{x}}_1)+{\mathbf{P}}_{\mathbf{A}}(h,{\mathbf{x}}_1)\mathbf{A}(h,{\mathbf{x}}_1)\hfill \\ =-{\mathbf{Q}}_{\mathbf{A}}(h,{\mathbf{z}}_1)\hfill \end{array}$

(13)

From the property P₃, the function f(x ₁ , x₂) is Lipschitz with regards to x₂, we can always find a bounded continuous function σ(e₂, t) and positive scalar μ satisfying the inequality:

$\Vert \mathbf{f}({\mathbf{x}}_1,{\mathbf{x}}_2)-{\stackrel{..}{\mathbf{q}}}_{ref}\Vert \le \sigma ({\mathbf{e}}_2,t)\le u\Vert \mathbf{e}\Vert .$

(14)

Now, we can state the following theorem.

Theorem 1: The equilibrium point of the nonlinear system (3) in closed loop with the feedback control (9) is asymptotically stable if the following inequality hold:

$\mu \<\frac{{\lambda }_{min}({\mathbf{Q}}_{\mathbf{A}}){\lambda }_{max}(\bar{\mathbf{P}})}{2r{\bar{M}}^2{\lambda }_{max}({\mathbf{P}}_{\mathbf{A}})}.$

Moreover, if r = 0, then the origin is asymptotically stable equilibrium point.

Proof:

Let V(e) = e ^T P_A e be a Lyapunov function candidate, the time derivative of this function along the trajectories (11) is:

$\begin{array}{l}\stackrel{.}{V}={\mathbf{e}}^T\mathbf{A}{(h,{\mathbf{x}}_1)}^T{\mathbf{P}}_A\mathbf{e}+{\mathbf{e}}^T{\mathbf{P}}_A\mathbf{A}(h,{\mathbf{x}}_1)\mathbf{e}+{\mathbf{e}}^T{\stackrel{.}{\mathbf{P}}}_A\mathbf{e}\hfill \\ +{\eta }^T{\mathbf{B}}^T{\mathbf{P}}_A\mathbf{e}+{\mathbf{e}}^T{\mathbf{P}}_A\mathbf{B}\mathbf{\eta }\hfill \end{array}$

By using the equality (13), we obtain:

$\stackrel{.}{V}=-{\mathbf{e}}^T{\mathbf{Q}}_A\mathbf{e}+{\mathbf{\eta }}^T{\mathbf{B}}^T{\mathbf{P}}_A\mathbf{e}+{\mathbf{e}}^T{\mathbf{P}}_A\mathbf{B}\mathbf{\eta }$

or $\stackrel{.}{V}\le -{\lambda }_{min}({\mathbf{Q}}_A){\Vert \mathbf{e}\Vert }^2+2\Vert \mathbf{B}\Vert \Vert {\mathbf{P}}_A\Vert \Vert \mathbf{e}\Vert$

From the inequality (14), we can write:

$\stackrel{.}{V}\le {\lambda }_{min}({\mathbf{Q}}_A){\Vert \mathbf{e}\Vert }^2+2\mu r{\bar{M}}^2{\lambda }_{max}({\bar{\mathbf{P}}}^{-1}){\lambda }_{max}({\mathbf{P}}_A){\Vert \mathbf{e}\Vert }^2$

Thus, V̇ ≤ –κ||e|| ² is negative definite if κ > 0, where $\kappa ={\lambda }_{min}({\mathbf{Q}}_A)-2\mu {\bar{M}}^{2}r\frac{{\lambda }_{max}({\mathbf{P}}_A)}{{\lambda }_{max}(\bar{\mathbf{P}})}$ . This ensures the asymptotic stability of the equilibrium point.

Note that a short steady state error will be observed in tracking position error when r ≠ 0. However, if r = 0 the time derivative of the Lyapunov function becomes: V̇ = –λ_min(Q_A)||e|| ² which is negative definite. Thus, we can conclude that the origin becomes the equilibrium point of the system (12) and is asymptotically stable, i.e:

$\underset{t\to \infty }{Lim}\mathbf{e}(t)=\underset{t\to \infty }{Lim}{\left|q-q_{ref}\stackrel{.}{\mathbf{q}}-{\stackrel{.}{\mathbf{q}}}_{ref}\right|}^T={\left|\begin{array}{cc}\mathit{0}& \mathit{0}\end{array}\right|}^T.$

2- Robustness

In order to incorporate modeling uncertainties into the model of the rigid robot (1), the matrices M(q), C(q, q̇) and the vector G(q) are split up into a nominal part (indicated by the subscript zero) and an uncertain part as:

$\begin{array}{l}({\mathbf{M}}_0(\mathbf{q})+\Delta \mathbf{M})\stackrel{..}{\mathbf{q}}+({\mathbf{C}}_0(\mathbf{q},\stackrel{.}{\mathbf{q}})+\Delta \mathbf{C})\stackrel{.}{\mathbf{q}}+{\mathbf{G}}_0(\mathbf{q})+\Delta \mathbf{G}+\mathbf{f}\hfill \\ =\mathbf{u}(t)\hfill \end{array}$

(15)

The friction torque f is included in the uncertain part given the difficulty to model it correctly. Obviously, only the nominal part of the model can be used by the nonlinear predictive control, given by:

$\begin{array}{l}\mathbf{u}(t)=-\frac{1}{h^2}{\mathbf{M}}_0(q)(\mathbf{e}+2h\stackrel{.}{\mathbf{e}})\hfill \\ +({\mathbf{C}}_0(\mathbf{q},\stackrel{.}{\mathbf{q}})\stackrel{.}{\mathbf{q}}+{\mathbf{G}}_0(\mathbf{q}))+{\mathbf{M}}_0(\mathbf{q}){\stackrel{..}{\mathbf{q}}}_{ref}\hfill \end{array}$

(16)

Where R is setting to zero. With the nonlinear control law (16), the closed loop system is:

$\begin{array}{l}\stackrel{..}{\mathbf{e}}+\frac{2}{h}\mathbf{M}{(\mathbf{q})}^{-1}{\mathbf{M}}_0(\mathbf{q})\stackrel{.}{\mathbf{e}}+\frac{1}{h^2}\mathbf{M}{(\mathbf{q})}^{-1}{\mathbf{M}}_0(\mathbf{q})\mathbf{e}\hfill \\ =-{\mathbf{M}}^{-1}(\mathbf{q})\left(\Delta \mathbf{M}{\stackrel{..}{\mathbf{q}}}_{ref}+\Delta \mathbf{C}\stackrel{.}{\mathbf{q}}+\Delta \mathbf{G}+\mathbf{f}\right)=\mathbf{\upsilon }(\mathbf{q},\stackrel{.}{\mathbf{q}},{\stackrel{..}{\mathbf{q}}}_{ref},t)\hfill \end{array}$

(17)

To estimate the worst case bound of the function υ, we make the following assumptions for all q ε ℜ ⁿ :

$\begin{array}{l}{\underset{\_}{M}}_0\le \Vert {\mathbf{M}}_0(\mathbf{q})\Vert \le {\bar{M}}_0;\Vert \Delta \mathbf{M}\Vert \le \bar{m}\le {\lambda }_{min}(M(\mathbf{q}));\hfill \\ \Vert \Delta \mathbf{C}\Vert \le \bar{c};\Vert \Delta \mathbf{G}\Vert \le \bar{g};\Vert f\Vert \le \bar{f};\hfill \end{array}$

Given these assumptions with the inequalities (2), we can find a bounded continuous vector function ρ(e, ė, t) satisfying the inequality (Spong M. W. & Vidyasagar M.; 1989):

$\Vert \upsilon \Vert \<\rho (\mathbf{e},\stackrel{.}{\mathbf{e}},t)\le \gamma \Vert \mathbf{e}\Vert \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{for all}\hspace{0.17em}\mathbf{q}\in {\Re }^n,$

where γ is a positive scalar and e = |q q̇| ^T = |e ₁ e ₂ | ^T .

In the state space representation, the system (17) can be transformed to:

$\stackrel{.}{\mathbf{e}}=\bar{\mathbf{B}}(h,{\mathbf{x}}_1)\mathbf{e}+\mathbf{B}\mathbf{\upsilon }$

(18)

where $\bar{\mathbf{B}}(h,{\mathbf{x}}_1)=\left|\begin{array}{ll}\hfill 0\hfill & \hfill {\mathbf{I}}_n\hfill \\ \hfill -\frac{\mathbf{b}({\mathbf{x}}_1)}{h^2}\hfill & \hfill -\frac{2\mathbf{b}({\mathbf{x}}_1)}{h}\hfill \end{array}\right|.$ .

Since both M⁻¹(x₁) and M₀(x₁) are symmetric definite positives, the matrix b(x₁) = M⁻¹ (x₁) M₀ (x₁) has all its eigenvalues reels and positives (Samson C.; 1983). Thus, from the Lemma-1, we can conclude that the matrix B̄(h, x₁) is Hurwitz, then for any symmetric positive definite matrix Q_B(h,x₁) there exists a positive definite matrix P_B(h,x₁) solution of the given Lyapunov equation:

$\begin{array}{l}{\stackrel{.}{\mathbf{P}}}_{\mathbf{B}}(h,{\mathbf{x}}_1)+\bar{\mathbf{B}}{(h,{\mathbf{x}}_1)}^T{\mathbf{P}}_{\mathbf{B}}(h,{\mathbf{x}}_1)+{\mathbf{P}}_{\mathbf{B}}(h,{\mathbf{x}}_1)\bar{\mathbf{B}}(h,{\mathbf{x}}_1)\hfill \\ =-{\mathbf{Q}}_{\mathbf{B}}(h,{\mathbf{x}}_1)\hfill \end{array}$

Theorem 2: Suppose that the inequality holds $\gamma \<\frac{{\lambda }_{min}({\mathbf{Q}}_B)}{\mathit{2}{\lambda }_{max}({\mathbf{P}}_B)}$ , then the equilibrium point of the nonlinear system with uncertainties (15) closed by the optimal control (16) is asymptotically stable.

Proof:

Let V = e ^T P _B e a Lyapunov function candidate. The differentiation of V along the trajectories (18) leads to:

$\begin{array}{l}\stackrel{.}{V}=-{\mathbf{e}}^T{\mathbf{Q}}_{\mathbf{B}}\mathbf{e}+{\mathbf{\upsilon }}^T{\mathbf{B}}^T{\mathbf{P}}_B\mathbf{e}+{\mathbf{e}}^T{\mathbf{P}}_B\mathbf{B}\mathbf{\upsilon }.\hfill \\ \stackrel{.}{V}\le -{\lambda }_{min}({\mathbf{Q}}_B){\Vert \mathbf{e}\Vert }^2+2\gamma {\lambda }_{max}({\mathbf{P}}_B){\Vert \mathbf{e}\Vert }^2\hfill \\ \le -({\lambda }_{min}({\mathbf{Q}}_B)-2\gamma {\lambda }_{max}({\mathbf{P}}_B)){\Vert \mathbf{e}\Vert }^2.\hfill \end{array}$

Which is definite negative if $\gamma \<\frac{{\lambda }_{min}({\mathbf{Q}}_{\mathbf{B}})}{\mathit{2}{\lambda }_{max}({\mathbf{P}}_B)}.$ . Therefore, By LaSalle's invariance theorem, the solution e(t) of (18) tends to the invariance set:

$S=\left\{\mathbf{e}/{\mathbf{e}}_{\mathit{1}}=h^{\mathit{2}}\mathbf{b}{({\bar{\mathbf{x}}}_1)}^{\mathit{-}1}\mathbf{\upsilon }(\bar{\mathbf{x}},\stackrel{..}{\mathbf{q}}ref),{\mathbf{e}}_2=\mathit{0}\right\}.$

We conclude that bounded uncertainties will introduce a steady state error on tracking position angular. Where x̄ is the equilibrium point of the system (3) in closed loop with the control law (16).

3- Integral action

It is known in the literature that the integral action increases the robustness of the closed loop system against the low frequency disturbances as long as the closed loop system is stable. In this part, we shall incorporate an integral action in the loop to eliminate the steady state error and enhance the robustness of the proposed control scheme with respect to model uncertainties and disturbances.

Thus, the cost function to minimize becomes:

$\begin{array}{l}J({\mathbf{e}}_0,\mathbf{u},t)=\frac{1}{2}{\int }_t^{t+T}\left({\mathbf{e}}_t^T(\tau )\mathbf{Q}{\mathbf{e}}_0(\tau )+\mathbf{u}{(\tau )}^T\mathbf{R}\mathbf{u}(\tau )\right)d\tau \hfill \\ =\frac{1}{2}\underset{t}{\overset{t+h}{\int }}L(\tau )d\tau \hfill \end{array}$

(19)

Where ė₀= e₁ = x₁ – q _ref = q – q _ref and L(τ) = e₀ ^T (τ) Q e₀ (τ) + u ^T (τ) R u (τ).

Also in this case, we use the Simpson's rule to approximate the integral in the cost function (19) by:

$\begin{array}{l}J({\mathbf{e}}_0,\mathbf{u},t)=\frac{h}{3}(L(t)+4L(t+h)+L(t+2h))\hfill \\ With:\hspace{0.17em}L(t)={\mathbf{e}}_0^T(t)\mathbf{Q}{\mathbf{e}}_0(t)+{\mathbf{u}}^T(t)\mathbf{R}\hspace{0.17em}\mathbf{u}(t);\hfill \\ L(t+h)={\mathbf{e}}_0^T(t+h)\mathbf{Q}{\mathbf{e}}_0(t+h)+{\mathbf{u}}^T(t+h)\mathbf{R}\hspace{0.17em}\mathbf{u}(t+h)\hfill \\ L(t+2h)={\mathbf{e}}_0^T(t+2h)\mathbf{Q}{\mathbf{e}}_0(t+2h).\hfill \\ +\mathbf{u}(t+2h)\mathbf{R}\hspace{0.17em}\mathbf{u}(t+2h)\hfill \end{array}$

Note that in this case, the Taylor approximation of the predicted vector e₀(t + h) is given by:

$\begin{array}{l}{\mathbf{e}}_0(t+h)={\mathbf{e}}_0(t)+h{\mathbf{e}}_1+\frac{h^2}{2}{\mathbf{e}}_2+\frac{h^3}{6}\left(\mathbf{f}(x)-{\stackrel{..}{\mathbf{q}}}_{ref}\right)\hfill \\ +\frac{h^3}{6}\mathbf{P}({\mathbf{x}}_1)\mathbf{u}(t)\hfill \end{array}$

Following same steps in paragraph 2, the optimal control vector U(t) that minimizes the new cost function is:

$\mathbf{U}(t)=-{\bar{\mathbf{H}}}^{-1}\bar{\mathbf{G}}(\mathbf{x})$

(20)

Where U(t) = |u(t) u(t + h) u(t + 2h)| ^T ; $\bar{\mathbf{H}}=diag(\mathbf{R}+\frac{5}{9}{h}^6{\mathbf{P}}^T({\mathbf{x}}_1)\mathbf{Q}\mathbf{P}({\mathbf{x}}_1),4\mathbf{R},\mathbf{R})$ and

$\begin{array}{l}\bar{\mathbf{G}}(\mathbf{x})=\hfill \\ \hfill \left|\begin{array}{l}\hfill \frac{2}{3}{h}^3\mathbf{P}({\mathbf{x}}_1)\mathbf{Q}\left(2{\mathbf{e}}_0+3h{\mathbf{e}}_1+2h^2{\mathbf{e}}_2+\frac{5}{6}{h}^3\left(\mathbf{f}(\mathbf{x})-{\stackrel{..}{\mathbf{q}}}_{ref}\right)\right)\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right|\hfill \end{array}$

The control signal to be applied to the nonlinear system at time t is:

$\begin{array}{l}\mathbf{u}(t)=-\frac{2}{3}{h}^3{\mathbf{M}}_0({\mathbf{x}}_1){\overline{\stackrel{̄}{\mathbf{P}}}}^{-1}(h,{\mathbf{x}}_1)\mathbf{Q}\times \hfill \\ \left(2{\mathbf{e}}_0+3h{\mathbf{e}}_1+2h^2{\mathbf{e}}_2+\frac{5}{6}{h}^3(\mathbf{f}(\mathbf{x})-{\stackrel{..}{\mathbf{q}}}_{ref})\right)\hfill \end{array}$

(21)

Where $\overline{\stackrel{̄}{\mathbf{P}}}(h,{\mathbf{x}}_1)=\frac{5}{9}{h}^6\mathbf{Q}+{\mathbf{M}}_0({\mathbf{x}}_1)\mathbf{R}{\mathbf{M}}_0({\mathbf{x}}_1).$ .

Let R = 0 in equation (21), the control signal becomes:

$\begin{array}{l}\mathbf{u}(t)=-\frac{9}{5}{\mathbf{M}}_0({\mathbf{x}}_1)\times \hfill \\ \left(\frac{4}{3h^3}{\mathbf{e}}_0+\frac{2}{h^2}{\mathbf{e}}_1+\frac{4}{3h}{\mathbf{e}}_2+\frac{5}{9}\left(\mathbf{f}(\mathbf{x})-{\stackrel{..}{\mathbf{q}}}_{ref}\right)\right)\hfill \end{array}$

The dynamic of the tracking error is given by the equations:

$\{\begin{array}{l}{\stackrel{.}{\mathbf{e}}}_0={\mathbf{e}}_1\hfill \\ {\stackrel{.}{\mathbf{e}}}_1={\mathbf{e}}_2\hfill \\ {\stackrel{.}{\mathbf{e}}}_2=-\frac{12}{5h^3}\mathbf{b}({\mathbf{x}}_1){\mathbf{e}}_0-\frac{18}{5h^2}\mathbf{b}({\mathbf{x}}_1){\mathbf{e}}_1-\frac{12}{5h}\mathbf{b}({\mathbf{x}}_1){\mathbf{e}}_2\hfill \\ +\mathbf{\upsilon }(\mathbf{x},{\stackrel{..}{\mathbf{q}}}_{ref})\hfill \end{array}$

(22)

Or in compact form:

$\stackrel{.}{\mathbf{e}}=\tilde{\mathbf{B}}(h,{\mathbf{x}}_1)\mathbf{e}+\mathbf{B}\hspace{0.17em}\mathbf{\upsilon }(\mathbf{x},{\stackrel{..}{\mathbf{q}}}_{ref})$

(23)

Where:

$\tilde{\mathbf{B}}(h,x_1)=\left|\begin{array}{lll}\hfill 0\hfill & \hfill I_n\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill I_n\hfill \\ -\frac{12}{5h^3}\mathbf{b}({\mathbf{x}}_1)\hfill & -\frac{18}{5h^2}\mathbf{b}({\mathbf{x}}_1)\hfill & -\frac{12}{5h}\mathbf{b}({\mathbf{x}}_1)\hfill \end{array}\right|$

Lemma 2: Suppose that λ_max(ε(x₁)) < 2.6, with $\varepsilon ({\mathbf{x}}_1)={\mathbf{M}}_0^{-1}\Delta \mathbf{M}$ , then the matrix $\overline{\stackrel{̄}{\mathbf{B}}}(h,{\mathbf{x}}_1)$ is Hurwitz.

Proof: Let x̄ ε ℜ ⁿ and $\bar{\lambda }\in {\Re }^{+}$ represent eigenvector and eigenvalue of the matrix b(x₁) respectively. Set $\mathbf{v}={\left|\begin{array}{lll}\bar{\mathbf{x}}\hfill & \tilde{\lambda }\bar{\mathbf{x}}\hfill & \tilde{\lambda }\bar{\mathbf{x}}\hfill \end{array}\right|}^T\in {\Re }^{3n}$ be a vector and $\tilde{\lambda }\in \Re$ a scalar. We have:

$\begin{array}{l}\hfill \tilde{\mathbf{B}}(h,{\mathbf{x}}_1)\mathbf{v}=\left|\begin{array}{l}\hfill \tilde{\lambda }\bar{\mathbf{x}}\hfill \\ \hfill {\tilde{\lambda }}^2\bar{\mathbf{x}}\hfill \\ -\frac{12}{5h^3}\bar{\lambda }\bar{\mathbf{x}}-\frac{18}{5h^2}\tilde{\lambda }\bar{\lambda }\bar{\mathbf{x}}-\frac{12}{5h}{\tilde{\lambda }}^2\bar{\lambda }\bar{\mathbf{x}}\hfill \end{array}\right|.\hfill \\ =\left|\begin{array}{l}\tilde{\lambda }\bar{\mathbf{x}}\hfill \\ {\tilde{\lambda }}^2\bar{\mathbf{x}}\hfill \\ {\tilde{\lambda }}^3\bar{\mathbf{x}}\hfill \end{array}\right|=\tilde{\lambda }\hspace{0.17em}\mathbf{v}\hfill \end{array}$

Then, $\tilde{\lambda }$ and v are the eigenvalue and eigenvector of the matrix B˜(h, x₁) respectively, where the eigenvalue $\tilde{\lambda }$ verify the equality:

${\tilde{\lambda }}^3+\frac{12}{5h}\bar{\lambda }{\tilde{\lambda }}^2+\frac{18}{5h^2}\bar{\lambda }\tilde{\lambda }+\frac{12}{5h^3}\bar{\lambda }=0.$

By using the Rooth-Hurwitz criterion, the solutions of the above equation lie in left half plane (stable domain) if: $\bar{\lambda }>\frac{5}{18}$ .

Since M^–1M ₀ = (I _n + ε(x ₁ ))^–1, the condition for stability becomes: ${\lambda }_{max}({\mathbf{I}}_n+\varepsilon ({\mathbf{x}}_1))\&lt;\frac{18}{5}$ or λ_max (ε(x₁)) < 2.6 which can be easily verified if the uncertainty ΔM is small with regards to the nominal value of the matrix M₀(x₁).

We conclude that the matrix B˜(h, x₁) is Hurwitz, then for any symmetric positive definite matrix Q˜(h, x₁) there exists a positive definite matrix P˜(h, x₁) solution of the given Lyapunov equation:

$\begin{array}{l}\stackrel{.}{\stackrel{~}{\mathbf{P}}}(h,{\mathbf{x}}_1)+\tilde{\mathbf{B}}{(h,{\mathbf{x}}_1)}^T\tilde{\mathbf{P}}(h,{\mathbf{x}}_1)+\tilde{\mathbf{P}}(h,{\mathbf{x}}_1)\tilde{\mathbf{B}}(h,{\mathbf{x}}_1)\hfill \\ =-\tilde{\mathbf{Q}}(h,{\mathbf{x}}_1)\hfill \end{array}$

Theorem 3: Suppose the inequality $\gamma \&lt;\frac{{\lambda }_{min}(\tilde{\mathbf{Q}})}{2{\lambda }_{max}(\tilde{\mathbf{P}})}$ hold, then the equilibrium point of the nonlinear system with uncertainties (23) is asymptotically stable.

The proof may be obtained in the same way as the proof of the theorem 2 and is therefore omitted. Therefore, the equilibrium point of (22) or (23) is asymptotically stable. The tracking error tends towards the set:

$S=\left\{\mathbf{e}/{\mathbf{e}}_0=\frac{5h^3}{12}\mathbf{b}{({\bar{\mathbf{x}}}_1)}^{-1}\mathbf{\upsilon }(\bar{\mathbf{x}},{\stackrel{..}{q}}_{ref}),{\mathbf{e}}_1=0,{\mathbf{e}}_2=0\right\}\mathrm{W}$

e conclude that the position and velocity tracking error converge to zero, therefore the integral action eliminates the position steady state error, i.e:

$\underset{t\to \infty }{Lim}{\mathbf{e}}_1=0\hspace{0.17em}\text{and}\hspace{0.17em}\underset{t\to \infty }{Lim}{\mathbf{e}}_2=0.$

The price to be paid by introducing an integral action in the loop is that the control signal will not vanish and this will increase the required energy to maintain the tracking performance as in the matched case.

5. Nonlinear observer

A drawback of the previous nonlinear predictive controller is that it requires at least the measurement of velocity on the link side. However, as pointed out in (Nicosia S. & Tomei P. 1990) and (Canudas W. C. et al 1992), in the practical robotic systems all the generalized coordinates can be precisely measured by the encoder for each joint, but the velocity measurements obtained through the tachometers are easily perturbed by noises. Therefore, in order to coincide with these physical constraints, a nonlinear observer proposed in (Bornard G. et al 1993) is used in this paper.

Define the state vector as:

$\mathbf{z}(t)=\mathbf{T}\hspace{0.17em}\mathbf{x}(t)=[\dots \dots q_i(t){\stackrel{.}{q}}_i(t)\dots \dots ]\in {\Re }^{2n}$

where q_i(t) and q̇_i(t) are the link position and the velocity of the i^th arm respectively. Tεℜ^2nx2n is the transformation matrix. The system (3) can be transformed to:

$\{\begin{array}{l}\stackrel{.}{\mathbf{z}}=\mathbf{A}\mathbf{z}+\mathbf{H}\mathbf{f}(\mathbf{z})+\mathbf{H}\mathbf{P}(\mathbf{q})\mathbf{u}(t)\hfill \\ \mathbf{y}=\mathbf{C}\mathbf{z}\hfill \end{array}$

(24)

where

$\begin{array}{l}\mathbf{A}=\text{diag}({\mathbf{A}}_{\mathrm{i}}),\hspace{0.17em}{\mathbf{A}}_i=\left|\begin{array}{ll}0\hfill & 1\hfill \\ 0\hfill & 0\hfill \end{array}\right|,\hspace{0.17em}\mathbf{C}=\text{diag}({\mathbf{C}}_{\mathrm{i}}),\hspace{0.17em}{\mathbf{C}}_i=\left|\begin{array}{ll}1\hfill & 0\hfill \end{array}\right|,\hfill \\ \mathbf{H}=\text{diag}({\mathbf{H}}_{\mathrm{i}}),\hspace{0.17em}{\mathbf{H}}_i={\left|\begin{array}{ll}0\hfill & 1\hfill \end{array}\right|}^T\hspace{0.17em}\text{for}\hspace{0.17em}\mathrm{i}=1,\dots \dots \dots ,\mathrm{n}.\hfill \end{array}$

With the assumption that the control torque u(t) is uniformly bounded, the high gain observer described in (Bornard G. et al, 1993) can be used to estimate angular positions and angular velocities of the n link rigid robot manipulator (24). The dynamic nonlinear observer is given by:

$\{\begin{array}{l}\stackrel{.}{\widehat{\mathbf{z}}}=\mathbf{A}\widehat{z}+\mathbf{H}\mathbf{f}(\widehat{\mathbf{z}})+\mathbf{H}\mathbf{P}(\mathbf{q})\mathbf{u}(t)+\mathbf{K}(\mathbf{y}-\widehat{\mathbf{y}}),\hfill \\ \widehat{\mathbf{y}}=\mathbf{C}\widehat{\mathbf{z}},\hfill \end{array}$

(25)

where K = Γ⁻¹(α) V is the gain of the observer with T (α) = diag(Γ _i (α)), ${\Gamma }_i(\alpha )=\left|\begin{array}{ll}\alpha \hfill & 0\hfill \\ 0\hfill & {\alpha }^2\hfill \end{array}\right|$ for any α >0 and due to the observability property of (A,C) the eigenvalues of (A-VC) can be assigned by V.

6. A Simulation example

To illustrate some of the conclusions of this paper, we have simulated the approximate receding-horizon control scheme on a two-link robot arm used in (Lee K. W. & Khalil H. K. 1997), (Spong M. W. 1992) with mechanic equations.

The arm is shown in Figure 1. The dynamic model is described in equation (1) with the following components:

$\begin{array}{l}{m}_{11}(\mathbf{q})=m_1{l}_{c1}^2+m_1{l}_{c2}^2+m_2{l}_{c1}^2;\hfill \\ +2m_2{l}_1{l}_{c2}cos(q_2)+I_1+I_2\hfill \\ m_{12}(\mathbf{q})=m_{21}(\mathbf{q})=I_2+m_2{l}_{c2}^2+m_2{l}_1{l}_{2}cos(q_2);\hfill \\ m_{22}(\mathbf{q})=I_2+m_2{l}_{c2}^2;C_{11}(\mathbf{q},\stackrel{.}{\mathbf{q}})=-{\stackrel{.}{q}}_2{m}_2{l}_1{l}_{c2}sin(q_2);\hfill \\ C_{12}(\mathbf{q},\stackrel{.}{\mathbf{q}})=-({\stackrel{.}{q}}_1+{\stackrel{.}{q}}_2)m_2{l}_1{l}_{c2}sin(q_2).\hfill \\ C_{21}(\mathbf{q},\stackrel{.}{\mathbf{q}})={\stackrel{.}{q}}_1{m}_2{l}_1{l}_{2}sin(q_2),\hspace{0.17em}{C}_{22}(\mathbf{q},\stackrel{.}{\mathbf{q}})=0.\hfill \\ G_1(\mathbf{q})=(m_1{l}_{c1}+m_2{l}_1)gcos(q_1)+m_2{l}_{c2}gcos(q_1+q_2)\hfill \\ G_2(\mathbf{q})=m_2{l}_{c2}gcos(q_1+q_2).\hfill \end{array}$

where 1_c1 is the mass center of gravity coordinate of the link 1, and 1_c2 is the mass center of gravity coordinate of the link 2. The values of the manipulator parameters are given in Table 1 (Lee K. W. & Khalil H. K. 1997) and (Spong M. W. 1992).

Fig.1

Two-linkage manipulators

Table 1

Physical parameters of two-link arm.

Link 1	m1=10kg   l1=1m   lc1=0.5m $I_1=\frac{10}{12}{\text{kgm}}^2$
Link 2	m2=5kg   l2 = 1m lc2=0.50m   $I_2=\frac{5}{12}{\text{kgm}}^2$

The reference models chosen in continuous time are:

$\begin{array}{lll}{\mathbf{q}}_{ref}=\left|\begin{array}{l}{q}_{ref1}\hfill \\ q_{ref2}\hfill \end{array}\right|\hfill & \text{with}\hfill & q_{ref1}(s)=\frac{{\omega }_1^2}{s^2+2\xi {\omega }_{1}s+{\omega }_1^2}\hfill \end{array}{r}_1(s)$

and $q_{ref2}(s)=\frac{{\omega }_2^2}{s^2+2\xi {\omega }_{2}s+{\omega }_2^2}{r}_2(s).$ .

The nonlinear predictive controller is used to track these desired trajectories with inputs (Lee K. W. & Khalil H. K. 1997):

$r_1(t)=r_2(t)=1.5(1-exp(-5t)(1+5t))rad.$

All simulations are carried with the nonlinear observer (25) with α=0.01 and the assigned eigenvalues σ(λ) = {–0.4, −0.8}. The initial displacements and velocities are chosen as:

$\begin{array}{l}{\mathrm{q}}_1(0)={\mathrm{q}}_2(0)=0^{\circ },\hspace{0.17em}{\stackrel{.}{q}}_1(0)={\stackrel{.}{q}}_2(0)=0,\hfill \\ {\widehat{q}}_1(0)={\widehat{q}}_2(0)=0.01\hspace{0.17em}\text{and}\hspace{0.17em}{\stackrel{.}{\widehat{q}}}_1(0)={\stackrel{.}{\widehat{q}}}_2(0)=0.\hfill \end{array}$

The parameter values of two-reference models are chosen as follows: ξ = 1, ω₁ = ω₂ = 10 rad/s.

The nonlinear controller (9), has been tested by simulation and the control parameters: Q=10⁷ I_n, R=10⁻¹⁴ I_n and h is set to 0.001. Simulation results are show in Figure 2. This Figure gives the angular position (q₁(t), q₂(t)) and the position tracking error. Although a very short steady state error is observed in the position tracking error and this was expected in the analysis part, a good tracking performance is achieved by this controller in matched case. Figure 3 illustrates the induced control torque applied to robot manipulator. Note that the control torque lie inside the saturation limits (Lee K. W. & Khalil H. K. 1997).

Fig.2

Position tracking error of joint q₁ and q₂. (matched case).

Fig.3

Induced torque control.

In mismatched case, the frictions are added to the joint of robot manipulator model in equation (1) and are modeled: F _r = f _s q + f_v sign(q̇), with the values f _s = f_v = diag(5,5).

If we regard an unknown load carried by the robot as part of the second link, then the parameters m₂, l_c2 and I₂ will change to: m₂ + Δm₂, l_c2 + Δl_c2 and I₂ + ΔI₂.

Let Δm₂ = 5Kg, Δl_c2 = 0.5m and and $\Delta I_2=\frac{1}{6}Kg{m}^2$ to be the maximum parameters variations of the second link due to unknown load. It is observed from Figure 4 that the output q₂(t) tracks tightly the reference trajectory with a steady state error. This results has been expected in the analysis did in section 4, i.e the uncertainties will introduce a steady state error in tracking error. Furthermore, the induced torque control lie outside the saturation limits. Figure 5 shows the results when the control law (21) is applied to robot, it is seen from this figure the error was eliminated and the torque control signal lie inside the saturation limits. These results prove the robustness of the rigid link manipulator under the approximate receding-horizon controller with integral action to payload uncertainties and viscous friction.

Fig.4

Position tracking performance in mismatched case.

Fig.5

Torque and tracking performance with integral action.

In addition, other simulations have been carried out and the following remarks have been observed:

Decreases in performances are obtained when we increase the step control parameter h.

High dynamics of the reference trajectories results in high increase of the control torque signal. To reduce this control torque amplitude, one should increase the predictive time increment h. It should be pointed that over a threshold value of h noted h_max, the performance decrease and instability mechanism will appear. This is due to the Taylor approximation used to derive the predictive controller, which becomes invalid.

7. Conclusion

In this paper, the approximation of receding-horizon controller of rigid link robot manipulator using output feedback via link position measurements were considered. Minimizing a quadratic function of the predicted tracking error and the predicted input over the fixed horizon, by using the Simpson's rule approximation, derives the control law. One of the main advantages of these control schemes is that it does not need to perform an online optimization and asymptotic tracking of the smooth reference signal is guaranteed.

To enhance the robustness property of the nonlinear predictive developed by Ping Lu, we proposed to incorporate an integral action in the loop. Simulation shown that payload uncertainties and friction have no effect on the robot manipulator closed by the proposed algorithm. Moreover, the obtained torque signal lie in saturations limit. The Lyapunov theory is used to prove the asymptotic stability of equilibrium point of both original and augmented system.

Finally, we expect that the results presented here can be explored and extended to discrete implementation of these continuous-time predictive controllers either through computers or by special purpose chips that can run at a higher speed.