Identifying Friction in a Nonlinear Chaotic System Using a Universal Adaptive Stabilizer

Ali Wadi, Shayok Mukhopadhyay, Lotfi Romdhane. A. Wadi is with the department of Mechanical Engineering, American University of Sharjah, P.O. Box 26666, Sharjah, United Arab Emirates.S. Mukhopadhyay is with the department of Electrical Engineering, American University of Sharjah, P.O. Box 26666, Sharjah, United Arab Emirates.L. Romdhane is with the department of Mechanical Engineering, American University of Sharjah, P.O. Box 26666, Sharjah, United Arab Emirates.Corresponding author S. Mukhopadhyay [email protected].

Abstract

This paper proposes a friction model parameter identification routine that can work with highly nonlinear and chaotic systems. The chosen system for this study is a passively-actuated tilted Furuta pendulum, which is known to have a highly non linear and coupled model. The pendulum is tilted to ensure the existence of a stable equilibrium configuration for all its degrees of freedom, and the link weights are the only external forces applied to the system. A nonlinear analytical model of the pendulum is derived, and a continuous friction model considering static friction, dynamic friction, viscous friction, and the stribeck effect is selected from the literature. A high-gain Universal Adaptive Stabilizer (UAS) observer is designed to identify friction model parameters using joint angle measurements. The methodology is tested in simulation and validated on an experimental setup. Despite the high nonlinearity of the system, the methodology is proven to converge to the exact parameter values, in simulation, and to yield qualitative parameter magnitudes in experiments where the goodness of fit was around 85% on average. The discrepancy between the simulation and the experimental results is attributed to the limitations of the friction model. The main advantage of the proposed method is the significant reduction in computational needs and the time required relative to conventional optimization-based identification routines. The proposed approach yielded more than 99% reduction in the estimation time while being considerably more accurate than the optimization approach in every test performed. One more advantage is that the approach can be easily adapted to fit other models to experimental data.

Index Terms:

Furuta pendulum; Rotary inverted pendulum; Parameter Identification; Universal Adaptive Stabilizer; Viscous Friction; Dry Friction

I Introduction

The Furuta pendulum (also known as a rotary inverted pendulum) is an under-actuated system, which is formed by two bodies in series. The first one is rotating around a vertical axis followed by a pendulum rotating around a horizontal axis. It was first proposed by Furuta [1] mainly to test different control laws [2]. The Furuta pendulum is mainly used as a testbed of nonlinear control strategies and is also of educational value. Many of the works in the literature use the Furuta pendulum to illustrate a proposed control law [2, 1, 3, 4]. However, these works did not elaborate on the dynamic model of the Furuta pendulum and they limited their work to a simplified model [5, 6].

Few papers dealt with the modeling of the Furuta pendulum. They highlighted the complexity of the dynamic model and the necessity to simplify it based on several assumptions. The authors of [1] and [7] proposed a linearized model based on a small angle assumption. Others neglected one or more of the cross coupling terms relating the two rotations [1, 3, 6, 8, 9]. The authors of [7] showed, through simulations, that the aforementioned assumption is not acceptable and it could have an important effect on the dynamic response. The authors of [8] derived a significantly simplified system and carried a brief three-second test. All the experimental works on the Furuta pendulum are directed at testing various control strategies [10, 11]. The nonlinear dynamic behavior of the system, which includes Coriolis and centrifugal forces in the two rotating frames of motion, is often overlooked or linearized. However, the knowledge of the dynamic behavior of a system is highly desirable to design advanced nonlinear controllers for a given nonlinear system. In this paper, the developed Furuta pendulum model, which encompasses the aforementioned dynamics, is validated experimentally.

The problem of identifying the different parameters of the Furuta pendulum is also treated in the literature [6, 12, 13]. Most of the authors dealing with this problem mentioned the high complexity of the model, especially the friction effects arising from the wide range of velocities the links go through. The parameters representing the dry and viscous frictions in the joints are one of the most difficult parameters to identify [14, 8, 15, 16, 17, 18, 19, 20]. Indeed, there is no universal model to represent damping, and the literature contains several models that attempt to capture the physics of this complex phenomenon [21, 22, 23].

When the Furuta pendulum is used in a passive mode, the gravity, which is the only external load, actuates only the second joint. In this case, the Furuta pendulum does not have a stable configuration, which makes any reproducible experimentation nearly impossible. Moreover, this will cause the passive Furuta pendulum to exhibit chaotic behavior such that the motion of both links become heavily dependent and sensitive to initial conditions. This fact could explain the reason of the nonexistence in the literature of dynamic models of a passively actuated Furuta pendulum. The closest work to passive excitation was in [24] where friction in a cart pendulum system was identified. Nevertheless, that work decoupled the pendulum from the cart and dealt with them separately in the identification phase. A significant portion of the literature work, attempting parameter estimation, chooses an active excitation of the setup, often alongside a controller that helps stabilize the system. However, a controller might mask system dynamics and thus hinder parameter identification results. Moreover, a dynamic model of the passive system could capture the physics of the system without the perturbation of the input torque. In this paper, to actuate both joints in the passive mode, we propose to tilt the Furuta pendulum. In this case, the Furuta pendulum has a stable equilibrium configuration for both links, which is defined by gravity. The oscillations of the pendulum around this equilibrium could then be studied analytically and experimentally. After tilting, the model reveals a contribution of the gravity in the two generalized forces applied to the two joints. The existence of the stable equilibrium configuration, around which the perturbed system oscillates, is also seen. The developed model includes the effect of the static and dynamic frictions, viscous friction, and the stribeck effect. We propose a novel Universal Adaptive Stabilizer-based model identification algorithm and apply it to estimate the parameters of the friction model in the joints of the Furuta pendulum. The proposed approach is tested in simulation and in experiment, and it is compared against grey box optimization-based parameter estimation. After identification, validation is performed by simulating the system model using the identified parameters and comparing that to the experimental time response of the system. The coefficient of determination, $R^{2}$ , is used as a benchmark statistic, and the computation time the algorithm takes to run is also noted for comparison with the optimization-based approach.

The rest of the paper is organized as follows: after the introduction, Section II introduces the model developed. The identification approach is presented in Section III. The comparison of the results, along with a discussion, are presented in Section IV. Some concluding remarks and future work are presented in Section V.

II Materials and Methods

Refer to caption — Figure 1: The Furuta pendulum in its tilted position.

TABLE I: Parameters of the Furuta Pendulum Model.

Parameter	Unit	Arm link	Pendulum link
Mass	$[Kg]$	$m_{1}=0.370$	$m_{2}=0.128$
Inertia	$[Kg\,m^{2}]$	$j_{1_{z}}\,=3.09\times 10^{-3}$	$\begin{matrix}j_{2_{x}}=5.25\times 10^{-3}\\ j_{2_{y}}=5.25\times 10^{-3}\\ j_{2_{z}}=2.91\times 10^{-6}\end{matrix}$
Pivot to $CG$ distance	$[m]$	$l_{1}\;=OG_{1}=0.0620$	$l_{2}=AG_{2}=0.0620$
Total length	$[m]$	$L_{1}\;=OA=0.216$	$L_{2}=AB=0.316$

II-A The Furuta Pendulum

Figure 1 shows the Furuta pendulum in its tilted position. The angle $\phi$ is taken around the fixed y-axis, which is pointing out of the page. The obtained tilted reference frame is $x_{0}-y_{0}-z_{0}$ , where $z_{0}$ is now the axis of the first joint between the ground and the arm. The joint is disconnected from the existing motor and the pendulum is subjected only to its weight. Figure 2 shows the schematic of the Furuta pendulum in its tilted position. The arm is rotating around the fixed axis $z_{0}$ by an angle $\theta_{0}$ and the pendulum is rotating with respect to the moving axis $x_{1}$ by an angle $\theta_{1}$ . The different parameters of the two links are defined in Table I.

II-B Dynamic Model of the Furuta Pendulum

Based on the schematic of the Furuta pendulum, shown in Figure 2, the dynamic model is elaborated. The Euler-Lagrange Formulation is used to derive the equations of the motion of the system. Two generalized coordinates are required to fully describe the dynamics of the two-degrees-of-freedom system, and they are the angular displacements of the Arm and Pendulum links given by

q=\begin{bmatrix}\theta_{0}\\ \theta_{1}\end{bmatrix}

(1)

The Euler-Lagrange equation, for the $i^{th}$ generalized coordinate of the system can be written as follows:

\frac{d}{dt}\left(\frac{\partial T}{\partial\dot{q}_{i}}\right)-\frac{\partial T}{\partial q_{i}}+\frac{\partial V}{\partial q_{i}}=Q_{i}\,

(2)

where the kinetic energy, $T$ , is written as the sum of the linear and angular kinetic energies of the two links, which is given by

T=\sum_{i=1}^{2}\frac{1}{2}m_{i}{v}_{G_{i}}^{2}+\sum_{i=1}^{2}\frac{1}{2}m_{i}{\omega}_{i-1}^{T}({J}_{i}{\omega}_{i-1}),

(3)

where ${\omega}_{0}=\dot{\theta}_{0}{z}_{0}$ , ${\omega}_{1}=\dot{\theta}_{0}{z}_{0}+\dot{\theta}_{1}{x}_{1}$ and ${v}_{G_{i}}$ is the velocity of $G_{i}$ ( $i=1$ for the Arm and $i=2$ for the Pendulum). The potential energy of gravity is obtained following the scalar product between the global vertical axis and the center of gravity position vector in the local tilted frame of reference.

V=-m_{1}g({z}\cdot{OG_{1}})-m_{2}g({z}\cdot{OG_{2}})

(4)

The generalized forces due to friction in the joints are written in (5) as the sum of friction torques acting on the links. This model, proposed in [25], is chosen for numerous reasons. It is versatile in modeling different friction phenomenon, its parameters are physically meaningful and not added on an ad hoc basis, and it is continuously differentiable. The latter two properties are required by the proposed algorithm as upper and lower limits on the estimated parameters steady-state values need to be placed in the algorithm, and the time derivative of the friction term is a necessary part in the algorithm formulation.

\begin{split}{f}_{i}(\dot{\theta}_{i})=&F_{n_{i}}\mu_{d_{i}}tanh\left(4\;\frac{\dot{\theta}_{i}}{\dot{\theta}_{t_{i}}}\right)+F_{n_{i}}\cfrac{(\mu_{s_{i}}-\mu_{d_{i}})\;\cfrac{\dot{\theta}_{i}}{\dot{\theta}_{t_{i}}}}{\cfrac{1}{4}\left(\cfrac{\dot{\theta}_{i}}{\dot{\theta}_{t_{i}}}\right)^{2}+\cfrac{3}{4}}\\ &+\mu_{v_{i}}\dot{\theta}_{i}tanh\left(4\;\frac{F_{n_{i}}}{F_{nt_{i}}}\right)\end{split}

(5)

where $F_{n_{i}}$ is the normal force in joint $i$ , which is a function of the system parameters as well as the state $q$ and its derivatives; $\mu_{d_{i}}$ , $\mu_{s_{i}}$ , and $\mu_{v_{i}}$ are the dynamic, static and viscous friction coefficients; $\dot{\theta}_{t_{i}}$ is the transition angular velocity responsible for shaping the stribeck curve; and $F_{nt_{i}}$ is the transition force that activates the viscous friction term.

The final model is given in the following form:

{H(q)}{\ddot{q}}+{B(q,\dot{q})}+{G(q)}={Q(\dot{q})}

(6)

where the generalized inertia matrix is given by

{H(q)}=\begin{bmatrix}\begin{pmatrix}j_{1_{z}}+m_{1}l_{1}^{2}\\ +(m_{2}l_{2}^{2}+j_{2_{y}})\;sin(\theta_{1})^{2}\\ +j_{2_{z}}\;cos(\theta_{1})^{2}\end{pmatrix}&-m_{2}l_{2}L_{1}\;cos(\theta_{1})\\ \\ -m_{2}l_{2}L_{1}\;cos(\theta_{1})&j_{2_{x}}+m_{2}l_{2}^{2}\end{bmatrix}

(7)

the centrifugal and coriolis forces are given by

{B(q,\dot{q})}=\begin{bmatrix}\begin{pmatrix}-(m_{2}l_{2}^{2}-j_{2_{y}}+j_{2_{x}})\;sin(2\theta_{1})\;\dot{\theta}_{0}\;\dot{\theta}_{1}\\ -m_{2}l_{2}L_{1}\;sin(\theta_{1})\;\dot{\theta}_{1}^{2}\end{pmatrix}\\ \\ (m_{2}l_{2}^{2}-j_{2_{y}}+j_{2_{z}})sin(\theta_{1})\;cos(\theta_{1})\;\dot{\theta}_{0}^{2}\end{bmatrix}

(8)

the effect of gravity is given by

{G(q)}=\begin{bmatrix}\begin{pmatrix}m_{1}gl_{1}\;sin(\theta_{0})\;sin(\phi)\\ \\ +m_{2}g\begin{pmatrix}L_{1}\;sin(\theta_{0})\\ -l_{2}sin(\theta_{1})\;cos(\theta_{0})\end{pmatrix}\;sin(\phi)\end{pmatrix}\\ \\ -m_{2}gl_{2}\begin{pmatrix}sin(\theta_{0})\;cos(\theta_{1})\;sin(\phi)\\ -sin(\theta_{1})\;cos(\phi)\end{pmatrix}\end{bmatrix}

(9)

and the generalized forces due to friction are given by

{Q(\dot{q})}=\begin{bmatrix}{f}_{0}(\dot{\theta}_{0})\\ {f}_{1}(\dot{\theta}_{1})\end{bmatrix}

(10)

where the full term is presented in Equation (5).

II-C Setup

The experimental setup is comprised of a Furuta Pendulum and a Dspace 1104 data acquisition system (Figure 1). The system is equipped with two Baumer rotary incremental encoders with a resolution of $40,000$ counts/rev to record the angular positions of the two joints as a function of time. A sampling frequency of $10$ kHz was used to capture the time response of the system. The acquisition time covers the full motion until the links stop moving.

Now that the setup is explained, the preliminaries required for Universal Adaptive Stabilizer (UAS) based parameters estimation are provided in Section III-A. This is followed by description of the proposed parameters identification framework in Section III-B, and the mathematical justification in Section III-C.

III UAS-based High Gain Adaptive Parameter Identification Routine

III-A Universal Adaptive Stabilizer (UAS)

Nussbaum functions are switching functions used in high-gain adaptive control such as with a UAS. They are employed where the control direction is generally unknown. Hence, they are selected as no assumptions and no restrictions on the parameter signs and magnitudes are made. This has been successfully employed in [26, 27] in estimating the parameters of various other systems. Equation (12) shows some examples of Nussbaum functions.

A Nussbaum function is a piecewise right continuous and locally Lipschitz function $N(.):[k^{\prime},\infty)\to\mathbb{R}$ that satisfies

		$\displaystyle\sup\limits_{k>k_{0}}\frac{1}{k-k_{0}}\int_{k_{0}}^{k}N(\tau)\,d\tau=+\infty$		(11)
		and
		$\displaystyle\inf\limits_{k>k_{0}}\frac{1}{k-k_{0}}\int_{k_{0}}^{k}N(\tau)\,d\tau=-\infty$

A Mittag-Leffler function was chosen to act as a Nussbaum function for the purpose of tuning the design parameters. It is of interest to note that a Mittag-Leffler function can act as a Nussbaum function under certain conditions that are documented in [28]. The Mittag-Leffler function depends on two positive real parameters, $\alpha$ and $\beta$ . Authors in [28] determined the conditions under which the Mittag-Leffler function acts as a Nussbaum function, which are if $\alpha\in(2,3]$ and $\beta=1$ . The Nussbaum function of Mittag-Leffler type is given in (12) as $N_{1}(-\lambda z^{\alpha})$ .

N_{1}(z)=\sum_{\gamma=0}^{\infty}\frac{z^{\gamma}}{\Gamma(\alpha\gamma+\beta)},\,\alpha\in(2,3],\,\beta=1

(12)

N_{2}(z)=z\;cos(\sqrt{|z|})

(13)

N_{3}(z)=z^{2}\;cos(|z|)

(14)

N_{4}(z)=cos(\frac{\pi z}{2})\;e^{z^{2}}

(15)

The above functions differ in how fast they vibrate. In other works [27, 26], it was advantageous to choose a rapidly vibrating Nussbaum functions in applications involving DC motor and Li-ion Battery Dynamics. For that reason, the Mittag-Leffler Nussbaum function is used in this work.

Figure 3: Algorithm Flowchart

III-B Parameter Identification

The models in Equation (6) and Equation (5) are rewritten here to setup the parameter identification procedure as the following:

{H(\hat{q})}{\ddot{\hat{q}}}+{B(\hat{q},\dot{\hat{q}})\dot{\hat{q}}}+{G(\hat{q})}={\hat{Q}(\dot{\hat{q}})}+{u}\\

(16)

here ${\ddot{\hat{q}},\dot{\hat{q}},\hat{q}}$ are the estimated acceleration, velocity and position, respectively, and ${u}$ is the UAS input given by Equation (22). Also note that the parameters of $H$ , $B$ , and $G$ in (16) are known, and ${\hat{Q}}=\begin{bmatrix}{\hat{f}_{0}(\dot{\hat{\theta}}_{0})}&{\hat{f}_{1}(\dot{\hat{\theta}}_{1})}\end{bmatrix}^{T}$ is the estimated friction vector whose entries are given in (17), and (18).

\begin{split}\hat{f}_{0}(\dot{\hat{\theta}}_{0})=&F_{n_{0}}\hat{z}_{1}tanh\left(4\;\frac{\dot{\hat{\theta}}_{0}}{\hat{z}_{4}}\right)+F_{n_{0}}\cfrac{(\hat{z}_{2}-\hat{z}_{1})\;\cfrac{\dot{\hat{\theta}}_{0}}{\hat{z}_{4}}}{\cfrac{1}{4}\left(\cfrac{\dot{\hat{\theta}}_{0}}{\hat{z}_{4}}\right)^{2}+\cfrac{3}{4}}\\ &+\hat{z}_{3}\dot{\hat{\theta}}_{0}tanh\left(4\;\frac{F_{n_{0}}}{\hat{z}_{5}}\right)\end{split}

(17)

\begin{split}\hat{f}_{1}(\dot{\hat{\theta}}_{1})=&F_{n_{1}}\hat{z}_{6}tanh\left(4\;\frac{\dot{\hat{\theta}}_{1}}{\hat{z}_{9}}\right)+F_{n_{1}}\cfrac{(\hat{z}_{7}-\hat{z}_{6})\;\cfrac{\dot{\hat{\theta}}_{1}}{\hat{z}_{9}}}{\cfrac{1}{4}\left(\cfrac{\dot{\hat{\theta}}_{1}}{\hat{z}_{9}}\right)^{2}+\cfrac{3}{4}}\\ &+\hat{z}_{8}\dot{\hat{\theta}}_{1}tanh\left(4\;\frac{F_{n_{1}}}{\hat{z}_{10}}\right)\end{split}

(18)

In these equations $\hat{z}_{n}$ represents the frictional model parameter estimates where $z_{1}=\mu_{d_{0}}$ , $z_{2}=\mu_{s_{0}}$ , $z_{3}=\mu_{v_{0}}$ , $z_{4}=\dot{\theta}_{t_{0}}$ , $z_{5}=F_{nt_{0}}$ , $z_{6}=\mu_{d_{1}}$ , $z_{7}=\mu_{s_{1}}$ , $z_{8}=\mu_{v_{1}}$ , $z_{9}=\dot{\theta}_{t_{1}}$ , $z_{10}=F_{nt_{1}}$ , the following equations complete the high-gain universal adaptive observer used for parameters estimation in this work.

e(t)=\dot{q}(t)-\dot{\hat{q}}(t)

(19)

\dot{k}(t)=\left\lVert e(t)\right\rVert^{2}_{2},\;k(t_{0})=k_{0}>0

(20)

N(k(t))=N_{1}(-\lambda k(t)^{\alpha})

(21)

u(t)=H(\hat{q})N(k(t))e(t)

(22)

In Equation (21), $\lambda=1$ , and $\alpha=3$ , thus resulting in a Nussbaum function of the Mittag-Leffler form. In Equation (22) $H$ is the generalized inertia matrix from (7). Equation (23) is used in identifying the friction model parameters $\hat{z}_{n}$ in Equation (17), and Equation (18).

\dot{\hat{z}}_{n}(t)=(\gamma+\lambda_{nu}(z_{nu}-\hat{z}_{n}(t))+\lambda_{nl}(z_{nl}-\hat{z}_{n}(t)))\left\lVert e(t)\right\rVert_{2}

(23)

where $\hat{z}_{n}(t)\in\mathbb{R}$ represents each parameter that needs to be identified, and $n\in{1,2,\dots,N}$ refers to the parameter number where $N=10$ is the number of unknown parameters that need to be identified. Also, $\gamma$ is a tuning parameter to control the rate at which the parameters adapt.

The proposed parameter evolution equation requires steady state upper and lower bounds $z_{nu}$ and $z_{nl}$ , as well as their confidence levels $\lambda_{nu}$ and $\lambda_{nl}$ , respectively. The parameters $\lambda_{nu},\lambda_{nl}>0$ and the parameters $z_{nu},z_{nl},\lambda_{nu},\lambda_{nl}\in\mathbb{R}$ . The upper and lower bounds $z_{nu}$ and $z_{nl}$ are constants which represent the limits within which $\hat{z}_{n}(t)$ is desired to settle, as $t\to\infty$ . The constants $\lambda_{nu}$ and $\lambda_{nl}$ represent the users confidence in their choice of steady state upper and lower bounds $z_{nu}$ and $z_{nl}$ , respectively. The flowchart of the proposed parameter identification routine is presented in Figure 3.

III-C Mathematical Justification

This section presents the mathematical justification behind the proposed friction model parameters identification routine. First, Lemma 1 shows that the adaptive observer’s error heading to zero leads to boundedness of the parameters. Second, Theorem 1 proves that the adaptive observer’s error is driven to zero by the designed UAS input. Finally, Theorem 2 establishes the convergence of the estimated parameters to the true parameters values.

Lemma 1.

Let $\hat{z}_{n}$ be given as in Equation (23). Given that $\lambda_{nu},\lambda_{nl}>0$ , if the parameter adaptation proposed in Equation (23) is used and $e(t)\to 0$ as $t\to\infty$ , then $\hat{z}_{n}$ is bounded $\forall t>t_{0}$ , for $n\in\{1,2,\dots,10\}$ .

Proof of Lemma 1.

Observe that (23) is a stable Linear Time Invariant (LTI) system driven by the input $(\gamma+\lambda_{nu}z_{nu}+\lambda_{nl}z_{nl})\left\lVert e(t)\right\rVert_{2}$ as shown in Equation (24).

\dot{\hat{z}}_{n}(t)+(\lambda_{nu}+\lambda_{nl})\hat{z}_{n}(t)=(\gamma+\lambda_{nu}z_{nu}+\lambda_{nl}z_{nl})\left\lVert e(t)\right\rVert_{2}

(24)

The solution of the above equation can be obtained by the properties of the convolution integral to be

\begin{split}&\hat{z}_{n}(t)=\hat{z}_{n}(t_{0})e^{-(\lambda_{nu}+\lambda_{nl})t}\\ &+(\gamma+\lambda_{nu}z_{nu}+\lambda_{nl}z_{nl})\int_{t_{0}}^{t}\left\lVert e(t-\tau)\right\rVert_{2}e^{-(\lambda_{nu}+\lambda_{nl})\tau}d\tau\end{split}

(25)

Because $\hat{z}_{n}(t),\lambda_{nl},\lambda_{nu},z_{nl},z_{nu}$ are all bounded for $n\in\{1,2,\dots,10\}$ , and because $\left\lVert e(t-\tau)\right\rVert_{2}>0$ for all $t$ , then Equation (25) yields $\hat{z}_{n}(t)\to\hat{z}_{n_{\infty}}$ that is bounded. This can be seen from the first term on the right hand of (25) growing smaller as $t\to\infty$ , and following the assumption that $e\to 0$ as $t\to\infty$ , the second term under the integral is bounded.

∎

With the boundedness of the parameters established in Lemma 1, the assumption that the estimation error grows to zero is now proven in Theorem 1.

Theorem 1.

Let $e(t)=\dot{q}(t)-\dot{\hat{q}}(t)$ as in Equation (19), where the dynamics of $q(t)$ are given by Equation (5) — Equation (10), and the dynamics of $\hat{q}(t)$ are given by Equation (16) — (18), $u$ is as given in Equation (22). If $H(q(t))$ and $H(\hat{q}(t))$ are invertible for all $t\geq t_{0}$ Then, $e(t)\to 0$ as $t\to\infty$ .

Proof of Theorem 1.

Let the system dynamics in Equation (6) and system model in Equation (16) be split into the dynamics and the friction counterparts and written as follows

\begin{split}{\ddot{{q}}}&={H({q})}^{-1}[{{Q}(\dot{{q}})}-{B({q},\dot{{q}})}\dot{{q}}-{G({q})}]\\ {\ddot{\hat{q}}}&={H(\hat{q})}^{-1}[{\hat{Q}(\dot{\hat{q}})}-{B(\hat{q},\dot{\hat{q}})}\dot{\hat{q}}-{G(\hat{q})}]+{H(\hat{q})}^{-1}{u}\\ \end{split}

(26)

where it is understood that all the estimated quantities $\hat{[\quad]}$ vary with time, ${u}$ is as defined in Equation (16), and ${H({q})}^{-1},{H(\hat{q})}^{-1}$ exist by assumption.

The error dynamics $\dot{e}={\ddot{q}}-{\ddot{\hat{q}}}$ can then be expanded and written as

\begin{split}\dot{e}&={H}^{-1}{Q}-\hat{H}^{-1}\hat{{Q}}+\hat{H}^{-1}\hat{{B}}\dot{\hat{q}}-{H}^{-1}{B}\dot{q}\\ &+\hat{H}^{-1}\hat{{G}}-{H}^{-1}{G}-\hat{H}^{-1}{u}\end{split}

(27)

where for ease of use we write $H\equiv H(q),\hat{H}\equiv H(\hat{q}),G\equiv G(q),\hat{G}\equiv G(\hat{q}),B\equiv B({q},\dot{{q}}),\hat{B}=B(\hat{q},\dot{\hat{q}}),Q\equiv Q(q)$ and $\hat{Q}=\hat{Q}(\hat{q})$ .

In Equation (27), we add and subtract $\Lambda e$ , where $\Lambda$ is a positive definite design matrix.

\begin{split}\dot{e}&={H}^{-1}{Q}-\hat{H}^{-1}\hat{{Q}}+\hat{H}^{-1}\hat{{B}}\dot{\hat{q}}-{H}^{-1}{B}\dot{q}\\ &+\hat{H}^{-1}\hat{{G}}-{H}^{-1}{G}-\hat{H}^{-1}{u}+\Lambda e-\Lambda e\\ \dot{e}&={H}^{-1}{Q}-\hat{H}^{-1}\hat{{Q}}+\hat{H}^{-1}\hat{{B}}\dot{\hat{q}}-{H}^{-1}{B}\dot{q}\\ &+\hat{H}^{-1}\hat{{G}}-{H}^{-1}{G}-\hat{H}^{-1}{u}+\Lambda\dot{q}-\Lambda\dot{\hat{q}}-\Lambda e\\ \dot{e}&={H}^{-1}{Q}-\hat{H}^{-1}\hat{{Q}}+\hat{H}^{-1}\hat{{G}}-{H}^{-1}{G}-\hat{H}^{-1}{u}-\Lambda e\\ &-(\Lambda-\hat{H}^{-1}\hat{{B}})\dot{\hat{q}}+(\Lambda-{H}^{-1}{B})\dot{q}\end{split}

(28)

The following notation simplifications are applied to Equation (28).

\begin{split}\dot{e}&=\underbrace{{H}^{-1}{Q}-\hat{H}^{-1}\hat{{Q}}}_{\mathcal{Q}}+\underbrace{\hat{H}^{-1}\hat{{G}}-{H}^{-1}{G}}_{\mathcal{G}}-\hat{H}^{-1}{u}-\Lambda e\\ &-\underbrace{(\Lambda-\hat{H}^{-1}\hat{{B}})}_{\mathcal{B}_{2}}\dot{\hat{q}}+\underbrace{(\Lambda-{H}^{-1}{B})}_{\mathcal{B}_{1}}\dot{q}\\ &=\mathcal{Q}+\mathcal{G}+\mathcal{B}_{1}\dot{q}-\mathcal{B}_{2}\dot{\hat{q}}-\hat{H}^{-1}{u}-\Lambda e\end{split}

(29)

In Equation (29), we add and subtract $\mathcal{B}_{2}\dot{q}$ .

\begin{split}\dot{e}&=\mathcal{Q}+\mathcal{G}+\mathcal{B}_{1}\dot{q}-\mathcal{B}_{2}\dot{\hat{q}}-\hat{H}^{-1}{u}-\Lambda e+\mathcal{B}_{2}\dot{q}-\mathcal{B}_{2}\dot{q}\\ \dot{e}&=\mathcal{Q}+\mathcal{G}+\mathcal{B}_{1}\dot{q}-\mathcal{B}_{2}\dot{q}+\mathcal{B}_{2}(\dot{q}-\dot{\hat{q}})-\hat{H}^{-1}{u}-\Lambda e\\ \dot{e}&=\mathcal{Q}+\mathcal{G}+(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}+\mathcal{B}_{2}e-\Lambda e-\hat{H}^{-1}{u}\\ \dot{e}&=\mathcal{Q}+\mathcal{G}+(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}-(\Lambda-\mathcal{B}_{2})e-\hat{H}^{-1}{u}\\ \end{split}

(30)

To prove the stability of the error dynamics above, pre-multiply (30) by $e^{T}$ to get

\begin{split}e^{T}\dot{e}&=e^{T}\mathcal{Q}+e^{T}\mathcal{G}+e^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}-e^{T}(\Lambda-\mathcal{B}_{2})e\\ &-e^{T}\hat{H}^{-1}{u}\\ \end{split}

(31)

Equation (31) can be rearranged and rewritten as the following inequality taking into account that terms like $e^{T}\mathcal{Q}$ can be used to form $2e^{T}\mathcal{Q}\leq e^{T}e+\mathcal{Q}^{T}\mathcal{Q}$ as explained in Appendix A. Also, The substitution $u=\hat{H}N(k)e$ is used.

\begin{split}e^{T}\dot{e}+e^{T}(\Lambda-\mathcal{B}_{2})e&\leq 3e^{T}e+\mathcal{Q}^{T}\mathcal{Q}+\mathcal{G}^{T}\mathcal{G}\\ &+\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}\\ &-e^{T}\hat{H}^{-1}\hat{H}N(k)e\\ \end{split}

(32)

which is further simplified considering $\dot{k}(t)=\left\lVert e(t)\right\rVert^{2}_{2}=e^{T}e$ and that $N(k)$ is scalar.

\begin{split}e^{T}\dot{e}+e^{T}(\Lambda-\mathcal{B}_{2})e&\leq 3\dot{k}+\mathcal{Q}^{T}\mathcal{Q}+\mathcal{G}^{T}\mathcal{G}-N(k)\dot{k}\\ &+\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}\\ \end{split}

(33)

Integrating both sides of the above system yields

\begin{split}&\int_{t_{0}}^{t}e^{T}\dot{e}\;d\tau+\int_{t_{0}}^{t}e^{T}(\Lambda-\mathcal{B}_{2})e\;d\tau\leq 3\int_{t_{0}}^{t}\dot{k}\;d\tau\\ &+\int_{t_{0}}^{t}\mathcal{Q}^{T}\mathcal{Q}\;d\tau+\int_{t_{0}}^{t}\mathcal{G}^{T}\mathcal{G}\;d\tau-\int_{t_{0}}^{t}N(k)\dot{k}\;d\tau\\ &+\int_{t_{0}}^{t}\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}\;d\tau\\ \end{split}

(34)

Applying a change in variable in the $\int_{t_{0}}^{t}N(k)\dot{k}\;d\tau$ term, we get

\begin{split}&\frac{1}{2}e^{T}e+\int_{t_{0}}^{t}e^{T}(\Lambda-\mathcal{B}_{2})e\;d\tau\leq 3(k(t)-k(t_{0}))\\ &+\int_{t_{0}}^{t}\mathcal{Q}^{T}\mathcal{Q}\;d\tau+\int_{t_{0}}^{t}\mathcal{G}^{T}\mathcal{G}\;d\tau-\int_{k_{0}}^{k}N(k)\;dk\\ &+\int_{t_{0}}^{t}\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}\;d\tau\\ \end{split}

(35)

Dividing both sides by $k(t)-k(t_{0})$ and rearranging yields

\begin{split}\frac{1}{2}\frac{e^{T}e}{k(t)-k(t_{0})}&+\frac{\int_{t_{0}}^{t}e^{T}(\Lambda-\mathcal{B}_{2})e\;d\tau}{k(t)-k(t_{0})}\leq 3\\ &-\frac{1}{k(t)-k(t_{0})}\int_{k_{0}}^{k}N(k)\;dk\\ &+\frac{\int_{t_{0}}^{t}\mathcal{Q}^{T}\mathcal{Q}\;d\tau}{k(t)-k(t_{0})}+\frac{\int_{t_{0}}^{t}\mathcal{G}^{T}\mathcal{G}\;d\tau}{k(t)-k(t_{0})}\\ &+\frac{\int_{t_{0}}^{t}\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}\;d\tau}{k(t)-k(t_{0})}\\ \end{split}

(36)

Dividing both sides by the terms on the R.H.S. except for $\frac{1}{k(t)-k(t_{0})}\int_{k_{0}}^{k}N(k)\;dk$ yields

\begin{split}&\begin{pmatrix}\frac{1}{2}\frac{e^{T}e}{k(t)-k(t_{0})}+\\ \\ \frac{\int_{t_{0}}^{t}e^{T}(\Lambda-\mathcal{B}_{2})e\;d\tau}{k(t)-k(t_{0})}\end{pmatrix}\begin{pmatrix}3+\frac{\int_{t_{0}}^{t}\mathcal{Q}^{T}\mathcal{Q}\;d\tau}{k(t)-k(t_{0})}+\frac{\int_{t_{0}}^{t}\mathcal{G}^{T}\mathcal{G}\;d\tau}{k(t)-k(t_{0})}\\ \\ +\frac{\int_{t_{0}}^{t}\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}\;d\tau}{k(t)-k(t_{0})}]\end{pmatrix}^{-1}\leq\\ &1-\frac{\int_{k_{0}}^{k}N(k)\;dk}{k(t)-k(t_{0})}\begin{pmatrix}3+\frac{\int_{t_{0}}^{t}\mathcal{Q}^{T}\mathcal{Q}\;d\tau}{k(t)-k(t_{0})}+\frac{\int_{t_{0}}^{t}\mathcal{G}^{T}\mathcal{G}\;d\tau}{k(t)-k(t_{0})}\\ \\ +\frac{\int_{t_{0}}^{t}\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}\;d\tau}{k(t)-k(t_{0})}]\end{pmatrix}^{-1}\\ \end{split}

(37)

It can be noted from the above that with the exception of the $\frac{1}{k(t)-k(t_{0})}\int_{k_{0}}^{k}N(k)\;dk$ term on the R.H.S. of Equation (36), all the terms are quadratic forms and positive. Now if the term $k(t)$ is monotonically increasing, i.e. $k(t)\to\infty$ as $t\to\infty$ , then the Nussbaum function term on the R.H.S. can, as a whole, take values approaching $\pm\infty$ by the definition of a Nussbaum function as presented in Equation (11). This violates the inequality above. Consequently, the assumption that $k(t)\to\infty$ as $t\to\infty$ is false and $k(t)$ must, therefore, be bounded. However, $\dot{k}(t)$ being a non-decreasing function by definition in (20) and $k(t)$ being bounded implies that $k(t)\to k_{\infty}$ as $t\to\infty$ . This means that $\dot{k}(t)\to 0$ as $t\to\infty$ i.e. $e^{T}e\to 0$ , or $e\to 0$ . This establishes the required result. ∎

With the estimation error going to zero established in Theorem 1, the convergence of the estimated parameters to the real parameters is now shown in Theorem 2. It is also worth noting that verifying the assumptions that $H(q(t))$ and $H(\hat{q}(t))$ are invertible is simple and can be explicitly checked by numerical computation. Further, before proceeding to the next result, some notation is defined, to facilitate the a condensed representation of the next result and its proof. Thus motivated let $i\in\{0,1\}$ , and $j,k\in\{1,2,\cdots,9,\ell\}$ . Here we use the set $\{1,2,\cdots,9,\ell\}$ and let $\ell=10$ so that in the statement of Theorem 2 the subscript value of 10 can be easily visualized. Now consider the following definitions.

a_{i}(t)=F_{n_{i}}\mu_{d_{i}}tanh\left(4\;\frac{\dot{\theta}_{i}(t)}{\dot{\theta}_{t_{i}}(t)}\right)\\

(38)

\hat{a}_{ijk}(t)=F_{n_{i}}\hat{z}_{j}(t)tanh\left(4\;\frac{\dot{{\theta}}_{i}(t)}{\hat{z}_{k}(t)}\right)

(39)

b_{i}(t)=\mu_{v_{i}}\dot{\theta}_{i}(t)tanh\left(4\;\frac{F_{n_{i}}}{F_{nt_{i}}}\right)\\

(40)

\hat{b}_{ijk}(t)=\hat{z}_{j}(t)\dot{{\theta}}_{i}(t)tanh\left(4\;\frac{F_{n_{i}}}{\hat{z}_{k}(t)}\right)\\

(41)

c_{i}(t)=\cfrac{F_{n_{i}}\mu_{s_{i}}\;\cfrac{\dot{\theta}_{i}(t)}{\dot{\theta}_{t_{i}}(t)}}{\cfrac{1}{4}\left(\cfrac{\dot{\theta}_{i}(t)}{\dot{\theta}_{t_{i}}(t)}\right)^{2}+\cfrac{3}{4}}\\

(42)

\hat{c}_{ijk}(t)=\cfrac{F_{n_{i}}\hat{z}_{j}(t)\;\cfrac{\dot{{\theta}}_{i}(t)}{\hat{z}_{k}(t)}}{\cfrac{1}{4}\left(\cfrac{\dot{{\theta}}_{i}(t)}{\hat{z}_{k}(t)}\right)^{2}+\cfrac{3}{4}}\\

(43)

d_{i}(t)=\cfrac{F_{n_{i}}\mu_{d_{i}}\;\cfrac{\dot{\theta}_{i}(t)}{\dot{\theta}_{t_{i}}(t)}}{\cfrac{1}{4}\left(\cfrac{\dot{\theta}_{i}(t)}{\dot{\theta}_{t_{i}}(t)}\right)^{2}+\cfrac{3}{4}}\\

(44)

\hat{d}_{ijk}(t)=\cfrac{F_{n_{i}}\hat{z}_{j}(t)\;\cfrac{\dot{{\theta}}_{i}(t)}{\hat{z}_{k}(t)}}{\cfrac{1}{4}\left(\cfrac{\dot{{\theta}}_{i}(t)}{\hat{z}_{k}(t)}\right)^{2}+\cfrac{3}{4}}

(45)

Theorem 2.

Considering the definitions in (38)-(45) let

A(t)=\begin{bmatrix}{a_{0}-\hat{a}_{014}}&{b_{0}-\hat{b}_{035}}&{c_{0}-\hat{c}_{024}}&{d_{0}-\hat{d}_{014}}\\ {a_{1}-\hat{a}_{169}}&{b_{1}-\hat{b}_{18\ell}}&{c_{1}-\hat{c}_{179}}&{d_{1}-\hat{d}_{169}}\end{bmatrix},

and $x=\begin{bmatrix}{1}&{1}&{1}&{1}\end{bmatrix}^{T}$ . For the parameter identification problem described in Section III-B, suppose the conditions required for Theorem 1 to hold are satisfied. Suppose there exists a time instant $t^{*}>t_{0}$ such that for $t\to\infty,t>t^{*}$ , $x$ does not belong to the nullspace of $A(t),A(t)\neq 0$ ,and all the following functions: $tanh\left(4\;\frac{\dot{\theta}_{i}(t)}{\dot{\theta}_{t_{i}}(t)}\right)$ , $tanh\left(4\;\frac{\dot{{\theta}}_{i}(t)}{\hat{z}_{k}(t)}\right)$ , $tanh\left(4\;\frac{F_{n_{i}}}{F_{nt_{i}}}\right)$ , $tanh\left(4\;\frac{F_{n_{i}}}{\hat{z}_{k}(t)}\right)$ tend to 1. Then the parameter adaptation law in Equation (23) leads to $\hat{z}_{n}\to z_{n_{true}}$ as $t\to\infty$ for $n\in\{1,2,\dots,10\}$ .

Proof of Theorem 2.

By the assumptions Theorem 1 is satisfied, this gives $e\to 0$ as $t\to\infty$ . As a result $k(t)\to k_{\infty}$ (a constant), which further leads ${u}\to 0$ following the definition of $u$ in (21)-(22). As $e\to 0$ , therefore by definition of $e$ in (19) we get that ${\dot{\hat{q}}}\to{\dot{q}}$ , and ${\ddot{\hat{q}}}\to{\ddot{q}}$ .

Considering this, the model of the system in Equation (6) and the estimator dynamics in Equation (16) can be written as

\begin{split}{H(q)}{\ddot{{q}}}+{B({q},\dot{{q}})}\dot{q}+{G({q})}&={{Q}(\dot{{q}})}\\ {H(\hat{q})}{\ddot{\hat{q}}}+{B(\hat{q},\dot{\hat{q}})\dot{\hat{q}}}+{G(\hat{q})}&={\hat{Q}(\dot{\hat{q}})}+{u}.\\ \end{split}

(46)

Because $u\to 0$ , ${\dot{\hat{q}}}\to{\dot{q}}$ , ${\ddot{\hat{q}}}\to{\ddot{q}}$ , subtracting the lower equation in (46) from the upper one, and substituting $\hat{q}=q$ provides

{Q}(\dot{{q}})-\hat{Q}(\dot{{q}})=\begin{bmatrix}{f}_{0}(\dot{{\theta}}_{0})\\ {f}_{1}(\dot{{\theta}}_{1})\end{bmatrix}-\begin{bmatrix}\hat{f}_{0}(\dot{{\theta}}_{0})\\ \hat{f}_{1}(\dot{{\theta}}_{1})\end{bmatrix}=0.

(47)

Now using the definitions of the friction force model in (5), and the concerned estimates as defined in (17), (18) one can rewrite (47) as (48) using the notation in (38)-(45).

\underbrace{\begin{bmatrix}{a_{0}-\hat{a}_{014}}&{b_{0}-\hat{b}_{035}}&{c_{0}-\hat{c}_{024}}&{d_{0}-\hat{d}_{014}}\\ {a_{1}-\hat{a}_{169}}&{b_{1}-\hat{b}_{18\ell}}&{c_{1}-\hat{c}_{179}}&{d_{1}-\hat{d}_{169}}\end{bmatrix}}_{A(t)}\underbrace{\begin{bmatrix}1\\ 1\\ 1\\ 1\end{bmatrix}}_{x}=0

(48)

We have arrived at (48) because as $t\to\infty$ we have $e\to 0$ . Recall that, by the assumptions of the theorem, there exists a time instant $t^{*}>t_{0}$ such that for $t\to\infty,t>t^{*}$ , $x$ does not belong to the nullspace of $A(t),A(t)\neq 0$ . So for such $t>t^{*}$ from (48) we can write ${\hat{a}_{014}=a_{0}}$ , ${\hat{b}_{035}=b_{0}}$ , ${\hat{c}_{024}=c_{0}}$ , ${\hat{d}_{014}=d_{0}}$ , ${\hat{a}_{169}=a_{1}}$ , ${\hat{b}_{18\ell}=b_{1}}$ , ${\hat{c}_{179}=c_{1}}$ , ${d_{1}=\hat{d}_{169}}$ . Now considering ${\hat{a}_{014}=a_{0}}$ , ${\hat{b}_{035}=b_{0}}$ , ${\hat{c}_{024}=c_{0}}$ , ${\hat{d}_{014}=d_{0}}$ along with their definitions in (38)-(45), the following are obtained.

F_{n_{0}}\hat{z}_{1}tanh\left(4\;\frac{\dot{{\theta}}_{0}}{\hat{z}_{4}}\right)=F_{n_{0}}\mu_{d_{0}}tanh\left(4\;\frac{\dot{\theta}_{0}}{\dot{\theta}_{t_{0}}}\right)

(49)

\hat{z}_{3}\dot{{\theta}}_{0}tanh\left(4\;\frac{F_{n_{0}}}{\hat{z}_{5}}\right)=\mu_{v_{0}}\dot{\theta}_{0}tanh\left(4\;\frac{F_{n_{0}}}{F_{nt_{0}}}\right)

(50)

\cfrac{F_{n_{0}}\hat{z}_{2}\;\cfrac{\dot{{\theta}}_{0}}{\hat{z}_{4}}}{\cfrac{1}{4}\left(\cfrac{\dot{{\theta}}_{0}}{\hat{z}_{4}}\right)^{2}+\cfrac{3}{4}}=\cfrac{F_{n_{0}}\mu_{s_{0}}\;\cfrac{\dot{\theta}_{0}}{\dot{\theta}_{t_{0}}}}{\cfrac{1}{4}\left(\cfrac{\dot{\theta}_{0}}{\dot{\theta}_{t_{0}}}\right)^{2}+\cfrac{3}{4}}\

(51)

As per the assumptions of the theorem, all the $tanh(\cdot)$ terms tend to unity. Considering the equations (49) and Equation (50) with this gives $\hat{z}_{1}=\mu_{d_{0}}$ , $\hat{z}_{3}=\mu_{v_{0}}$ . Now considering (49) again with $\hat{z}_{1}=\mu_{d_{0}}$ and explcitly writing the $tanh(\cdot)$ terms gives

tanh\left(4\;\frac{\dot{{\theta}}_{0}}{\hat{z}_{4}}\right)=tanh\left(4\;\frac{\dot{\theta}_{0}}{\dot{\theta}_{t_{0}}}\right).

(52)

Now from the definition of the hyperbolic tangent function, and from its graph, it is known that its inverse can be found, i.e. so let $p=tanh\left(4\;\frac{\dot{{\theta}}_{0}}{\hat{z}_{4}}\right)=tanh\left(4\;\frac{\dot{\theta}_{0}}{\dot{\theta}_{t_{0}}}\right)$ . Then, $tanh^{-1}(p)=4\;\frac{\dot{{\theta}}_{0}}{\hat{z}_{4}}=4\;\frac{\dot{\theta}_{0}}{\dot{\theta}_{t_{0}}}$ , leading to $\hat{z}_{4}=\dot{\theta}_{t_{0}}$ . Now using this exact same process with (50), considering $\hat{z}_{3}=\mu_{v_{0}}$ and inverting the $\tanh(\cdot)$ terms provides $\hat{z}_{5}=F_{nt_{0}}$ . Further using the fact that $\hat{z}_{4}=\dot{\theta}_{t_{0}}$ and considering (51) gives us $\hat{z}_{2}=\mu_{s_{0}}$ .

The exact same procedure as followed from after (48) to the line above, can be used to explicitly analyze ${\hat{a}_{169}=a_{1}}$ , ${\hat{b}_{18\ell}=b_{1}}$ , ${\hat{c}_{179}=c_{1}}$ , ${d_{1}=\hat{d}_{169}}$ . Using this process provides the following to $z_{6}=\mu_{d_{1}}$ , $z_{7}=\mu_{s_{1}}$ , $z_{8}=\mu_{v_{1}}$ , $z_{9}=\dot{\theta}_{t_{1}}$ , $z_{10}=F_{nt_{1}}$ . This completes the proof. ∎

It is worth mentioning that the assumptions of the above theorem are not necessarily restrictive. This is because during experimentation it has been routinely observed that all the $tanh(\cdot)$ terms frequently approach unity. So it is easy to pick time intervals at which $e\approx 0$ , and $tanh(\cdot)\approx 1$ so that the the values of $\hat{z}_{i}(t),i\in\{1,2,\cdots,10\}$ in such an interval of time can be used as candidates for parameters estimates. Also, once such estimated parameters values are available, they can be fed into (38)-(45) to reconstruct matrix $A$ mentioned in Theorem 2. This will easily allow ascertaining if the nullspace related assumptions in Theorem 2 are met. If they are not met then estimates from another interval in time where $e\approx 0$ , and $tanh(\cdot)\approx 1$ , can be checked for satisfaction of the nullspace assumption.

TABLE II: Friction model parameter identification problem setup.

Parameter	$\mu_{d_{0}}$	$\mu_{s_{0}}$	$\mu_{v_{0}}$	$\dot{\theta}_{t_{0}}$	$F_{nt_{0}}$
$Z_{nl}$	$2.22\times 10^{-16}$	$2.22\times 10^{-16}$	$2.22\times 10^{-16}$	$2.22\times 10^{-16}$	$2.22\times 10^{-16}$
$Z_{nu}$	$0.0750$	$0.0750$	$0.0100$	$0.0100$	$0.1$
$\lambda_{l}$	$50$	$50$	$50$	$50$	$50$
$\lambda_{u}$	$1$	$1$	$1$	$1$	$1$
Parameter	$\mu_{d_{1}}$	$\mu_{s_{1}}$	$\mu_{v_{1}}$	$\dot{\theta}_{t_{1}}$	$F_{nt_{1}}$
$Z_{nl}$	$2.22\times 10^{-16}$	$2.22\times 10^{-16}$	$2.22\times 10^{-16}$	$2.22\times 10^{-16}$	$2.22\times 10{-16}$
$Z_{nu}$	$0.150$	$0.151$	$0.0100$	$0.0100$	$0.100$
$\lambda_{l}$	$50$	$50$	$50$	$50$	$50$
$\lambda_{u}$	$1$	$1$	$1$	$1$	$1$

IV Results and Discussion

In this section, the proposed parameter identification algorithm is tested in simulation, then experimental validation is performed. In all experiments, the links start with a zero initial angular velocity and a specified initial angular position. The motion is monitored until the system returns to the equilibrium position.

IV-A Simulation

Here, the performance of the proposed algorithm is evaluated in simulation. The dynamic equations of motion for the Furuta pendulum were simulated from some initial conditions, and the angular position and velocity states are fed into the algorithm. The initial conditions were $\theta_{0}(0)=0^{\circ}$ , $\theta_{1}(0)=120^{\circ}$ , $\dot{\theta}_{0}(0)=0^{\circ}/sec$ , and $\dot{\theta}_{1}(0)=0^{\circ}/sec$ . The Runge-Kutta method of order $4$ with a sampling time of $0.001s$ was used to propagate the equations in time and solve the ODEs. Normally distributed Gaussian random white noise characterized by $v\sim\mathcal{N}(0,\,\sigma=0.1)\,$ was injected into the simulated states to serve as measurements noise. Results achieved despite this, show the resilience of the proposed approach to noise. The initial conditions fed into the algorithm were also injected with similar noise to account for the difficulty associated with acquiring highly accurate and precise initial conditions. Additionally, this serves to simulate the ramifications of initial conditions mismatch in chaotic systems such as the Furuta pendulum. The adaptation parameters used to test the algorithm are detailed in Table II. The upper/lower parameters limits in Table II were used to form a normal distribution centered around the mean of each respective parameter bounds. The initial friction model parameters were sampled from this distribution to start the algorithm. The actual parameters and the initial guesses used to simulate the response are presented in Table III.

TABLE III: Friction model parameter simulation identification results with the UAS method.

Parameter	$\mu_{d_{0}}$	$\mu_{s_{0}}$	$\mu_{v_{0}}$	$\dot{\theta}_{t_{0}}$	$F_{nt_{0}}$
Initial Guess	$5.135\times 10^{-3}$	$5.875\times 10^{-3}$	$2.531\times 10^{-3}$	$4.853\times 10^{-2}$	$1.029\times 10^{-1}$
Estimate	$4.793\times 10^{-4}$	$5.740\times 10^{-4}$	$2.424\times 10^{-4}$	$4.742\times 10^{-3}$	$9.479\times 10^{-3}$
Actual Value	$5\times 10^{-4}$	$6\times 10^{-4}$	$2.5\times 10^{-4}$	$5\times 10^{-3}$	$10\times 10^{-3}$
Parameter	$\mu_{d_{1}}$	$\mu_{s_{1}}$	$\mu_{v_{1}}$	$\dot{\theta}_{t_{1}}$	$F_{nt_{1}}$
Initial Guess	$5.705\times 10^{-3}$	$6.683\times 10^{-3}$	$2.399\times 10^{-3}$	$4.820\times 10^{-2}$	$9.720\times 10^{-2}$
Estimate	$5.740\times 10^{-4}$	$6.688\times 10^{-4}$	$2.424\times 10^{-4}$	$4.742\times 10^{-3}$	$9.479\times 10^{-3}$
Actual Value	$6\times 10^{-4}$	$7\times 10^{-4}$	$2.5\times 10^{-4}$	$5\times 10^{-3}$	$10\times 10^{-3}$

Figure 4 shows the UAS tracking results where it is noticed that the estimated states quickly approach and converge to the true simulated response of the system. Further, as the estimated states approach the simulated states, the norm of the error is observed to subside in magnitude and approaches zero, and the UAS gain is observed to also grow smaller as a consequence of the reduction in estimation error. Figure 5 shows the time trajectories of the parameters of the friction model as the UAS converges to the correct states of the real system. It is seen that the parameters do indeed stabilize around some final values. The values considered to be the parameter estimates were taken just after the error norm reduced below some user defined threshold, which here is recommended to be below $\left\lVert e(t)\right\rVert_{2}=0.01$ . The identified values are presented in Table III. Figure 6 showcases a validation test of the chosen parameter estimates. The experimental response is plotted against the response realized from both the UAS estimated parameters as well as the initial parameters. The UAS-provided parameter estimates generate a response very similar to the simulated response, which establishes that the proposed approach was successful at identifying high quality parameter estimates in simulation.

Further analysis is presented in frequency domain in Figure 7 where both the phase and spectrum of the time responses of the angular positions of the arm and pendulum links are presented. It is evident that the algorithm is successful at moving from the initial estimates to ones that are more representative of the system. The portrait is quite complex considering the nonlinearity of the system, and it helps explain the compression and expansion effect present in the time responses presented as a wide array of harmonics is present in the response. The spectrum plots show good matching between estimated and actual parameters in this test.

IV-B Experiment

Here, the performance of the proposed algorithm is validated experimentally. The initial conditions were approximately $\theta_{0}(0)=0.2^{\circ}$ , $\theta_{1}(0)=126^{\circ}$ , $\dot{\theta}_{0}(0)=0^{\circ}/sec$ , and $\dot{\theta}_{1}(0)=0^{\circ}/sec$ . As with the test in simulation, the initial friction model parameter estimates were randomly generated also considering the upper/lower limits from Table II. The Runge-Kutta method of order $4$ with a sampling time of $0.0001s$ was also used here to propagate the equations in time and solve the ODEs. The adaptation parameters used to test the algorithm are detailed in Table II.

Figure 8 shows the UAS tracking results where it is noticed that the estimated states quickly approach and converge to the true response of the system. Throughout that process, the norm of the error is observed to subside in magnitude and approach zero, and the UAS gain is observed to follow a similar trend to the estimation error norm. Unlike the simulated case in the previous section, the UAS gain does not completely reduce to zero and intermittent spikes are observed in the response. Figure 9 shows the time trajectories of the parameters of the friction model as the UAS tracks the states of the real system. It is seen that the parameters settle and oscillate within some final range of values, which makes sense following the behavior of the UAS gain. The values considered to be the parameter estimates were taken just after the error norm reduced below some user defined threshold, which here is recommended to be below $\left\lVert e(t)\right\rVert_{2}=0.05$ . Due to the extra uncertainty present in the experimental data, it is recommended to average the parameters estimates during adaptation between the error subsiding and the system motion seizing. This provided tighter bounds on the parameters estimates in the experimental test. The identified values are presented in Table IV. Figure 10 showcases a validation test of the chosen parameter estimates. The experimental response is plotted against the response realized from both the UAS estimated parameters as well as the initial parameter guesses. It is clear that the UAS-provided parameter estimates generate a response that is qualitatively similar to the experiment. However, unlike the simulation test case, an exact match is not realized between the experimental response and the estimated one. We reason this to be due to the limitations of friction models available in the literature in that they only approximate and not fully characterize the damping phenomenon. Nevertheless, the proposed approach proved to be computationally highly efficient at identifying quality parameter estimates using experimental data.

The same analysis is performed in the frequency domain. In Figure 11, the phase and the spectrum of the time responses of the angular positions of the arm and pendulum links are presented. Similar to the experimental results, the algorithm is successful at moving from the initial estimates to the ones that are more representative of the system. However, the spectrum plots do not show as good matching between the estimated and the actual parameters in this test as it is the case with simulation. This is due to the difficulty concerning the description of the actual friction present in the system. Still, very good qualitative parameter estimates were realized using the proposed approach. In the next section, we explore alternative methods to address this problem.

TABLE IV: Friction model parameter experimental identification results with the UAS method.

Parameter	$\mu_{d_{0}}$	$\mu_{s_{0}}$	$\mu_{v_{0}}$	$\dot{\theta}_{t_{0}}$	$F_{nt_{0}}$
Initial Guess	$1.090\times 10^{-2}$	$8.239\times 10^{-3}$	$4.971\times 10^{-3}$	$2.726\times 10^{-3}$	$1.297\times 10^{-2}$
Estimate	$1.575\times 10^{-3}$	$1.575\times 10^{-3}$	$3.006\times 10^{-4}$	$3.006\times 10^{-4}$	$2.065\times 10^{-3}$
Parameter	$\mu_{d_{1}}$	$\mu_{s_{1}}$	$\mu_{v_{1}}$	$\dot{\theta}_{t_{1}}$	$F_{nt_{1}}$
Initial Guess	$1.957\times 10^{-2}$	$1.796\times 10^{-2}$	$1.607\times 10^{-2}$	$1.342\times 10^{-2}$	$9.213\times 10^{-3}$
Estimate	$3.046\times 10^{-3}$	$3.065\times 10^{-3}$	$3.006\times 10^{-4}$	$3.006\times 10^{-4}$	$2.065\times 10^{-3}$

IV-C Comparison with conventional methods

Here, we compare the proposed approach to a conventional optimization-based one. Specifically, a nonlinear grey-box identification in the MATLAB environment is selected to exploit the known structure of the dynamic model and generate estimates for its parameters. The results for all the cases are shown in Figure 12, where the quality of fit is measured through the coefficient of determination $R^{2}$ metric and is presented in Equation 53.

{fit=100\times\left(1-\frac{\sum(\hat{x}-x_{ref})^{2}}{\sum(x_{ref}-\bar{x}_{ref})^{2}}\right)}

(53)

where $x_{ref}$ is the reference signal, and $\hat{x}$ is the estimate of $x_{ref}$ .

Aditionally, the computation time is computed using

\begin{split}&\text{Normalized Computation Time}=\\ &\frac{\text{Computation Time}}{\text{Experiment Time}\times\text{Sampling Frequency}}\end{split}

(54)

Here, the same simulation parameters and initial conditions presented in the previous sections are used. Due to the high nonlinearity of the problem, optimization alone is reportedly not always a viable solution to address this problem. Testing shows that in simulation, optimization alone is not capable of reaching acceptable parameter estimates with a goodness of fit of only $13\%$ and $3\%$ for the Arm Joint and the Pendulum Joint responses, respectively. The proposed algorithm, in contrast, provided quick estimates that provided much better goodness of fit metrics of $98\%$ for both responses. We conclude that the proposed approach is more robust to noisy data than the optimization approach. In the experimental test, the proposed method, while it is able to quickly generate quality estimates, it provides estimates, which can be improved on with a goodness of fit of $77\%$ and $91\%$ for the responses of the Arm Joint and the Pendulum Joint, respectively. The optimization approach, however, exhibited a low performance in terms of goodness of fit with only $3\%$ and $24\%$ for the responses of the Arm Joint and the Pendulum Joint, respectively. That result was also achieved at a considerably higher computational cost with optimization taking approximately $500$ and $1300$ times more computation time than the UAS approach for the simulation and the experiment, respectively. Table V summarizes these results.

TABLE V: Computation time and quality of estimates.
Computation time is presented along with experiment runtime.

		Proposed		Optimization
		Arm Joint	Pendulum Joint	Arm Joint	Pendulum Joint
Sim.	Test Duration [s]	$35$		$35$
Sim.	Computation Time [s]	$3.43$		$4504.6$
	Goodness of fit [%]	$98.97$	$99.53$	$12.77$	$2.80$
Exp.	Test Duration [s]	$20.28$		$20.28$
Exp.	Computation Time [s]	$19.49$		$10373.7$
	Goodness of fit [%]	$77.84$	$91.99$	$3.40$	$24.11$

V Conclusion

This paper presented a dynamic analysis of a passive tilted Furuta pendulum. The pendulum was tilted to ensure the existence of a stable equilibrium configuration with the weight of each link as the only external force applied to the system. A complete analytical model was derived where full inertia of the pendulum is taken into account, and a comprehensive friction model was selected form the literature. The problem of parameter estimation of the friction in the setup was addressed. A UAS-based high gain observer was designed to estimate the friction model parameters, and it was tested in simulation and experimentally. The obtained results showed that the developed algorithm is effective at identifying parameters in this class of systems. Nonetheless, it was also concluded that the friction phenomenon is not fully described by the 5-parameter per joint friction model. Simulation results, even with a high level of injected and propagated dynamic noise, were close to the dynamic behavior of the system. Experimental results were qualitatively good when compared to conventional optimization based parameter estimation strategies. A frequency analysis was performed and it showed that the system is highly nonlinear with multiple harmonics detected. Differences between the simulation and the experiments are mainly attributed to the difficulty associated with characterizing friction in the joints. Future work includes the implementation of nonlinear control strategy using the developed model to showcase the performance a more accurate model of the system would permit.

Appendix A Proof Addendum

Remark 1.

Here we show how the inequalities in the mathematical justification are formed.

Take the difference of quadratic forms below where the terms are the same as those in Equation (32)

(e-(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q})^{T}(e-(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q})\geq 0

(55)

Expanding the above difference gives

\begin{split}e^{T}e-e^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}&-((\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q})^{T}e\\ &+((\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q})^{T}((\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q})\geq 0\end{split}

(56)

Simplifying

\begin{split}e^{T}e-e^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}&-\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}e\\ &+\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}\geq 0\end{split}

(57)

Rearranging

\begin{split}e^{T}e+\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}\geq 2e^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}\end{split}

(58)

where the substitution $e^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}+\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}e=2e^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}$ was made following the fact that $e^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}$ and $\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}e$ are scalar quantities implying $e^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})\dot{q}=(\dot{q}^{T}(\mathcal{B}_{1}-\mathcal{B}_{2})^{T}e)^{T}$ .

The same procedure can be applied on the following quadratic difference form

(e-\mathcal{Q})^{T}(e-\mathcal{Q})\geq 0

(59)

which is modified by multiplying $\mathcal{Q}$ by the $2\times 2$ identity matrix $I_{2\times 2}$ .

(e-I_{2\times 2}\mathcal{Q})^{T}(e-I_{2\times 2}\mathcal{Q})\geq 0

(60)

which leads to the following after applying the same procedure above

e^{T}e+\mathcal{Q}^{T}\mathcal{Q}\geq 2e^{T}\mathcal{Q}

(61)

References

[1] K. Furuta, M. Yamakita, and S. Kobayashi, “Swing-up Control of Inverted Pendulum Using Pseudo-State Feedback,” Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, vol. 206, no. 4, pp. 263–269, nov 1992.
[2] M. Yamakita, T. Hoshino, and K. Furuta, “Control practice using pendulum,” in Proceedings of the American Control Conference, vol. 1. IEEE, 1999, pp. 490–494.
[3] M. Gafvert, “Dynamic model based friction compensation on the Furuta pendulum,” in IEEE Conference on Control Applications - Proceedings, vol. 2. IEEE, 1999, pp. 1260–1265.
[4] A. Wadi, J. H. Lee, and L. Romdhane, “Nonlinear sliding mode control of the Furuta pendulum,” in 11th International Symposium on Mechatronics and its Applications, ISMA 2018, vol. 2018-Janua, mar 2018, pp. 1–5.
[5] S. Jadlovská and J. Sarnovský, “Modelling of classical and rotary inverted pendulum systems - A generalized approach,” Journal of Electrical Engineering, vol. 64, no. 1, pp. 12–19, jan 2013.
[6] O. García-Alarcon, S. Puga-Guzman, and J. Moreno-Valenzuela, “On parameter identification of the Furuta pendulum,” Procedia Engineering, vol. 35, pp. 77–84, 2012.
[7] B. S. Cazzolato and Z. Prime, “On the dynamics of the Furuta pendulum,” Journal of Control Science and Engineering, vol. 2011, pp. 1–8, 2011.
[8] M. Antonio-Cruz, R. Silva-Ortigoza, J. Sandoval-Gutiérrez, C. A. Merlo-Zapata, H. Taud, C. Márquez-Sánchez, and V. M. Hernández-Guzman, “Modeling, simulation, and construction of a furuta pendulum test-bed,” in 25th International Conference on Electronics, Communications and Computers, CONIELECOMP 2015. IEEE, feb 2015, pp. 72–79.
[9] M. Antonio-Cruz, R. Silva-Ortigoza, C. A. Merlo-Zapata, M. G. Villarreal-Cervantes, D. Muñoz-Carrillo, and V. M. Hernández-Guzmán, “Modeling and Construction of a Furuta Pendulum Prototype,” in Proceedings - 2014 IEEE International Conference on Mechatronics, Electronics, and Automotive Engineering, ICMEAE 2014. IEEE, nov 2014, pp. 98–103.
[10] L. T. Aguilar, I. Boiko, L. Fridman, and A. Ferreira, “Identification based generation of self-excited oscillations for underactuated mechanical systems via Two-relay algorithm,” in IEEE 10th International Workshop on Variable Structure Systems, VSS’08. IEEE, jun 2008, pp. 41–46.
[11] L. Freidovich, A. Shiriaev, F. Gordillot, F. Gomez-Estern, and J. Aracil, “Partial-energy-shaping control for orbital stabilization of high frequency oscillations of the Furuta pendulum,” in Proceedings of the IEEE Conference on Decision and Control. IEEE, 2007, pp. 4637–4642.
[12] K. Schlee, S. Gangadharan, J. Ristow, C. Hubert, J. Sudermann, and C. Walker, “Modeling and parameter estimation of spacecraft fuel slosh mode,” in Proceedings - Winter Simulation Conference, vol. 2005. IEEE, 2005, pp. 1265–1273.
[13] G. Habib, A. Miklos, E. T. Enikov, G. Stepan, and G. Rega, “Nonlinear model-based parameter estimation and stability analysis of an aero-pendulum subject to digital delayed control,” International Journal of Dynamics and Control, vol. 5, no. 3, pp. 629–643, aug 2017.
[14] A. Wadi, J. H. Lee, and L. Romdhane, “Dynamic Analysis of the Tilted Furuta Pendulum,” in MATEC Web of Conferences, vol. 104, 2017.
[15] E. I. Drewniak, G. D. Jay, B. C. Fleming, and J. J. Crisco, “Comparison of two methods for calculating the frictional properties of articular cartilage using a simple pendulum and intact mouse knee joints,” Journal of Biomechanics, vol. 42, no. 12, pp. 1996–1999, aug 2009.
[16] C. Makkar, W. E. Dixon, W. G. Sawyer, and G. Hu, “A new continuously differentiable friction model for control systems design,” in IEEE/ASME International Conference on Advanced Intelligent Mechatronics, AIM, vol. 1. IEEE, 2005, pp. 600–605.
[17] J. J. Crisco, J. Blume, E. Teeple, B. C. Fleming, and G. D. Jay, “Assuming exponential decay by incorporating viscous damping improves the prediction of the coefficient of friction in pendulum tests of whole articular joints,” Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, vol. 221, no. 3, pp. 325–333, mar 2007.
[18] H. T. Teixeira, V. S. De Mattos Siqueira, and C. J. Munaro, “Closed-loop quantification and compensation of friction in an inverted pendulum,” Journal of Control, Automation and Electrical Systems, vol. 24, no. 6, pp. 794–805, sep 2013.
[19] H. B. Fang and J. Xu, “Dynamics of a three-module vibration-driven system with non-symmetric Coulomb’s dry friction,” Multibody System Dynamics, vol. 27, no. 4, pp. 455–485, feb 2012.
[20] L. Freidovich, A. Robertssont, A. Shiriaev, and R. Johansson, “Friction compensation based on LuGre model,” in Proceedings of the IEEE Conference on Decision and Control
[21] E. Pennestrì, V. Rossi, P. Salvini, and P. P. Valentini, “Review and comparison of dry friction force models,” Nonlinear Dynamics, vol. 83, no. 4, pp. 1785–1801, nov 2016.
[22] A. Keck, J. Zimmermann, and O. Sawodny, “Friction parameter identification and compensation using the ElastoPlastic friction model,” Mechatronics, vol. 47, pp. 168–182, nov 2017.
[23] F. Marques, P. Flores, J. C. Claro, and H. M. Lankarani, “Modeling and analysis of friction including rolling effects in multibody dynamics: a review,” Multibody System Dynamics, vol. 45, no. 2, pp. 223–244, aug 2019.
[24] C. Hu and F. Wan, “Parameter identification of a model with coulomb friction for a real inverted pendulum system,” in 2009 Chinese Control and Decision Conference, CCDC 2009. IEEE, jun 2009, pp. 2869–2874.
[25] P. Brown and J. McPhee, “A Continuous Velocity-Based Friction Model for Dynamics and Control with Physically Meaningful Parameters,” Journal of Computational and Nonlinear Dynamics, vol. 11, no. 5, jun 2016.
[26] D. Ali, S. Mukhopadhyay, H. Rehman, and A. Khurram, “UAS based Li-ion battery model parameters estimation,” Control Engineering Practice, vol. 66, pp. 126–145, 2017.
[27] H. M. Usman, S. Mukhopadhyay, and H. Rehman, “Permanent magnet DC motor parameters estimation via universal adaptive stabilization,” Control Engineering Practice, vol. 90, pp. 50–62, sep 2019.
[28] Y. Li and Y. Chen, “When is a Mittag–Leffler function a Nussbaum function?” Automatica, vol. 45, no. 8, pp. 1957–1959, 2009.
[29] L. Ljung, System Identification. Boston, MA: Birkhäuser Boston, 1998.