Min–Max Optimal Control of Robot Manipulators Affected by Sensor Faults

This paper is concerned with the control law synthesis for robot manipulators, which guarantees that the effect of the sensor faults is kept under a permissible level, and ensures the stability of the closed-loop system. Based on Lyapunov’s stability analysis, the conditions that enable the application of the simple bisection method in the optimization procedure were derived. The control law, with certain properties that make the construction of the Lyapunov function much easier—and, thus, the determination of stability conditions—was considered. Furthermore, the optimization problem was formulated as a class of problem in which minimization and maximization of the same performance criterion were simultaneously carried out. The algorithm proposed to solve the related zero-sum differential game was based on Newton’s method with recursive matrix relations, in which the first- and second-order derivatives of the objective function are calculated using hyper-dual numbers. The results of this paper were evaluated in simulation on a robot manipulator with three degrees of freedom.

Furthermore, an optimal and robust control in the presence of sensor and actuator faults has, for the last decade, been a topic of import for researchers on the control of dynamical systems. Many authors consider the L 2 -gain of a nonlinear system (often called nonlinear H ∞ , because, according to [25,26], the L 2 -gain is equivalent to the H ∞ norm of a linear system) to be a measure of the influence of the faults. A fixed L 2 -gain fault-tolerant controller, for a class of Lipschitz nonlinear system with actuator saturation, was designed in [27]. In [28], a data-driven output-feedback approach to the fault-tolerant control (FTC) problem, considering L 2 -gain properties, was proposed. A robust H ∞ FTC scheme to regulate the quadrotor system was adopted in [29]: the approach was based on linear matrix inequalities (LMI), so the overall dynamics of the quadrotor were linearized. In [30], backstepping and adaptive control methods based on L 2 -gain were combined together, to provide a passive fault-tolerant attitude controller of a nonlinear dynamic model of morphing aircraft.
All the above-mentioned works were based on the formulation of the problem in the form of LMI, or on the determination of the approximate solution of the associated Hamilton-Jacobi-Isaacs (HJI) equation. The application of LMI requires the linearization of system dynamics, and therefore robustness cannot be guaranteed in all operating points; on the other hand, solutions of the HJI equation can be complex, and therefore difficult to apply in real control tasks.
The difficulties described above motivated us to conduct the research presented here. We formulated the control problem of a robot manipulator affected by sensor faults, as a typical representative of nonlinear dynamic systems, in the form of optimal L 2 -gain control and the related min-max optimization, and we approached its solution without LMI formalism or the need to approximate the solution of the HJI equation.
The main idea of the approach presented in this paper was based on the algorithmic procedure described in the authors' previous works [31,32]. The idea of using a Newton-like algorithm with recursive matrix relations to calculate exact gradients and Hessians was applied to the synthesis of a control system which aimed to overcome the sensor malfunctions whilst maintaining desirable stability and optimal L 2 -gain performance properties.
This paper considered the synthesis of the PID control law of a robot manipulator, which kept the influence of sensor faults below the permissible limits, and ensured the robust stability of the overall closed-loop system. The problem was presented as a twoplayer zero-sum differential game, with the objective function including a parameter in such a way that the control vector represented the "player" that minimized the objective function, while the vector containing the sensor faults represented the "player" that maximized the same objective function.
The algorithm proposed in this paper, in one part, was based on Lyapunov's stability analysis. The derived stability conditions were used in the optimization process for the appropriate initial setting of the algorithm parameters: these conditions were crucial in the application of the bisection method by which we determined the minimum L 2 -gain. Also, by using Lyapunov method, the stopping criteria of the algorithm was improved, by deriving inequalities dependent on PID controller gains that bound the distance of the current point from the saddle point, and thus added bound on the number of iterations.
The main contributions of this work are summarized as follows: • A suitable mathematical tool and systematic algorithmic procedure for the synthesis of an optimal robust PID controller for robot manipulators affected by sensor faults was developed. The presented approach admitted a min-max formulation, and hence provided a guarantee of robustness, which was achieved by optimizing the worst-case performance. To the best of the authors' knowledge, such an algorithmic PID controller synthesis procedure had not previously been investigated. • In many other similar optimal control-based algorithms that have been proposed in the literature (see, for example, [33,34] and references therein), the nonlinear system dynamics were treated as equality constraints, and were included in the optimization process, using the method of Lagrange multipliers: the results were HJI equations that were very difficult or almost impossible to solve. For this reason, many approximation methods have been developed (see, for example, [35][36][37][38] and references therein) in which actual computational complexity increased with the number of system states which needed to be estimated. In contrast to such approaches, the algorithm proposed in this paper had no high-dimensional structure: this followed from the fact that, instead of incorporating the robot dynamics in the closed loop with the controller directly into an objective function, and solving the corresponding HJI equation, the state variables and PID controller gains were coupled by recursive matrix relations, used to calculate the first-and second-order derivatives that appear in the Newtonlike method.
Based on the contributions noted above, the importance of the research presented in this paper lies in the fact that the proposed control architecture does not rely directly on sensor fault estimation, in order to improve the desired positioning of the robot manipulator, and is closely related to nonlinear L 2 -gain robust control, where a PID-type controller, with its well-known simplicity of structure, is designed to be robust against a sensor fault in the system. Such a control strategy has the potential for practical applications, due to its simplicity of design, it having less lag between sensor fault appearance and accommodation, and it not requiring significant computing resources.
The rest of the paper is organized as follows. In Section 2, the robot manipulator dynamics and their main properties are presented. In Section 3, the min-max optimal control problem, of a robot manipulator affected by sensor faults and in a closed loop with a PID controller, is defined. The main results are presented in Sections 4 and 5. Based on the theory of L 2 stability and Lyapunov's approach, the stability conditions of the entire control system are established in Section 4. In Section 5, the bisection method, the Newton-like algorithm, the Adams discretization method and the recursive matrix calculation of gradients and Hessians are integrated into an efficient algorithmic procedure. The closed-loop response of the proposed control strategy and a robot manipulator with three degrees of freedom (3DOF) are evaluated in computer simulations, and the results are given and discussed in Section 6. Section 7 concludes the paper. There are three appendices included. Appendix A describes the notation used throughout the paper. In Appendix B, a detailed derivation of expressions for the kinetic and potential energy of the cylindrical robot manipulator is given. Appendix C contains the expressions for the evaluation of coefficients resulting from the mechanical properties of a cylindrical robot manipulator with 3DOF.

Dynamic Model and Properties
Consider a dynamic model of an N-DOF robot manipulator system of the form where q ∈ R N is the vector of generalized coordinates, u ∈ R N is the vector of control forces/torques applied to the system, M(q) ∈ R N×N is the inertia matrix, C(q,q)q ∈ R N is the vector of centrifugal and Coriolis forces/torques and g(q) ∈ R N is the vector of gravitational forces/torques. It is well-known that Equation (1) can be obtained through modeling a robot manipulator system by Euler-Lagrange equations. For this paper, we considered a robot manipulator that had both rotational and translational generalized coordinates, and hence had the following properties (see, for example, [39][40][41][42][43]): is skew-symmetric for all q,q ∈ R N : this implies thaṫ M(q) = C(q,q) + C(q,q) T ; Property 2. The inertia matrix M(q) is a positive-definite symmetric matrix which satisfies for all ξ and q ∈ R N , where a 1 , a 2 > 0, and c 2 , d 2 ≥ 0.

Property 3.
There exist constants, c 1 and d 1 , such that the Coriolis and centrifugal term, C(q,q)q, satisfies for all q andq ∈ R N , where c 1 > 0, d 1 ≥ 0.
As is well-known, the vector of gravitational forces/torques is a gradient of the potential energy. The considered class of robot manipulator was such that its potential energy depended linearly on translational generalized coordinates, while rotational generalized coordinates appeared as a trigonometric function, with period 2π. Hence, the following property was imposed: Property 4. There exist positive constants, k g1 and k g2 , such that the Jacobian of the gravity vector satisfies ∂g(q) ∂q ≤ k g1 + k g2 q , (6) for all q ∈ R N . Remark 1. If the system (1) has no translational generalized coordinates, then c 2 , d 2 , d 1 , k g2 = 0. See, for example, [44][45][46][47], and references therein.

Optimal Control Problem
As robot manipulators often work in extreme environments, their sensors are prone to faults during operation: for example, due to environmental noise, like a magnetic field. Faults in robot manipulator sensors often occur due to vibrations that cause disruptions in communication between the control unit and the sensor, or even a short circuit. Furthermore, intermittent sensor connection, bias in sensor measurement and sensor gain drop can also be presented as faults.
For this paper, we were interested in designing a decentralized system controller (1) that guaranteed the stability around a desired constant of generalized coordinate positions q d ∈ R N , while the influence of sensor faults was kept under a permissible level. In this context, we considered the least upper bound (supremum) of the ratio of the L 2 -norm of the vector of the output signals to be controlled, and the L 2 -norm of the vector representing sensor faults. It should be noted that considering the control problem in this way was actually equivalent to the L 2 -gain (or nonlinear H ∞ ) optimal control defined in [25,26].
In a sensor fault situation, the internal robot dynamic properties are not affected, but the controller of the robot manipulator (1) receives generalized coordinates that are not exact, and are given as follows [2,8,48,49]: where ϕ p ∈ L 2 t 0 , t f , R S and ϕ v ∈ L 2 t 0 , t f , R S are the vector functions representing sensor faults, and F p and F v are the matrices with appropriate dimensions. The assumptions about the sensor fault terms were employed as follows: Assumption 1. There exist positive-definite time-dependent functions f p and f v , such that the sensor faults are bounded by Remark 2. In many other works-for example, in [2,7,8,50]-sensor fault bounds are considered as known positive real constants. In this work, we did not treat these bounds as known specific constants, but as unknown functions that needed to be determined in such a way that they maximized the objective function, so that the proposed control strategy might be available not only for robotic manipulators but also for a wide class of nonlinear dynamical systems.
A PID-type control law was proposed, with gravity vector compensation in the following form: whereq =q − q d was a generalized coordinates errors vector, and K P , K D and K I were N × N positive-definite diagonal matrices. Note that, as K P , K D and K I were diagonal matrices, the proposed control law was formed in a fully decentralized fashion. The robot manipulator system (1) with mixed rotational and translational generalized coordinates, affected by sensor faults (7) and (8) in a closed loop with a PID-type controller (11) and (12), can be represented by a block diagram, as shown in Figure 1. Each DOF of the robot manipulator was driven by servo motors with gears or pneumatic and hydraulic cylinders. The decentralized PID control law was implemented in the microcontroller, and point-to-point communication was established by electrical interface and physical layers. The received encoder data were used to close the feedback loop in the microcontroller.
Regarding the block diagram shown in Figure 1, in our approach the emphasis was not on the structure of the controller itself, but on the method for adjusting the gains of the controller. The PID controller structure was considered because it is the most common structure used in industrial applications.

Remark 3.
It is important to emphasize the main difference between a conventional PID controller and our proposed one. In conventional applications of a PID controller, the gain values are independently and freely chosen by the designer, often using a trial-and-error method; however, in our approach, a proportional gain matrix K P , a derivative gain matrix K D and an integral gain matrix K I were tuned automatically and correlatively by a numerical algorithm that will be presented in detail in the following sections. Inserting (7) and (8) into (11) and (12) gave Finally, inserting (13) and (14) into (1) gave closed-loop equations in the following form: Based on the definition of finite L 2 -gain and the definition of L 2 stability (see, for example, [25,51]), we defined the min-max optimal control problem of a robot manipulator affected by sensor faults, as follows: Problem 1. Determine the PID controller gains, which are the elements of the matrices K P , K D and K I in (11) and (12), and determine the "worst case" of sensor fault functions ϕ p and ϕ v in (7) and (8), such that the system (15) and (16) is stable around the desired constant generalized coordinates positions, the ratio of the L 2 -norm of the vector of the output signals is controlled and the L 2 -norm of the vector representing sensor faults is minimized. In other words, solve the following zero-sum differential game: whereq, K P , K D , K I and ϕ = ϕ T p ϕ T v T are coupled via closed-loop system dynamics (15) and (16), and µ > 0 is a finite L 2 -gain; furthermore, q 0 is an a priori known vector of initial states of generalized coordinates.

L 2 Stability Conditions
It is well-known from the dissipativity theory [52][53][54] that a nonlinear dynamic system is L 2 -stable if, for all initial conditions, all uncertainties w and all t f ≥ t 0 there exists a storage function V, which is also a Lyapunov function, such that the following inequality holds: where z is the performance variable, and γ > 0 is the finite L 2 -gain. Note that the right-hand side of the inequality (18) in our case actually corresponded to the argument of the min-max operator in (17). Based on this observation, we give the following proposition, by which we set the stability conditions of the control system proposed in this paper: Proposition 1. If the following conditions are satisfied: for some α > 0, where µ > 0 is a finite L 2 -gain, and where a 1 , c 1 and d 1 are defined in Properties 2 and 3, then the system (15) and (16) is locally L 2 -stable.
Proof. The first step was to transform the system (15) and (16) into a form with a zero steady state. The stationary state of the system (15) and (16) isq = 0,q = 0 ⇒ q = q d , so it was obtained thus: where υ * , ϕ * p and ϕ * v were steady state values. Subtracting (16) and (25) gaveq = −F p ϕ p − ϕ * p ; then, subtracting (15) and (24), we obtained the error equations in the following form: Following the methodology from [39,46], as the first step in the construction of the Lyapunov function, the error Equations (26) were multiplied on the left side by αq T +q T with some α > 0: Some of the terms in the expression (27) can be written as follows: Inserting (28)- (32) in (27) and, asṀ(q) − 2C(q,q) was skew-symmetric (see Property 1), which impliedq T Ṁ (q) − 2C(q,q) q = 0, we obtained: Based on the nonlinear differential form (33), it could be concluded that on the lefthand side we had the Lyapunov function candidate V(q,q,υ), while on the right-hand side we had a candidate for the time derivative of the Lyapunov functionV(q,υ,φ v ).
The Lyapunov function candidate could be rewritten as follows: In order for the above expression to be positive-definite, it was necessary to determine the conditions under which Using definitions of vector and matrix norms, bound of quadratic forms (see Notation, Appendix A) and, finally, applying Property 2 we obtained the following: which was positive-definite if Next, according to the theory of L 2 stability and the passivity properties of Euler-Lagrange systems, the following inequality needed to be satisfied: where µ > 0 was a finite L 2 -gain. Comparing (38) to the right-hand side of (33), using Cauchy-Schwarz inequality, triangle inequality, definitions of vector and matrix norms, bound of quadratic forms (see Notation, Appendix A) and, finally, applying Properties 2 and 3, we got the following inequality: The conditions that satisfied the above inequality were as follows: Remark 4. As the condition (40) depended on the generalized coordinates of the system, it was concluded that the proposed controller guaranteed only local stability. For a more detailed analysis of stability, it would be necessary to determine the domain of attraction that guarantees asymptotic stability; furthermore, in order to ensure global stability, it would be necessary to introduce a nonlinear term in the control law (see, for example, [39]). Determining the domain of attraction and deriving the global stability conditions were beyond the scope of the research presented in this paper. The stability conditions derived in this work were calculated numerically, and served only for the appropriate initialization of the algorithm, thereby ensuring satisfactory convergence properties.

Algorithm for the Controller Synthesis
In order to derive the algorithm for the controller synthesis, i.e., the algorithm to solve Problem 1, the PID control law (13) had to be written in parametrized form first. To do this, a vectorization operation was performed: Then, by substituting (15) and (16), we got the following first-order nonlinear differential equations: where Note that u = R(x, ϕ)k.
In accordance with the previous substitution, expression (17) then became where x, k and ϕ were coupled via nonlinear differential Equations (45) and (46).

L 2 -Gain Minimization
The stability conditions from Proposition 1 suggested that a simple and well-known schematic bisection algorithm [55] could be applied, to minimize the L 2 -gain: the steps of this method were as follows.
Step initialization: Choose the initial elements of the vector k, i.e., the initial PID controller gains, such that the lower µ l 0 and upper µ u 0 bound of the L 2 -gain satisfy the conditions in expression (19). Choose a small enough positive constant, , as the stopping criteria of the bisection method.
Step 2: Set then calculate the value of function J * µ k (x 0 ) by solving the zero-sum differential game (47).

Zero-Sum Differential Game Solution
In the second step of the algorithmic procedure described in the previous subsection, it was necessary to solve the problem defined by expression (47). To solve this zerosum differential game, we implemented the Newton method, which is described in the following steps.
Step initialization: Choose a small enough positive constant, ε, as the stopping criterion of the Newton algorithm.
Step 2: Determine the search direction vector s j by solving a system of linear equations using Cholesky factorization: where J is the argument of the min max operator in (47), ∇ k J, ∇ ϕ J, ∇ 2 k J, ∇ 2 ϕ J, and ∇ 2 k,ϕ J are gradients and Hessians with respect to the vectors k and ϕ, respectively. Note that the maximization with respect to ϕ is achieved simply by the minus sign in front of the gradient and Hessians.
Step 3: Use the line search strategy satisfying the Wolfe conditions [56] to compute the step-size η j > 0.
As is well-known, Newton's method has a locally quadratic convergence to the saddle point, assuming that the initial point is close enough to the saddle point. In the approach we propose here, the proximity to the saddle point is ensured by the L 2 stability conditions set in Proposition 1: this further implies that ∇ 2 k J > 0 and −∇ 2 ϕ J > 0. As these Hessians are positive-definite submatrices on the main diagonal of the matrix on the left-hand side of (49), and submatrices ∇ 2 k,ϕ J and −∇ 2 ϕ,k J are negative transpositions of each other, the matrix on the left-hand side of Equation (49) is asymmetric positive-definite and, therefore, the linear system (49) is well-defined; therefore, the previously described algorithmic procedure, which is based on the Newton method, produces: in the k-th iteration of the bisection algorithm proposed in Section 5.1.
To perform the steps of the proposed Newton-like algorithm, we needed expressions for the gradients and Hessians that appear in (49). A detailed procedure for deriving recursive matrix relations for computing these gradients and Hessians is given in the next subsection.

Matrix Relations for Recursive Calculation of Gradients and Hessians
In order to derive recursive relations for calculating gradients and Hessians, we used the fact that derivatives are a measure of the sensitivity of functions to small changes in their variables. This suggested that, for their calculation, we could perform a time discretization of the system dynamics (45). In the approach presented in this paper, the discretization of system dynamics using the fourth-order Adams approximation with a small time step was applied. In [32] it is shown that an explicit Adams method can be conveniently transformed into discrete-time state space.
The explicit Adams approximation of system (45) can be simply written in the following discrete-time state-space form: such that the time grid consists of points t i = t 0 + iτ for i = 0, 1, 2, . . . , N − 1, where τ = (t f − t 0 )/N is the time step length, andx is the extended 4n-dimensional state vector (n is a dimension of state vector x in (45), and 4 is the order of the Adams method). More details on the Adams method can be found in [57]. The discrete-time form of the objective function resulted in where F was the sub-integral function of the argument of the min max operator in (47) The gradient of the objective function (55), with respect to the k in the j-th iteration of the Newton algorithm (Section 5.2), was given by for l = 0, 1, 2, . . . , N − 1. Note that, because of the causality principle, the terms for i < l were equal to zero. As the terms under the sum on the right-hand side of (57) depended on k(l) implicitly throughx(i) for i > l, this gave Using the chain rule for ordered derivatives, from (54) it followed that for i = l + 2, . . . , N − 1.
To simplify the following expressions, we introduced Starting from l = N − 2; i = N − 1, it followed that then, substituting (58) in (61) and substituting (59) in (62), and taking into account (54), we obtained The above procedure could be further continued for l = N − 3; i = N − 1, N − 2; l = N − 4; i = N − 3, N − 2, N − 1, etc., and the final recursive expressions for calculation of gradient ∇ k J had the following form: for l = N − 2, N − 3, . . . , 0. Before giving recursive expressions for computing the Hessians, we adopted the usual convention such that, for some function z = f (v), the Hessian is defined by Then, in order to calculate ∇ 2 k J in the j-th iteration of the Newton algorithm (Section 5.2), the derivatives of Equations (64)-(67), with respect to vectorsx and k, were taken. We obtained the following recursive matrix relation: for l = N − 2, N − 3, . . . , 0.
In the previous expressions, (64)-(71), the derivatives of the system dynamics (54) and sub-integral function (56), with respect to the vectorsx and k, were calculated by the hyper-dual number method [58,59], which enabled us to get the first-and second-order derivatives in one step, accurately to machine epsilon, which cannot be achieved by the finite difference method. Compared to the widely used and popular method of adjoint automatic differentiation, the hyper-dual number method is more robust and easier to apply. In the previous considerations from expression (54)-(71) the derivation of recursive matrix relations for the calculation of gradients and Hessians, with respect to the vector k containing the gains of the PID controller, are shown. Matrix relations for calculating gradients and Hessians, with respect to vector ϕ containing the sensor faults, can be obtained by the same procedure, with obvious changes in notation, and there is no need to give them here.
A flowchart of the entire algorithmic procedure described in Sections 5.1-5.3 is shown in Figure 2.

Numerical Simulations
In this subsection, the simulation results for control of the robot manipulator, using the controller synthesis algorithm that is described in the previous sections are presented. Figure 3 shows the considered structure of a robot manipulator. The structure of a robot manipulator of this type is also the subject of consideration in [60,61]. However, for the sake of completeness, in this paper a more detailed derivation of mathematical model is presented. Detailed expressions for kinetic and potential energy and parameters that define the bounds of the robot's dynamic properties are given (see Appendices B and C).
The first joint is rotary, which enables the rotation of the robot around the vertical axis, while the second and third are prismatic joints, which enable movement in the vertical and horizontal directions, respectively. Such robot manipulators have a workspace in the shape of a cylinder. Robot manipulators with a cylindrical structure, due to their compact design, are most often used for simple assembly tasks, handling at machine tools, applying coatings. Also, such robots can be fast, so a compromise should be sought because speeds imply problems with rotational inertia, which can affect repeatability if the entire system is not configured in accordance with its capabilities. It is usual for cylindrical robots that the up-and-down motion is achieved by pneumatic or hydraulic cylinders, while the rotation is usually achieved by electric motors and gears.
To model the kinematics and dynamics we denoted by q 1 (rad) the rotational coordinate of the first joint, q 2 (m) and q 3 (m) the translational coordinates of second and third joint, respectively.

Numerical Simulations
In this subsection, the simulation results for control of the robot manipulator, using the controller synthesis algorithm that is described in the previous sections, are presented. Figure 3 shows the considered structure of a robot manipulator. The structure of a robot manipulator of this type was also the subject of consideration in [60,61]; however, for the sake of completeness, in this paper a more detailed derivation of the mathematical model is presented. Detailed expressions for kinetic and potential energy, and parameters that define the bounds of the robot's dynamic properties, are given (see Appendices B and C).
The first joint is rotary, which enables the rotation of the robot around the vertical axis, while the second and third are prismatic joints, which enable movement in the vertical and horizontal directions, respectively. Such robot manipulators have a workspace in the shape of a cylinder. Robot manipulators with a cylindrical structure, due to their compact design, are most often used for simple assembly tasks, handling machine tools, applying coatings, etc. In addition, such robots can be fast, so a compromise should be sought, because speed implies problems with rotational inertia, which can affect repeatability if the entire system is not configured in accordance with its capabilities. It is usual for cylindrical robots that the up-and-down motion is achieved by pneumatic or hydraulic cylinders, while the rotation is usually achieved by electric motors and gears.
To model the kinematics and dynamics we denoted by q 1 (rad) the rotational coordinate of the first joint, and by q 2 (m) and q 3 (m) the translational coordinates of the second and third joints, respectively. First, the forward kinematic equations for the cylindrical structure of a robot manipulator with 3DOF, using the Denavit-Hartenberg convention were developed. The Denavit-Hartenberg parameters are given in Table 1, where θ = q 1 , d 1 = 0, d 2 = L 1 + q 2 and d 3 = L 2 + q 3 . The homogeneous transformation matrices were as follows: The transformation matrices were as follows: where 0 T 3 was the corresponding Denavit-Hartenberg matrix and, based on the last column, we had the position vector of the robot end-effector in the Cartesian space, as follows: The solution of the inverse kinematic problem could be obtained, based on the following equation:     p x cos(q 1 ) + p y sin(q 1 ) −p x sin(q 1 ) + p y cos(q 1 ) Based on the last columns of the previous matrices on the left and right sides, we obtained: −p x sin(q 1 ) + p y cos(q 1 ) = L 2 + q 3 =⇒ q 3 = −L 2 − p x sin(q 1 ) + p y cos(q 1 ).
In order to derive the corresponding dynamic equations of motion, the kinetic energy of the robot was first determined. The overall kinetic energy was calculated as the sum of kinetic energies from corresponding links, as follows: where Next, the overall potential energy was calculated, as the sum of potential energies from corresponding links, as follows: where The Euler-Lagrange equations of motion for a considered robot manipulator were given by where T was the total kinetic energy defined by (85)-(88), U was the total potential energy defined by (89)-(92), and τ i was the torque/force applied to the i-th robot link. Energy dissipation in rotational joints, and viscous friction in translational joints were neglected. The Equation (93) can be written in matrix form (1), such that q = [q 1 q 2 q 3 ] T , u = [τ 1 τ 2 τ 3 ] T and A detailed derivation of expressions for the kinetic and potential energy of the considered robot manipulator is given in Appendix B.
The numerical values of the parameters relevant to the derivation of the dynamic model of the considered robot, using the Euler-Lagrange formalism, are given in Table 2. Furthermore, for the appropriate initialization and numerical efficiency of the proposed algorithmic procedure for PID controller synthesis, we needed the constant values of the coefficients defined in properties (4)- (6). The numerical values of these parameters are given in Table 3, and the expressions for determining them for the considered cylindrical robot are given in Appendix C. The parameters of the algorithm proposed in this paper, for this particular robot manipulator, were set as follows. The vector of initial conditions of the generalized coordinates was q 0 = 0. The initial values of the PID controller gaining satisfying stability conditions were chosen as The stop criterion of the bisection method was chosen as = 10 −3 . The criterion for stopping Newton's method was set as ε = 10 −4 . The final time was t f = 2 sec., and the number of optimization time intervals was N = 2000, so that the sampling interval was τ = 0.002 sec.
The desired positions of the generalized coordinates were selected as follows: q d1 = π/2; q d2 = 0.2 m; q d3 = 0.1 m. By running the algorithm, we obtained the following PID controller gains: and the minimum L 2 -gain, µ = 532.3096. The time responses of the rotational coordinates of the first joint q 1 and the translational coordinates of the second and third joints, q 2 and q 3 , respectively, are presented in Figure 4. It can be seen from the figures that the transient response lasted approximately 0.8 seconds. The duration of the transient response could be further improved by adding weighting matrices into the objective function (47). When selecting these matrices, one should strive to reach a compromise between the time required to achieve the desired position, and the maximum value of the control torque/force, such that the effect of the sensor faults is kept under a permissible level, and ensures the stability of the closed-loop system. Figure 5 illustrates the time dependence of the positioning errors. It can be seen that the errors from the desired positions were approximately 10 −3 . These errors could be further reduced by reducing the parameters and ε of the bisection and Newton methods. Reducing these parameters would mean an increase in the number of iterations and, thus, the execution time of the algorithm. Furthermore, the positioning errors could also be reduced by choosing a smaller sampling time τ of the Adams method; however, as is well-known, there is a minimal value of τ that guarantees numerical stability.   Figure 6 shows the time dependence of the force and torques applied to the robot links. Furthermore, Figure 7 illustrates the simulation results for the time dependence of the sensor faults calculated by the algorithm proposed in Sections 5.1 and 5.2. Note that f p1 , f p2 , f p3 , f v1 , f v2 and f v3 are elements of vector ϕ, which is incorporated in the objective function of the zero-sum differential game (17), i.e., (47) and, since ϕ is the maximizing player, the sensor fault functions shown in Figure 7 affected the robot manipulator system in the worst possible manner. The obtained results can be interpreted as a sensor failure due to a deviation belonging to the class of bounded L 2 functions.
It is obvious from the figures that the robot manipulator reached the desired positions at an acceptable settling time, with a negligible steady state error and with acceptable amounts of applied forces and torques to achieve the desired motion: therefore, it can be concluded that the proposed control strategy is efficient in the presence of sensor faults.

Discussion of Comparison with Other Methods
Here, we discuss the features of our approach, in comparison to other similar approaches.
In fault-tolerant control approaches based on LMI formalism (see, for example, [6,7,23,29] and references therein), the nonlinear dynamics of the robot manipulator must be linearized, which means that robustness cannot be guaranteed in all operating points. Then, the decentralized PID controller must be transformed into a state feedback controller or an output feedback controller, and the optimization problem is presented in the form of an LMI, in which it is necessary to determine the elements of the positive-definite Lyapunov matrix and the transformation matrix. In the case of a robot with 3DOF, the Lyapunov matrix is 6 × 6 and the transformation matrix is 3 × 6: this means that the problem has at least 21 + 9 optimization parameters. In our approach, the PID controller gains were directly optimized, which meant that we only had 9 optimization parameters. Furthermore, the Mehrotra-type predictor-corrector variant of interior-point method algorithms is commonly used to solve the LMI problem: in our approach, we use Newton's method, which-near the saddle point of the associated differential game-produces a sequence that quadratically converges to a desired solution, while the Mehrotra-type predictor-corrector has linear convergence. Thus, compared to methods based on solving LMI, the approach proposed in this paper has less optimization variables, and requires less iterations.
Many of the approaches to fault-tolerant control (see, for example, [2,6,12] and references therein), in addition to the synthesis of the controller, also imply a certain detection or estimation of the faults that appear in the nonlinear dynamic system: this actually means additional computational requirements that increase with the number of state variables that need to be estimated due to the synthesis of an additional subsystem. Compared to these approaches, our approach does not require an additional subsystem to determine how to modify the controller structure and gains, but focuses on the L 2 -gain robust stability of the robot manipulator, considering the worst-case scenario of sensor faults rather than the desired performance for each fault occurrence scenario.
The algorithmic procedure presented in this paper is intended for the offline solution of the zero-sum differential game related to the L 2 -gain optimal control problem, with the explicit calculation of the PID controller gains: therefore, the computational complexity of our approach is not as much of a limitation as for approaches intended for online execution, such as model predictive control (see, for example, [15,16] and references therein).

Conclusions
In this paper, the tuning of PID-type controller gains for robot manipulators affected by sensor faults, using an algorithmic procedure that gives an explicit solution to the L 2 -gain optimality criterion and related zero-sum differential game, is presented. The main contribution of the paper includes the integration of the simple bisection method, the Newton method for solution of related zero-sum differential games without solving the HJI equation, the Adams method for time discretization, and hyper-dual numbers, to provide an effective method for control law synthesis. By applying Lyapunov stability analysis, and the dissipativity theory and passivity properties of systems described by Euler-Lagrange equations, we derived local L 2 stability conditions that we used for appropriate initialization of the bisection method, and for ensuring and controlling the convergence of Newton's method. The simulation results for control of the 3DOF cylindrical robot manipulator showed that the proposed algorithmic procedure could efficiently calculate the controller gains, in order to position the system affected by sensor faults, as desired.
The extension of the proposed approach could be continued in the following directions: • Although the case of sensor faults is considered in this paper, the proposed algorithm for control law synthesis could easily be extended to the case of dynamic systems affected by actuator faults, without significant increase in its complexity. • Improvements in the control strategy proposed in this paper could also go towards the synthesis of a complete model-free control law, i.e., control law without gravity vector compensation: this would complicate the derivation of the stability conditions, and the controller gains should be presented as nonlinear functions of generalized coordinates. • Instead of the initial conditions being known in advance, one could consider a case where the initial conditions were treated as an unknown uncertainty, i.e., the variables of the min-max optimization problem that maximizes the objective function. • From the numerical optimization algorithms point of view, to perform a detailed analysis and comparison between a proposed algorithm and other existing methods, such as genetic algorithm or particle swarm optimization, in terms of convergence, accuracy and computational efficiency.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: DOF degrees of freedom FTC fault-tolerant control HJI Hamilton-Jacobi-Isaacs LMI linear matrix inequalities PID proportional-integral-derivative

Appendix A. Notation
Matrices and vectors are represented in bold upper and bold lower case, respectively. Scalars are represented in italic lower case. I is an identity matrix, and 0 is a null matrix. The dimensions of the matrices and vectors can generally be determined trivially by the context. The symbols ∇ and ∇ 2 stand for the gradient and Hessian, respectively. The symbol T denotes transposition.
The vec(·) is an operator that stacks the columns of a matrix one underneath the other. The Kronecker product of the two matrices A (m × n) and B (p × q), denoted by A ⊗ B, is an mp × nq matrix defined by A ⊗ B = (a ij B) ij . The definitions of the matrix differentials calculus, and the algebras related to Kronecker products, can be found in [62,63].
The Euclidean norm of vector is defined as x = √ x T x. λ min {A} and λ max {A} represent the smallest and the largest eigenvalues, respectively, of the symmetric positive-definite matrix A. The induced norm of matrix A is defined as where σ max {A} represents the largest singular value of matrix A. If A is a symmetric positive-definite matrix, it implies A = λ max {A}; furthermore, for a symmetric positivedefinite matrix, it holds that L 2 (I, R n ) stands for the standard Lebesgue spaces of vector-valued square-integrable and essentially bounded functions mapping an interval I ⊂ R to R n : this space is equipped with an L 2 norm defined by · L 2 = t f t 0 · 2 dt. We avoided explicitly showing the dependence of the variables from the time when not needed.

Appendix B. Detailed Derivation of Expressions for the Kinetic and Potential Energy of the Cylindrical Robot
The movement of the first link of the robot can be viewed as a rotation around an axis, so the kinetic energy can be determined as follows: where I 1 is dynamic moment of inertia I 1 = 1 2 m 1 R 2 . The potential energy of the first link is defined by where T is the vector of the position of the center of gravity of the infinitesimal mass of the first link in relation to the fixed coordinate system, and g = 0 0 −g 0 T is the vector of gravitational acceleration. By inserting them into the Equation (A2), we get The kinetic energy of the second link is If we assume that the second link is homogeneous, then it follows: where dm 2 is infinitesimal mass and du 2 is infinitesimal displacement. Substituting expression (A5) into (A4) and changing the bounds of integration, it follows: The position of the dm 2 in the robot coordinate system is determined by the position vector The position vector of the mass dm 2 , with respect to the fixed coordinate system, is By inserting expressions (77) and (A7) into expression (A8), we get The velocity of the mass dm 2 is and the square of the velocity is Inserting Equation (A11) into (A6), and integrating from 0 to L A , it follows: A 3q 2 1 cos 2 (q 1 ) + sin 2 (q 1 ) + L Aq The potential energy of the second link is calculated as follows: where the vector p 2 is determined by (A9) for u 2 = L A /2, so the potential energy is equal to The same assumptions about the homogeneity of the second link and the linearity of the mass distribution are applied to the calculation of the kinetic and potential energy of the third link.
The kinetic energy of the third link is and, with the dm 3 = m 3 L 3 du 3 , it follows: The position vector of the infinitesimal mass dm 3 in the robot coordinate system is The position vector of the mass dm 3 , with respect to the fixed coordinate system, is equal to p 3 = 0 T 3 R 3 =     u 3 sin(q 1 ) − (L 2 + q 3 ) sin (q 1 ) (L 2 + q 3 ) cos(q 1 ) − u 3 cos(q 1 ) where the transfer matrix from the third coordinate system to the fixed coordinate system is defined by (78). The velocity of the mass dm 3 is v 3 = dp 3 dt =     u 3q1 cos(q 1 ) − (L 2 + q 3 )q 1 cos(q 1 ) −q 3 sin(q 1 ) u 3q1 sin(q 1 ) − (L 2 + q 3 )q 1 sin(q 1 ) +q 3 cos(q 1 ) while the square of the velocity is v 2 3 = v 3 · v 3 =q 2 1 (L 2 − u 3 + q 3 ) 2 +q 2 2 +q 2 3 . (A20) By inserting Equation (A20) into (A16), and after integrating within the limits from 0 to −L 3 , we get kinetic energy of the third link as follows: The vector of the position of the center of gravity of the third link with respect to the fixed coordinate system is obtained by inserting u 3 = L 3 /2 into (A18): The potential energy of the third link is calculated as follows: according to which, we get

Appendix C. Properties of the Dynamic Model of the Cylindrical Robot
Based on the inertia matrix defined by (94) and (95), using expression (4) with q = [q 1 q 2 q 3 ] T and ξ = [ξ 1 ξ 2 ξ 3 ] T , the parameters a 1 , a 2 , c 2 and d 2 are determined as follows.
We can write the above expression as follows: from which it follows: that is, we have Based on the Coriolis vector (96), using expression (5), the parameters c 1 and d 1 are determined as follows.
Inserting (A33) into (A32) gives the following inequality: so that, in the end, we get The parameter k g2 from (6) is equal to zero, because the considered cylindrical robot structure does not have translational degrees of freedom of motion at an angle that changes in relation to the gravitational field. As the gravitational vector defined by (97) is not a function of the controlled coordinates, it also follows from (6) that the parameter k g1 is equal to zero.