\documentclass[10pt,openany]{book}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{euscript}
\usepackage{amssymb}
\usepackage{enumitem}
\usepackage{mdwlist}
\usepackage{float}
\usepackage{graphicx}
\usepackage[english,activeacute]{babel}
\decimalpoint
\usepackage{fancyhdr}
\usepackage{multicol}
\usepackage{vmargin}
\usepackage{mathrsfs}
\usepackage[titletoc]{appendix}

\renewcommand{\baselinestretch}{1.15}

\newtheorem{theorem}{Theorem}[section]
\newtheorem{assumption}{Assumption}[section]
\newtheorem{corollary}{Corollary}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{proposition}{Proposition}[section]

\newtheorem{definition}{Definition}[section]
\newtheorem{property}{Property}[section]
    \newenvironment{proof}{\textbf{Proof.}}{\hfill$\Box$\\}
\newtheorem{rem}{Remark}[section]
    \newenvironment{remark}{\begin{rem}\rm}{\hfill$\triangle$\end{rem}}
    
\begin{document}
\begin{titlepage}
\centering
{\includegraphics[width=0.7\textwidth]{IPICYT.png}\par}
\vspace{1cm}
{\bfseries\LARGE Instituto Potosino de Investigación Científica y Tecnológica \par}
\vspace{1cm}
{\scshape\Large Posgrado en Control y Sistemas Dinámicos \par}
\vspace{3cm}
{\scshape\LARGE PID-type continuous regulation with non-Lipschitz control actions for mechanical systems with bounded inputs \par}
\vspace{3cm}
{\large Tesis que presenta \par}

{\Large \textbf{Mariana Barrera Velázquez} \par}
\vspace{1cm}
{\large para obtener el grado de \par}
\vspace{0.5cm}
{\Large \textbf{Maestra en Control y Sistemas Dinámicos }}
\vfill
{\large Director de Tesis \par}
{\Large \textbf{Dr. Arturo Zavala Río} \par}
\vfill
\raggedleft
{San Luis Potosí, S. L. P. Agosto 2022 \par}
\end{titlepage}
\newpage
\begin{center}
\Large\textbf{Créditos Institucionales}
\end{center}

Esta tesis fue elaborada en la División de Control y Sistemas Dinámicos del Instituto Potosino de Investigación Científica y Tecnológica, A.C., bajo la dirección del Dr. Arturo Zavala Río.\\
\par
Durante la realización del trabajo el autor recibió una beca académica del Consejo Nacional de Ciencia y Tecnología con número de registro 1079858.

\newpage
\begin{center}
    \Large\textbf{Acknowledgments}\par
\end{center}
    
    I would like to thank my parents and family, for their trust and unconditional support, who always encouraged me to do my best.\par
    To Dr. Arturo Zavala, for his help, motivation, guidance and entirely dedication throughout the project.\par
    And finally, to IPICYT and CONACYT, for the grant student provided during the development of this thesis.
    \newpage

\begin{center}
    \Large\textbf{Abstract}
\end{center}
This work proposes a generalized Proportional-Integral-Derivative (PID) like continuous control law which solves the regulation problem for mechanical systems with constrained inputs. The proposal is inspired by finite-time stabilization schemes by permitting the proportional and derivative type actions to loose Lipschitz-continuity at their respective zero-error values. The resulting generalized PID type controller, whose implementation does not need to a priori know parameter values or the structure of the system model, turns out to give rise to closed-loop performance improvements, such as faster responses with smaller overshoot, when non-Lipschitz proportional and derivative type actions are involved. Global asymptotic stability is proven by means of Lyapunov's direct method and invariance theory, and an exhaustive description of the design requirements is explicitly presented. Stability results, performance improvements and input saturation avoidance are corroborated through simulations using a model of a 2-degree-of-freedom manipulator robot.
\par

\begin{center}
    \Large\textbf{Resumen}
\end{center}
Este trabajo propone una ley de control continua generalizada tipo Proporcional-Integral-Derivativo (PID) que resuelve el problema de regulación para sistemas mecánicos con entradas acotadas. La propuesta está inspirada en esquemas de estabilización en tiempo finito permitiendo que las acciones tipo proporcional y derivativa pierdan Lipschitz-continuidad cuando las respectivas variables de error son cero. El controlador tipo PID generalizado resultante, cuya implementación no necesita conocer previamente los valores de los parámetros o la estructura del modelo del sistema, resulta en mejoras en el desempeño en lazo cerrado, tales como respuestas más rápidas con menor sobretiro, cuando se involucran acciones tipo proporcional y derivativa que no son Lipschitz-continuas. Se demuestra estabilidad asintótica global por medio del método directo de Lyapunov y la teoría de invarianza, además de que se presenta una descripción exhaustiva de los requerimientos de diseño. Los resultados de estabilidad y las mejoras en desempeño se corroboran a través de simulaciones usando el modelo de un robot manipulador de 2 grados de libertad; donde además se corrobora que se evita la saturación en las entradas.
\tableofcontents
    \chapter{Introduction}

   This thesis deals with the position control problem of mechanical systems. Such a problem aims at the regulation of the system coordinates at desired (pre-defined) values, which shall be reached from any initial (position and velocity) conditions. Within the context of robot manipulators, such a regulation problem is described in \cite{ControlOfRobots} as that of driving the manipulator’s end-effector reach a desired position with a desired orientation,
    regardless of the initial posture. This approach deals only with constant  desired (generalized) positions, unlike the tracking problem, which involves a time-varying desired configuration trajectory.\par
    On the other hand, by considering the actual physical capabilities of the actuators used to move the mechanical system, constrained control inputs are taken into account. Moreover, the use of a control law without a saturation condition may lead to instability or performance degradation \cite{Kapasouris}, or in the case of non-linear actuator saturation, difficulties such as large overshoot, the appearance of limit cycle or an unstable output response  (as it is mentioned in \cite{Chen88}).  Therefore, the control objective shall be achieved avoiding input saturation bounds in order to provide a feasible solution to the considered regulation problem. \par
   
    \section{Previous works}
     With the goal of providing a new solution with its own benefits, it is worth considering the existing control schemes. Those considered in the context of the present thesis are Proportional-Integral-Derivative (PID) type controllers and finite-time regulators. As for the former case, a successful but quite complex way of solving the regulation problem is presented in \cite{Gorez}, and consists on a saturating PID-like control law for global regulation. In that article, the author is able to achieve the result either with velocity measurement or with its estimation by differentiation of the position measurement.
     \par
     Another solution is a recent approach presented in \cite{IJC15}, where a PID-type control is proposed. Departing from the dynamic model of mechanical systems, the authors present the PID-type control law along with its integral-action dynamics, which achieve the stabilization objective in the formulated analytical framework. In order to avoid the input bound, the PID-type control structure adopts a generalized form that incorporates generalized saturation functions. Some conditions on the proportional gain and on some parameters of the generalized saturation functions are obtained as a result of the analysis made.
    It is worth mentioning that given the generalized form of the control law, the work provides different cases for the election of its saturating structure. Therefore, the article covers the SP-SI-SD, SPD-SI, SPID, and SP-SID structures. Moreover, experiments on a manipulator robot with 3 degrees of freedom (DOF) are implemented, where the stabilization objective is corroborated to be achieved, by a set of tests involving the mentioned structures.\par
    
    Nevertheless, other authors suggest that the use of the integral action might cause trouble in the case of input constraints, as suggested in \cite{Krikelis}. This work proposes a design procedure for a PID-like control law, given a linear plant and a tracking problem with actuator saturation. Then, a general control law is generated from a plant model in state space representation, within the consideration of a control limiter (to avoid input saturation) and an \textit{intelligent integrator}. The last one is presented in order to improve the performance of the linear system that suffers from the effects of the so-called ``integrator wind-up'', which basically causes actuator saturation and large overshoot. Therefore a practical design of tracking systems is provided considering varying setpoint, disturbances and parameter changes.\par
    
    The passivity theory is also useful to solve the regulation problem, as shown in \cite{Meza05}. This article considers the existence of maximum torques in the dynamics of a mechanical system with friction, for which a saturated nonlinear PID control law is proposed.  This way, considering that the closed loop dynamics can be rewritten as a two-block interconnected feedback system, global asymptotic stabilization towards the desired equilibrium is proven to be achieved by means of passivity conditions on each block. Furthermore, the actuator constraints are avoided, regardless of initial conditions, which is corroborated by simulation results.
    
    Another approach including a PID-like controller is presented in \cite{Su10}, which is also addressed to mechanical systems with constrained inputs. The control law includes saturation functions (more specifically, the hyperbolic tangent function, $\tanh(\cdot)$), in each one of the three actions, in order to guarantee that the actuators avoid to reach their physical saturation level. In this case, the advantage of the proposed control law is the possibility to get the result either with the case of velocity measurement or the case when this value is dynamically approximated. Therefore, two analyses are provided in which, by following the Lyapunov direct method along with the invariance theory, the global asymptotic stability is obtained. The efficiency of the proposed control against other approaches is shown by two examples with and without velocity measurement, where the proposed approach is shown to be faster than the rest (although this is achieved by forcing a specific initial value on the auxiliary variable related to the integral action). Nevertheless, not all the tuning conditions are exhaustively derived, and some of them are left implicit, rendering difficult their application or verification which, in the best of the cases, shall be carried out numerically.\par
    
    On the other hand, other schemes include the finite-time stabilization of the origin as the main objective. Briefly, such a stabilization objective is achieved when the origin of the closed-loop system is rendered Lyapunov stable and finite-time attractive. A formal definition can be found in \cite{TAC14}.\par
    
    Finite-time stabilization is achieved, for instance, in the work presented in \cite{IET19} involving the notion of local homogeneity. Departing from the dynamic model of mechanical systems without damping effects, this work designs an PD-type with desired gravity compensation control law in order to achieve stabilization of the origin. Such a controller does not only include the saturating structure that involves both the proportional and derivative type actions within a single saturation function (at each link), but also that where each one is subjected to its own saturation function. In this case, the involved functions are strictly increasing, strictly passive, and strongly passive ones. From the definitions of such type of functions (which can be recalled from \cite{IET19}), specifically in the strongly passive case, it can be noted that there is an exponential weight involved (the parameter denoted by $\beta$ in the experimental results section), which, when less than unity, generates an infinite slope on the function around the origin. This effect, which is illustrated in the article through examples, is used in order to obtain finite-time stabilization.
    Furthermore, exponential stabilization is achievable using the same scheme with exponential weights equal to 1, so a comparison of both types of stabilization can be made. Even though this control scheme solves the problem in a satisfactory way, one of the disadvantages is the dependence on the exact parameters of the model through the inclusion of the desired conservative open-loop terms in the control law, $g(q_d)$; whereby, an error in those values would prevent the closed-loop system to achieve exact stabilization at the origin. The authors also mention that the closed-loop analysis results in a larger amount of requirements in the election of parameters, with respect to the online gravity compensation scheme.
    \par
   As for the finite-time tracking problem, a solution is given in \cite{TAC20}. The control law, in this case, corresponds to an SP-SD type scheme, where the saturation functions are bounded strongly passive ones. The functions involved in the SP and SD type actions also use exponential weights (denoted in \cite{TAC20} as $a_1$ and $a_2$) taking values in $(0,1]$. As in the previous mentioned work, the election of those parameters define the type of stabilization achieved (finite-time or exponential). Let us note that, the finite-time tracking result is proven in a different way than in \cite{IET19}, because the developed analysis does not involve local-homogeneity on the closed-loop dynamics (vector field). The authors are able to get conditions on the bounds of the defined functions in the proportional and derivative actions in order to avoid input saturation. As for the gain matrices, both related to the P and D type actions are independent of the latter conditions. The controller is tested by simulation, where it can be seen some advantages of the proposed approach compared to other schemes, like the less effort required by the actuators in order to achieve the stabilization.\par
    
    
    Both works involving finite-time stability (\cite{IET19} and \cite{TAC20}) show the advantages of the use of functions with fractional exponential weights, like: less overshoot, less effort from the control signal to get the stabilization, and a faster convergence time. Furthermore, robustness is achieved through finite-time controllers in the case of bounded perturbations, as it is presented in \cite{RNC20}. This paper presents a revisited closed-loop analysis involving the scheme proposed in \cite{TAC20} under perturbation. After a robustness analysis, the main result shows that whenever a sufficiently small perturbation term is introduced, the finite-time approach (with fractional exponential weights) outperforms its exponential counterpart (with unitary exponential weights), by giving rise to closed-loop responses with lower post-transient variations.
    
    \section{Motivation and objective}
    Based on the works presented above, one can see that there are various schemes that are able to solve the regulation problem, either using a PID-like scheme or a PD with (desired) conservative-force compensation. In fact, by analyzing the articles focusing on finite-time stabilization, the idea of achieving that result by using a PID-like control law results interesting. This was, in view of the nonexistence of such result, the original motivation of the thesis. That idea contemplated using that kind of control law to get finite-time stabilization by utilizing functions with fractional weights considering the condition of input constraints.\par
    
    Unfortunately, after our initial analyses in this direction, we realized that as the idea involves the use of functions with an infinite slope around the origin, it is necessary to know the value of stabilization in all the variables (including the one related to the integral action). Nevertheless, that idea opposes the general purpose of a PID controller, namely that such control law (particularly the integral action) must be able to compensate for the open-loop conservative force term (at equilibrium) independently of the desired position, and certainly without the need to a priori know any information on the structure or parameter values of those conservative forces.\par
    
    Therefore, the objective of this thesis had to be redefined, considering the benefits of the use of fractional exponential weights included in the functions to which the proportional and derivative actions would be subjected, as observed in the previous finite-time control approaches. Then, the motivation of the new objective is to design a more general PID type scheme for mechanical systems with bounded inputs, that incorporates exponential weights that may take fractional values in the P and D type actions, and explore the advantages that those definitions of functions can bring to the design in comparison to previous approaches, specially to the PID-like one, in which exponential weights are equal to unity. In fact, effects like the ones described in the previous section (benefits in the finite-time stabilization schemes), such as improved closed-loop performance, are expected to be corroborated through simulations.\par
    
    Thus, the objective of the thesis is to develop a research focused on designing a generalized PID-like control law that solves the regulation problem for mechanical systems with bounded inputs. The generalized scheme shall permit the use of exponential weights in the P and D type actions. Asymptotic and exponential stabilization are aimed to be achieved through less-than-or-equal-to-unity values of the incorporated exponential weights. The impact of such an incorporation is to be evaluated through simulation results.
    
   
    \section{Notation}
    Denoting $\mathbb{R}^n$ the set of $n$-dimensional vectors whose entries are real numbers, $0_n$ as its origin, and $\mathbb{R}^{m\times n}$ the set of $m\times n$-dimensional matrices with the same type of entries, let $x\in\mathbb{R}^n$ and $X, Y\in \mathbb{R}^{m\times n}$ (where $m$ and $n$ are positive integers). Thereby, throughout this document, $X_{ij}$ stands for the the element of $X$ at its $i$-th row and $j$-th column, $X_i$ refers to the $i$-th row of $X$, and $x_i$ denotes the $i$-th element of the vector $x$. In the case where $m=n$, $I_n$ stands for the $n\times n$ identity matrix, $X>0$, resp. $X\geq0$, symbolizes that $X$ is a positive definite, resp. semi-definite, matrix, and $X>Y$, resp. $X\geq Y$, indicates that $X-Y$ is a positive definite, resp. semi-definite, matrix. Moreover, $\lambda_m(X)$ and $\lambda_M(X)$ represent the minimum and maximum eigenvalue of a symmetric positive semi-definite matrix $X$, respectively.\par
    Consider also $\mathbb{R}_{>0}=\{x\in\mathbb{R}: x>0\}$, $\mathbb{R}_{\geq0}=\{x\in\mathbb{R}: x\geq0\}$, $\mathbb{R}_{>0}^n=\{x\in\mathbb{R}^n: x_i>0, \forall i=1,...,n\}$, and $\mathbb{R}_{\geq0}^n=\{x\in\mathbb{R}^n: x\geq0, \forall i=1,...,n\}$. For vectors, $||\cdot||$ denotes the standard Euclidean norm, i.e. $||x||=\sqrt{\sum_{i=0}^n|x_i|^2}$, while for matrices, the same notation is used for the induced 2-norm, i.e. $||X||=\sqrt{\lambda_M(X^TX)}$. More generally, $||\cdot||_p$ represents the $p$-norm for vectors with $p\geq1$, which is defined as $||x||_p=\left[\begin{aligned}
    \sum_{i=0}^n|x_i|^p
    \end{aligned}\right]^{1/p}$. 
    Furthermore, $\mathcal{B}_c^n$ symbolizes an $n$-dimensional ball with radius $c>0$, i.e. $\mathcal{B}_c^n=\{x\in\mathbb{R}^n: ||x||\leq c\}$, while $\mathcal{S}_c^{n-1}$ symbolizes an $(n-1)$-dimensional sphere with radius $c>0$, i.e. $\mathcal{S}_c^{n-1}=\{x\in\mathbb{R}^n: ||x||=c\}$.
    \par
Let $\mathcal{A}$ and $\mathcal{E}$ be subsets of the vector spaces $\mathbb{A}$ and $\mathbb{E}$ respectively. Then, for any $m\geq0$,  $\mathcal{C}^m(\mathcal{A},\mathcal{E})$ stands for the set of continuous functions from $\mathcal{A}$ to $\mathcal{E}$ being $m$ times continuously differentiable when $m>0$ (with differentiability at any point on the boundary of $\mathcal A$ meant as the limit from the interior of $\mathcal{A}$).  Furthermore, with a continuously differentiable scalar function $V\in\mathcal{C}^1(\mathbb{R}^n;\mathbb{R})$ and a vector function $f:\mathbb{R}^n\longrightarrow\mathbb{R}^n$, $D_fV$ denotes the directional derivative of $V$ along $f$, i.e. $D_f V(x) = \frac{\partial V}{\partial x}(x) f(x)$. In particular, if $f$ turns out to be the representation of a vector field, $\dot V$ will denote the derivative of $V$ along $f$, i.e. $\dot V(x) = \frac{\partial V}{\partial x}(x) f(x)$. Consider the next definitions
\begin{itemize}
    \item Sign function:
    \begin{equation*}
        \text{sign}(\varsigma)=\begin{cases}
         \frac{\varsigma}{|\varsigma|} & \text{if } \varsigma\neq0
             \\
             0 &\text{if } \varsigma=0
        \end{cases}
    \end{equation*}
    \item Unitary saturation (scalar) function
    \begin{equation*}
        \text{sat}(\varsigma)=\text{sign}(\varsigma)\min\{|\varsigma|,1\}
    \end{equation*}
\end{itemize}
    
Other facts that are used through the document are \cite{FTTOFRNC09}:
\begin{itemize}
    \item Young's inequality\\
    \begin{equation}
        \forall \phi,\psi\in(1,\infty) \text{ such that } \frac{1}{\phi}+\frac{1}{\psi}=1 \text{ and } \forall a,b\in\mathbb{R}_{\geq0}: ab\leq\frac{a^\phi}{\phi}+ \frac{b^\psi}{\psi}
    \end{equation}
    \item Hölder inequality
    \begin{equation}
        \forall \phi,\psi\in[1,\infty) \text{ such that } \frac{1}{\phi}+\frac{1}{\psi}=1 \text{ and } \forall x,y\in\mathbb{R}^n: |x^Ty|\leq||x||_\phi||y||_\psi
    \end{equation}
\end{itemize}
    \begin{lemma}
    \cite{TAC20} For any $x\in\mathbb{R}^n, ||x||_p$ is non-decreasing in $p$. In other words $\forall x\in\mathbb{R}^n, ||x||_\phi\geq||x||_\psi, \forall\phi,\psi$ such that $\phi\leq\psi$. 
    \label{NormaNoDecreciente}
    \end{lemma}
    
    \begin{remark}
    \cite{FTTOFRNC09} By equivalence of $p$-norms, for any $||\cdot||_{\phi}$ and $||\cdot||_\psi$, with $\phi\neq\psi$, and $\forall x\in\mathbb{R}^n$,  there are positive constants $\bar{c}_{\phi,\psi}>c_{\phi,\psi}$ such that:
    \begin{equation}
        c_{\phi,\psi}||x||_\psi\leq||x||_\phi\leq \bar{c}_{\phi,\psi}||x||_\psi
        \label{CotaNormaC}
    \end{equation}
    Indeed, by \ref{NormaNoDecreciente}, and considering that $||x||_\phi=[\sum_{i=1}^n|x_i|^\phi]^{1/\phi}\leq [\sum_{i=1}^n||x||_\psi^\phi]^{1/\phi}=n^{1/\phi}||x||_\psi$ one can see that $||x||_\psi\leq||x||_\phi\leq n^{1/\phi}||x||_\psi$ if $\phi\leq\psi$ and $n^{-1/\psi}||x||_\psi\leq||x||_\phi\leq||x||_\psi$ if $\phi\geq\psi$. Therefore, \eqref{CotaNormaC} is satisfied with:
    \begin{equation*}
        \begin{matrix}
             c_{\phi,\psi}=n^{[\text{sign}(\psi-\phi)-1]/2\psi}
            , &\bar{c}_{\phi,\psi}=n^{[\text{sign}(\psi-\phi)+1]/2\phi}
        \end{matrix}
    \end{equation*}
    \end{remark}


    \section{Mechanical systems dynamics}
     The dynamical model of mechanical systems is obtained through the same methodology followed to get that of robot manipulators. As it is explained in \cite{ControlOfRobots}, robot manipulators are \textit{``articulated mechanical systems composed of links
    connected by joints''}, where the joints can be either \textit{prismatic}
    or \textit{revolute}, or a combination of both types. As a matter of fact, robot manipulators formed by an open chain will be taken into account, as it is seen in Fig. \ref{ManipulatorScheme}
    \begin{figure}[H]
        \centering
        \includegraphics[width=.8
       \textwidth]{diagrama.png}
        \caption{Diagram of a n-DOF manipulator}
        \label{ManipulatorScheme}
    \end{figure}
    \par
    When a mathematical model of the robot is required, it is necessary to locate a 3-dimensional reference frame at the base of the robot, which is generally denoted as $[x_0,y_0,z_0]$. Then, the angular (revolute case), resp. translational (prismatic case), displacement of the \textit{j}- joint around, resp. along, the axis $z_{j}$, denoted as $q_j$, corresponds to the generalized joint coordinate. Notice that in most cases, the number of joints determines the number of \textit{degrees of freedom (DOF)} of the robot. Therefore, for robots with $n$-DOF the vector $q$ generally has $n$ elements.\par
    
    The dynamic model, as it consists in an ordinary differential equation, can be found using at least 2 methods: Newton's equations of motion or Lagrange's equations. However, as the first one becomes more complex when the number of DOF increases, the second method is generally considered. This method begins by obtaining the kinetic and potential energy of a manipulator robot with $n$-DOF (like the one in Fig. \ref{ManipulatorScheme}), whose sum is equal to the total energy of the robot, i.e.
    \begin{equation*}
        \mathcal{E}(q,\dot{q})=\mathcal{K}(q,\dot{q})+\mathcal{U}(q)
    \end{equation*} where
    $q=[q_1,...,q_n]^T$. Then, the \textit{Lagrangian} $\mathcal{L}(q,\dot{q})$ of the robot is calculated, which is equal to the difference between the kinetic and potential energies of the robot, i.e.
    \begin{equation*}
        \mathcal{L}(q,\dot{q})=\mathcal{K}(q,\dot{q})-\mathcal{U}(q)
    \end{equation*}
    Finally, the Lagrange equations of the robot 
are expressed as
    \begin{equation*}
        \frac{d}{dt}\left[\begin{aligned}
        \frac{\partial\mathcal{L}}{\partial \dot{q}}(q,\dot{q})
        \end{aligned}\right]- \frac{\partial\mathcal{L}}{\partial q}(q,\dot{q})=\tau
    \end{equation*}
    where $\tau$ represents the vector containing the external input forces at each joint. A more exhaustive explanation of the method can be found in \cite{ControlOfRobots}. In fact, the model of any $n$-DOF robot may be obtained in compact form using a method that involves a developed form of the Lagrangian equations. Let us mention that some cases include frictional effects, for which there are approximate models of friction forces, as the one that is involved in the present work.
    
    
    \section*{Dynamic model and properties}\par
       Consider the $n$-degree-of-freedom (DOF) fully actuated mechanical system dynamics with linear damping effects:
        \begin{equation}
            H(q)\Ddot{q}+C(q,\dot{q})\dot{q}+F\dot{q}+g(q)=\tau
            \label{Modelo}
        \end{equation}
    where $q, \dot{q},\Ddot{q}\in\mathbb{R}^{n}$ are the position, velocity and acceleration vectors, respectively. Furthermore, $H(q)\in\mathbb{R}^{n\times n}$ is the inertia matrix; $C(q,\dot{q})\in\mathbb{R}^{n\times n}$ is the Coriolis and centrifugal effect matrix defined through the Christoffel symbols of the first kind; $F \in \mathbb R^{n \times n}$ is the (a priori symmetric positive semi-definite) damping effect matrix; $g(q)=\nabla
    \mathscr{U}_{ol}(q)$, with $\mathscr{U}_{ol}:\mathbb{R}^n\to\mathbb{R}$ being the potential energy function of the open-loop system or equivalently: $\mathscr{U}_{ol}(q)=\mathscr{U}_{ol}(q_0)+\int_{q_0}^{q}g^T(z)dz$ for any $q,q_0\in\mathbb{R}^n$; and $\tau$ is the external input force vector. Recalling \cite[Chapter 4]{ControlOfRobots}, \cite{IJC15} and \cite{FTTOFRNC09}, some of the properties of the enlisted terms of the model are:
    \begin{property}
    $H(q)$ is a continuously differentiable symmetric matrix function satisfying $H(q)\geq\mu_mI_n$, which implies that $\|H(q)\| \geq \mu_m$, for some positive constant $\mu_m$, $\forall q\in\mathbb{R}^n$.
    \label{Cotas|H|}
    \end{property}
    \begin{property}
    The Coriolis and centrifugal effect matrix defined through the Christoffel symbols of the first kind satisfies:
        \begin{enumerate}
            \item $\dot{H}(x,y)=C(x,y)+C^T(x,y) $, and consequently $ z^T[\frac{1}{2}\dot{H}(x,y)-C(x,y)]z=0$, $\forall x,y,z\in\mathbb{R}^n$
            \item $C(x, y)z = C(x, z)y, \forall x, y, z \in\mathbb{R}^n.$
            \item $||C(x,y)||\leq\psi(x)||y||$, $\forall x,y\in\mathbb{R}^n$, for some $\psi:\mathbb{R}^n\longrightarrow\mathbb{R}_{\geq0}$
        \end{enumerate}
        \label{CoriolisProperty}
    \end{property}
    Moreover, consider the next assumptions, which according to \cite{ControlOfRobots} are valid, for instance, in the case of robot manipulator with only revolute joints.
    \begin{assumption}
    The inertia matrix is bounded, i.e.
    \begin{equation*}
        ||H(q)||\leq\mu_M
    \end{equation*} for a positive constant $\mu_M\geq\mu_m$, $\forall q\in\mathbb{R}^n$.
    \label{CotaH}
    \end{assumption}
    \begin{assumption}
    $\psi(x)$ in Property \ref{CoriolisProperty} is bounded by a non-negative constant $k_C$, and consequently $||C(x,y)||\leq k_C||y||$.
    \end{assumption}
    \begin{assumption}
     The conservative force vector $g(q)$ along with its Jacobian matrix $\frac{\partial
    g}{\partial q}(q)$ satisfy:
    \begin{enumerate}
        \item $|g_j(q)|<B_{g_j}$ for some $B_{gj} > 0$, $j = 1,\dots,n$, $\forall q\in\mathbb{R}^n$
        \item $||\frac{\partial
    g}{\partial q}(q)||\leq k_g \Longrightarrow ||g(x)-g(y)||\leq k_g||x-y||$, $\forall q,x,y\in\mathbb{R}^n$, for some $k_g\geq0$
    \end{enumerate}
    \label{Cotas|g|}
    \end{assumption}
    \begin{assumption}
     The damping effect matrix $F$ is symmetric positive definite, and consequently $f_m\|x\|^2 \leq x^TFx \leq f_M\|x\|^2$, $\forall x \in \mathbb R^n$,with $f_M \geq \lambda_{\max}(F) \geq \lambda_{\min}(F) \geq f_m > 0$.
    \end{assumption}

In this thesis, it is considered that the absolute value of each input $\tau_j$ is constrained to be smaller
than a given saturation bound $T_j>0$, i.e. $|\tau_j|\leq T_j, j=1,...,n$. In fact, by taking $u_j$ as the control variable relative to the $j$-th DOF we have:
\begin{equation}
    \tau_j=T_j\text{sat}(u_j/T_j)
    \label{Cond.tau}
\end{equation}
    Observe that from \eqref{Modelo} and \eqref{Cond.tau}, a necessary condition for the manipulator to be stabilisable at any desired equilibrium configuration $q_d\in\mathbb{R}^n$ is $T_j>B_{gj}, \forall j\in\{1,...,n\}$ \cite{IJC15}. 
    Thus, the next assumption is considered;
    \begin{assumption}
     $T_j>\alpha B_{gj}, \forall j\in\{1,...,n\}$, for some $\alpha\geq1$.
    \end{assumption}
    \section{Thesis structure }
    The rest of the document is divided and organized as follows.    Chapter 2 includes some definitions, Lemmas  and Theorems related to stability and the characteristics of the proposed controller. This theory becomes important specially in the stability analysis chapter.
    The third chapter presents the control law, as well as the design requirements. The conditions for (asymptotic and exponential) stabilization are given, along with the respective proofs.\par
    
    Simulation results are shown in chapter 4, where the stabilization objective is corroborated to be achieved and the closed-loop performance is evaluated, considering the model of a 2-DOF robot presented in \cite{IET19}. In order to carry out the implementations, a set of parameters are established as well as the type of functions that are involved in the control law.  Furthermore, several tests with the proposed scheme in different conditions are implemented, as well as comparison using alternative schemes. Finally, the conclusions are presented in Chapter 5.
    
    \chapter{Mathematical Background}
    This chapter provides some definitions, Lemmas and Theorems which are used throughout the document, related to Lipschitz continuity, Lyapunov stability, Invariance theory and certain types of functions. References \cite{Khalil} and \cite{FTTOFRNC09} support the contents presented here. 
    
    \section{Lipschitz continuity}
   Recalling \cite{Khalil} and \cite{Marquez}, Lipschitz continuity is a property of some functions, which one refers to either as continuous functions satisfying the Lipschitz condition or as Lipschitz-continuous functions. Formal definitions are stated next \cite{Khalil}.
   \begin{definition}
   A function $f : \mathbb R^n \to \mathbb R^m$ is said to be:
    \begin{itemize}
        \item \textit{locally Lipschitz} on a domain (open and  connected set) $D\subset\mathbb{R}^n$ if each point of $D$ has a neighborhood $D_0$ such that $f$ satisfies the Lipschitz condition:
        \begin{equation*}
            ||f(x)-f(y)||\leq L_0||x-y||
            \label{LipschitzCond}
        \end{equation*}
        for all points $x$ and $y$ in $D_0$ with some positive constant $L_0$;
        \item \textit{Lipschitz} on a set $W$, if it satisfies the Lipschitz condition
        \[
          \|f(x) - f(y)\| \leq L\|x-y\|
        \] 
        for all points $x$ and $y$ in $W$, with the same positive constant $L$, usually called the Lipschitz constant;
        \item \textit{globally Lipschitz} if it is Lipschitz on $\mathbb{R}^n$.
        \end{itemize}
   \end{definition}
   As a way to determine if a function fullfills the latter definition, consider the next Lemmas, whose proofs can be found in \cite{Khalil}.
   \begin{lemma}
   Let $f: D\longrightarrow\mathbb{R}^n$ be continuous for some domain $D\subset\mathbb{R}^n$. Suppose that $[\partial f/\partial x]$ exists and is continuous on $D$. If, for a convex subset $W\subset D$, there is a constant $L\geq 0$ such that
        \begin{equation*}
            \left\|\begin{aligned}
            \frac{\partial f}{\partial x}(x)
            \end{aligned}\right\|\leq L
        \end{equation*} on $W$, then
        \begin{equation*}
            \|f(x)-f(y)\|\leq L\|x-y\|
        \end{equation*} for every  $x$ and $y$ in $W$.
   \end{lemma}
   \begin{lemma}
   If $f(x)$ and $[\partial f/\partial x](x)$ are continuous on $D$, for some domain $D\subset\mathbb{R}^n$, then $f$ is \textit{locally Lipschitz} in $x$ on $D$.
   \end{lemma}
   \begin{lemma}
   If $f(x)$ and $[\partial f/\partial x](x)$ are continuous on $\mathbb{R}^n$, then $f$ is \textit{globally Lipschitz} in $x$ on $\mathbb{R}^n$ if and only if $[\partial f/\partial x]$ is \textit{uniformly bounded} on $\mathbb{R}^n$.
   \end{lemma}
   Versions of these definitions and Lemmas for $f : [0,\infty) \times \mathbb R^n \to \mathbb R^m$ are presented in \cite{Khalil}.
   
    \section{Lyapunov stability}
    The following definitions and theorems which refer to stability in the sense of Lyapunov, are taken from \cite{Bacciotti} and \cite{Khalil}, where also the proofs can be found. The analytical context of autonomous systems is considered due to the features of the system that will be taken into account in the following chapter.\par
    Consider the $n$-th order autonomous system:
    \begin{equation}
        \dot{x}=f(x)
        \label{AutSystem}
    \end{equation}
   where $f : D \to \mathbb R^n$ is a continuous vector field defined on a domain $D \subset \mathbb R^n$ containing the origin, which guarantees that for each $x_0 \in D$, there exists at least one (classical) solution $x : [0,\delta) \to D$, such that $x(0) = x_0$, for some $\delta \in (0,\infty]$. In fact, let $\mathcal{S}_{x_0}$ denote the set of all the solutions with $x(0)=x_0$. Let us further consider that the system has an equilibrium point $\bar x \in D$, and that this is taken to be at the origin, i.e. $\bar{x}=0_n$, considering that there is no loss of generality in doing so, due to the possibility of a change of variables.
   \begin{definition}
    The equilibrium point $x=0_n$ of \eqref{AutSystem} is:
    \begin{itemize}
        \item \textit{stable} if, for every $\varepsilon>0$, there is $\delta>0$ such that for each $x_0$ with $\|x_0\|<\delta$ and for all the solutions $x(\cdot)\in\mathcal{S}_{x_0}$: $x(t)$ exists for $t\in[0,\infty)$ and 
        \begin{equation*}
                ||x(t)||<\varepsilon \quad \forall t\geq0
        \end{equation*}
        \item \textit{unstable} if it is not stable
        \item\textit{asymptotically stable} if it is stable and if there exists $\delta_0>0$ such that for each $x_0$ satisfying $\|x_0\|<\delta_0$ and  for every $ x(\cdot)\in\mathcal{S}_{x_0}$:
        \begin{equation*}
           \lim_{t\longrightarrow\infty}\|x(t)\|=0
        \end{equation*}
        The origin is said to be \textit{globally asymptotically stable} if $\delta_0$ can be taken as large as desired.
    \end{itemize}
   \end{definition}
    Consider the next results as a way to know if the equilibrium point satisfies the last definition.
    \begin{theorem}
    Let $x=0_n$ be an equilibrium point of \eqref{AutSystem}, and $D\in\mathbb{R}^n$ be a domain containing the origin. Let $V:D\longrightarrow\mathbb{R}$ be a continuously differentiable function such that:
    \begin{equation}
        V(0_n)=0 \text{ and } V(x)>0 \text{ in } D\setminus\{0_n\}
    \end{equation}
    \begin{equation}
        \dot{V}(x)\leq0 \text{ in } D
    \end{equation}
    Then, $x=0_n$ is \textit{stable}. Moreover, if 
    \begin{equation}
        \dot{V}(x)<0 \text{ in } D\setminus\{0_n\}
    \end{equation}
    then $x=0_n$ is \textit{asymptotically stable}.
    \end{theorem}
    \begin{theorem}
    Let $x=0_n$ be an equilibrium point of \eqref{AutSystem}. Let $V:D\longrightarrow\mathbb{R}$ be a continuously differentiable function such that:
    \begin{equation}
        V(0_n)=0 \text{ and } V(x)>0\quad\forall x\neq0_n
    \end{equation}
    \begin{equation}
        ||x||\longrightarrow\infty \Longrightarrow V(x)\longrightarrow\infty
    \end{equation}
    \begin{equation}
        \dot{V}(x)<0 \quad \forall x\neq0_n
    \end{equation}
     then $x=0_n$ is \textit{globally asymptotically stable}.
    \end{theorem}
    
    \section{Invariance principle}
    The idea of \textit{LaSalle's invariance principle} appears from the situation in which asymptotic stability cannot be concluded due to a negative semi-definite derivative of the Lyapunov function. As it is explained in \cite{Khalil}, ``the argument shows, formally, that if in a domain about the origin we can find a Lyapunov function whose derivative along the trajectories of the system is negative semi-definite, and if we can establish that no trajectory can stay identically at points where $\dot{V}(x)=0$, except at the origin, then the origin is asymptotically stable.'' Versions of the invariance principle are given in \cite[Section 4.2]{Khalil}(under the consideration of a locally Lipschitz-continuous vector field in \eqref{AutSystem}) and [11, Section 7.2] (under the consideration of a continuous vector field in \eqref{AutSystem}). A Corollary of the invariance principle that will be applied within the analytical context of this thesis is reproduced next from \cite[Corollary 7.2.1]{Michel}. We begin by giving a useful definition.
    \begin{definition}
     \textbf{Invariant set}\\
        Given \eqref{AutSystem} and $x(t)$ as a solution of the latter, a set $M$ is set to be an \textit{invariant set} with respect to \eqref{AutSystem} if
        \begin{equation}
            x(0)\in M\Longrightarrow x(t)\in M \quad \forall t\in\mathbb{R}
        \end{equation}
    \end{definition}
    \begin{corollary}
    Let $x=0_n$ be an equilibrium point for \eqref{AutSystem}. Assume that there exists a continuously differentiable positive definite radially unbounded function $V: \mathbb R^n\longrightarrow\mathbb{R}$ such that $\dot{V}(x)\leq0,  \forall x\in\mathbb{R}^n$. Suppose that the origin $x = 0_n$ is the only invariant subset of the set $Z = \{x \in \mathbb R^n : \dot V(x) = 0\}$. Then, the equilibrium $x = 0_n$ of \eqref{AutSystem} is \textit{globally asymptotically stable}.
    \label{LaSalle}
    \end{corollary}
    
    \section{Exponential Stability}
    Considering an unforced system (which may be autonomous or non-autonomous)
    \begin{equation}
        \dot{x}=f(t,x)
        \label{NonAutSystem}
    \end{equation}
    a specific case of asymptotic stability is the \textit{exponential stability}. We recall from \cite{Khalil} the next definition and the subsequent theorem.
    \begin{definition}
     The equilibrium point $x=0_n$ of \eqref{NonAutSystem} is exponentially stable if there exist positive constants $c, k$, and $\lambda$ such that
     \begin{equation}
         \|x(t)\|\leq k\|x(t_0)\|e^{-\lambda(t-t_0)} \quad\forall\|x(t_0)\|<c
         \label{CondExpStability}
     \end{equation}
     and globally exponentially stable if \eqref{CondExpStability} is satisfied for any initial state $x(t_0)$.
    \end{definition}
    Exponential stability can be verified by means of the next theorem.
    \begin{theorem}
       Let $x=0_n$ be an equilibrium point for \eqref{NonAutSystem} and $D\subset\mathbb{R}^n$ be a domain containing $x=0_n$. Let $V: [0,\infty)\times D\longrightarrow\mathbb{R}$ be a continuously differentiable function such that
    \begin{equation}
        k_1||x||^a\leq V(t,x)\leq k_2||x||^a
    \end{equation}
    \begin{equation}
        \frac{\partial V}{\partial t}+\frac{\partial V}{\partial x}f(t,x)\leq -k_3||x||^a
    \end{equation}
    $\forall t\geq0$ and $\forall x\in D$, where $k_1, k_2, k_3$ and $a$ are positive constants. Then, $x=0_n$ is \textit{exponentially stable}. If the assumptions hold globally, then $x=0_n$ is \textit{globally exponentially stable}.
    \label{ExponentialStab}
    \end{theorem}
   
    \section{Scalar functions with particular properties}
    From \cite{FTTSFgRNC03a}, we recall the next definition and Lemmas.
   \begin{definition}
    A continuous scalar function $\sigma:\mathbb{R}\longrightarrow\mathbb{R}$ will be said to:\\
\begin{enumerate}
    \item be \textit{strictly passive} if $\varsigma\sigma(\varsigma)>0, \forall \varsigma\neq0$
    \item be \textit{strongly passive}---for ($\kappa,a,b$)---if $\varsigma\sigma(\varsigma)>0, \forall\varsigma\neq0$, and satisfies $|\sigma(\varsigma)|\geq\kappa(\min\{|\varsigma|,b\})^a, \forall\varsigma\in\mathbb{R}$ for positive constants $\kappa, a$ and $b$.
    \item have a \textit{local ($\bar{\kappa}, \bar{a}, \bar{b}$)-growth restriction} ---or satisfy a \textit{local ($\bar{\kappa}. \bar{a},\bar{b}$)-growth condition}--- if $|\sigma(\varsigma)|\leq \bar{\kappa}|\varsigma|^{\bar{a}}, \forall|\varsigma|\leq \bar{b}$, for positive constants $\bar{\kappa}, \bar{a}$ and $\bar{b}$.
    \item \textit{have a bounded} 
    ($\bar{\kappa},\bar{a},\bar{b}$)-\textit{growth restriction} ---or satisfy a \textit{bounded ($\bar{\kappa}. \bar{a},\bar{b}$)-growth condition}---
    if $|\sigma(\varsigma)|\leq\bar{\kappa}(\min\{|\varsigma|,\bar{b}\})^{\bar{a}}, \forall\varsigma\in\mathbb{R}$, for positive constants $\bar{\kappa}, \bar{a}$ and $\bar{b}$.
    \item be \textit{bounded} --- by $M$--- if $|\sigma(\varsigma)|\leq M, \forall \varsigma\in\mathbb{R}$, for a positive constant $M$.
    \end{enumerate}
    \label{FunctionTypes}
   \end{definition}
    \begin{remark}
    Let $\sigma$ be a strictly passive function. Then, from Item 1 of Definition \ref{FunctionTypes}, there exists positive constants $\kappa,b,a\geq \bar{a}, \bar{\kappa},\bar{b}$ such that
        \begin{equation*}
            \begin{matrix}
                 |\sigma(\varsigma)|\geq\kappa|\varsigma|^a \quad\forall|\varsigma|\leq b & 
                \text{and}& |\sigma(\varsigma)|\leq\bar{\kappa}|\varsigma|^{\bar{a}}\quad\forall|\varsigma|\leq \bar{b}
            \end{matrix}
        \end{equation*}
        With $\sigma$ being additionally nondecreasing, one can notice that:
        \begin{equation*}
            |\sigma(\varsigma)|\geq\kappa(\min\{|\varsigma|,b\})^a \quad \forall\varsigma\in\mathbb{R}
        \end{equation*}
        i.e. a nondecreasing strictly passive function is strongly passive for $\kappa,a,b$ and satisfies a local ($\bar{\kappa}, \bar{a}, \bar{b}$)-growth condition. Furthermore, if $\sigma(\varsigma)$ is Lipschitz continuous at $\varsigma=0$ then, $a\geq1$; while if $\min\{D^+\varsigma(0),D^-\varsigma(0)\}>0$ then, $\bar{a}\leq 1$.
        \label{sigma_StroPassive}
    \end{remark}
    \begin{remark}
    Consider a function $\sigma$ with a local ($\bar{\kappa}, \bar{a}, \bar{b}$)-growth condition and being bounded by $M$ i.e. $|\sigma(\varsigma)|\leq M, \forall\varsigma\in\mathbb{R}$, then 
        \begin{equation*}
            |\sigma(\varsigma)|\leq\bar{\kappa}|\varsigma|^{\bar{a}}\leq\max\{\bar{\kappa},M/\bar{b}^{\bar{a}}\}|\varsigma|^{\bar{a}}\leq \Bar{M}, \forall|\varsigma|\leq\bar{b}
        \end{equation*}
        with $\Bar{M}=\max\{\bar{\kappa},M/\bar{b}^{\bar{a}}\}\bar{b}^{\bar{a}}=\max\{\bar{\kappa}\bar{b}^{\bar{a}},M\}\geq M$ and in consequence
        \begin{equation*}
            |\sigma(\varsigma)|\leq\min\left\{\max\{\bar{\kappa},M/\bar{b}^{\bar{a}}\}|\varsigma|^{\bar{a}},\max\{\bar{\kappa},M/\bar{b}^{\bar{a}}\}\bar{b}^{\bar{a}}\right\}=\max\{\bar{\kappa},M/\bar{b}^{\bar{a}}\}(\min\{|\varsigma|,\bar{b}\})^{\bar{a}}, \forall\varsigma\in\mathbb{R}
        \end{equation*}
        which proves that $\sigma$ satisfies a bounded ($\max\{\bar{\kappa},M/\bar{b}^{\bar{a}}\},\bar{a}, \bar{b}$)-growth condition.
        \label{sigma_BoundedRestr}
    \end{remark}
    \begin{remark}
    For an increasing continuous scalar function $\sigma_0$, the function $\sigma(\varsigma_1,\varsigma_2)=\sigma_0(\varsigma_1+\varsigma_2)-\sigma_0(\varsigma_1)$ is increasing strictly passive with respect to $\varsigma_2$, uniformly in $\varsigma_1$. Indeed, from the increasing character of $\sigma_0$, for any $\varsigma_1$ one can see that
        \begin{equation*}
            \sigma(\varsigma_1,\varsigma_2)=\sigma_0(\varsigma_1+\varsigma_2)-\sigma_0(\varsigma_1)>0 \Longleftrightarrow \varsigma_2>0
        \end{equation*} and 
        \begin{equation*}
             \sigma(\varsigma_1,\varsigma_2)=\sigma_0(\varsigma_1+\varsigma_2)-\sigma_0(\varsigma_1)<0 \Longleftrightarrow \varsigma_2<0
        \end{equation*}
        Consequently, $\varsigma_2 \sigma(\varsigma_1,\varsigma_2)=\varsigma_2[\sigma_0(\varsigma_1+\varsigma_2)-\sigma_0(\varsigma_1)]>0, \forall\varsigma_2\neq0, \forall\varsigma_1\in\mathbb{R}$. Therefore, the enunciated result is obtained.
        \label{sigma_FunctionSum}
    \end{remark}
    
    An equivalent version of the next lemma was stated and proven in \cite[Lemma 4]{FTTOFRNC09}.
    \begin{lemma}
    For every $j\in\{1,...,n\}$, let $\sigma_j$ be a strongly passive function for ($\kappa,a,b$), $k_j$ be a positive constant, $k_m=\min_j\{k_j\}$, $k_M=\max_j\{k_j\}$ and, for any $x\in\mathbb{R}^n$ and $c>0$, let
    \begin{equation}
        S_0(x;a,c)=\begin{cases}
             ||x||^{1+a} & \text{ if } ||x||\leq c\\ 
             c^a||x|| & \text{ if } ||x||> c
        \end{cases}
    \end{equation}
    Then, letting $\varpi_a=n^{[\text{sign}(1-a)-1](1+a)/4}$, we have that:
    \begin{equation}
    \int_0^{x}s^T(Kz)dz=\sum_{j=1}^n\int_0^{x_{j}}\sigma_{j}(k_{j}z_j)dz_j\geq\frac{\kappa k_{m}^{a} \varpi_a}{1+a}S_0(x;a,b/k_{M})
    \label{CotaIntSigma}
\end{equation}
\begin{equation}
    \sum_{j=1}^nx_{j}\sigma_{j}(k_{j}x_{j})\geq\kappa k_m^a \varpi_aS_0(x;a,b/k_M), \forall x\in\mathbb{R}^n
    \label{CotaSumXSigma}
\end{equation}
    If $\forall j\in\{1,...,n\}$, $\sigma_j$ satisfies  a bounded ($\bar{\kappa},\bar{a},\bar{b}$)---growth condition then, in addition to the items above:
    \begin{equation}
    \sum_{j=1}^n\int_0^{x_{j}}\sigma_{j}(k_{j}z_j)dz_j\leq\bar{\kappa} k_{M}^{\bar{a}}nS_0(x;\bar{a},\bar{b}/k_{M})
    \label{Cota2IntSigma}
\end{equation}
and
\begin{equation}
    \sum_{j=1}^nx_{j}\sigma_{j}(k_{j}x_{j})\leq\bar{\kappa} k_{M}^{\bar{a}}nS_0(x;\bar{a},\bar{b}/k_{M})
    \label{Cota2SumXSigma}
\end{equation}
    \label{CotasSigmas}
    \end{lemma}
    \begin{remark}
    Note that if $a\leq1$ then---from the expression defining $\varpi$ (for $a<1$) and the proof of Lemma \ref{CotasSigmas} (for $a=1$)---$\varpi$ in expressions \eqref{CotaIntSigma} and \eqref{CotaSumXSigma} of Lemma \ref{CotasSigmas} is equal to unity.
    \end{remark}
    \chapter{Proposed Control scheme}
    This chapter presents the proposed control law with the design requirements in accordance to the formulated objective. The control law corresponds to a PID-like scheme, which involves a saturation function in each one of the control actions. 
        Consider the control law along with its auxiliary dynamics:
        \begin{equation}
    u(\Bar{q},\dot{q},\varphi)=-s_1(K_1\Bar{q})-s_2(K_2\dot{q})+s_3(K_3\varphi)
\end{equation}
    \begin{equation}
       \dot{\varphi}=-\dot{q}-\varepsilon_1\rho(\Bar{q})
       \label{AuxDynamics}
    \end{equation}
where:
\begin{itemize}
    \item $\Bar{q}=q-q_d$, with $q_d$ standing for the desired position vector
    \item $s_1(x)=[\sigma_{1_1}(x_1), ..., \sigma_{1_n}(x_n)]^T$
    \item $s_2(x)=[\sigma_{2_1}(x_1), ..., \sigma_{2_n}(x_n)]^T$
    \item $s_3(x)=[\sigma_{3_1}(x_1), ..., \sigma_{3_n}(x_n)]^T$
\end{itemize}\\
Throughout the rest of the chapter, we denote $k_{im} = \min_j\{k_{ij}\}$,  $k_{iM}=\max_{j}\left\{k_{ij}\right\}$.
\\[3ex]
 \textbf{Design requirements}\\
$K_1, K_2 $ and $K_3$ are positive definite diagonal matrices, $\varepsilon_1$ is a small enough positive constant, and for every $j\in\{1,..,n\}$ 
    $\sigma_{ij}$, $i=1,2$, are strongly passive functions, for ($\kappa_i$, $a_i$,  $b_i$), satisfying a bounded ($\Bar{\kappa}_i$, $a_i$,  $b_i$)- growth restriction
   such that $0<a_1\leq a_2\leq 1$, both ($i=1,2$) being locally Lipschitz-continuous on $\mathbb R \setminus \{0\}$;
    $\sigma_{3j}$ being bounded Lipschitz-continuous strictly increasing strictly passive functions, and all of them ($i = 1,2,3$) being such that
\begin{equation}
    \sup_{\forall(\varsigma_1,\varsigma_2,\varsigma_3)\in\mathbb{R}^{3}}\left|\begin{aligned}
    \sigma_{1j}(\varsigma_1)+\sigma_{2j}(\varsigma_2)+\sigma_{3j}(\varsigma_3)
    \end{aligned}\right|\triangleq B_j<T_j \quad \forall j\in\{1,...,n\}
    \label{Cond.u<T}
\end{equation}
     Furthermore, to ensure that the input is going to be able to compensate the conservative force steady-state value, it is required
\begin{equation}
    \min\left\{\begin{aligned}
\lim_{\varsigma\to\infty}\sigma_{3j}(\varsigma),-\lim_{\varsigma\to-\infty}\sigma_{3j}(\varsigma)
\end{aligned}\right\}>B_{gj} \quad \forall j\in\{1,...,n\}
\label{Cond.s3>Bg}
\end{equation}
Additionally,  for each $j\in\{1,...,n\}$, $k_{1j}$ and $\sigma_{1j}$ must be such that 
\begin{equation}
    |\sigma_{1j(k_{1j}\varsigma)}|> \min\{k_g|\varsigma|, 2B_{gj}\} \quad \forall \varsigma\neq 0
    \label{Orig.Cond.s1>min}
\end{equation}
with $k_g$ and $B_{gj}$ as defined through Assumption \ref{Cotas|g|}.\\
Moreover
 \begin{equation*}
     \rho(\bar{q})=h\Big(\bar{q},\frac{b_1}{k_{1M}}\Big)\bar{q}
 \end{equation*}
 where $b_1$ is the common parameter involved in the definition of $\sigma_{1j}$, $h\in C^0(\mathbb{R}^n\times\mathbb{R}_{>0};(0,1])$, being continuously differentiable on $\mathbb R^n \setminus 0_n$, uniformly in $\mathbb{R}_{>0}$, and such that $\forall c>0$, $\rho$ is a continuously
differentiable function satisfying:
\begin{equation}
    ||\rho(x)||=h(x;c)||x||\leq \min\{||x||,c \}, \forall x\in\mathbb{R}^n
    \label{NormaRho}
\end{equation}
and
\begin{equation}
    -h(x;c)<D_xh(x;c)<0, \forall x\neq 0_n
    \label{DerH}
\end{equation}
\begin{remark}
\cite{TAC20} By \eqref{DerH}, \textit{h} is decreasing on any radial direction, and consequently (as $h:\mathbb{R}^n\times\mathbb{R}_{>0}\longrightarrow(0,1]$) $h(x;c)\longrightarrow \omega$ as $||x||\longrightarrow\infty$ for some nonnegative constant $\omega$, whereas, on any compact connected neighborhood of the origin $\Upsilon\subset\mathbb{R}^n$, \textit{h} is lower bounded by a positive value $h_{m,\Upsilon}$ or more precisely: $1\geq h(0_n;c)\geq h(x;c)\geq \inf_{x\in\Upsilon}h(x;c)\triangleq h_{m,\Upsilon}=\inf_{x\in\partial\Upsilon}h(x;c)>\omega\geq0, \forall x\in\Upsilon$.
\label{Rem.h_1^*}
\end{remark}
    \begin{remark}
    As it is demonstrated in \cite{TAC20}, a family of functions which fullfill the properties required by $h$, is:
    \begin{equation*}
        h(x;c)=\frac{c}{[c^\varpi+\|x\|^\varpi]^{1/\varpi}} 
    \end{equation*}
    \label{h's}
    \end{remark}
    for any positive constant $\varpi$.
    \begin{remark}
    From \cite{TAC20}, and considering \eqref{DerH}:
     \begin{equation*}
         x^T\frac{\partial\rho}{\partial x}(x)x=x^T\left[\begin{aligned}
         h(x;c)I_n+x\frac{\partial h}{\partial x}(x;c)
         \end{aligned}\right]x=||x||^2[h(x;c)+D_xh(x;c)] \quad \forall x\neq0_n
     \end{equation*}
    whence, in view of (3.7), we have that $0 < x^T\frac{\partial \rho}{\partial x}(x)x < \|x\|^2$, $\forall x \neq 0$. Consequently, $0<\frac{\partial\rho}{\partial x}(x)\leq I_n$, which implies $||\frac{\partial\rho}{\partial x}(x)||\leq1$.
    \label{CotaJacobianoRho}
    \end{remark}

\par
For the rest of the analysis let:
\begin{equation*}
    h_1(\bar{q})\triangleq h\left(\bar{q},\frac{b_1}{k_{1M}}\right)
\end{equation*}
 \begin{remark}
 From \cite[Remark 12]{FTTOFRNC09} we have that for every $\nu\geq 1+a_1$:
 \begin{equation}
     ||\rho(\bar{q})||^\nu\leq \left(\frac{b_1}{k_{1M}}\right)^{\nu-1-a_1}h_1(\bar{q})S_1(\bar{q})\leq\left(\frac{b_1}{k_{1M}}\right)^{\nu-1-a_1}S_1(\bar{q})
     \label{NormaRhoNu}
 \end{equation}
 \end{remark}
\begin{remark}
\eqref{Orig.Cond.s1>min} implies the existence of positive constants $\hat{k}_{1j}>k_g$ and $b_j>2B_{gj}$ such that:
\begin{equation}
    |\sigma_{1j}(k_{1j}\varsigma)|\geq \min\{\hat{k}_{1j}|\varsigma|,b_j\} \quad \forall j\in\{1,...,n\}
    \label{Cond.s1>min}
\end{equation}
\label{Cond.s1>min2}
\end{remark}
By the satisfaction of \eqref{Cond.u<T}, the closed loop system becomes:
\begin{align*}
    H(q)\Ddot{q}+C(q,\dot{q})\dot{q}+F\dot{q}+g(q)&=-s_1(K_1\Bar{q})-s_2(K_2\dot{q})+s_3(K_3\varphi)
    \\
    \dot{\varphi}&=-\dot{q}-\varepsilon_1\rho(\Bar{q})
\end{align*}
Let $\Bar{\varphi}=\varphi-\varphi^*$, and $\Bar{s}_3(\Bar{\varphi})=s_3(K_3\Bar{\varphi}+K_3\varphi^*)-s_3(K_3\varphi^*)$ 
with $\varphi^*=[\varphi^*_1,...,\varphi^*_n]^T$ such that $s_3(K_3\varphi^*)=g(q_d)$, or equivalently, $\varphi_j^* = \sigma_{3j}^{-1}(g_j(q_d))/k_{3j}$, $j=1,\dots,n$ (the strictly increasing character of $\sigma_{3j}$ and \eqref{Cond.s3>Bg} ensure invertibility of $\sigma_{3j}$ on $[-B_{gj} , B_{gj}]$).
Thus, the closed-loop dynamics can be rewritten as:
\begin{equation}
    \begin{aligned}
        H(q)\Ddot{q}+C(q,\dot{q})\dot{q}+F\dot{q}+g(q)&=-s_1(K_1\Bar{q})-s_2(K_2\dot{q})+\Bar{s}_3(\Bar{\varphi})+g(q_d)
        \label{CLsystem}
        \\
        \dot{\bar\varphi}&=-\dot{q}-\varepsilon_1\rho(\Bar{q})
    \end{aligned}
\end{equation}
    \begin{remark}
    Notice that by the definition of $\sigma_{3j}(\varsigma)$, $j = 1 , \dots , n$, and Remarks \ref{sigma_StroPassive}, \ref{sigma_BoundedRestr} and \ref{sigma_FunctionSum}, there exist positive constants $\kappa_3$, $\bar\kappa_3$ and $b_3$ such that, in $\bar{s}_{3}(\sigma)=[\bar{\sigma}_{31}(\sigma_1),..., \bar{\sigma}_{3n}(\sigma_n))]$, $\bar{\sigma}_{3j}$ is a strongly passive function, for $(\kappa_3,1,b_3)$, satisfying a bounded ($\bar{\kappa}_3,1,b_3$)-growth condition.
    \label{PropBarSigma3}
    \end{remark}
   Let $x_1 = \bar q$, $x_2 = \dot q$, $x_3 = \bar\varphi$, and consider the consequent state-space representation of the closed-loop system, which takes the form:
   \begin{equation}
        \begin{aligned}
        \dot{x}_1&=x_2\\
        \dot{x}_2&=H^{-1}(x_1+q_d)[-C(x_1+q_d,x_2)x_2-Fx_2-g(x_1+qd)+g(q_d)-s_1(K_1x_1)-s_2(K_2x_2)+\Bar{s}_3(x_3)]\\
        \dot{x}_3&=-x_2-\varepsilon_1\rho(x_1)
    \end{aligned}
    \label{CLsystemSS}
   \end{equation}

    \begin{proposition}
            Consider the closed loop system \eqref{CLsystemSS} along with the design conditions \eqref{Cond.u<T},\eqref{Cond.s3>Bg}, and Remark \ref{Cond.s1>min2}. Then, $|\tau_j(t)|=|u_j(\bar q(t),\dot q(t),\varphi(t))|<T_j$ for all $j\in\{1,...,n\}$, $\forall t\geq0$, and the origin of the closed-loop system $(x_1,x_2,x_3)=(0_n,0_n,0_n)$ is globally asymptotically stable. In particular, the origin of the closed-loop system $(x_1,x_2,x_3)=(0_n,0_n,0_n)$ is additionally (locally) exponentially stable if $a_1 = a_2 = 1$.
    \end{proposition}
            
        \textbf{Proof}\par
        Considering \eqref{Cond.tau} and \eqref{Cond.u<T}, it is possible to notice that, along the system trajectories, \\
        $|u_j(\bar q(t),\dot q(t),\varphi(t))|<T_j$, thus $|u_j(\bar q(t),\dot q(t),\varphi(t))|=|\tau_j(t)|$, which leads to the first result. The rest of the proof is divided in 2 stages: global asymptotic stability and local exponential stability. 
    \paragraph{First stage: Global asymptotic stabilization\\} 
 Under the consideration of the direct Lyapunov's method, define the continuously
differentiable scalar function:
\begin{multline}
    V_1(x_1,x_2,x_3)=\frac{1}{2}x_2^TH(x_1+q_d)x_2 +\int_{0_n}^{x_1}s_1^T(K_1z)dz + \mathscr{U}(x_1)+\varepsilon_1\rho^T(x_1)H(x_1+q_d)x_2\\
    +\int_{0_n}^{x_3}\Bar{s}_3^T(z)dz
    \label{Lyapunov1}
\end{multline}
where
\begin{equation*}
    \int_{0_n}^{x_1}s_1^T(K_1z)dz=\sum_{j=1}^n\int_0^{x_{1j}}\sigma_{1j}(k_{1j}z_j)dz_j
\end{equation*}
\begin{equation*}
     \int_{0_n}^{x_3}\Bar{s}_3^T(z)dz=\sum_{j=1}^n\int_0^{x_{3j}}\Bar{\sigma}_{3j}(z_i)dz_i
 \end{equation*}
\begin{align*}
    \mathscr{U}(x_1)&=\mathscr{U}_{ol}(x_1+q_d)-\mathscr{U}_{ol}(q_d)-g^T(q_d)x_1\\
    &=\int_{0_n}^{x_1}\left[g(z+q_d)-g(q_d)\right]^Tdz\\
    &=\int_{0_n}^{x_1}\left[\int_{0_n}^z\frac{\partial g}{\partial q}(\Bar{z}+q_d)d\Bar{z}\right]^Tdz\\
\end{align*}
and $\varepsilon_1$ is a positive constant such that:
\begin{equation}
      \varepsilon_1<\min\{\varepsilon_{1M1},\varepsilon_{1M2},\varepsilon_{1M3}\}
      \label{Epsilon1Cotas's1}
\end{equation} where
\begin{equation*}
    \begin{matrix}
     \varepsilon_{1M1}=\begin{aligned}
         \frac{1}{\mu_M}\left(\begin{aligned}
     \frac{k_{1M}}{b_1}
    \end{aligned}\right)^{\frac{1-a_1}{2}}\sqrt{
     \frac{2\kappa_1k_{1m}^{a_1}\mu_m \gamma_{1M}}{1+a_1}}
     \end{aligned}\\
     \\
     \varepsilon_{1M2}=\begin{aligned}
         \frac{\gamma_{1M}\kappa_1k_{1m}^{a_1}f_m}{f_M^2(b_1/k_{1M})^{1-a_1}+2\gamma_{1M}\kappa_1k_{1m}^{a_1}(k_C\frac{b_1}{k_{1M}}+\mu_M)}
     \end{aligned}\\
      \\
      \varepsilon_{1M3}=\begin{aligned}
          \frac{\eta(1+a_2)\gamma_{1M}^{1/a_2}}{a_2 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}}\left(
       \frac{b_1}{k_{1M}}
      \right)^{\frac{a_1-a_2}{a_2} }\left(
     \frac{\kappa_1k_{1m}^{a_1}(1+a_2)}{2n^{\frac{2-a_2}{2}}\Bar{\kappa}_2k_{2M}^{a_2}}
    \right)^{1/a_2}
      \end{aligned}\\
    \end{matrix}
\end{equation*}
with 
\begin{equation}
    \gamma_{1M}&=1-\max\left\{
     \frac{k_g}{\hat{k}_{1m}},\frac{2B_{gM}}{b_m}
   \right\}
   \label{Gamma1M}
\end{equation}
   and
  \begin{equation}
      \eta \triangleq \min\left\{
\frac{\kappa_2 k_{2m}^{a_2}}{2},\frac{f_m}{2}\left(
 \frac{b_2}{k_{2M}} \right)^{1-a_2}
\right\}
\label{Eta}
  \end{equation}
(observe that in view of \eqref{Orig.Cond.s1>min} and Remark \ref{Cond.s1>min2},  $0<\gamma_{1M}<1$). Also consider the definitions
 \begin{equation*}
     S_0(x;a,c)\triangleq||x||(\min\{||x||,c\})^{a}
 \end{equation*} 
 \begin{equation*}
     S_1(x_1)\triangleq S_0\left(x_1;a_1,\frac{b_1}{k_{1M}}\right)
 \end{equation*}
 \begin{equation*}
     S_2(x_2)\triangleq S_0\left(x_2;a_2,\frac{b_2}{k_{2M}}\right)
 \end{equation*}
 and
 \begin{equation*}
     S_3(x_3)\triangleq S_0\left(x_3;1,\frac{b_3}{k_{3M}}\right)
 \end{equation*}\par
Notice that the scalar function $V_1$ can be rewritten as:
\begin{equation*}
    \begin{aligned}
         V_1(x_1,x_2,x_3)&=\frac{1}{2}x_2^TH(x_1+q_d)x_2
    +\gamma_1\int_{0_n}^{x_1}s_1^T(K_1z)dz+\varepsilon_1\rho^T(x_1)H(x_1+q_d)x_2+\int_{0_n}^{x_3}\Bar{s}_3^T(z)dz 
    \\
    &+(1-\gamma_1)\int_{0_n}^{x_1}s_1^T(K_1z)dz+ \mathscr{U}(x_1)
    \end{aligned}
\end{equation*}
with $\gamma_1$ satisfying:
\begin{equation}
    \max\{\gamma_{1m_1},\gamma_{1m_2},\gamma_{1m_3}\}\triangleq\gamma_{1m}<\gamma_1<\gamma_{1M}
    \label{CotasGamma1}
\end{equation}
\begin{align*}
      \gamma_{1m_1}&=\frac{(\varepsilon_1\mu_M)^2(1+a_1)}{2\kappa_1 k_{1m}^{a_1}\mu_m}\left(
     \frac{b_1}{k_{1M}}
    \right)^{1-a_1}\\
     \gamma_{1m_2}&=\frac{\varepsilon_1 f_M^2}{(\kappa_1k_{1m}^{a_1})(f_m-2\varepsilon_1(k_C\frac{b_1}{k_{1M}}+\mu_M))}\left(
       \frac{b_1}{k_{1M}}
      \right)^{1-a_1}
    \\\gamma_{1m_3}&=\frac{2n^{\frac{2-a_2}{2}}\Bar{\kappa}_2k_{2M}^{a_2}}{\kappa_1k_{1m}^{a_1}(1+a_2)}\left(
       \frac{b_1}{k_{1M}}
      \right)^{a_2-a_1}\left(
    \frac{\varepsilon_1 a_2 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}}{\eta(1+a_2)}
    \right)^{a_2}\\
\end{align*}
(one can see, from \eqref{Epsilon1Cotas's1}, that $\varepsilon_1 < \varepsilon_{1M\ell} \implies \gamma_{1m\ell} < \gamma_{1M}$, $\ell = 1,2,3$).
The terms of $V_1(x_1,x_2,x_3)$ are analyzed as follows:
\begin{enumerate}
    \item $\frac{1}{2}x_2^TH(x_1+q_d)x_2$\\
    Considering the lower bound property of the inertia matrix (Property \ref{Cotas|H|}):
    \begin{equation*}
        \frac{1}{2}x_2^TH(x_1+q_d)x_2\geq\frac{\mu_m}{2}||x_2||^2
    \end{equation*}
    \item $\gamma_1\int_{0_n}^{x_1}s_1^T(K_1z)dz$\\
    Recalling \eqref{CotaIntSigma} (in Lemma \ref{CotasSigmas})
    \begin{equation}
        \gamma_1\int_{0_n}^{x_1}s_1^T(K_1z)dz=\gamma_1\sum_{j=1}^n\int_0^{x_{1j}}\sigma_{1j}(k_{1j}z_j)dz_j\geq\frac{\gamma_1\kappa_1 k_{1m}^{a_1}}{1+a_1}S_1(x_1)
        \label{V1_c2}
    \end{equation}\par
    \item $\varepsilon_1\rho^T(x_1)H(x_1+q_d)x_2$\par
    Recalling \eqref{NormaRhoNu} for $\nu=2$:
    \begin{equation*}
        \begin{matrix}
         ||\rho(x_1)||^2\leq \left(\frac{b_1}{k_{1M}}\right)^{1-a_1}S_1(x_1)\Longleftrightarrow ||\rho(x_1)||\leq \left(\frac{b_1}{k_{1M}}\right)^{\frac{1-a_1}{2}}S_1^{1/2}(x_1)
         \\ \\
         \Longrightarrow
         \varepsilon_1\rho^T(x_1)H(x_1+q_d)x_2\geq-\varepsilon_1\mu_M\left(\frac{b_1}{k_{1M}}\right)^{\frac{1-a_1}{2}}S_1^{1/2}(x_1)||x_2||
        \end{matrix}
    \end{equation*}
    \item $\int_{0_n}^{x_3}\Bar{s}_3^T(z)dz$\\
Considering Lemma \ref{CotasSigmas} (specifically Eq.\eqref{CotaIntSigma}) and Remark \ref{PropBarSigma3}:
\begin{equation*}
    \int_{0_n}^{x_3}\Bar{s}_3^T(z)dz\geq\frac{\kappa_3 k_{3m}}{2}S_3(x_3)
\end{equation*}

\item $(1-\gamma_1)\int_{0_n}^{x_1}s_1^T(K_1z)dz+ \mathscr{U}(x_1)$\\
From \cite{IET19}
\begin{equation}
     \mathscr{U}(x_1)\leq\sum_{j=1}^n\int_0^{x_{1j}}k_gz_jdz_j
     \label{Cota1U(x_1)}
 \end{equation}
\begin{equation} \mathscr{U}(x_1)\leq\sum_{j=1}^n\int_0^{x_{1j}}\text{sign}(z_j)2B_{gj}dz_j
\label{Cota2U(x_1)}
 \end{equation}
Following a procedure analogous to that
shown in \cite{IET19}, under the consideration of Remark \ref{Cond.s1>min2}:
\begin{align*}
    (1-\gamma_1)\int_{0_n}^{x_1}s_1^T(K_1z)dz+ \mathscr{U}(x_1)&\geq \sum_{j=1}^n\int_0^{x_{1j}}\text{sign}(z_j)\min\{((1-\gamma_1)\hat{k}_{1j}-k_g)|z_j|,((1-\gamma_1)b_j-2B{gj})\}dz_j\\
    &\geq \sum_{j=1}^n\int_0^{x_{1j}}\text{sign}(z_j)\min\{\Bar{k}_{1j}|z_j|,\Bar{b}_{1j}\}dz_j =\sum_{j=1}^n\omega_{1j}(x_{1_j})
\end{align*}
where 
\begin{equation}
    \omega_{1j}(x_{1_j})=
   \begin{cases}
    \frac{\Bar{k}_{1j}}{2}x_{1j}^2 & \text{if } |x_{1j}|\leq\Bar{b}_{1j}/\Bar{k}_{1j} \\
    \Bar{b}_{1j}[|x_{1j}|-\Bar{b}_{1j}/2\Bar{k}_{1j}] & \text{if } |x_{1j}|>\Bar{b}_{1j}/\Bar{k}_{1j}
   \end{cases}
    \label{omega}
\end{equation}
and considering
\begin{equation*}
    0 <\Bar{k}_{1j}\leq(1-\gamma_1)\hat{k}_{1j}-k_g 
\end{equation*} and
\begin{equation*}
    0 <\Bar{b}_{1j}\leq (1-\gamma_1)b_j-2B{gj} 
\end{equation*}
(keeping in mind \eqref{Gamma1M} and the positive character of $\gamma_{1m}$ in \eqref{CotasGamma1}, notice that $\eqref{CotasGamma1} \implies [(1-\gamma_1)\hat k_{1j} - k_g > 0] \wedge [(1-\gamma_1)b_j - 2B_{gj} > 0]$). Recalling Lemma \ref{Lema2.3TAC20} which appears in \cite{TAC20}, and which is included in the Appendix \ref{ApendiceB}, notice that $\sum_{j=1}^n\omega_{1j}(x_{1_j})$ is analogous to $S(\varsigma)$ in the referred Lemma, i.e.
\begin{equation}
\omega_{1j}(x_{1j})=\hat{S}(x_{1j})=
    \begin{cases}
     \frac{\kappa k_j^a}{1+a}|x_{1j}|^2 & \forall |x_{1j}|\leq b_j/k_j \\
    \kappa b_j[|x_{1j}|-\frac{ab_j}{k_j(1+a)}] & \forall |x_{1j}|>b_j/k_j
    \end{cases}
    \label{S}
\end{equation}
with
$\kappa=1$, $a=1$,  $b_j=\Bar{b}_{1j}$ and $k_j=\Bar{k}_{1j}$.\par
Then, by Lemma 2.4 in \cite{TAC20} $\sum_{j=1}^n\omega_{1j}(x_{1_j})$ is lower bounded as follows:
\begin{equation}
    (1-\gamma_1)\int_{0_n}^{x_1}s_1^T(K_1z)dz+ \mathscr{U}(x_1)\geq \sum_{j=1}^n\hat{S}(x_{1j})\geq\frac{\Bar{k}_{1m}}{2}S_0\left(x_1;1,\frac{\Bar{b}_{1m}}{\Bar{k}_{1M}}\right)
    \label{Cota(1-gamma)w}
\end{equation}
where $\Bar{k}_{1m}=\min_j\{\Bar{k}_{1j}\}$, $\Bar{k}_{1M}=\max_j\{\Bar{k}_{1j}\}$ and      $\Bar{b}_{1m}=\min_j\{\Bar{b}_{1j}\}$. 
\end{enumerate}
Thus, from the bound gotten in items 1--5 above, $V_1(x_1,x_2,x_3)$ is lower bounded by:
\begin{align*}
     V_1(x_1,x_2,x_3)&\geq \frac{\mu_m}{2}||x_2||^2+\frac{\gamma_1\kappa_1 k_{1m}^{a_1}}{1+a_1}S_1(x_1)-\varepsilon_1\mu_M\left(\frac{b_1}{k_{1M}}\right)^{\frac{1-a_1}{2}}S_1^{1/2}(x_1)||x_2||+\frac{\kappa_3 k_{3m}}{2}S_3(x_3)\\
     &+\frac{\Bar{k}_{1m}}{2}S_0\left(x_1;1,\frac{\Bar{b}_{1m}}{\Bar{k}_{1M}}\right)
     \\
     &\geq\frac{\gamma_1\kappa_1 k_{1m}^{a_1}}{1+a_1}S_1(x_1)-\varepsilon_1\mu_M\left(\frac{b_1}{k_{1M}}\right)^{\frac{1-a_1}{2}}S_1^{1/2}(x_1)||x_2||+\frac{\mu_m}{2}||x_2||^2+\frac{\kappa_3 k_{3m}}{2}S_3(x_3)
\end{align*}
\begin{equation}
     =\frac{1}{2}\begin{pmatrix}
      S_1^{1/2}(x_1)\\
      ||x_2||
     \end{pmatrix}^TQ_{11}\begin{pmatrix}
      S_1^{1/2}(x_1)\\
      ||x_2||
     \end{pmatrix}+\frac{\kappa_3 k_{3m}}{2}S_3(x_3)\triangleq W_{11}(x_1,x_2,x_3)
     \label{Lyapunov1CotaInf}
\end{equation}
where 
\begin{equation}
    Q_{11}=\begin{pmatrix}
     \frac{2\gamma_1\kappa_1 k_{1m}^{a_1}}{1+a_1} & -\varepsilon_1\mu_M(\frac{b_1}{k_{1M}})^{\frac{1-a_1}{2}}\\
     -\varepsilon_1\mu_M(\frac{b_1}{k_{1M}})^{\frac{1-a_1}{2}} & \mu_m
    \end{pmatrix}
    \label{Q_11}
\end{equation}

and where it has been considered that $S_0$ in \eqref{Cota(1-gamma)w} is positive definite (with respect to $x_1$). Notice that $W_{11}$ in \eqref{Lyapunov1CotaInf} is positive definite if and only if $Q_{11}>0$, which is satisfied since $\gamma_1 > \gamma_{1m1} \implies Q_{11} > 0$ (see \eqref{Epsilon1Cotas's1}); moreover, $W_{11}$ in \eqref{Lyapunov1CotaInf} is radially unbounded, and consequently $V_{1}(x_1,x_2,x_3)$ is positive definite and radially unbounded.\par

Furthermore, the derivative of $V_1$ along the system trajectories is obtained as:
\begin{align*}
    \Dot{V}_1(x_1,x_2,x_3)&=x_2^TH(x_1+q_d)\dot{x}_2+\frac{1}{2}x_2^T\Dot{H}(x_1+q_d,x_2)x_2+[s_1(K_1x_1)+g(x_1+q_d)-g(q_d)]^T\dot{x}_1\\
      &+\varepsilon_1\rho^T(x_1)H(x_1+q_d)\dot{x}_2+\varepsilon_1\rho^T(x_1)\Dot{H}(x_1+q_d,x_2)x_2+\varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}\dot x_1\\
      &+\Bar{s}_3^T(x_3)\dot{x}_3\\ \\
      &=x_2^T[-(s_1(K_1x_1)+s_2(K_2x_2)-\Bar{s}_3(x_3))-C(x_1+q_d,x_2)x_2-Fx_2]\\
      &-x_2^T[g(x_1+q_d)-g(q_d)] +\varepsilon_1\rho^T(x_1)[-(s_1(K_1x_1)+s_2(K_2x_2)-\Bar{s}_3(x_3))]\\
      &+\varepsilon_1\rho^T(x_1)[-(C(x_1+q_d,x_2)x_2+Fx_2+g(x_1+q_d)-g(q_d))]+\frac{1}{2}x_2^T\Dot{H}(x_1+q_d,x_2)x_2\\
      &+[s_1(K_1x_1)+g(x_1+q_d)-g(q_d)]^Tx_2+\varepsilon_1\rho^T(x_1)\Dot{H}(x_1+q_d,x_2)x_2\\
      &+\varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}x_2+\Bar{s}_3^T(x_3)[-x_2-\varepsilon_1\rho(x_1)]\\ \\
      &=x_2^T[s_1(K_1x_1)-(s_1(K_1x_1)+s_2(K_2x_2)-\Bar{s}_3(x_3))]-x_2^TFx_2      -\varepsilon_1\rho^T(x_1)s_1(K_1x_1)\\
     &-\varepsilon_1\rho^T(x_1)[s_2(K_2x_2)-\Bar{s}_3(x_3)+g(x_1+q_d)-g(q_d)]-\varepsilon_1\rho^T(x_1)Fx_2\\&+\varepsilon_1 x_2^TC(x_1+q_d,x_2)\rho(x_1)+\varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}x_2+\Bar{s}_3^T(x_3)[-x_2-\varepsilon_1\rho(x_1)]
\end{align*}
\begin{equation}
    \begin{aligned}
     \Dot{V}_1(x_1,x_2,x_3)&=-x_2^Ts_2(K_2x_2)-x_2^TFx_2-\gamma_1\varepsilon_1\rho^T(x_1)s_1(K_1x_1)-\varepsilon_1\rho^T(x_1)s_2(K_2x_2)\\
      &-\varepsilon_1\rho^T(x_1)[(1-\gamma_1)s_1(K_1x_1)+(g(x_1+q_d)-g(q_d))]-\varepsilon_1\rho^T(x_1)Fx_2\\
      &+\varepsilon_1 x_2^TC(x_1+q_d,x_2)\rho(x_1)
      +\varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}x_2
          \label{V.}
\end{aligned}
\end{equation}

    where the system dynamics has been replaced and Property \ref{CoriolisProperty} has been considered.\par
As in the previous analysis, the terms of $\dot V_1$ are analyzed by considering a positive constant $\gamma$ satisfying
\begin{equation*}
    \gamma_m<\gamma<\gamma_M
\end{equation*}
\begin{equation*}
    \gamma_m=\begin{matrix}
      \frac{\varepsilon_1 a_2 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}}{\eta(1+a_2)}
    && \gamma_M=\left(
     \frac{\gamma_1\kappa_1k_{1m}^{a_1}(1+a_2)}{2n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}}
    \right)^{1/a_2}\left(
       \frac{b_1}{k_{1M}}
      \right)^{\frac{a_1-a_2}{a_2}}
    \end{matrix}
\end{equation*}
where $\hat{k}_{1m}\triangleq \min_j\{\hat{k}_{1j}\}$, $B_{gM}\triangleq \max_j\{B_{gj}\}$, $b_m\triangleq \min_j\{b_j\}$, and $\eta$ as defined in \eqref{Eta}. 
Then
\begin{enumerate}
    \item $-x_2^Ts_2(K_2x_2)-x_2^TFx_2$\\
    Considering the definition of the function $s_2(\cdot)$ whose elements are bounded strongly passive functions, from \eqref{CotaSumXSigma} (Lemma \ref{CotasSigmas}), notice that:
    \begin{equation*}
        -x_2^Ts_2(K_2x_2)=-\sum_{j=1}^nx_{2j}\sigma_{2j}(k_{2j}x_{2j})\leq-\kappa_2 k_{2m}^{a_2}S_0(x_2;a_2,b_2/k_{2M}
        )=-\kappa_2 k_{2m}^{a_2}S_2(x_2)
    \end{equation*}
On the other hand, by the properties of $F$:
     \begin{equation*}
 -x_2^TFx_2\leq-f_m||x_2||^2\\
    \end{equation*}
    Therefore:
    \begin{equation*}
        -x_2^Ts_2(K_2x_2)-x_2^TFx_2\leq -\kappa_2 k_{2m}^{a_2}S_2(x_2)-f_m||x_2||^2
    \end{equation*}
    Then, considering the definition and properties of $S_2(x_2)$ observe that for all $||x_2||\leq b_2/k_{2M}$:
\begin{align*}
 -\kappa_2 k_{2m}^{a_2}S_2(x_2)-f_m||x_2||^2&=-\kappa_2 k_{2m}^{a_2}||x_2||^{1+a_2}-f_m||x_2||^2\\
 &=-\frac{\kappa_2 k_{2m}^{a_2}}{2}||x_2||^{1+a_2}-\frac{\kappa_2 k_{2m}^{a_2}}{2}||x_2||^{a_2-1}||x_2||^2-f_m||x_2||^2\\
 &\leq -\frac{\kappa_2 k_{2m}^{a_2}}{2}||x_2||^{1+a_2}-\frac{\kappa_2 k_{2m}^{a_2}}{2}\left(
 \frac{b_2}{k_{2M}}
 \right)^{a_2-1}||x_2||^2-f_m||x_2||^2\\
 &=-\frac{\kappa_2 k_{2m}^{a_2}}{2}||x_2||^{1+a_2}-\left(
 \frac{\kappa_2 k_{2m}^{a_2}}{2}\left(
 \frac{b_2}{k_{2M}}
 \right)^{a_2-1}+f_m
 \right)||x_2||^2\\
  &\leq -\eta||x_2||^{1+a_2}-\frac{f_m}{2}||x_2||^2
\end{align*}
where $\eta$ is as defined in \eqref{Eta}; moreover, for all $||x_2||> b_2/k_{2M}$ note that:
\begin{align*}
     -\kappa_2k_{2m}^{a_2}S_2(x_2)-f_m||x_2||^2&\leq -f_m||x_2||^2\\
     &= -\frac{f_m}{2}||x_2||^{1-a_2}||x_2||^{1+a_2}-\frac{f_m}{2}||x_2||^2\\
     &\leq-\frac{f_m}{2}\left(\begin{aligned}
 \frac{b_2}{k_{2M}}
 \end{aligned}\right)^{1-a_2}||x_2||^{1+a_2}-\frac{f_m}{2}||x_2||^2\\
 &\leq -\eta||x_2||^{1+a_2}-\frac{f_m}{2}||x_2||^2
\end{align*}\par
whence
\begin{equation*}
   \begin{matrix}
     -\kappa_2 k_{2m}^{a_2}S_2(x_2)-f_m||x_2||^2\leq-\eta||x_2||^{1+a_2}-\frac{f_m}{2}||x_2||^2
     \\
     \Longrightarrow -x_2^Ts_2(K_2x_2)-x_2^TFx_2\leq-\eta||x_2||^{1+a_2}-\frac{f_m}{2}||x_2||^2
   \end{matrix}
\end{equation*} $\forall x_2\in\mathbb{R}^n$
     \item $-\gamma_1\varepsilon_1\rho^T(x_1)s_1(K_1x_1)$\\
    First of all, consider the definition of $\rho$, i.e.
    \begin{equation*}
        -\gamma_1\varepsilon_1\rho^T(x_1)s_1(K_1x_1)= -\gamma_1\varepsilon_1 h_1(x_1)x_1^Ts_1(K_1x_1)
    \end{equation*}
    Now, from \eqref{CotaSumXSigma} with $(\kappa_1,a_1,b_1)$:
    \begin{equation}
        \begin{matrix}
   \sum_{j=1}^nx_{1j}\sigma_{1j}(k_{1j}x_{1j})\geq \kappa_1k_{1m}^{a_1}S_0(x_1,a_1,b_1/k_{1M})\\
     \Longrightarrow -\gamma_1\varepsilon_1 h_1(x_1)x_1^Ts_1(K_1x_1)\leq -\gamma_1\varepsilon_1\kappa_1k_{1m}^{a_1} h_1(x_1)S_1(x_1)
 \end{matrix}
 \label{V1.c3}
    \end{equation}
    \item $-\varepsilon_1\rho^T(x_1)s_2(K_2x_2)$\\
     Considering Hölder inequality (with $\phi=\frac{2}{2-a_2}, \psi=\frac{2}{a_2}$), the definition of a strongly passive function, Young's inequality (with $\phi=1+a_2, \psi=\frac{1+a_2}{a_2}$), 
     $\gamma>0$ and the properties of $\rho(x_1)$ through \eqref{NormaRho} and \eqref{NormaRhoNu} with $\nu=1+a_2$, we have that:
     \begin{equation}
     \begin{aligned}[b]
     -\varepsilon_1\rho^T(x_1)s_2(K_2x_2)&\leq\varepsilon_1|\rho^T(x_1)s_2(K_2x_2)|\\
     &\leq \varepsilon_1||\rho(x_1)||_{\frac{2}{2-a_2}}||s_2(K_2x_2)||_{\frac{2}{a_2}}\\
     &\leq\varepsilon_1 n^{\frac{2-a_2}{2}}||\rho(x_1)||\Bar{\kappa}_2||K_2x_2||^{a_2}\\
     &\leq\varepsilon_1 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}(\gamma^{\frac{a_2}{1+a_2}}||\rho(x_1)||)(\gamma^{-\frac{a_2}{1+a_2}}||x_2||^{a_2})\\
     &\leq\varepsilon_1 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}\left(\frac{\gamma^{a_2}}{1+a_2}||\rho(x_1)||^{1+a_2}+\frac{a_2\gamma^{-1}}{1+a_2}||x_2||^{1+a_2}\right)\\
     &\leq\frac{\varepsilon_1 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}}{1+a_2}\left[\begin{aligned}
     \gamma^{a_2}\left(\begin{aligned}
     \frac{b_1}{k_{1M}}
     \end{aligned}\right)^{a_2-a_1}h_1(x_1)S_1(x_1)+a_2\gamma^{-1}||x_2||^{1+a_2}
     \end{aligned}\right]
     \label{V1.c4}
 \end{aligned}
 \end{equation}
     
  where the condition $a_1\leq a_2$ has been considered.
     \item $-\varepsilon_1\rho^T(x_1)[(1-\gamma_1)s_1(K_1x_1)+g(x_1+q_d)-g(q_d)]$\\
 Considering the definition of $\rho$ note that:
 \begin{equation*}
     \varepsilon_1\rho^T(x_1)[(1-\gamma_1)s_1(K_1x_1)+g(x_1+q_d)-g(q_d)]=\varepsilon_1 h_1(x_1)x_1^T[(1-\gamma_1)s_1(K_1x_1)+g(x_1+q_d)-g(q_d)]
 \end{equation*}
 Recalling \cite{IET19}, under the consideration of the properties of strictly passive functions (in accordance to Definition \ref{FunctionTypes}), and by taking $0 <\Bar{k}_{1j}\leq(1-\gamma_1)\hat{k}_{1j}-k_g$ and $0 <\Bar{b}_{1j}\leq(1-\gamma_1)b_j-2B_{gj}$:
 \begin{equation}
     \begin{aligned}[b]
       x_1^T[(1-\gamma_1)&s_1(K_1x_1)+(g(x_1+q_d)-g(q_d))]\\
      &=\sum_{j=1}^n|x_{1j}|[(1-\gamma_1)|\sigma_{1j}(k_{1j}x_{1j})|+\text{sign}(x_{1j})(g_j(x_1+q_d)-g_j(q_d))]\\
      &\geq\sum_{j=1}^n|x_{1j}|[(1-\gamma_1)|\sigma_{1j}(k_{1j}x_{1j})|-|g_j(x_1+q_d)-g_j(q_d)|]\\
  &\geq \sum_{j=1}^n|x_{1j}|\min\{((1-\gamma_1)\hat{k}_{1j}-k_g)|x_{1j}|,(1-\gamma_1)b_j-2B_{gj}\}\\
    &\geq \sum_{j=1}^n|x_{1j}|\min\{\Bar{k}_{1j}|x_{1j}|,\Bar{b}_{1j}\}>0
 \label{V1.c5}
 \end{aligned}
 \end{equation}
 
  where the right-hand-side inequality of \eqref{CotasGamma1} has been taken into account.
  Therefore:
  \begin{equation*}
      \begin{matrix}
        -\varepsilon_1\rho^T(x_1)[(1-\gamma_1)s_1(K_1x_1)+g(x_1+q_d)-g(q_d)]\leq-\varepsilon_1 \sum_{j=1}^nh_1(x_{1j})|x_{1j}|\min\{\Bar{k}_{1j}|x_{1j}|,\Bar{b}_{1j}\}
      \end{matrix}
  \end{equation*}
  Observe that the upper bound is negative definite with respect to $x_1$. Therefore, it is upper bounded by 0.
  \item $-\varepsilon_1\rho^T(x_1)Fx_2$\\
 From \eqref{NormaRhoNu} with $\nu=2$ and considering the properties of $F$ observe that:
 \begin{equation}
     \begin{aligned}[b]
     -\varepsilon_1\rho^T(x_1)Fx_2&\leq\varepsilon_1 f_M||\rho(x_1)||\cdot||x_2||\\
     &\leq \varepsilon_1 f_M\left[
      \begin{aligned}
      \left(\begin{aligned}
       \frac{b_1}{k_{1M}}
      \end{aligned}\right)^{1-a_1}h_1(x_1)S_1(x_1)
      \end{aligned}\right]^{1/2}||x_2||
       \label{V1.c6}
     \end{aligned}
 \end{equation}
     \item $\varepsilon_1 x_2^TC(x_1+q_d,x_2)\rho(x_1)
      +\varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}x_2$\\
      For the first term, taking into account the properties of $\rho$ (particularly that stated through \eqref{NormaRho}):
      \begin{align*}
     \varepsilon_1 x_2^TC(x_1+q_d,x_2)\rho(x_1)&\leq\varepsilon_1 k_C||\rho(x_1)||\cdot||x_2||^2\\
     &\leq\varepsilon_1 k_C\cdot \frac{b_1}{k_{1M}}||x_2||^2
    \end{align*}
For the second term, recalling Assumption \ref{CotaH} and Remark \ref{CotaJacobianoRho}:
\begin{align*}
     \varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}x_2&\leq\varepsilon_1\mu_M\left\|\frac{\partial \rho(x_1)}{\partial x_1}\right\|\cdot||x_2||^2\\
     &\leq \varepsilon_1\mu_M||x_2||^2
\end{align*}
Whence
\begin{equation}
    \varepsilon_1 x_2^TC(x_1+q_d,x_2)\rho(x_1)
      +\varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}x_2\leq \varepsilon_1\left[
      \begin{aligned}
      \mu_M+k_C\frac{b_1}{k_{1M}}
      \end{aligned}\right]||x_2||^2
      \label{V1.c7}
\end{equation}
\end{enumerate}
    
     Therefore, the derivative of $V_1$ is upper bounded as follows:
 \begin{align*}
     \Dot{V}_1(x_1,x_2,x_3)
     &=-x_2^Ts_2(K_2x_2)-x_2^TFx_2-\gamma_1\varepsilon_1\rho^T(x_1)s_1(K_1x_1)-\varepsilon_1\rho^T(x_1)s_2(K_2x_2)\\
      &-\varepsilon_1\rho^T(x_1)[(1-\gamma_1)s_1(K_1x_1)+(g(x_1+q_d)-g(q_d))]-\varepsilon_1\rho^T(x_1)Fx_2\\
      &+\varepsilon_1 x_2^TC(x_1+q_d,x_2)\rho(x_1)
      +\varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}x_2\\
      \\
     &\leq -\eta||x_2||^{1+a_2}-\frac{f_m}{2}||x_2||^2 -\gamma_1\varepsilon_1\kappa_1k_{1m}^{a_1} h_1(x_1)S_1(x_1)\\
     &+\frac{\varepsilon_1 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}}{1+a_2}\left[\begin{aligned}
     \gamma^{a_2}\left(\begin{aligned}
     \frac{b_1}{k_{1M}}
     \end{aligned}\right)^{a_2-a_1}h_1(x_1)S_1(x_1)+a_2\gamma^{-1}||x_2||^{1+a_2}
     \end{aligned}\right]\\
     &-\varepsilon_1\rho^T(x_1)[(1-\gamma_1)s_1(K_1x_1)+(g(x_1+q_d)-g(q_d))]\\
     &+\varepsilon_1 f_M\left[\left(
       \frac{b_1}{k_{1M}}
      \right)^{1-a_1}h_1(x_1)S_1(x_1)\right]^{1/2}||x_2||+\varepsilon_1\left[\mu_M+k_C\frac{b_1}{k_{1M}}\right]||x_2||^2\\ \\
     &\leq -\eta||x_2||^{1+a_2}-\frac{f_m}{2}||x_2||^2 -\gamma_1\varepsilon_1\kappa_1k_{1m}^{a_1} h_1(x_1)S_1(x_1)\\
     &+\frac{\varepsilon_1 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}}{1+a_2}\left[\begin{aligned}
     \gamma^{a_2}\left(\begin{aligned}
     \frac{b_1}{k_{1M}}
     \end{aligned}\right)^{a_2-a_1}h_1(x_1)S_1(x_1)+a_2\gamma^{-1}||x_2||^{1+a_2}
     \end{aligned}\right]
     \\
     &+\varepsilon_1 f_M\left[\left(\begin{aligned}
       \frac{b_1}{k_{1M}}
      \end{aligned}\right)^{1-a_1}h_1(x_1)S_1(x_1)\right]^{1/2}||x_2||+\varepsilon_1\left[\mu_M+k_C\frac{b_1}{k_{1M}}\right]||x_2||^2\\ \\
      &=-\varepsilon_1\left[\begin{aligned}
     \gamma_1\kappa_1k_{1m}^{a_1}-\frac{ n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}(b_1/k_{1M})^{a_2-a_1}}{1+a_2}\gamma^{a_2}
     \end{aligned}\right]h_1(x_1)S_1(x_1)\\
     &+\varepsilon_1 f_M\left(\begin{aligned}
       \frac{b_1}{k_{1M}}
      \end{aligned}\right)^{\frac{1-a_1}{2}}[h_1(x_1)S_1(x_1)]^{1/2}||x_2||\\
    &-\left(\frac{f_m}{2}-\varepsilon_1\left[\mu_M+k_C\frac{b_1}{k_{1M}}\right]\right)||x_2||^2-\left[\eta-\frac{\varepsilon_1 a_2 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}}{1+a_2}\gamma^{-1}\right]||x_2||^{1+a_2}
\end{align*}

    The last expression can be rewritten as
\begin{align*}
     \Dot{V}_1(x_1,x_2,x_3)&\leq -\frac{1}{2}\gamma_1\varepsilon_1\kappa_1k_{1m}^{a_1}h_1(x_1)S_1(x_1)+\varepsilon_1 f_M\left(
       \frac{b_1}{k_{1M}}
      \right)^{\frac{1-a_1}{2}}[h_1(x_1)S_1(x_1)]^{1/2}||x_2||\\
     &-\left[\frac{f_m}{2}-\varepsilon_1\left(\mu_M+k_C\frac{b_1}{k_{1M}}\right)\right]||x_2||^2
     \\
     &-\varepsilon_1\left[
     \frac{\gamma_1\kappa_1k_{1m}^{a_1}}{2}-\frac{ n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}(b_1/k_{1M})^{a_2-a_1}}{1+a_2}\gamma^{a_2}
     \right]h_1(x_1)S_1(x_1)\\
     &-\left[\eta-\frac{\varepsilon_1 a_2 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}}{1+a_2}\gamma^{-1}\right]||x_2||^{1+a_2}\\
     \end{align*}
     \begin{equation}
         \begin{aligned}[b]
         \dot{V}_1(x_1,x_2)&\leq-\frac{1}{2}\begin{pmatrix}
          [h_1(x_1)S_1(x_1)]^{1/2}\\
          ||x_2||
         \end{pmatrix}^TQ_{12}\begin{pmatrix}
          [h_1(x_1)S_1(x_1)]^{1/2}\\
          ||x_2||
         \end{pmatrix}\\
         &-\varepsilon_1p^{14}_{1}h_1(x_1)S_1(x_1)
         -p^{14}_2||x_2||^{1+a_2}
       \triangleq W_{14}(x_1,x_2)
         \end{aligned}
        \label{DLyapunov1CotaSup}
    \end{equation}
where \begin{equation}
    Q_{12}=\begin{pmatrix}
     \gamma_1\varepsilon_1\kappa_1k_{1m}^{a_1} & -\varepsilon_1 f_M\left(\begin{aligned}
       \frac{b_1}{k_{1M}}
      \end{aligned}\right)^{\frac{1-a_1}{2}}\\
     -\varepsilon_1 f_M\left(\begin{aligned}
       \frac{b_1}{k_{1M}}
      \end{aligned}\right)^{\frac{1-a_1}{2}} & f_m-2\varepsilon_1(\mu_M+k_C\frac{b_1}{k_{1M}})
    \end{pmatrix}
\end{equation}
\begin{equation*}
    p^{14}_1= \frac{\gamma_1\kappa_1k_{1m}^{a_1}}{2}-\frac{ n^{\frac{2-a_2}{2}}     \Bar{\kappa}_2k_{2M}^{a_2}(b_1/k_{1M})^{a_2-a_1}}{1+a_2}\gamma^{a_2}
\end{equation*}
\begin{equation*}
    p^{14}_2=\eta-\frac{\varepsilon_1 a_2 n^{\frac{2-a_2}{2}}
     \Bar{\kappa}_2k_{2M}^{a_2}}{1+a_2}\gamma^{-1}
    \end{equation*}
    Let us note that $\gamma_{1m2} <\gamma_1 < \gamma_{1M} \implies Q_{12} > 0
$. Moreover, one can corroborate that $ \gamma_{1m3}<\gamma_1<\gamma_{1M}$ and $ \gamma_m<\gamma<\gamma_M \Longrightarrow p_1^{14}>0$ and $p_2^{14}>0$, whence, $W_{14}$ in \eqref{DLyapunov1CotaSup} is concluded to be negative definite. Now, let:
\begin{equation*}
    \begin{matrix}
     \Omega\triangleq\{(x_1,x_2,x_3) \in \mathbb R^n\times\mathbb R^n\times\mathbb R^n: \dot{V}_1=0\}
     =\{ (x_1,x_2,x_3) \in \mathbb R^n\times\mathbb R^n\times\mathbb R^n : x_1 = x_2 = 0_n \}
    \end{matrix}
\end{equation*}
Moreover, notice that $x_1(t)\equiv x_2(t)\equiv 0_n \Longrightarrow \dot{x}_2(t)\equiv0_n$. Hence, for a solution to remain in $\Omega, \forall t\geq0$, it would be necessary that $x_1(t) \equiv x_2(t) \equiv \dot{x}_2(t) \equiv 0_n$ and by \eqref{CLsystem}:
\begin{equation*}
    H(0_n+q_d)0_n+C(0_n+q_d,0_n)0_n+F0_n+g(0_n+q_d)=-s_1(K_10_n)-s_2(K_20_n)+\Bar{s}_3(x_3)+g(q_d)
\end{equation*}
\begin{equation}
    \Longrightarrow \Bar{s}_3(x_3)=0_n \Longleftrightarrow x_3=0_n \quad \forall t\geq0
    \label{InvarianceP}
\end{equation}
Thus, $\{0_{3n}\}$ is the only invariant set in $\Omega$. Hence, since  $V_1(x_1,x_2,x_3)$ is positive definite and radially unbounded, $\dot{V}_1(x_1,x_2,x_3)$ is negative semi-definite, by the invariance theory (Corollary \ref{LaSalle}), we conclude that the origin of the closed-loop system \eqref{CLsystemSS} is globally asymptotically stable.
    \paragraph{2nd Stage: Exponential stabilization\\}
    Since the first stage of the proof includes the case $a_1=a_2=1$, global asymptotic stability of the trivial solution of the closed loop system is already proven; consequently, exponential stability is left to be proven. First consider a $3n$ dimensional ball $\mathcal{B}^{3n}_{\varrho}$ of radius $\varrho$ for any positive $\varrho\leq\min_{i=1,2,3}\left\{\begin{aligned}
\frac{b_i}{k_{iM}}
\end{aligned}\right\}$.\par

Notice that:
\begin{itemize}
    \item $(x_1^T \ x_2^T \ x_3^T)^T\in\mathcal{B}^{3n}_{\varrho}\Longrightarrow \max\{||x_1||,||x_2||,||x_3||\}\leq\varrho $
    \item $\forall i\in\{1,2,3\}, \varrho\leq\frac{b_i}{k_{iM}}\leq\frac{b_i}{k_{ij}} \quad\forall j\in\{1,\ldots,n\}$
\end{itemize}
Whence $(x_1^T \ x_2^T\ x_3^T)^T\in\mathcal{B}^{3n}_{\varrho}\Longrightarrow \forall j\in\{1,...,n\}$,
\begin{align*}
    |x_{1_j}|&\leq||x_1||\leq\frac{b_1}{k_{1M}}\leq\frac{b_1}{k_{1j}}, \\
 |x_{2_j}|&\leq||x_2||\leq\frac{b_2}{k_{2M}}\leq\frac{b_2}{k_{2j}}, \text{and}\\ 
 |x_{3_j}|&\leq||x_3||\leq\frac{b_3}{k_{3M}}\leq\frac{b_3}{k_{3j}}\\
\end{align*}
Consider also $a_1=a_2=1$ for the rest of the analysis.\par
Now, let:
\begin{equation*}
    V_2(x_1,x_2,x_3)=V_1(x_1,x_2,x_3)-\varepsilon_2x_3^TH(x_1+q_d)x_2
\end{equation*}
i.e. 
\begin{equation}
    \begin{aligned}
        V_2(x_1,x_2,x_3)&=\frac{1}{2}x_2^TH(x_1+q_d)x_2 +\int_{0_n}^{x_1}s_1^T(K_1z)dz + \mathscr{U}(x_1)+\varepsilon_1\rho^T(x_1)H(x_1+q_d)x_2\\
        &+\int_{0_n}^{x_3}\Bar{s}_3^T(z)dz
        -\varepsilon_2x_3^TH(x_1+q_d)x_2
     \label{Lyapunov2}
    \end{aligned}
\end{equation}

with $\varepsilon_1$ satisfying
\begin{equation}
    \varepsilon_1<\min\{\varepsilon_{1M1}, \varepsilon_{1M2}, \varepsilon_{1M3}, \varepsilon_{1M4},\varepsilon_{1M5}\}
    \label{Epsilon1Cotas's2}
    \end{equation} where $\varepsilon_{1M1}, \varepsilon_{1M2}$ and $\varepsilon_{1M3}$ are defined as in \eqref{Epsilon1Cotas's1}, and
\begin{align*}
        \varepsilon_{1M4}&=\frac{1}{\mu_M}\left[\begin{aligned}
        \frac{\kappa_1 k_{1m}\mu_m\gamma_{2M}}{2}
        \end{aligned}\right]^{1/2}\\
        \\
    \varepsilon_{1M5}&=\begin{aligned}
     \frac{\gamma_{2M}\kappa_1k_{1m}h_1^*(f_m+\kappa_2 k_{2m})}{2(f_M+\bar{\kappa}_2k_{2M})^2+4\gamma_{2M}\kappa_1k_{1m}h_1^*(k_c\varrho +\mu_M)}
     \end{aligned}
\end{align*}
with 
\[h_1^*\triangleq h_{m,\mathcal{B}^{n}_{\varrho}}=\inf_{\bar{q}\in\partial\mathcal{B}^{n}_{\varrho}}h_1(\bar{q})= \inf_{\|\bar q\| = \varrho} h_1(\bar q)


\] $\varepsilon_2$ being a positive constant such that:
\begin{equation}
    \varepsilon_2<\min\{\varepsilon_{2M1},\varepsilon_{2M2},\varepsilon_{2M3},\varepsilon_{2M4}\}
    \label{Epsilon2Cotas's}
\end{equation} where
\begin{align*}
         \varepsilon_{2M1}&=\frac{1}{\mu_M}\left[\begin{aligned}
       \frac{\kappa_3k_{3m}\mu_m}{2}
        \end{aligned}\right]^{1/2} &
        \varepsilon_{2M2}&=\frac{1}{2\mu_M}\left[\begin{aligned}
            \frac{\gamma_{2M}\kappa_1k_{1m}h_1^*d_m}{\varepsilon_1}
        \end{aligned}\right]^{1/2}\\
         \\\varepsilon_{2M3}&=\begin{aligned}
         \frac{\gamma_{2M}\varepsilon_1\kappa_1k_{1m}h_1^*\kappa_3k_{3m}d_m}{d_m(k_g+\bar{\kappa_1}k_{1M})^2+\gamma_{2M}\varepsilon_1\kappa_1k_{1m}h_1^*\theta^2}
         \end{aligned} &
         \varepsilon_{2M4}&=\frac{d_m}{4(k_c\varrho+\mu_M)}
\end{align*}
with 
\begin{subequations}
\begin{equation}
    \theta=f_M+\bar{\kappa}_2k_{2M}
    \label{Theta}
\end{equation}
\begin{equation}
    d_m=f_m+\kappa_2 k_{2m}
    \label{dm}
\end{equation}
\label{ThetayDm}
\end{subequations}and with $\gamma_2<1$ being a positive constant satisfying:
\begin{equation}
    \max\{\gamma_{2m1},\gamma_{2m2},\gamma_{2m3},\gamma_{2m4}\}<\gamma_2<\gamma_{2M}
    \label{Gama2Cotas}
\end{equation} where
\begin{equation*}
    \begin{matrix}
         \gamma_{2m1}=\begin{aligned}
             \frac{2\varepsilon_1^2\mu_M^2}{\kappa_1 k_{1m}\mu_m}
         \end{aligned} && \gamma_{2m2}=\begin{aligned}
             \frac{\varepsilon_1\theta^2}{\kappa_1k_{1m}h_1^*[\frac{d_m}{2}-2\varepsilon_1(k_c\varrho +\mu_M)]}
         \end{aligned}\\
         \\ 
         \gamma_{2m3}=\begin{aligned}
        \frac{4\varepsilon_1\varepsilon_2^2\mu_M^2}{\kappa_1k_{1m}h_1^*d_m}
        \end{aligned} && \gamma_{2m4}=\begin{aligned}
         \frac{\varepsilon_2d_m(k_g+\bar{\kappa_1}k_{1M})^2}{\varepsilon_1\kappa_1k_{1m}h_1^*[\kappa_3k_{3m}d_m-\varepsilon_2\theta^2]}
         \end{aligned}\\
         \\
        \gamma_{2M}= 1-\max\left\{\begin{aligned}
     \frac{k_g}{\hat{k}_{1m}},\frac{2B_{gM}}{b_m}
    \end{aligned}\right\}
    \end{matrix}
\end{equation*}
where it can be corroborated that $\varepsilon_{1}<\varepsilon_{1M4}$, $\quad \varepsilon_1<\varepsilon_{1M5}$, $\quad \varepsilon_2<\varepsilon_{2M2}$, and $\varepsilon_{2}<\varepsilon_{2M3}$ imply $\gamma_{2mj}<\gamma_{2M},j=1,...,4$, respectively. 
Consider also $0 <\Bar{k}_{1j}\leq(1-\gamma_2)\hat{k}_{1j}-k_g$ , $0 <\Bar{b}_{1j}\leq (1-\gamma_2)b_j-2B{gj}$.

    
    On the one hand, note that on $\mathcal{B}^{3n}_{\varrho}$:
   \begin{equation}
        \begin{aligned}[b]
        V_2(x_1,x_2,x_3)&\geq \frac{\mu_m}{2}||x_2||^2+ \frac{\gamma_2\kappa_1 k_{1m}}{2}||x_1||^2+\frac{\Bar{k}_{1m}}{2}S_0\left(x_1;1,\frac{\Bar{b}_{1m}}{\Bar{k}_{1M}}\right)-\varepsilon_1\mu_M||x_1||\cdot||x_2||\\
        &+\frac{\kappa_3k_{3m}}{2}||x_3||^2-\varepsilon_2\mu_M||x_2||\cdot||x_3||
        \\
        &\geq \frac{\mu_m}{2}||x_2||^2+ \frac{\gamma_2\kappa_1 k_{1m}}{2}||x_1||^2-\varepsilon_1\mu_M||x_1||\cdot||x_2||+\frac{\kappa_3k_{3m}}{2}||x_3||^2\\
        & -\varepsilon_2\mu_M||x_2||\cdot||x_3||\\
        &=
        \frac{1}{2}\begin{pmatrix}
        ||x_1||\\||x_2||
        \end{pmatrix}^TQ_{211}
        \begin{pmatrix}
        ||x_1||\\||x_2||
        \end{pmatrix}+\frac{1}{2}\begin{pmatrix}
        ||x_2||\\||x_3||
        \end{pmatrix}^TQ_{212}
        \begin{pmatrix}
        ||x_2||\\||x_3||
        \end{pmatrix}\triangleq W_{21}(x_1,x_2,x_3)
        \label{CotaInfLyapunov2}
\end{aligned}
   \end{equation}

where
\begin{equation}
    Q_{211}=\begin{pmatrix}
    \gamma_2\kappa_1 k_{1m} && -\varepsilon_1\mu_M\\
    -\varepsilon_1\mu_M && \frac{\mu_m}{2}
    \end{pmatrix}
\end{equation}
\begin{equation}
    Q_{212}=\begin{pmatrix}
    \frac{\mu_m}{2} && -\varepsilon_2\mu_M\\
    -\varepsilon_2\mu_M && \kappa_3k_{3m}
    \end{pmatrix}
\end{equation}
    Notice that $\gamma_{2m1}<\gamma_2<\gamma_{2M} \implies Q_{211}>0$. Moreover $\varepsilon_2<\varepsilon_{2M1}\implies Q_{212}>0$. Thus, $W_{21}$ in \eqref{CotaInfLyapunov2} is positive definite.
    \par
    On the other hand, recalling property \ref{Cotas|H|} (on the inertia matrix), \eqref{Cota2IntSigma}, \eqref{Cota1U(x_1)}, \eqref{Cota2U(x_1)}, \eqref{Cond.s1>min}, the properties of $\rho$, and Young's Inequality with $\phi=\psi=2$,  $V_2(x_1,x_2,x_3)$ is upper bounded in $\mathcal{B}^{3n}_\varrho$ as follows:
    \begin{equation}
    \begin{aligned}[b]
    V_2(x_1,x_2,x_3)&\leq \frac{1}{2}\mu_M||x_2||^2+\bar{\kappa}_1k_{1M}n||x_1||^2+\sum_{j=1}^n\int_{0}^{x_1j}\text{sign}(z_j)\min\left\{\begin{aligned}
    k_g|z_j|,2B_{gj}
    \end{aligned}\right\}\\
    &+\varepsilon_1\mu_m||\rho(x_1)||\cdot||x_2||+\bar{\kappa}_3k_{3M}n||x_3||^2+\varepsilon_2\mu_M||x_2||\cdot||x_3||\\
    &\leq \frac{1}{2}\mu_M||x_2||^2+\bar{\kappa}_1k_{1M}n||x_1||^2+\sum_{j=1}^n\int_{0}^{x_1j}\sigma_1(k_{1j}z_j)
    +\varepsilon_1\mu_m||x_1||\cdot||x_2||\\
    &+\bar{\kappa}_3k_{3M}n||x_3||^2+\varepsilon_2\mu_M||x_2||\cdot||x_3||
    \\
    &\leq \frac{1}{2}\mu_M||x_2||^2+2\bar{\kappa}_1k_{1M}n||x_1||^2
    +\frac{\varepsilon_1\mu_m}{2}[||x_1||^2+||x_2||^2]+\bar{\kappa}_3k_{3M}n||x_3||^2\\
    &+\frac{\varepsilon_2\mu_M}{2}[||x_2||^2+||x_3||^2]\\
    &= p_{1}^{22}||x_1||^2+p_{2}^{22}||x_2||^2+p_3^{22}||x_3||^2\triangleq W_{22}(x_1,x_2,x_3)
    \label{CotaSupLyapunov2}
    \end{aligned}
\end{equation}
where
\begin{align*}
    p_{1}^{22}&=2\bar{\kappa}_1k_{1M}n+\frac{\varepsilon_1\mu_m}{2}
    \\
    p_{2}^{22}&=\frac{1}{2}\mu_M+\frac{\varepsilon_1\mu_m}{2}+\frac{\varepsilon_2\mu_M}{2}
    \\
    p_{3}^{22}&=\bar{\kappa}_3k_{3M}n+\frac{\varepsilon_2\mu_M}{2}
\end{align*}
whence positive  definiteness of $W_{22}$ in \eqref{CotaSupLyapunov2} can be corroborated. Then, from the conclusions gotten for both $W_{21}$ and $W_{22}$, $V_2(x_1,x_2,x_3)$ is concluded to be positive definite.\par

    The derivative of $V_2(x_1,x_2,x_3)$ along the system trajectories is
    \begin{align*}
     \dot{V}_2(x_1,x_2,x_3)&=x_2^TH(x_1+q_d)\dot{x}_2+\frac{1}{2}x_2^T\Dot{H}(x_1+q_d,x_2)x_2+s_1^T(K_1x_1)\dot x_1+[g(x_1+q_d)-g(q_d)]^T\dot x_1\\
      &+\varepsilon_1\rho^T(x_1)H(x_1+q_d)\dot{x}_2+\varepsilon_1\rho^T(x_1)\Dot{H}(x_1+q_d,x_2)x_2+\varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}\dot{x}_1\\
      &+\Bar{s}_3^T(x_3)\dot{x}_3-\varepsilon_2x_3^TH(x_1+q_d)\dot{x}_2-\varepsilon_2x_3^T\dot{H}(x_1+q_d,x_2)x_2-\varepsilon_2x_2^TH(x_1+q_d)\dot{x}_3\\ \\
      &=x_2^TH(x_1+q_d)[-C(x_1+q_d,x_2)x_2-Fx_2-(g(x_1+q_d)-g(q_d))-s_1(K_1x_1)]\\
    &+x_2^TH(x_1+q_d)[-s_2(K_2x_2)+\bar{s}_3(x_3)]
     +\frac{1}{2}x_2^T\Dot{H}(x_1+q_d,x_2)x_2+s_1^T(K_1x_1)x_2\\
     &+[g(x_1+q_d)-g(q_d)]^Tx_2
      +\varepsilon_1\rho^T(x_1)H(x_1+q_d)[-C(x_1+q_d,x_2)x_2-Fx_2]\\
      &+\varepsilon_1\rho^T(x_1)H(x_1+q_d)[-(g(x_1+q_d)-g(q_d))-s_1(K_1x_1)-s_2(K_2x_2)+\bar{s}_3(x_3)]\\
      &+\varepsilon_1\rho^T(x_1)[C(x_1+q_d,x_2)+C^T(x_1+q_d,x_2)]x_2+\varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}x_2+\Bar{s}_3^T(x_3)\dot{x}_3\\
     &-\varepsilon_2x_3^T[-C(x_1+q_d,x_2)x_2-Fx_2-(g(x_1+q_d)-g(q_d))-s_1(K_1x_1)-s_2(K_2x_2)+\bar{s}_3(x_3)]\\
     &-\varepsilon_2x_3^T[C(x_1+q_d,x_2)+C^T(x_1+q_d,x_2)]x_2-\varepsilon_2x_2^TH(x_1+q_d)[-x_2-\varepsilon_1\rho(x_1)]\\ \\
      &=-x_2^TFx_2-x_2^Ts_2(K_2x_2)-\varepsilon_1\rho^T(x_1)Fx_2-\gamma_2\varepsilon_1\rho^T(x_1)s_1(K_1x_1)-\varepsilon_1\rho^T(x_1)s_2(K_2x_2)\\
      &-\varepsilon_1\rho^T(x_1)[(1-\gamma_2)s_1(K_1x_1)+(g(x_1+q_d)-g(q_d))]+\varepsilon_1 x_2^TC(x_1+q_d,x_2)\rho(x_1)\\
      &
      +\varepsilon_1 x_2^TH(x_1+q_d)\frac{\partial \rho(x_1)}{\partial x_1}x_2+\varepsilon_2x_3^TFx_2+\varepsilon_2x_3^T(g(x_1+q_d)-g(q_d))+\varepsilon_2x_3^Ts_1(K_1x_1)\\
      &+\varepsilon_2x_3^Ts_2(K_2x_2)
      -\varepsilon_2x_3^T\bar{s}_3(x_3)-\varepsilon_2x_2^TC(x_1+q_d,x_2)x_3+\varepsilon_2x_2^TH(x_1+q_d)x_2\\
     &+\varepsilon_1\varepsilon_2x_2^TH(x_1+q_d)\rho(x_1)
\end{align*}
    According to the properties of each element, and recalling the bound of the 4th element in $\dot{V}_1$ in \eqref{V1.c5} and Remark 1, it can be seen that in $\mathcal{B}^{3n}_\varrho$:
    \begin{align*}
     \dot{V}_2(x_1,x_2,x_3)&\leq -(f_m+\kappa_2 k_{2m})||x_2||^2+\varepsilon_1 f_M||x_1||\cdot||x_2||-\gamma_2\varepsilon_1\kappa_1k_{1m} h_1^*||x_1||^2+\varepsilon_1|\rho^T(x_1)s_2(K_2x_2)|\\
     &+\varepsilon_1k_c\varrho||x_2||^2+\varepsilon_1\mu_M||x_2||^2
     +\varepsilon_2f_M||x_2||\cdot||x_3||+\varepsilon_2k_g||x_1||\cdot||x_3||+\varepsilon_2|x_3^Ts_1(K_1x_1)|\\
     &+\varepsilon_2|x_3^Ts_2(K_2x_2)|-\varepsilon_2x_3^T\bar{s}_3(x_3)+\varepsilon_2\varrho k_c||x_2||^2+\varepsilon_2\mu_M||x_2||^2+\varepsilon_1\varepsilon_2\mu_M||x_1||\cdot||x_2||\\
     \\
     &\leq -(f_m+\kappa_2 k_{2m})||x_2||^2+\varepsilon_1 f_M||x_1||\cdot||x_2||-\gamma_2\varepsilon_1\kappa_1k_{1m} h_1^*||x_1||^2+\varepsilon_1\bar{\kappa}_2k_{2M}||x_1||\cdot||x_2||\\
     &+\varepsilon_1k_c\varrho||x_2||^2+\varepsilon_1\mu_M||x_2||^2
     +\varepsilon_2f_M||x_2||\cdot||x_3||+\varepsilon_2k_g||x_1||\cdot||x_3||+\varepsilon_2\bar{\kappa_1}k_{1M}||x_1||\cdot||x_3||\\
     &+\varepsilon_2\bar{\kappa}_2k_{2M}||x_2||\cdot||x_3||-\varepsilon_2\kappa_3k_{3m}||x_3||^2+\varepsilon_2\varrho k_c||x_2||^2+\varepsilon_2\mu_M||x_2||^2+\varepsilon_1\varepsilon_2\mu_M||x_1||\cdot||x_2||\\
     \\
     &=-\gamma_2\varepsilon_1\kappa_1k_{1m}h_1^*||x_1||^2+\varepsilon_1[(f_M+\bar{\kappa}_2k_{2M})+\varepsilon_2\mu_M]||x_1||\cdot||x_2||\\
     &-[(f_m+\kappa_2 k_{2m})-\varepsilon_1(k_c\varrho +\mu_M)+\varepsilon_2(k_c\varrho+\mu_M)]||x_2||^2+\varepsilon_2(f_M+\bar{\kappa}_2k_{2M})\\
     &-\varepsilon_2\kappa_3k_{3m}||x_3||^2+\varepsilon_2(k_g+\bar{\kappa_1}k_{1M})||x_1||\cdot||x_3||
\end{align*}
By taking $\theta$ and $d_m$ as in \eqref{ThetayDm}, the last expression can be rewritten as:
\begin{equation}
\begin{aligned}[b]
    \dot{V}_2(x_1,x_2,x_3)&\leq-\frac{1}{2}\begin{pmatrix}
          ||x_1||\\||x_2||   
    \end{pmatrix}^TQ_{241}\begin{pmatrix}
          ||x_1||\\||x_2||   
    \end{pmatrix}
    -\frac{1}{2}\begin{pmatrix}
          ||x_1||\\||x_2||   
    \end{pmatrix}^TQ_{242}\begin{pmatrix}
          ||x_1||\\||x_2||   
    \end{pmatrix}
    -\frac{1}{2}\begin{pmatrix}
          ||x_1||\\||x_2||\\||x_3||   
    \end{pmatrix}^TQ_{243}\begin{pmatrix}
          ||x_1||\\||x_2||\\||x_3|| 
    \end{pmatrix}\\
    &-p_1^{24}||x_2||^2\triangleq W_{24}(x_1,x_2,x_3)
        \label{CotaSupDLyapunov2}
    \end{aligned}
    \end{equation}
where
\begin{equation*}
Q_{241}=\begin{pmatrix}
        \gamma_2\varepsilon_1\kappa_1k_{1m}h_1^* && -\varepsilon_1\theta\\
        
        -\varepsilon_1\theta && \frac{d_m}{2}-2\varepsilon_1(k_c\varrho +\mu_M)
        \end{pmatrix}
    \end{equation*}
\begin{equation*}
    Q_{242}=\begin{pmatrix}
    \frac{1}{2}\gamma_2\varepsilon_1\kappa_1k_{1m}h_1^*&& -\varepsilon_1\varepsilon_2\mu_M\\
    \\
    -\varepsilon_1\varepsilon_2\mu_M && \frac{d_m}{2}
    \end{pmatrix}
\end{equation*}
    \begin{equation*}
        Q_{243}=\begin{pmatrix}
        \frac{1}{2}\gamma_2\varepsilon_1\kappa_1k_{1m}h_1^* && 0 && -\varepsilon_2(k_g+\bar{\kappa_1}k_{1M})
        \\
        0 && \frac{d_m}{2}
        && -\varepsilon_2\theta
        \\
        -\varepsilon_2(k_g+\bar{\kappa_1}k_{1M})
        &&-\varepsilon_2\theta
        && 2\varepsilon_2\kappa_3k_{3m}
        \end{pmatrix}
    \end{equation*}
    \begin{equation*}
        p_1^{24}=\frac{d_m}{4}-\varepsilon_2(k_c\varrho+\mu_M)
    \end{equation*}
Notice that $\gamma_{2m2}<\gamma_2<\gamma_{2M}\implies Q_{241}>0$, $\gamma_{2m3}<\gamma_2<\gamma_{2M}\implies Q_{242}>0$, and $\gamma_{2m4}<\gamma_2<\gamma_{2M}\implies Q_{243}>0$. Also, since $\varepsilon_2<\varepsilon_{2M4}$, $p_1^{24}$ is positive. Thus, $W_{24}$ in \eqref{CotaSupDLyapunov2} is negative definite, and so is $\dot{V}_2(x_1,x_2,x_3)$.
\par
Recalling \eqref{CotaInfLyapunov2}, \eqref{CotaSupLyapunov2} and \eqref{CotaSupDLyapunov2}, and defining $\zeta=(x_1^T \ x_2^T\   x_3^T)^T$ the established bounds can be considered as:
\begin{align*}
         V_2(\zeta)&\geq \frac{1}{2}\begin{pmatrix}
        ||x_1||\\||x_2||
        \end{pmatrix}^TQ_{211}
        \begin{pmatrix}
        ||x_1||\\||x_2||
        \end{pmatrix}+\frac{1}{2}\begin{pmatrix}
        ||x_2||\\||x_3||
        \end{pmatrix}^TQ_{212}
        \begin{pmatrix}
        ||x_2||\\||x_3||
        \end{pmatrix}\\
        \\
        &\geq \bar{p}_1^{21}||x_1||^2+ \bar{p}_2^{21}||x_2||^2+\bar{p}_3^{21}||x_3||^2\\
        &\geq\bar{p}^{21}\left[\begin{aligned}
        ||x_1||^2+||x_2||^2+||x_3||^2
        \end{aligned}\right]
        \\
        &=\bar{p}^{21}||\zeta||^2
\end{align*}
where 
\begin{equation*}
    \begin{matrix}
         \bar{p}_1^{21}=\frac{\lambda_m(Q_{211})}{2} && , \bar{p}_2^{21}=\frac{1}{2}[\lambda_m(Q_{211})+\lambda_m(Q_{212})] , && \bar{p}_3^{21}=\frac{\lambda_m(Q_{212})}{2} \\
         \\
        && \bar{p}^{21}=\min\{\bar{p}_1^{21}, \bar{p}_2^{21}, \bar{p}_3^{21}\}
    \end{matrix}
\end{equation*}
\par
\begin{align*}
         V_2(\zeta)&\leq p_{1}^{22}||x_1||^2+p_{2}^{22}||x_2||^2+p_3^{22}||x_3||^2\\
         &\leq \bar{p}^{22}\left[\begin{aligned}
        ||x_1||^2+||x_2||^2+||x_3||^2
        \end{aligned}\right]\\
       &=\bar{p}^{22}||\zeta||^2
\end{align*}
where
\begin{equation*}
    \bar{p}^{22}=\max\{p_1^{22}, p_2^{22}, p_3^{22}\}
\end{equation*}
\par
\begin{align*}
         \dot{V}_2(\zeta)&\leq
         -\bar{p}_1^{24}||x_1||^2- \bar{p}_2^{24}||x_2||^2-\bar{p}_3^{24}||x_3||^2\\
        &\leq-\bar{p}^{24}\left[\begin{aligned}
        ||x_1||^2+||x_2||^2+||x_3||^2
        \end{aligned}\right]
        \\
        &=-\bar{p}^{24}||\zeta||^2
\end{align*}
where
\begin{equation*}
    \begin{matrix}
         \bar{p}_1^{24}=\frac{1}{2}[\lambda_m(Q_{241})+\lambda_m(Q_{242})+\lambda_m(Q_{243})] && \bar{p}_2^{24}=\bar{p}_1^{24}+p_1^{24} && \bar{p}_3^{24}=\frac{1}{2}\lambda_m(Q_{243})\\
         \\
         &&\bar{p}^{24}=\min\{\bar{p}_1^{24},\bar{p}_2^{24},\bar{p}_3^{24}\}
    \end{matrix}
\end{equation*}
\par
Therefore, by Theorem \ref{ExponentialStab}, the origin of the system \eqref{CLsystemSS} is (locally) exponentially stable.

    \chapter{Simulation Results}
    This section has the objective to provide some simulation results that were implemented in order to show the efficiency of the proposed controller and the benefits in comparison to previous schemes. For the simulation, the model of a 2-DOF mechanical system has been considered. Such a model has been involved, for instance, in \cite{IET19}.
    
    \par
    The various matrices and vectors involved in the open-loop system dynamics \eqref{Modelo} are\footnote{For the sake of simplicity, units will be omitted.}
    \begin{align*}
             H(q)&=\begin{pmatrix}
             2.351+0.168\cos{q_2} && 0.102+0.084\cos{q_2}\\
             0.102+0.084\cos{q_2} && 0.102\\
             \end{pmatrix} 
            \\ \\ C(q,\dot{q})&=\begin{pmatrix}
             -0.084\dot{q}_2\sin{q_2} && -0.084(\dot{q}_1+\dot{q}_2)\sin{q_2}\\
              0.084\dot{q}_1\sin{q_2} && 0
             \end{pmatrix}\\
             \\
             F&=\begin{pmatrix}
             2.288 && 0\\
             0 && 0.175
             \end{pmatrix} \\ \\ g(q)&=\begin{pmatrix}
             38.465\sin{q_1}+1.825\sin{q_1+q_2}\\
             1.825\sin{q_1+q_2}
             \end{pmatrix}
    \end{align*}
    whence Assumption \ref{Cotas|g|} can be corroborated to be satisfied with $B_{g_1}=40.29$, $B_{g_2}=1.825$ and $k_g=40.37$. Consider also the input saturation values as $T_1=150$, $T_2=15$. The desired position will be considered to be $q_d=[\pi/4,\pi/2]$ at every one of the implemented simulations tests.\par
    In order to satisfy the design requirements on the functions involved in the proposed control scheme, let us define the next ones based on \cite{IJC15} and \cite{IET19}:
   \begin{subequations}
    \begin{equation}
        \sigma_{bh}(\varsigma;a,M)=\text{sign}(\varsigma)\min\{|\varsigma|^a,M\}
    \end{equation}
    \begin{equation}
        \sigma_{bs}(\varsigma;L,M)=
        \begin{cases}
             \varsigma & \text{if } |\varsigma|\leq L\\
             \text{sign}(\varsigma)L+(M-L)\tanh\left(\begin{aligned}
             \frac{\varsigma- \text{sign}(\varsigma)L}{M-L}
             \end{aligned}\right) & \text{if } |\varsigma|>L
        \end{cases}
    \end{equation}
    \label{SatFunctions}
   \end{subequations}
    with $L$ and $a$ as positive constants such that $0<L<M$ and $ a\leq1$. Then, the controller saturation functions are: 
   \begin{subequations}
    \begin{align}
         \sigma_{ij}(\varsigma)&=\sigma_{bh}(\varsigma;a_i,M_{ij}), \quad \forall i\in\{1,2\}, \forall j\in\{1,2\} \\[1ex] \sigma_{3j}(\varsigma)&=\sigma_{bs}(\varsigma;L_{3j},M_{3j}), \quad\forall j\in\{1,2\}
    \end{align}
    \label{FunctionChoice}
   \end{subequations}
    As for the auxiliary dynamics, the function $h$ is taken as shown in Remark \ref{h's}, with $\varpi=2$. Thus
    \begin{equation*}
        \rho(x_1;b_1/k_{1M})=\frac{(b_1/k_{1M})x_1}{[(b_1/k_{1M})^2+\|x_1\|^2]^{1/2}}
    \end{equation*}
    
    Conditions under which \eqref{Cond.s1>min} and \eqref{Cond.s3>Bg} are satisfied respectively are:
    \begin{subequations}
    \begin{equation}
        k_{1j}>k_g(2B_{gj})^{\frac{1-a_1}{a_1}}
        \label{k1>kg}
    \end{equation}
    \begin{equation}
        M_{1j}>2B_{gj}
        \label{CondM1}
    \end{equation}
    \label{CondParaCond.s1>min}
    \end{subequations}
    \begin{equation}
        M_{3j}>B_{gj}
    \end{equation}
    Inequalities \eqref{CondParaCond.s1>min} are derived in Appendix \ref{ApendiceA}, which in fact constitutes a special case of the derivation presented in \cite[Appendix 8.2]{IET19}.\par
    Therefore, $M_{11}=85, M_{12}=4$, $M_{31}=45, M_{32}=2$. With respect to $M_{2j}$, as \eqref{Cond.u<T} must be satisfied, $M_{21}=19, M_{22}=5$. Then, given the definition of $\sigma_{1j}(\varsigma)$, $b_{1}$ can be calculated by:
    \begin{equation}
        b_1=\min_j\{b_{1j}\}=\min_j\{M_{1j}^{1/a_1}\}
    \end{equation}
    \textbf{Note:} The control law using $0 < a_1 \leq a_2 < 1$, will be referred to as ``FEW-controller" (for fractional exponential weights) or for simplicity ``FEW", while the case where $a_1=a_2=1$ will be denoted as ``UEW- controller" (for unitary exponential weights) or for simplicity ``UEW''.
    \section{FEW vs UEW}
    A first test compares the asymptotic behavior of the system using the generalized controller with $0 < a_1 < a_2 < 1$ versus the use of $a_1 = a_2 = 1$. In the first case, the exponential weights were taken as $a_1=0.8$ and $a_2=0.9$. In both cases $\varepsilon_1=100$ and the gain matrices are 
    \begin{equation}
        \begin{matrix}
             K_1=\begin{pmatrix}
             300 && 0\\0 && 250
             \end{pmatrix}
             &&
             K_2=\begin{pmatrix}
             40 && 0\\0 && 3
             \end{pmatrix}
             &&
             K_3=\begin{pmatrix}
             25 && 0\\0 && 35
             \end{pmatrix}
        \end{matrix}
        \label{GainsTest1}
    \end{equation}\\
    Fig.\ref{GvsEx} shows the position error and the control signal in the first and second DOF for both controllers. Observe that even when in both cases the objective is achieved, the FEW controller is faster than the UEW one (which is clearer from the graph corresponding to link 2).  Moreover, the last one presents a more important overshoot than the FEW (particularly seen in the graph corresponding to link 1). About the input signal $u(t)$, one can see that input saturation avoidance is accomplished in both DOF. \par
   \begin{figure}[H]
       \centering
       \includegraphics[width=.9
       \textwidth]{graficas_matlab.jpg}
        \caption{FEW vs UEW controller}
         \label{GvsEx}
   \end{figure}
   
     In order to measure the approximate stabilization time, we obtain the ``$\rho$-stabilization time" as $t_{\rho} = \inf\{t_s \geq 0 : \|\bar{q}^T(t),\ \dot{q}^T(t)\| \leq \rho,  \forall t \geq t_s\}$. Then, for $\rho=0.001$ the values that are shown in Table \ref{ts}, are obtained from the simulation data. Notice that the $\rho$-stabilization time for the FEW controller is smaller than when using the UEW controller. Furthermore, the graph of the position error signal with approximately 1000 times zoom shows that the UEW controller keeps oscillating around the origin longer than the FEW one does, which is clearer in the first link.\par
     Fig. \ref{s3vsgqd} shows the integral action reaching the desired conservative-force
    compensation for both FEW and UEW controllers. In this case, the stabilization time does not differ too much; nevertheless, the fact that both controllers are able to stabilize at $g(q_d)$ shows that by including the integral action, the steady position error is able to reach $0$ without using the exact value of $g(q_d)$ as in other schemes.
   \begin{figure}[H]
       \centering
       \includegraphics[width=0.9
       \textwidth]{s3vsg_qd.jpg}
        \caption{$s_3(K_3\varphi)$ vs $g(q_d)$}
         \label{s3vsgqd}
   \end{figure}
    \begin{table}[H]
       \centering
        \begin{tabular}{c|c|c}
           &  FEW & UEW \\
        $t_{0.001}$    &  3.2846 s        & 3.959 s
    \end{tabular}
       \caption{$\rho$-stabilization time}
       \label{ts}
   \end{table}
   For further comparisons among both FEW and UEW controllers, and following the procedure shown in \cite{IJC15}, two performance indices are obtained. A first one refers to the sum of the integral of the square of the position error in each DOF, which will be referred to as ``ISE'': $\sum_{j=1}^2\int_{t_{0_j}}^{t_{f_j}}\bar{q_j}(t)^2dt$, where $t_{0_j}$ represents the moment in which the position error signal reaches the equilibrium value for the first time for the $j_{th}$ link, and $t_{f_j}=t_{0_j}+\delta$ with $\delta$ being a positive constant. This index measures the effect of the overshoot as well as the existent oscillations around the equilibrium point. The second index is the integral of the square of the input torques (ISI): $\int_{0}^{t_f}[\sum_{j=1}^2u_j(t)^2]dt$ with $t_f$ as the final simulation time; in this case, the objective is to measure the effort made by the control signal. These indices for both FEW and UEW controllers are shown in the table \ref{PerformanceIndices}. Observe that the ISE shows a considerably smaller value in the FEW case than in the UEW one, which corroborates the presented overshoot in the first link, as well as the oscillation phenomenon. Notice also that the ISI index corroborates that the regulation  objective is achieved with less effort from the control signal in the FEW controller with respect to the UEW one.
   \begin{table}[H]
       \centering
       \begin{tabular}{c|c|c}
        Index  &  FEW. & UEW \\
        ISE    &  4.871$\exp{-05}$         & 0.0247865 \\
        
        ISI    &  4036          & 5279  
       \end{tabular}
       \caption{Performance index evaluation for FEW and UEW controllers}
       \label{PerformanceIndices}
   \end{table}
   \par
      \section{FEW vs Desired conservative force compensation}
      Fig. \ref{GvsExsesgo} shows the comparison between the proposed PID-like generalized scheme presented here and the SP-SD with desired conservative force compensation approach \cite{IET19} under the consideration of a biased estimation of such a compensation term, i.e.
\[ u(q,\dot q) = s_1(K_1\bar q)+s_2(K_2\dot q)+ \hat g(q_d)
\]
From \cite{IET19} one can corroborate that the saturation functions, involved in the proportional and derivative actions in the control law, are analogous to the ones that are expressed in \eqref{FunctionChoice}. Considering the input saturation bounds $T_{1}=150$ and $T_{2}=15$ and the requirement \cite[Eq. (12)]{IET19}, we obtain $B_{gj}$ for $j=1,2$ as $B_{g1}=37.77$ and $B_{g2}=1.825$; then, in order to satisfy the mentioned requirement, the bounds $M_{ij}$ for $i=1,2$ and $j=1,2$ are kept as $M_{11}=85, M_{12}=4$, $M_{21}=19, M_{22}=6$. As for the control gain values, these are considered as in \eqref{GainsTest1} for both controllers. This test simulates an error in the parameter estimated values, more precisely taking $\hat{g}(q_d)=1.2*g(q_d)$. Parameter estimation bias is a possibility when using desired gravity compensation. As it can be seen, when a biased value is entered, the steady state error cannot reach $0$, while the FEW controller reaches that value independently of the parameter values. In fact this is the main advantage of the PID-like controllers.
    \begin{figure}[H]
       \centering
       \includegraphics[width=.9
       \textwidth]{g_qd sesgada.jpg}
        \caption{FEW vs Biased desired conservative force compensation}
         \label{GvsExsesgo}
   \end{figure}
   
   \section{Changing desired position value}
   
   The next test has the objective to show the speed at which the controller is able to respond when a change in the reference occurs. In this case, the change is made at 5 seconds, from $q_d=[\pi/4,\pi/2]$ to $q_d=[\pi/6,\pi/3]$, and with bounds and gains for each action keeping the same values than in the first test. One can see in Fig. \ref{ReferenceChange} that the position error signal presents a better behavior when using the FEW controller, than when using unitary $a$'s, e.g. less overshoot, faster stabilization (particularly seen in the graphs corresponding to link 1).
   \begin{figure}[H]
       \centering
       \includegraphics[width=.9
       \textwidth]{CambioReferencia.jpg}
        \caption{FEW controller vs UEW controller: change of desired position}
         \label{ReferenceChange}
   \end{figure}
       Moreover, the $\rho$-stabilization time from 5 seconds (the instant when the reference change is introduced) for $\rho=0.0001$ is as shown in Table \ref{tsChangeOfRef}, where it can be corroborated that the FEW controller still achieves stabilization in less time, than when using unitary $a$'s. 
       \begin{table}[H]
       \centering
      \begin{tabular}{c|c|c}
           &  FEW & UEW \\
        $t_{0.0001}$    &  2.7962 s        & 3.2969 s
    \end{tabular}
     \caption{$\rho$-stabilization time at the reference change}
       \label{tsChangeOfRef}
   \end{table}
   \section{FEW vs S10}
   For comparison purposes (involving the controller developed here), a simulation of another PID-like scheme is worthy; specifically the one proposed in \cite{Su10}, which will be referred to as \textbf{S10}, and is expressed as
   \begin{subequations}
   \begin{equation}
          u=-K_P\tanh(\bar q)-K_I\tanh(\phi)-K_D\tanh(\dot q)
   \end{equation}
   \begin{equation}
       \dot \phi =\alpha^2\dot q+\alpha\tanh(\bar q)
   \end{equation}
   \label{S10}
   \end{subequations}
   where $\alpha$ is a large enough positive constant and with $\tanh(x)=[\tanh x_1,..., \tanh x_n]$ for any  $x\in\mathbb{R}^n$. As it is explained in the experimental results section in \cite{IJC15}, the initial condition for the variable $\phi$ is forced in \cite{Su10} to be taken as $\phi(0)=\alpha^2\bar q(0)$. Moreover, the design requirements derived in \cite{Su10} are taken into account. Then, for the FEW controller, the gain matrices are taken as in \eqref{GainsTest1} with $a_1=0.8$ and $a_2=0.9$. On the other hand, for the S10 controller, the gain matrices are\footnote{Notice that in the case of the S10 controller, in Eqs. \eqref{S10}, the control gains simultaneously play the role of the saturation value of the corresponding (P, I, D) control action, which imposes a restriction on their size through which the control inputs be guaranteed to avoid their corresponding input saturation value (in view of which they cannot be arbitrarily large).}
   \begin{align}
       K_1=\begin{pmatrix}
       25 & 0 \\ 0 & 9
       \end{pmatrix} 
       && K_2=\begin{pmatrix}
       30 & 0 \\ 0 & 3.1
       \end{pmatrix} 
       && K_3=\begin{pmatrix}
       45 & 0 \\ 0 & 2
       \end{pmatrix} 
       \label{GainS10}
   \end{align}
   A first simulation using $\alpha=23$ returns the results shown in Fig. \ref{S10graficas23}, where for both controllers, the stabilization objective is achieved avoiding  input saturation. Notice that the S10 controller seems to be faster than the FEW, but when zoomed in (about 100 times), it can be seen that the S10 does not get as close to the origin throughout the duration of the simulation as the other one. Nevertheless, a tendency of decrease is noted in the S10 case.
    \begin{figure}[H]
       \centering
       \includegraphics[width=.9
       \textwidth]{Su10_Ks5_alpha23.jpg}
        \caption{Position error and control signal for FEW and S10 controllers, $\alpha=23$}
         \label{S10graficas23}
   \end{figure}

   \begin{figure}[H]
       \centering
       \includegraphics[width=.9
       \textwidth]{Su10_Ks5_alpha23_Norm.jpg}
        \caption{Position and velocity errors norm for FEW and S10 controllers, $\alpha=23$}
         \label{S10graficas23norm}
   \end{figure}
   On the other hand,  the norm on the position and velocity errors for the two cases are illustrated in Fig. \ref{S10graficas23norm}, where once again it can be seen that even when the norm in S10 is very close to the origin, the FEW signal converge considerably faster.\par
   In order to apply a more specific criterion for comparison purposes, the $\rho$-stabilization time is obtained with $\rho=0.01$ and $\rho=0.0065$ (Table \ref{tsSu10_23}). These data make it easier to notice that the FEW error norm converge faster than the S10 one. Observe that for both values of $\rho$ the stabilization time values remain close in the case of the FEW controller, while for the S10 one, the difference between its stabilization time values is quite large (about 5.6 seconds). It is further worth adding that when using a smaller value of $\rho$, the simulation time (10 seconds) was not enough to get a $\rho$-stabilization time for the S10 case.
   \begin{table}[H]
       \centering
      \begin{tabular}{c|c|c|c}
           &  FEW  & S10\\
        $t_{0.01}$    &  2.8921 s         & 3.2301 s
        \\
        $t_{0.0065}$    &  2.9212 s         & 8.8477 s
    \end{tabular}
     \caption{$\rho$-stabilization time for FEW and S10 controllers}
       \label{tsSu10_23}
   \end{table}
   Let us emphasize that with $\alpha=23$ the design requirements derived in \cite{Su10} are satisfied giving rise to a simulation with a response that is not optimal. Nevertheless a second simulation with the same gain matrices and using an $\alpha$ value for which all the requirements are not fulfilled was made, i.e. $\alpha=5$, obtaining faster responses and consequently a smaller $\rho$-stabilization time. It is worth mentioning that the design requirement that was not satisfied is not related to the input saturation avoidance condition. Position error and control signal are shown in Fig. \ref{S10graficas5}, where in each DOF, the position error gets closer to the origin quicker than when using $\alpha=23$. 
    \begin{figure}[H]
       \centering
       \includegraphics[width=.9
       \textwidth]{Su10_Ks5_alpha5.jpg}
        \caption{Position error and control signal for FEW and S10 controllers, $\alpha=5$}
         \label{S10graficas5}
   \end{figure}
   As in the previous $\alpha$ choice, the norm on the position and velocity errors is obtained for both cases (Fig \ref{S10graficas5norm}), which corroborates a faster convergence to the origin. However, even with the convergence improvement, the FEW controller shows a faster response behavior. In addition, the $\rho$-stabilization time is derived with $\rho=0.002$ for the two control laws as it is presented in Table \ref{tsSu10_5}. Notice that as happened with $\alpha=23$, the S10 control law gives rise to the largest stabilization time.
      \begin{figure}[H]
       \centering
       \includegraphics[width=.9
       \textwidth]{Su10_Ks5_alpha5_Norm.jpg}
        \caption{Position and velocity errors norm for FEW and S10 controllers, $\alpha=5$}
         \label{S10graficas5norm}
   \end{figure}
    \begin{table}[H]
       \centering
      \begin{tabular}{c|c|c|c}
           &  FEW  & S10\\
        $t_{0.002}$    &  3.3820 s        & 6.3886 s
    \end{tabular}
     \caption{$0.002$-stabilization time for FEW and S10 controllers}
       \label{tsSu10_5}
   \end{table}
   \section{Noise in the feedback signal}
   Let us note that the simulation results do not need to deal with the way of getting the position or velocity, either by direct measurement or by an approximation, since this was not considered in the problem formulation and those vectors are consequently obtained directly from the simulation. Nevertheless in order to include a possible difficulty during the measurement of those vectors, a test of the controller under noise in the feedback (position or velocity) signal is made. In order to simulate the noise, a random number signal was used, for which the variance and sample time can be changed; those values were chosen according to the test situation. As for the tests, three cases are developed, for which FEW and UEW controllers are compared. For all the cases the gain matrices were kept as:
   \begin{align}
       K_1=\begin{pmatrix}
       200 & 0 \\ 0 & 150
       \end{pmatrix} 
       && K_2=\begin{pmatrix}
       40 & 0 \\ 0 & 3
       \end{pmatrix} 
       && K_3=\begin{pmatrix}
       15 & 0 \\ 0 & 25
       \end{pmatrix} 
       \label{GainNoise}
   \end{align}
additionally taking $\varepsilon_1=90$, and $a_1=0.7$ and $a_2=0.8$ for the FEW controller; the desired position is kept in $q_d=[\pi/4 ,\pi/2]$, and $T_1=150, T_2=15$. A first case includes noise just in the feedback position vector (see Fig. \ref{RuidoPGráficas}), with a sample time of $0.05$ and a variance of $1\times10^{-5}$. Despite the presence of noise, the control signal (for both links) is kept below the saturation bound $T_j$; nevertheless, a significant detail is the important effect of the noise in the control signal for the 2nd link compared to the 1st one. On the other hand, it can be noted that the difference between the position error response of both controllers is minimum, as the noise affects both of them almost in the same way.  \par
 \begin{figure}[H]
       \centering
       \includegraphics[width=.9
       \textwidth]{RuidoPGráficas.jpg}
        \caption{Position error and control signal with noise in the position feedback}
         \label{RuidoPGráficas}
   \end{figure}
As a more exact way to measure the effect of the noise in both controllers responses, the ISE is obtained as $\int_{t_0}^{t_f}[\sum_{j=1}^2\bar q_j(t)^2]dt$, where $t_0=t_\rho$ (with $\rho$ corresponding to the $\rho$-stabilization time); and $t_{f}=t_{0}+\delta$ with $\delta=5$. For this measurement, $\rho$ took the value $\rho=0.6$, which unlike previous cases, was higher due to the presence of the noise, and even when the magnitude of the noise is relatively small. The ISE results (see Table \ref{ISENoiseP}) corroborate the first sight about the difference between the responses from both controllers. Notice that in this case, the calculation of the ISE index indicates a higher impact of the noise for the UEW controller.  
   
   \begin{table}[H]
       \centering
      \begin{tabular}{c|c|c|c}
           &  FEW  & UEW\\
        ISE    &  0.006902        & 0.007732
    \end{tabular}
     \caption{ISE for FEW and UEW controllers with noise in the position feedback}
       \label{ISENoiseP}
   \end{table}
  A second test, contemplates noise involved in the velocity feedback vector. The results for position error and control signal are presented in Fig. \ref{RuidoVGráficas}, where  the sample time is $0.05$ and the variance was taken as $1\times10^{-3}$. Notice that the variance value for the velocity is considerably smaller than the position one, due to the major impact that the noise signal has in the development of the response by being involved in the position feedback. As in the previous case, the position error responses for both controllers are almost the same, which is corroborated to the calculation of the ISE (Table \ref{ISENoiseV}) with an initial time of $t_\rho$ and $t_{f}=t_{0}+\delta$ with $\delta=5$. In this case, as the impact of the noise is less than in the previous one, the $\rho$ value was taken as $\rho=0.3$. Notice that in this case, the ISE index calculation indicates a higher impact for the FEW controller, although this was only around 40\% higher than the that calculated for the UEW controller.
   \begin{figure}[H]
       \centering
       \includegraphics[width=.9
       \textwidth]{RuidoVGráficas.jpg}
        \caption{Position error and control signal with noise in the velocity feedback}
         \label{RuidoVGráficas}
   \end{figure}
   \begin{table}[H]
       \centering
      \begin{tabular}{c|c|c|c}
           &  FEW  & UEW\\
        ISE    &  0.0003544        & 0.0002526
    \end{tabular}
     \caption{ISE for FEW and UEW controllers with noise in the velocity feedback}
       \label{ISENoiseV}
   \end{table}
   The last test contemplates the noise in both position and velocity vectors simultaneously, for which the same values of variance were taken respectively, i.e. $1\times10^{-5}$ for position and $1\times10^{-3}$ for velocity, with a sample time of 0.05 for both cases. When these two signals of noise are involved, the position error response is affected in a larger scale than when the noise is acting separately, as it can be seen in Fig. \ref{RuidoPVGráficas}. The control signal stays below the saturation bound, while as for the position error, the similarity between FEW and UEW is kept. Moreover, by involving the ISE with an initial time of $t_\rho$ for $\rho=0.7$ and  $t_{f}=t_{0}+\delta$ with $\delta=5$ (see Table \ref{ISENoisePV}), we notice that as in the previous case, despite the similarity between the values of ISE, the FEW controller gets a slightly larger value than the UEW one. This fact, as expected, is explained due to the infinite slope around the origin in the case of the FEW, since fractional exponential weights are involved. However, as previously observed, the relation among both ISE calculations reflects that the impact of noise is not considerably higher in the case of the FEW controller as one might expect (around 39\% higher than in the UEW case).
    \begin{figure}[H]
       \centering
       \includegraphics[width=.9
       \textwidth]{RuidoPVGráficas.jpg}
        \caption{Position error and control signal with noise in both position and velocity feedback}
         \label{RuidoPVGráficas}
   \end{figure}
   \begin{table}[H]
       \centering
      \begin{tabular}{c|c|c|c}
           &  FEW  & UEW\\
        ISE    &  0.01156        & 0.008328
    \end{tabular}
     \caption{ISE for FEW and UEW controllers with noise in both position and velocity feedback}
       \label{ISENoisePV}
   \end{table}
   \chapter{Conclusions}
   A generalized PID-like scheme was proposed in this work, through which it is possible to solve the regulation problem of mechanical systems with bounded inputs. Based on previous approaches, including finite-time schemes, PD-like ones with desired gravity compensation and other type of PID-like control laws, the proposed scheme's structure involves a special type of functions which do not only prove to be instrumental to achieve the avoidance of input saturation but also bring to the approach certain advantages over previous schemes. The adopted structure of the control law includes exponential weights acting on these special type of functions, whose election gives rise to the so-called ``non-Lipschitz'' characteristic of the control actions. Global asymptotic stability and (local) exponential stability for a certain choice of the exponential weights were thoroughly proven to be obtained through the Lyapunov method and invariance principle. This and the avoidance of the saturation bound are corroborated through both the analytical proof and simulation tests. As for the simulation results, these were gotten involving a model of a 2-DOF robot manipulator, considering bounded inputs and all the necessary conditions that were derived in the stability analysis. The efficiency of the proposed controller was corroborated through several tests, from which some advantages like less overshoot on the position error response and less effort from the control signal were observed. These simulations were made in order to evaluate the controller with different choices of the exponential weights, and also to compare the proposed scheme against other control laws.  Furthermore, some tests using a simulated noise in the position and velocity feedback were developed. The simulation results encourage the use of the proposed scheme with fractional exponential weights by mostly showing performance improvements over the use of unitary exponential weights and over other controllers.
   \par
   The proposed control law  contains the proportional, derivative and integral like actions, where the last one is calculated by means of an auxiliary dynamics, which involves the position error and velocity vectors. Let us note that this auxiliary dynamics has such a structure that makes possible to develop the stability analysis giving rise to the expected result, so as to guarantee the global stabilization objective. As mentioned earlier, a special type of functions are part of the structure, specifically, strongly passive functions with a bounded growth restriction on the proportional and derivative actions, where exponential weights are incorporated; and bounded strictly passive functions for the integral action. The choice of the mentioned exponential weights can be made between $(0,1]$, which makes possible to get a more generalized form of the PID-like control law, and which also resembles to structures used in finite-time approaches. It is worth mentioning that through the involved definition of functions, the control law is not forced to adopt a specific shape, like $\tanh(\cdot)$, as in other approaches. Moreover, the bounds of the involved type of functions are directly responsible for satisfying the constrained input condition. Other design requirements were presented, so that the conditions for stability could be achieved. 
   \par
   As for the stability analysis, a first proof was developed by proposing a Lyapunov function, by means of which global asymptotic stability was demonstrated through Lyapunov's direct method and invariance theory. On the other hand, through a local analysis on a ball with a suitable radius, (local) exponential stability was proven to be achieved involving an alternative Lyapunov's function depending on the first one. An important detail of these two analyses, is the fact that by changing the value of the exponential weights from a fractional value to the unity, the origin is demonstrated to go from just having global asymptotic stability to have also (local) exponential stability. However, through fractional exponential weights, the simulation tests showed improvement on several closed-loop performance aspects, compared to results obtained with unitary exponential weights.
   \par
   In order to corroborate the results, several simulation tests were made involving a model of a 2-DOF manipulator robot, where the availability of the position and velocity data was considered without the need for extra calculations. Beginning with the satisfaction of the conditions obtained from the closed-loop analysis, it was possible to see that in every case, the control signal remains always within the saturation limits for each link (which is one of the initial requirements). Moreover, as it was seen through the graphs, the closed-loop performance is improved when using fractional exponential weights (FEW), in the sense of less overshoot on the position error response, a faster stabilization and less effort from the control signal. Furthermore, it was possible to see that the UEW (for \emph{unitary exponential weights}) controller presents more oscillation around the origin. By incorporating the ISE and ISI performance indices, the differences between the FEW and UEW controllers were clearer, considering that in almost every test, smaller values of ISE and ISI were obtained for the FEW case, which can be interpreted as a better closed-loop performance. As for the $\rho$-stabilization time, the use of this index was useful to notice that a larger stabilization time was gotten when unitary exponential weights were used. All these results were also corroborated through the test in which a change of reference was made.\par
   
   Regarding the noise test, in addition to what was observed in the graphs, the ISE index was very useful in order to determine which one of the error signals of the controllers stays closer to the reference. A first thought was that the FEW controller was going to have a worse behavior because of the fractional exponential weights, which have a direct impact on the form of the functions around the origin (the closer to the origin, the larger the slope), and which in consequence would cause the noise to be amplified. Nevertheless, as observed on the graph, despite the advantage that the UEW has over the FEW controller, it is worth noting that the behavior in the FEW case was relatively similar to the UEW one; in fact, the difference on the ISE index was relatively small. This is encouraging for the proposed scheme with fractional exponential weights considering that the closed-loop performance is better and with no significant difference on the noise impact.
   \par
   All in all, through the proposed scheme, numerous advantages have been evidenced, including the more generalized form of the control law, the attributes of the response behavior and the ease with which the parameters can be chosen in an implementation. Recalling the initial motivation of the work, despite the apparent impossibility to achieve finite-time stabilization, let us note that the incorporation of a similar structure to the ones that are provided in approaches like \cite{IET19} and \cite{FTTOFRNC09}, was useful in order to improve the closed-loop performance, which was indeed the new main goal.
   Future work might be addressed to formally conclude if finite-time stabilization is indeed impossible to be achieved through PID-type controllers (particularly under input constraints), or to eventually find an alternative (probably more complex) PID-type structure through which finite-time stabilization could be achieved.
   
\begin{thebibliography}{99}
\bibitem{Bacciotti}
A. Bacciotti \& L. Rosier. \textit{Liapunov Functions and Stability in Control Theory} . The Netherlands. Springer. 2005.
\bibitem{Chen88}
B. S. Chen \& S. S. Wang. The stability of feedback control with nonlinear
saturating actuator: Time domain approach.  \textit{IEEE Transactions on Automatic Control}. 33(5): 483-488. May 1988.
\bibitem{Gorez} 
R. Gorez. Globally stable PID-like control of mechanical systems. \textit{Systems} \&  \textit{Control Letters}. 38: 61-72, May 1999.
\bibitem{Kapasouris} 
P. Kapasouris \& M. Athans. Control systems with rate and magnitude saturation
for neutrally stable open loop systems. \textit{29th IEEE Conference on Decision and Control}. 6: 3404-3409. Dec. 1990.
\bibitem{ControlOfRobots}
R. Kelly, V. Santibáñez, A. Loría. \textit{Control of Robot Manipulators in Joint Space}. Springer. London. 2005.
\bibitem{Khalil}
H. Khalil. \textit{Nonlinear Systems Third Edition}. Prentice Hall. United States of America. 2002.
\bibitem{Krikelis}
N. J. Krikelis \& S. K. Barkas. Design of tracking systems subject to actuator saturation an integrator wind-up. \textit{International Journal of Control}. 39(4): 667-682. 1984.
\bibitem{Marquez}
H. Márquez. \textit{Nonlinear control system}.Wiley. Hoboken, New Jersey. pp.52-54. 2003.
\bibitem{IJC15} 
M. Mendoza, A. Zavala, V. Santibáñez \& F. Reyes. A generalised PID-type
control scheme with simple tuning for the global regulation of robot manipulators with constrained inputs. \textit{International
Journal of Control}. 88(10):1995-2012. 2015.
\bibitem{Meza05}
J. L. Meza, V. Santibáñez \& V. M. Hernández. Saturated nonlinear PID global regulator for robot manipulators passivity based analysis. \textit{In P. Zitek (Ed.), Proceedings of the
16th IFAC World Congress}. Prague: IFAC. 2005.
\bibitem{Michel}
A. Michel, L. Hou \& D. Liu. \textit{Stability of Dynamical Systems}. Birkhäuser. Boston. 2008.
\bibitem{Su10}
Y. Su, P. Müller \& C. Zheng. Global asymptotic saturated PID control
for robot manipulators. \textit{IEEE Transactions on Control Systems Technology}. 18(6): 1280-1288. Nov. 2010.
\bibitem{TAC20} 
G. I. Zamora , A. Zavala, D. J. López, E. Cruz \& E. Nuño. Continuous control for fully damped mechanical systems with input
constraints: Finite-time and exponential tracking. \textit{IEEE Transactions on Automatic Control}. 65(2): 882-889. Feb. 2020.
\bibitem{RNC20}
G. I. Zamora, A. Zavala, E. Vázquez, F. Reyes \& V. Santibáñez. A robustness study of a finite-time/exponential tracking
continuous control scheme for constrained-input
mechanical systems: Analysis and experiments. \textit{International Journal of Robust Non-linear Control}. 30: 3923–3944. Nov. 2020
\bibitem{IET19} 
G. I. Zamora, A. Zavala, D. J. López \& V. Santibáñez. Further results on the global continuous
control for finite-time and exponential
stabilization of constrained-input mechanical
systems: desired conservative-force
compensation and experiments. \textit{IET Control Theory \& Applications}. 13(2): 159-170. 2019.

\bibitem{TAC14} 
A. Zavala \& I. Fantoni. Global finite-time stability characterized through
a local notion of homogeneity. \textit{IEEE Transactions on Automatic Control}. 59(2): 471-477. Feb. 2014.
\bibitem{TCST06} 
A. Zavala \& V. Santibáñez. Simple extensions of the PD-with-gravity-compensation control law
for robot manipulators with bounded inputs. \textit{IEEE Transactions on Control Systems Technology}. 14(5): 958-965. Sep. 2006.
\bibitem{IJC17} 
A. Zavala \& G. I. Zamora. Local-homogeneity based
global continuous control for mechanical systems with constrained inputs: finite time and exponential stabilization. \textit{International Journal of Control}. 90(5): 1037-1051. 2017.
\bibitem{FTTSFgRNC03a}
A. Zavala-Rio, G. I. Zamora-Gomez, T. Sanchez, F. Reyes & D.J. Lopez-Araujo. A more exhaustive study on the finite-time and exponential tracking continuous control for constrained-input mechanical systems: improved design and experiments. International Journal of Robust Nonlinear Control. 2022. Accepted.

\bibitem{FTTOFRNC09} 
A. Zavala, G. I. Zamora, T. Sánchez, F. Reyes \& V. Santibáñez. Output-feedback finite-time and exponential tracking
continuous control for mechanical systems with
constrained inputs. \textit{International Journal of Robust Nonlinear Control}. 3(32): 1393–1424. 2022. 

\end{thebibliography}
\begin{appendices}
    \chapter{}
On inequalities
 \eqref{CondParaCond.s1>min}:\par 
Observe that on $\{\varsigma\in\mathbb{R}: 0<|\varsigma|\leq\frac{2B_{gj}}{k_g}\}$:
\begin{equation*}
    \begin{matrix}
         |\varsigma|\leq\begin{aligned}
         \frac{2B_{gj}}{k_g}
         \end{aligned}\Longleftrightarrow |\varsigma|^{1-a_1}\leq\left(\begin{aligned}
         \frac{2B_{gj}}{k_g}
         \end{aligned}\right)^{1-a_1}
         \Longleftrightarrow k_{1j}^{a_1}\begin{aligned}
         \frac{|\varsigma|}{|\varsigma|^{a_1}}
         \end{aligned}\leq \left(\begin{aligned}
         \frac{2B_{gj}}{k_g}
         \end{aligned}\right)^{1-a_1}k_{1j}^{a_1}
         \\
         \Longleftrightarrow k_{1j}\left(\begin{aligned}
         \frac{2B_{gj}}{k_g}
         \end{aligned}\right)^{a_1-1}|\varsigma|\leq|k_{1j}\varsigma|^{a_1}
    \end{matrix}
\end{equation*}
Then, from \eqref{k1>kg}, it can be noticed that $\forall \varsigma\neq0$:
\begin{equation*}
    \begin{matrix}
         \eqref{k1>kg} \Longrightarrow k_{1j}>k_g(2B_{gj})^{\frac{1-a_1}{a_1}}
         \Longleftrightarrow
         k_g(2B_{gj})^{\frac{1-a_1}{a_1}}|\varsigma|^{1/a_1}<k_{1j}|\varsigma|^{1/a_1}
         \\ \\
         \Longleftrightarrow k_g^{a_1}(2B_{gj})^{1-a_1}|\varsigma|<k_{1j}^{a_1}|\varsigma|\\ \\
         \Longleftrightarrow k_g|\varsigma|\left(\begin{aligned}
         \frac{2B_{gj}}{k_g}
         \end{aligned}\right)^{1-a_1}<k_{1j}^{a_1}|\varsigma|\\ \\
         \Longleftrightarrow
         k_g|\varsigma|<k_{1j}^{a_1}\left(\begin{aligned}
         \frac{2B_{gj}}{k_g}
         \end{aligned}\right)^{a_1-1}|\varsigma|
         \end{matrix}
\end{equation*}
    Therefore, on $\{\varsigma\in\mathbb{R}: 0<|\varsigma|\leq\frac{2B_{gj}}{k_g}\}$, \eqref{k1>kg} $\Longrightarrow k_g|\varsigma|<|k_{1j}\varsigma|{a_1}$. Moreover, note that $\forall \varsigma\neq0$, $\eqref{k1>kg}\Longrightarrow\min\{k_g|\varsigma|,2B_{gj}\}<|k_{1j}\varsigma|^{a_1}$. Consequently, additionally considering \eqref{CondM1} it can be seen that:
    \begin{equation*}
        \eqref{CondParaCond.s1>min} \Longrightarrow\min\{k_g|\varsigma|,2B_{gj}\}<\min\{|k_{1j}\varsigma|^{a_1},M_{1j}\}=|\sigma_{1j}(k_{ij})|, \forall\varsigma\neq0
    \end{equation*}
    \label{ApendiceA}
\end{appendices}
\begin{appendices}
    \chapter{}
    From \cite{TAC20} we recall the next Lemma.
    \begin{lemma}
    Let $\sigma:\mathbb{R}\to\mathbb{R}$ be a strongly passive function for $(\kappa,a,b)$ and $k$ be a positive constant. Then, for all $\varsigma\in\mathbb{R}$
    \begin{equation}
        \int_0^{\varsigma}\sigma(kz)dz\geq S(\varsigma)=\begin{cases}
             \frac{\kappa k^a}{1+a}|\varsigma|^{1+a} & \forall|\varsigma|\leq\frac{b}{k}\\
             \kappa b\left(|\varsigma|-\frac{ab}{k(1+a)}\right) & \forall |\varsigma|>\frac{b}{k}
        \end{cases}
    \end{equation}
    \label{Lema2.3TAC20}
    \end{lemma}
    \label{ApendiceB}
\end{appendices}
\end{document}