ACC2021.tex 58.3 KB
Newer Older
Fnadi Mohamed's avatar
Fnadi Mohamed committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2345678901234567890123456789012345678901234567890123456789012345678901234567890
%        1         2         3         4         5         6         7         8

\documentclass[letterpaper, 10 pt, conference]{ieeeconf}  % Comment this line out if you need a4paper

%\documentclass[a4paper, 10pt, conference]{ieeeconf}      % Use this line for a4 paper

\IEEEoverridecommandlockouts                              % This command is only needed if 
                                                          % you want to use the \thanks command

\overrideIEEEmargins                                      % Needed to meet printer requirements.

%In case you encounter the following error:
%Error 1010 The PDF file may be corrupt (unable to open PDF file) OR
%Error 1000 An error occurred while parsing a contents stream. Unable to analyze the PDF file.
%This is a known problem with pdfLaTeX conversion filter. The file cannot be opened with acrobat reader
%Please use one of the alternatives below to circumvent this error by uncommenting one or the other
%\pdfobjcompresslevel=0
%\pdfminorversion=4

% See the \addtolength command later in the file to balance the column lengths
% on the last page of the document

% The following packages can be found on http:\\www.ctan.org
% The following packages can be found on http:\\www.ctan.org
\usepackage{graphics} % for pdf, bitmapped graphics files
\usepackage{epsfig} % for postscript graphics files
\usepackage{mathptmx} % assumes new font selection scheme installed
\usepackage{times} % assumes new font selection scheme installed
\usepackage{amsmath} % assumes amsmath package installed
\usepackage{amssymb}  % assumes amsmath package installed
\usepackage{amsmath,bm}
\usepackage{blindtext}
\usepackage{lipsum}
\usepackage{mathtools}
\usepackage{relsize}
\usepackage{cuted}

\usepackage{stmaryrd}
\usepackage[mathscr]{eucal}
%\usepackage{arevmath}     % For math symbols

\usepackage{graphics} % for pdf, bitmapped graphics files
\newcommand{\dd}{
	\mathop{}\mathopen{}\mathrm{d}
}
\usepackage{epsfig} % for postscript graphics files
\usepackage{mathptmx} % assumes new font selection scheme installed
\usepackage{times} % assumes new font selection scheme installed
\usepackage{amsmath} % assumes amsmath package installed
\usepackage{amssymb}  % assumes amsmath package installed
\usepackage{amsfonts,bm}
\usepackage{stackengine}
\usepackage{stfloats}

\usepackage{algorithm}
\usepackage{algorithmic}



%\usepackage{graphics} % for pdf, bitmapped graphics files
%\usepackage{epsfig} % for postscript graphics files
%\usepackage{mathptmx} % assumes new font selection scheme installed
%\usepackage{times} % assumes new font selection scheme installed
%\usepackage{amsmath} % assumes amsmath package installed
%\usepackage{amssymb}  % assumes amsmath package installed

\title{\LARGE \bf
%Dynamic Identification of Nonlinear Inverted Pendulum with Quantified Uncertainties via Interval Analysis
Guaranteed Identification of Viscous Friction for a Nonlinear Inverted Pendulum  Through Interval Analysis and Set Inversion
}


\author{Mohamed Fnadi$^{1}$,  Julien Alexandre dit Sandretto$^{1}$, Gabriel Ballet$^{1}$,  and Laurent Fribourg$^{2}$% <-this % stops a space
%\thanks{*This work was not supported by any organization}% <-this % stops a space
\thanks{$^{1}$ENSTA Paris, Institut Polytechnique de Paris, U2IS laboratory, 828 bd des mar\'echaux, 91762
	Palaiseau Cedex, France
        {\tt\small \{firstname.lastname\}@ensta-paris.fr}}%
\thanks{$^{2}$LSV,  Ecole Normale Sup\'erieure (ENS)  \& CNRS, Paris-Saclay, France
        {\tt\small fribourg@lsv.ens-cachan.fr}}%
}


\begin{document}



\maketitle
\thispagestyle{empty}
\pagestyle{empty}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}

This paper focuses on the guaranteed identification  of viscous friction parameters for a nonlinear inverted pendulum. The method is based on the interval analysis (IA) and set-inversion tools to determine the set of all the feasible friction parameters from a prior domain of interest, i.e. initial interval vector or box, that are consistent with all the experimental and theoretical datasets including their uncertainties. The capabilities of our proposed guaranteed identification are  compared with the more  commonly used   approach based on the least square  method identification (LSMI), which is used especially to adjust the inertial and geometric parameters of our experimental plant. Both of them have been investigated through several  experiments on a real inverted pendulum and simulations with uncertain ODEs  via the DynIbex library.

%The validation and evaluation of our identification strategies is performed via simulations and experiments using a nonlinear inverted pendulum.

\end{abstract}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{INTRODUCTION}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Control of  autonomous systems is often designed relying on their physical or mathematical models (i.e. dynamic and/or kinematic models), which are more specifically  derived as a nonlinear  ordinary differential equations (ODEs) related to some uncertain/unknown parameters, such as inertial and frictional contact coefficients. On the whole, the control of such kind of systems should particularly be guaranteed and robust towards all the uncertainties coming from the system modeling and manufacturing as well. That's why the guaranteed estimation of the model parameters is highly required  such that the coverage rate between the physical reality and the model is as high as possible. On the other hand, friction phenomenon  occurs in several  industrial and mechanical  systems (e.g. transmission, bearings,  soil/tire contacts, etc.), leading these controllers to be less accurate and less reliable.  Once again, it is becoming increasingly obvious  to characterize  this phenomenon in order to minimize its undesirable influence (e.g. heat, instability, energy waste). As at now, however, the friction cannot be modeled perfectly, it is   described in general by an empirical model which is globally based on experiments and observations \cite{c2x2}\cite{c20}. 







%Control of  autonomous systems is often designed relying on their physical or mathematical models (i.e. dynamic and/or kinematic models), which are more specifically  derived as a nonlinear  ordinary differential equations (ODEs). On the whole, most of them contain many uncertain/unknown parameters, such as inertial and frictional contact coefficients.  This makes  the efficiency of these control approaches highly dependent on their estimation. Basically,  friction phenomenon  occurs in several  industrial and mechanical  systems (e.g. transmission, bearings,  soil/tire contacts, etc.), leading their controllers to be less accurate and less reliable.  Once again, it is becoming increasingly obvious  to characterize  this phenomenon in order to minimize its undesirable influence (e.g. heat, energy consumption, instability). As at now, however, the friction cannot be modeled perfectly, it is   described in general by an empirical model which is globally based on experiments and observations \cite{c2x2}\cite{c20}. 


%Control of  autonomous systems is often designed using either a dynamic and/or a kinematic models, most of them use frictional contact parameters. This makes the efficiency of these control approaches highly dependent on the estimation of these parameters. Indeed, friction phenomenon  occurs in several  industrial and mechanical  systems (e.g. gearboxes, transmission, bearings,  soil/tire contacts, etc.), leading their controllers to be less accurate and less reliable.  Therefore, it is becoming increasingly obvious  to characterize  this phenomenon in order to minimize its undesirable influence (e.g. heat, energy consumption, instability). As at now, however, the friction cannot be modeled perfectly, it is   described in general by an empirical model which is globally based on experiments and observations \cite{c2x2}\cite{c20}. 

To properly control robotic devices, the inertial and friction parameters have to be adequately identified  in order to adapt their controllers and achieve the high accuracy and stability. However, many studies consider that the friction parameters are neglected and/or  defined   arbitrarily \cite{c23}. These strategies cannot guarantee the system safety and stability if these parameters are not well-defined. On the other hand,  the friction compensation    can  enhance  the  effectiveness  of the  control  laws. In the literature, many works on the identification procedures were proposed. For instance, authors in \cite{c10} introduce a method to adjust the friction in the chain of transmission from the motor to the robot axis. This method is based on the minimization of the input error, i.e. error between the measured and the computed torques. Moreover, \cite{c11} compares the performances of two kinds of friction models for a linear drive with ball-screw (LuGre model against   the static, viscous and Stribeck friction with hysteresis) which are identified by applying the recursive identification method as well as the least square method identification (LSMI).  In addition, a frequency domain identification technique for  the dynamic LuGre friction model is presented in \cite{c100}. Nevertheless, the quality of the estimation by these deterministic algorithms cannot be ascertained, i.e. the estimated variables cannot reflect the real friction, since they assume that uncertainties, related to the system's design and modeling as well as the degree of correctness of measuring systems, are omitted. Consequently, an estimation in a guaranteed way is needed to account for all the uncertainties involved in the hardware and software parts. 


% do not account for uncertainties related to the system's design and modelling, as well as the on-board sensors' errors. Consequently, an estimation in a guaranteed way is needed to account for uncertainties related to the hardware and software parts.


 More recently, many researchers have explored the possibility of using the interval analysis  strategies (IA) for parameter estimation in order to overcome the limits of the existing approaches. In fact,  IA  tools turn out to be more promising and reliable to manage all the noises and uncertainties contained in the experimental and theoretical data,  contrary to these classical approaches that are generally based on the global minimization of some cost functions. The most well-known example of such kind of application is the bounded error estimation method, which mainly based on the IA and set-inversion techniques \cite{c7}. This method is extremely popular in several research fields due to its  efficiency and reliability, such as in the robot localization from noisy landmarks \cite{c777}, tuning the gains of a controller, and kinetic parameters estimation in elector-chemistry \cite{c15}, etc. The algorithm, named SIVIA (Set Inversion Via Interval Analysis), is much more beneficial and widely used to determine the set of all acceptable parameters which match with the collected data and their uncertainties, by minimizing the error between the available output data and model output \cite{c155}\cite{c14}. As another example, two methodologies (set-inversion and Taylor series)   were used and combined in \cite{l1} and \cite{l2} to improve the algorithm's computation time for  uncertain dynamic systems in the context of the parameter estimation  via bounded measurement error.  The SIVIA algorithm is also, by the way, applied in the framework of this paper   (brief resume is given in the Appendix B).
 
In the short term, we are looking forward to developing new guaranteed and robust Nonlinear Model Predictive Control (NMPC) to stabilize the inverted pendulum in the upright position,  which ultimately  need an accurate parameters estimation. To address these objectives, the novelty of this paper comes from the guaranteed identification of viscous friction parameters of a new inverted pendulum (manufactured in our lab) via IA and set inversion techniques, and its application to uncertain nonlinear ODEs to evaluate the coverage ratios   in such a way their guaranteed solutions (given as tubes at each time-step) fit to the  physical reality. Besides that, the LSMI is initially deployed for multiples experiments in  order to approximate with uncertainties, as closely as possible, the mechanical parameters of our experimental device (geometric and inertial). Then, the SIVIA algorithm is used  to find the set of all possible viscous friction coefficients from an initial research domain of parameters that are consistent with the predicted dynamic model output and  the system output.  

%All the results of this presented paper will serve as a solid  basis for designing  a new guaranteed and robust Nonlinear Model Predictive Control (NMPC), inspired from \cite{c2}\cite{c3}.
 
%  while considering all the uncertainties included in the system's design and the embedded sensors (represented by intervals)
  
% Typically,  the performance of our guaranteed identification is evaluated and analyzed using the uncertain higher-order ODEs describing the system's dynamic such that their guaranteed solutions (given as tubes at each time-step) match with the measurements provided by the real system. This constitutes also one of the main contributions of this presented paper. All the results of this paper will serve as a solid  basis for designing of a new guaranteed and robust Nonlinear Model Predictive Control (NMPC) \cite{c2}\cite{c3}.
 
 
% In order to approximate with uncertainties, as closely as possible, the mechanical parameters of our experimental device (geometric and inertial), . Then, the SIVIA algorithm is used  to find all the guaranteed viscous friction parameters (solutions of the set-inversion problem) from an initial research domain of parameters. Our approach takes into account all the uncertainties included in the system's design and the embedded sensors (represented by intervals). Typically,  the performance of our guaranteed identification is evaluated and analyzed using the uncertain higher-order ODEs describing the system's dynamic such that their guaranteed solutions (given as tubes at each time-step) match with the measurements provided by the real system. This constitutes also one of the main contributions of this presented paper. All the results of this paper will serve as a solid  basis for designing of a new guaranteed and robust Nonlinear Model Predictive Control (NMPC) \cite{c2}\cite{c3}.
 
%  which is based on the IA and set inversion tools
%
%(solutions of the set-inversion problem) from an initial research domain of parameters.

% The main contribution of this paper lies on the guaranteed identification of all possible viscous friction coefficients in such a way the higher-order ODEs solution, with respect to the uncertain parameters, remain consistent with their corresponding measurements. 
% 
%Typically,  the performances of our guaranteed identification is evidenced by the uncertain higher-order ODEs describing the system's dynamic such that their guaranteed solutions (presented as tubes at each time-step) match with the measurements provided by the real system.
%
%by the comparisons between the nonlinear ODE solutions and their real measurements.
% 
%% the set-inversion algorithm is employed in \cite{l1} and \cite{l2} for the guaranteed parameter estimation of an uncertain nonlinear dynamic systems in a bounded measurement error framework.
% 
%
%The main contribution of this paper lies on the identification of viscous friction parameters of an ordinary differential equations (ODEs), which describes the dynamic model of our experimental plant (nonlinear inverted pendulum). This method makes use of the set-inversion algorithm based on IA tools to determine the feasible parameters compatible with the ODEs and the available measurements datasets. The results of this paper will serve as the basis for the design of a new guaranteed and robust Nonlinear Model Predictive Control (NMPC). This controller will be synthesized based on our previous works \cite{c2}\cite{c3}. 
%
%The results of this paper  will serve as a basis for the achievement of these goals.
%
%Our motivation is to find an accurate identification of viscous friction parameters for our experimental plant (nonlinear inverted pendulum) in a guaranteed way via IA and set-inversion approaches. While there are significant uncertainties in the system's design and modeling, the friction variables should be estimated precisely and accurately to compensate the friction effects. In fact, the precise estimation of these parameters is very useful for the design of robust and guaranteed controllers and/or observers that deserves further study. Shortly we plan to synthesize a new guaranteed Nonlinear Model Predictive Control (NMPC) relying on \cite{c2}\cite{c3}. To do this, we are going to need  highly the present  work.


%This point will be treated in our future works using a new guaranteed Nonlinear Model Predictive Control strategy (NMPC) \cite{c2}\cite{c3}.

%In other words, it can  actually   guarantee the system stability and achieve high accuracy during the control process.

The remainder of this paper is structured as follows: Firstly, an overview of the dynamic model of our inverted pendulum is presented. Secondly, section \ref{Sec2} introduces the identification approaches. Thirdly, experimental results and discussions are reported in section \ref{Sec3}.  Finally, section \ref{Sec4} closes the paper with conclusion and future works.


%the  efficiency  of  the proposed guaranteed  identification compared to the traditional LSMI  is illustrated by experiments. Finally, section \ref{Sec4} closes the paper with conclusion and future works.
% using a nonlinear inverted pendulum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{RECALL ON THE DYNAMIC MODELING OF AN INVERTED PENDULUM}
\label{Sec1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Robotic Device}

 
The inverted pendulum is a  very simple system with two serial rotational joints (active and passive). We have designed and manufactured one in our laboratory (see Fig. \ref{Pendule}). Our experimental device  is equipped with a brushless DC motor to actuate the first joint,  and two incremental encoders to measure the rotational angles of both the motor shaft and the passive joint. Its kinematics is sketched in Fig. \ref{Pendule}. This device has to be identified taking into consideration  all the errors and uncertainties related to the system modeling and data measurements. 

\begin{figure}[!h]
	\vspace{-0.42cm}
	\begin{center}
		\includegraphics[scale=0.8]{Figures/SchemCine.pdf}
			\vspace{-0.45cm}
		\caption{[Left] Definition of links frames (configuration $q_1=q_2=0^\circ$). [Right] Representation of the real rotary inverted pendulum ($S_0$: Pendulum housing, $S_1$: horizontal arm, $S_2$: pendulum arm, $S_3$: balancing mass).}
		\label{Pendule} 
	\end{center}
\end{figure}

\vspace{-0.71cm}
\subsection{Dynamic Modeling of the Inverted Pendulum}
  As described in \cite{c4}, the equations of motion are given by,
  \begin{equation}\label{DynEqu}
  \Gamma={M} \ddot{{q}}+{B}({q}, \dot{{q}}) +{Q}(q),
  \end{equation}
  where  $q=[q_1,q_2]$ is the joints' positions of the inverted pendulum, $\Gamma= [\Gamma_1,\Gamma_2]^T$ is the vector of torques applied at the joints. Since the second joint is passive,  $\Gamma_2$ depends mainly on friction variables.  ${M}\in \mathbb{R}^{2 \times 2}$ is the inertia matrix, ${B}\in\mathbb{R}^{2 \times 1}$ is the Coriolis and centrifugal matrix computed applying the Christoffel symbol, ${Q}\in \mathbb{R}^{2\times 1}$ is the gravity forces, and all these matrices are defined as, 
  \begin{equation*}
  \begin{aligned}
  M(q) =\left[\begin{array}{ccc} \mu_1\sin(q_2)^2+\mu_2&& \mu_3 \cos(q_2)  \\\\
  
  \mu_3 \cos(q_2)  && \mu_4\end{array}\right],
  \end{aligned}
  \end{equation*}
   \begin{equation*}
  \begin{aligned}
{B}({q}, \dot{{q}})=
  \left[\begin{array}{c}   - \mu_3\sin(q_2)\dot{q}_2^2 + 2\mu_1\cos(q_2)\sin(q_2)\dot{q}_1\dot{q}_2
  \\\\ 
   - \mu_1\cos(q_2)\sin(q_2)\dot{q}_1^2  \end{array}\right],
  \end{aligned}
  \end{equation*}
       \begin{equation*}
  \mu_1 = l_p^2\left[ \dfrac{m_p}{4}+M\right],  \hspace{0.4cm}
  \mu_2=l_a^2[m_p+M]+\dfrac{J_m}{N_g^2}+J_a,
  \end{equation*}
  \begin{equation*}
  \mu_3 = l_pl_a\left[ \dfrac{m_p}{2}+ M,\right], \hspace{0.4cm}
  \mu_4 =  l_p^2 \left[ \dfrac{m_p}{4}+M\right]+J_p,\hspace{0.15cm} 
  \end{equation*}
     \begin{equation*}
  \begin{aligned}
  {Q}({q}, \dot{{q}})=
  \left[\begin{array}{c} 0
  \\\\ 
   \mu_g\sin(q_2) \end{array}\right],
  \end{aligned}
  ~~~\mu_g=\left[\dfrac{m_p}{2}+M\right]l_pg,
  \end{equation*}
  where $m_a$ and $J_a$ denote the horizontal arm's mass and its inertia, respectively. Similarly, for the pendulum arm with the mass $m_p$ and inertia $J_p$. $M$ is the mass of the load attached to the pendulum arm, $J_m$ is the DC motor inertia and $N_g$ its gear ratio. $l_a$, $l_p$ and $r_0$ are the device's dimensional parameters (as defined in Fig. \ref{Pendule}) and $g$ is the gravitational acceleration.
  
  The external torque $\Gamma$ should be expressed taking into account friction conditions. The relationship between friction and velocity is well modeled in the literature \cite{c2x2}\cite{c20}. In this paper, we use a simple friction model, in which it is related to the friction viscous terms which are  proportional to the joint speed $\dot{q}$,
  \begin{equation}
  \begin{aligned}
  \Gamma =\left[\begin{array}{c} \tau -f_{v_1} \dot{q}_1  \\\\ 
  -f_{v_2} \dot{q}_2\end{array}\right],
  \end{aligned}
  \end{equation}
  where $\tau$ is the motor's torque and $f_{v_i}$ is the  viscous Colombe friction coefficient of joint $i$. These coefficients are often considered as constant values in the majority of literature works. Here the main goal is to identify them in a guaranteed way by reducing the error between the theoretical model  output and measured data output, i.e. they must belong to a certain subsets that are consistent with the model and sensors outputs. Finally, the dynamic model (\ref{DynEqu}) can simply be written as follows, 
   \begin{equation}\label{DynEqu2}
 y_m=f(q,\dot{q}, \ddot{q}, p),
  \end{equation}
 where $y_m=[\tau,0]^T$ is the model output and ${p}=[\mu_1,\mu_2,\mu_3,\mu_4,\mu_g,f_{v_1},f_{v_2}]\in  \mathbb{R}^{n_p=7}$ are the model parameters to be identified. %\mathbb{P} \subset
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{IDENTIFICATION APPROACHES}
\label{Sec2}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The purpose of this section is to propose an identification strategy for the unknown vector ${p}$ embedding all the model   parameters.
%$${p}=[\mu_1,\mu_2,\mu_3,\mu_4,\mu_g,f_{v1},f_{c1},f_{v2},f_{c2}]\in \mathbb{R}^{n_y=9}$$
%where $n_p$ is the number of the identified parameters. 
Assuming the availability of an experimental datasets, measured at each  sampling time, which contains the joints' position, velocity and acceleration as well as the  torque applied by the DC motor. 

%All these data provided by the embedded sensors are gathered  as  interval vectors (e.g, boxes or Cartesian product of intervals)) taking into account the sensors' accuracy, i.e, $[q]\in\left[[0,2\pi]\times[-\pi,\pi]\right]^{n_y}$, $[\dot{q}]\in\mathbb{IR}^{2n_y}$, $[\ddot{q}]\in\mathbb{IR}^{2n_y}$, $[y]\in\mathbb{IR}^{2n_y}$, where $n_y$ is the size of measurements performed on the dynamic system and $\mathbb{IR}^n$ denote the set of all boxes of $\mathbb{R}^n$. The identification approach is  focused on the set inversion and interval analysis techniques to find  the friction parameters that satisfy the equation (\ref{DynEqu2}).

\subsection{Torque, speed and acceleration computation }

Because our experimental datasets are noised, the signals are processed by applying an eighth order  low-pass  Butterworth filter. To identify the needed parameters,  the input of the system (torque signal) and the joints' velocity and acceleration are required. The velocity and the acceleration are computed using a non causal derivative filter. The function \textit{firpm} of Matlab allows us to synthesize this kind of filters to determine the suitable filter's parameters (impulse response coefficients). The second stage is the actual application of the designed filter to the  input data ($q_1$ and $q_2$), using  \textit{filtfilt} command, enabling to process data  in both the forward and reverse directions (zero-phase distortion).






On the other hand, the torque generated at the base of the rotary arm (i.e. at the load gear) can be expressed as a function of the angular velocity $\omega$ as below,
\begin{eqnarray}\label{torque}
\tau=\dfrac{N_mN_g\eta_m\eta_g(\omega_0-\omega)}{R_v}
\end{eqnarray}
where $N_m$ and $N_g$ are respectively the motor and gear ratios, $\eta_m$ and $\eta_g$ are the motor and gearbox efficiencies, $\omega$ and $\omega_0$ are the instantaneous and no-load angular velocities, respectively, and $R_v$ is the motor speed-torque constant.

\subsection{Identification with classical least square method}

%Because of its simplicity and effectiveness to optimize the error between the model estimated output and measured output data 

The least square method identification (LSMI) is commonly used in robotics and control engineering  for parameter identification  due to its simplicity  to optimize the error between the model estimated output and measured output data (e.g. \cite{c10} and \cite{c11}). Since the inverse dynamic model (\ref{DynEqu2})  can be written in a linear way according  to the needed parameters $p$,  LSMI can be used to adjust and approximate the inertial parameters and friction coefficients. Then, the dynamic model (\ref{DynEqu2})  can be rewritten as follows:

\begin{equation}
\begin{aligned}
y_m=D_{}(q, \dot{q}, \ddot{q}) p
\end{aligned}
\end{equation}
where   $D \in\mathbb{R}^{2\times n_p} $ is the matrix of the measured values, expressed as
a function of joint positions $q$, velocities $\dot{q}$  and joint acceleration vector $\ddot{q}$. It is defined as, 
\small
$$D_{}(q, \dot{q}, \ddot{q})=\left[\begin{array}{ccccccc} \sin(q_2)^2\ddot{q}_1+2cos(q_2)sin(q_2)\dot{q}_1\dot{q}_2 & \ddot{q}_1 & \ldots \\\\ 
-\cos(q_2)\sin(q_2)\dot{q}_1^2 & 0 & \ldots  \end{array}\right.$$
$$\left.\begin{array}{cccccc} \cos(q_2)\ddot{q}_2 - \sin(q_2)\dot{q}_2^2 & 0& 0 & \dot{{q}}_1  & 0 \\\\ 
\cos(q_2)\ddot{q}_1 & \ddot{q}_2 & \sin(q_2) & 0 &\dot{{q}}_2  \end{array}\right],$$

\normalsize
To identify the inertial and geometric parameters, we apply the LSMI method to the augmented form relying wholly on measurement data, 
\begin{equation}
\begin{aligned}
Y=W p+\rho
\end{aligned}
\end{equation}
where   $W=[D(1),D(2),\ldots, D(n_y)]^T \in\mathbb{R}^{(2\times n_y) \times n_p} $ is the global measured matrix, i.e. observability  matrix, $n_y$ is the size of measurements performed on the dynamic system.  $Y= [y(1),\ldots, y(n_y)]^T \in\mathbb{R}^{(2\times n_y) \times 1} $ is the output signals, and $\rho$ denotes the error vector between the real system data $Y$ and the output of the dynamic model $Y_m=[y_m(1),\ldots, y_m(n_y)]^T$.

The optimal solution $\widehat{p}$ can be expressed as follows, 
\begin{equation}
\begin{aligned}
\widehat{p}=\underset{p}{\operatorname{\arg\min}}\|Y-Wp\|_{_2}=\left(W^{T} W\right)^{-1} W^{T} Y=W^{+} Y
\end{aligned}
\end{equation}
where $\|{a}\|_{_2}$ denotes the second euclidean norm of the vector $a$, and  $W^{+}$ denote the Moore-Penrose inverse matrix of $W$.
In practice, it is recommended to use at least $500\times n_p$  measures to identify correctly the needed parameters by LSMI.  


\subsection{Guaranteed identification with uncertainty quantification}

Our motivation is to adjust in a better way the viscous friction parameters of the dynamic model to fit the datasets with certain uncertainties. That is why we make use of the most effective tools of  interval analysis (IA) and set inversion. Our guaranteed identification method is inspired from the work presented in \cite{c7} which is based on the bounded-error estimation with uncertain independent variables (a short review of IA and set inversion techniques is presented in Appendices A and B).

All available data are gathered as  interval vectors (i.e. boxes or Cartesian product of intervals) taking into account the sensors' accuracy, i.e. $[\mathbf{q}]\in\left[[0,2\pi]\times[-\pi,\pi]\right]^{\times n_y}$, $[\mathbf{\dot{q}}]\in\mathbb{IR}^{n_y\times 2}$, $[\mathbf{\ddot{q}}]\in\mathbb{IR}^{n_y\times 2}$, $[\mathbf{y}]\in\mathbb{IR}^{ n_y\times 2}$, where  $\mathbb{IR}^n$ denotes the set of all boxes of $\mathbb{R}^n$. The numerical values of the corresponding intervals have been computed easily by  adding a bounded uniform error interval to each collected component data, this error is assumed to belong to a known box.

As already said, the guaranteed identification approach is  focused on the set-inversion and IA techniques to find the feasible sets of   the viscous friction coefficients that satisfy the equation (\ref{DynEqu2}). As a result, the set of all parameters $p\in\left[\mathbf{p}\right]$ that are coherent with the $i^{th}$ measurements is,
%\in \mathbb{R}^{n_p} 
\begin{eqnarray}\label{P_set}\small
\begin{split}
\mathbb{P}_{i}&=\Big \{ \left. p\in[\mathbf{p}]~\mid ~\exists q{(i)} \in[\mathbf{q}](i),~ \exists \dot{q}(i) \in[\mathbf{\dot{q}}](i),~\exists \ddot{q}(i) \in[\mathbf{\ddot{q}}](i)\right.\\
&~~~~~~~~~~~~~~~~~~\left.\text { s.t. }~~ f(q(i),\dot{q}(i),\ddot{q}(i),p)~\in~ [\mathbf{y}](i) ~ \right. \Big \}
\end{split}
\end{eqnarray}
where $[\mathbf{q}](i)$, $[\mathbf{\dot{q}}](i)$, $[\mathbf{\ddot{q}}](i)$ and $[\mathbf{y}](i)$ correspond to the measurement components expressed as intervals at time-step $t_i$.

Then, the feasible set of accepted parameters $\mathbb{P}$ can be derived as,
\begin{eqnarray}\label{P_sett}
\small
\begin{split}
\mathbb{P}
%&=\Big\{p\in[\mathbf{p}] ~\mid ~\exists q(1) \in[q](1),~ \exists \dot{q}(1) \in[\dot{q}](1),~\exists \ddot{q}(1) \in[\ddot{q}](1)\\
%&\left.\text { s.t. }~~ f(q(1),\dot{q}(1),\ddot{q}(1),p)-y(1)=0,~ y(1)\in [y](1)\right. ~~\&\\
%%& \&\\
%&~~~~~~~~~~~~~~ \vdots ~~~~~~~~~\vdots~~~~~~~~~\vdots ~~~~~~~~~\vdots~~~~~~~~~\vdots~ \\
%&~~~\&~\exists q(n_y) \in[q](n_y),~ \exists \dot{q}(n_y) \in[\dot{q}](n_y),~\exists \ddot{q}(n_y) \in[\ddot{q}](n_y)\\
%&\text { s.t. }~~ f(q(n_y),\dot{q}(n_y),\ddot{q}(n_y),p)-y(n_y)=0,~ y(n_y)\in [y](n_y) \Big \},\\
%\vspace{0.1cm}
& = \mathbb{P}_{1}\cap \mathbb{P}_{2}\ldots \cap \mathbb{P}_{n_y}
 = \bigcap_{i=1}^{n_y} \mathbb{P}_i
\end{split}
\end{eqnarray}

The constraints (\ref{P_sett}) can equivalently be written in a matrix format depending on the measurement data (expressed here as interval vectors) and model output,
 \begin{equation}\small
\begin{split}
&~~~~~~~\underbrace{\left[\begin{array}{c}  y_m(1)\\\\ \vdots \\\\ 
	y_m(n_y)\end{array}\right]}_{Y_m}
 \in
 \underbrace{\left[\begin{array}{c}[\mathbf{y}](1) \\\\ \vdots \\\\ \left[\mathbf{y}\right](n_y)\end{array}\right]}_{[\mathbf{Y}]},\\\\
& ~\text{with} ~~y_m(i)=f(q(i),\dot{q}(i),\ddot{q}(i),p),~~\forall i\in\llbracket 1,n_y\rrbracket
\end{split}
\label{inclusion}
\end{equation}

The feasible set $\widehat{\mathbb{P}}$ is a set-inversion problem since the reciprocal image of the set $[\mathbf{Y}]$ must be computed. The SIVIA algorithm, introduced in \cite{c7} and recalled in the Appendix B, is employed to find this desired domain $\widehat{\mathbb{P}}$.
\begin{eqnarray}\label{P_settt}\small
\widehat{\mathbb{P}}=f^{-1}([\mathbf{Y}])\cap[\mathbf{p}]
\end{eqnarray}
% \begin{equation}\small
%\begin{split}
%\left[\begin{array}{c}f(q(1),\dot{q}(1),\ddot{q}(1),p) \\\\ \vdots \\\\ 
%f(q(n_y),\dot{q}(n_y),\ddot{q}(n_y),p)-y(n_y)\end{array}\right] \in
%\left[\begin{array}{c}[y](1) \\\\ \vdots \\\\ \left[y\right](n_y)\end{array}\right],
%\end{split}
%\end{equation}



%\begin{eqnarray}\label{P_settt}
%\begin{split}
%\mathbb{P}&=\big\{p\in \mathbb{R}^{9} ~\mid ~\exists \textbf{q} \in[\textbf{q}],~ \exists \dot{\textbf{q}} \in[\dot{\textbf{q}}],~\exists \ddot{\textbf{q}} \in[\ddot{\textbf{q}}]\\
%&~~~~~~\left.\text { s.t. }~~ \textbf{y}_m-\textbf{y}=0,~~ \text{with}~ \textbf{y}\in [\textbf{y}]\right. \big \},
%\end{split}
%\end{eqnarray}

%Finally, the contracted feasible set of adequate variable can be obtained by interesting all the sets $\mathbb{P}_i$. The algorithm \ref{algo1} shows the mean steps to compute the paving of needed parameters.

\begin{figure*}[b]
	\vspace{-0.6cm}
	\begin{minipage}[c]{.46\linewidth}
		\begin{center}
			\centering
			\includegraphics[scale=0.21]{Figures/q.png}\hspace{-0.5cm}
			\includegraphics[scale=0.21]{Figures/dq.png}\\\vspace{-0.1cm}
			\includegraphics[scale=0.21]{Figures/ddq.png}\hspace{-0.5cm}
			\includegraphics[scale=0.21]{Figures/tau.png}\\\vspace{-0.1cm}
			\caption{ Uncertainty rectangles associated to the measurements of acquisition \#1: (a) angular positions; (b) Angular velocities; (c) Angular accelerations; (d) DC motor torque.}
			\label{Experiments1}
		\end{center}
	\end{minipage} \hfill
	\vspace{-0.6cm}
	\begin{minipage}[c]{.46\linewidth}
		\begin{center}
			\hspace{-0.3cm}	\centering
			\includegraphics[scale=0.21]{Figures/q2.png}\hspace{-0.5cm}
			\includegraphics[scale=0.21]{Figures/dq2.png}\\\vspace{-0.1cm}
			\includegraphics[scale=0.21]{Figures/ddq2.png}\hspace{-0.5cm}
			\includegraphics[scale=0.21]{Figures/tau2.png}\\\vspace{-0.1cm}
			\caption{ Uncertainty rectangles associated to the measurements of acquisition \#2: (a) angular positions; (b) Angular velocities; (c) Angular accelerations; (d) DC motor torque.}
			\label{Experiments2}
		\end{center}
	\end{minipage}
\end{figure*}


%\begin{algorithm}
%	\caption{Compute the paving of the friction Parameters}
%	\begin{algorithmic}[1]
%		\REQUIRE $\mathbb{P}_0=[-3,3]\times[-3,3],~ \mathbb{P}=[]\times[]$, $~[q],~[\dot{q}],~[\ddot{q}],~[y_m]$%, \mathcal{P}:=\emptyset .$
%%		\ENSURE $y = x^n$ %\COMMENT{put some comments here}
%		\FOR{$i = 1$ to $n_y$}		
%%		\STATE {\COMMENT{//$m\leq n^2$}}
%		\STATE $\mathbb{P}_{[if]} \gets f_{Fwd}([q](i),[\dot{q}](i),[\ddot{q}](i),p)$ 
%		\STATE $\mathbb{P}_{[ib]} \gets f_{Bwd}([q](i),[\dot{q}](i),[\ddot{q}](i),p)$
%		\STATE $\mathbb{P}_{i} \gets \mathbb{P}_{[if]}\cap \mathbb{P}_{[ib]}$
%		\STATE Push $\mathbb{P}_i$ onto $\mathbb{P}$
%		\ENDFOR
%		\STATE $\mathbb{P} \gets \mathbb{P}_{1}\cap \mathbb{P}_{2}\ldots \cap \mathbb{P}_{n_y}$
%	\end{algorithmic}
%\label{algo1}
%\end{algorithm}


%This contractor is typically used where we have a set of contractors that result from measurements (each measurement enforces a constraint), some of which can be incorrect.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{EXPERIMENTS}
\label{Sec3}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Identification  results}
All the experiments were performed with our nonlinear inverted pendulum shown in Fig. \ref{Pendule}. It is controlled in an opened-loop to identify its dynamic variables using a Phidget Motor Control.  This control card is connected to our laptop via a USB cable to control the motor's velocity, using C language. To allow high accurate identification, two data acquisition were recorded,   with significant dynamics, at the sampling period of $16ms$. Thus,  we have implemented a bang-bang strategy to send an exciting impulse to the motor's controller. In fact, by appropriately switching   the sign of the pivot velocity values, the bang-bang motion is obtained.  All the recorded information during these tested scenarios ($q(t)$, $\dot{q}(t)$, $\ddot{q}(t)$ and $\tau(t)$) are drawn in Fig. \ref{Experiments1} (acquisition \#1) and Fig. \ref{Experiments2} (acquisition \#2). 

However, from a practical point of view, data collection must take into account the measurement inaccuracy. Indeed, the embedded and low-cost sensors are  unfortunately not always realistic, i.e. they may fail along data collection and send an erroneous measurements. Therefore, the uncertainty associated to each input and output data is characterized as an interval by adding a bounded and known error boxes to each of them. After several trials and attempts on the real system, these  error boxes were quantified and determined for each measured component. The resulting prior intervals for both tested scenarios are depicted in  Fig. \ref{Experiments1} and \ref{Experiments2}.



% Many trails on the real system were performed to  quantify these error boxes.  
%In order to excite the dynamics of the inverted pendulum,

%($[q_1]=[q_1-\delta_1,q_1-\delta_1]$, $[\dot{q}_1]=[\dot{q}_1-\delta_2,\dot{q}_1-\delta_2]$, $[\ddot{q}_1]=[\ddot{q}_1-\delta_3,\ddot{q}_1-\delta_3]$).   

%Since the first acquisition 


Moreover, the motor and gearbox are of type \textit{MDP-Maxon}, their characteristics are provided by the manufacturer. Here is a recap of their properties: the motor inertia  $J_m=9.06e^{-7}kg.m^2$, its no-load speed $w_0=11700~ tr.min^{-1}$, the slope speed-torque $Rv=34000~tr.min^{-1}.N^{-1}.m^{-1}$, the motor and gearbox efficiencies $\eta_m=90.6\%$ and $\eta_g=90\%$ and gear ratios $N_m=5.2$ and $N_g=5$.

\begin{center}
	\begin{table}[h!]
		\vspace{-0.05cm}
		\caption{Identification results by LSMI method.}
		\centering
		\begin{tabular}{lccc}
			Symbol & Values with  \#1 & Values with  \#2 & CAD Values \\
			\hline
			$\mu_1[kg.m^2]$ & $1.04082e^{-3}$ & $1.06171e^{-3}$ & $1.18152e^{-3}$\\ 
			$\mu_2[kg.m^2]$ &$2.56211e^{-3}$ &  $2.37309e^{-3}$ & $2.55064e^{-3}$\\ 
			$\mu_3[kg.m^2]$ & $8.21e^{-4}$ &  $8.29e^{-4}$ & $7.59e^{-4}$\\
			$\mu_4[kg.m^2]$ & $1.31781e^{-3}$  &   $1.19171e^{-3}$ & $1.18152e^{-3}$ \\ 
			$\mu_g[kg.m^2.s^{-2}]$ &  $7.4175e^{-2}$ & $7.2591e^{-2}$ &  $7.3575e^{-2}$\\ 
			$f_{v_1}[N.m.s]$ & $8.13e^{-2}$  & $7.02e^{-2}$ & $-$\\
			$f_{v_2}[N.m.s]$ & $6.42e^{-4}$  &  $5.93e^{-4}$ &  $-$\\  \hline
		\end{tabular}
		\label{Table_LSIM}
	\end{table}
\end{center}

\vspace{-0.45cm}
The  parameters identified via LSMI are reported in Table \ref{Table_LSIM}. Otherwise, some of these physical parameters can be deduced from the CAD  model  (Computer-Aided Design - \textit{SolidWorks}) which are  also summarized in the same Table \ref{Table_LSIM}. As it can be noticed, most of the theoretical values calculated  in both tested cases are very close, which is a good indication of the model consistency. Furthermore, the  identified inertial and geometric parameters are almost similar to those provided by the CAD tool with regard to the mechanical parts of the system. To quantify the error between the measured and identified parameters ($\mu_1$, $\mu_2$, $\mu_3$, $\mu_4$ and $\mu_g$), we can add some uncertainties to these variables which are approximated through the identification errors, as displayed in Table \ref{Tab1_uncetainties}.  From that perspective, the main aim is to identify  the sets of adequate friction parameters  $[f_{v_1}]$ and $[f_{v_2}]$, which are consistent with these uncertainties and satisfy the inclusion  given in (\ref{inclusion}).



\begin{center}
	\begin{table}[h!]
		\vspace{-0.3cm}
		\caption{Inertial and geometric parameters with uncertainties.}
		\centering
		\begin{tabular}{l c c}
			Symbol & Mean value & Uncertainty \\
			\hline
			$\mu_1[kg.m^2]$ & $1.09468e^{-3}$ & $\pm 12 \%$ \\ 
			$\mu_2[kg.m^2]$ & $2.49528e^{-3}$ & $\pm 12\%$\\ 
			$\mu_3[kg.m^2]$ & $8.03e^{-4}$ &  $\pm 10\%$ \\ 
			$\mu_4[kg.m^2]$ & $1.23035e^{-3}$  &  $\pm 12\%$ \\  
			$\mu_g [kg.m^2.s^{-2}]$ & $7.3447e^{-2}$  &  $\pm 15\%$\\ \hline
		\end{tabular}
		\label{Tab1_uncetainties}
	\end{table}
\end{center}

		\vspace{-0.3cm}
From the initial  search domain $[\mathbf{p}]=[-1.5,1.5]^{\times 2}$, the SIVIA algorithm generates for each studied cases the subpavings depicted in Fig. \ref{Vibes}, that enable to bracket the posterior acceptable set for viscous friction coefficients between the inner and outer approximation. The SIVIA accuracy was taken as $\epsilon=0.01$. The red boxes (earth color) have been proved to belong to the feasible sets $\widehat{\mathbb{P}}$ and the blue ones (ocean color) have been considered to create an empty intersection with $\widehat{\mathbb{P}}$. finally, the undermined sets are materialized by the yellow boxes (beach color). As far as can be determined from Fig. \ref{Vibes}, the centered feasible box, in both cases,  is with bounds  $[f_{v1}]\times [f_{v2}]=[0.043012,0.13002]\times[0.000454,0.001174]$.

% These subpavings, by the way, have been easily drawn using the Visualizer for Intervals and Boxes (VIBEs) \cite{c16}
% (red boxes are feasible sets, yellow boxes are undetermined,  and the blues ones are the infeasible sets)
\vspace{-0.4cm}
\begin{figure}[!h]
	\centering
	\includegraphics[scale=0.56]{Figures/Siva1.pdf}\hspace{0.cm}
	\includegraphics[scale=0.56]{Figures/Siva.pdf}\hspace{-0.3cm}
	\caption{ Pavings generated by SIVIA algorithm to show the feasible set $\widehat{\mathbb{P}}$  for the friction parameters. [Left] : scenario \#1. [Right] : scenario \#2.}
	\label{Vibes}
\end{figure}

\vspace{-0.55cm}
\subsection{Identification results validation}

To demonstrate the usefulness of our identification, the cross-validation method is initially applied on an additional trajectory not used in the identification procedure. From the parameters obtained by the LSMI and IA methods, the torques applied at the horizontal arm  are computed using (\ref{DynEqu2}) (red and blue lines shown in Fig. \ref{Torques}),  and compared to the measured one (black line Fig. \ref{Torques}). As can be expected, there are a good similarity and an agreement between these curves. One way to evaluate this is by computing the root-mean-square error (RMSE) that can be calculated by the equation (\ref{RMSE}).  The RMSE is around  $6.4\%$ with LSMI against $2.6\%$ with IA method, which indicates a  reliable fit and a good coherence when the variables uncertainties are accounted. 

  \begin{equation}\small
 \mathrm{RMSE}=\sqrt{\dfrac{1}{n} \sum_{i=1}^{n}\left(\tau_{i}-\widehat{\tau}_{i}\right)^{2}}
 \label{RMSE}
 \end{equation}
 

 
 



 Furthermore, we can solve the Ordinary Differential Equations (ODE) so as to compute the angular positions using the estimated parameters. The dynamic model (\ref{DynEqu2}) can be formulated as a nonlinear  ODE given by,
\begin{equation}\label{ODE}
\left\{\begin{array}{l}
\dot{{x}}(t)=g(t, {x}(t),u(t)) \\
{x}(0) \in\left[\mathbf{{x}}_{0}\right] \subseteq {\mathbb{IR}}^{4}
\\
{u}(0) \in\left[\mathbf{{u}}_{0}\right] \subseteq {\mathbb{IR}}
\end{array}\right.
\end{equation}
where $x=[x_1,x_2,x_3,x_4]^T=[q_1,\dot{q}_1,q_2,\dot{q}_2]^T$ is the state vector  and $u(t)=\tau$ is the input variable. The function $g: \mathbb{R} \times \mathbb{R}^{4}\times \mathbb{R} \rightarrow \mathbb{R}^{4}$ is computed from the inverse dynamic model (\ref{DynEqu}), i.e. $\ddot{{q}}={M}^{-1}[\Gamma-[{B}({q}, \dot{{q}})+{Q}(q)]]$. The sets $\left[\mathbf{{x}}_{0}\right]$ and $\left[\mathbf{{u}}_{0}\right]$, both expressed as boxes, are the initial conditions to model uncertainties related to the state and input vectors. Since the function $g$ is continuous in time and globally Lipschitz in state and input vectors, the equation (\ref{ODE}) has a unique solution at each time-step. 



To further illustrate the quality of our identification,  the sets of the numerical solutions of (\ref{ODE}) are calculated, i.e. the set of possible tight enclosures $[\mathbf{x}_1],\ldots, [\mathbf{x}_n]$ at a sequence of time-instants  $t_1,\ldots,t_n$ with respect to the initial conditions $[\mathbf{x}_0]$ and $[\mathbf{u}_0]$. To do so, our library DynIbex is employed to solve this ODE equation in a guaranteed way, which is mainly based on IA and Runge-Kutta schemes \cite{c12}. 
\vspace{-0.4cm}
\begin{figure}[!h]
	\begin{center}
		\includegraphics[scale=0.37]{Figures/validation/Torques.png}\hspace{-0.3cm}
		\vspace{-0.3cm}
		\caption{Comparison between measured and estimated torques via IA and LSMI methods.}
		\label{Torques} 
	\end{center}
\end{figure}

\vspace{-0.3cm}
Fig. \ref{Experiments_valid1} shows the measured   pendulum angle $q_2$ (black lines), as well as the ones calculated with the identified parameters with LSMI  (red lines) and IA  (blue lines) approaches.  All tests are conducted under the same settings and conditions (e.g. same  initial conditions, same inputs, etc.). The simulation duration is fixed to $4s$ and the precision of $10^{-7}$. We can notice that the simulated $q_2$ by the sophisticated solver DynIbex, using IA and LSMI parameters, are quite similar to the real measured pendulum angle. Nonetheless, the simulated $q_2$ with IA parameters comply as closely as possible with  the measured one at different initial conditions. This is due to the fact that all the system variables are considered as intervals characterizing uncertainties instead of constant values.
\begin{figure}[!h]
	\vspace{-0.4cm}
	\begin{center}
		\centering
		\includegraphics[scale=0.2245]{Figures/validation/q2_comp100.png}\hspace{-0.6cm}
		\includegraphics[scale=0.2245]{Figures/validation/q2_comp.png}\\\vspace{-0.1cm}
		\includegraphics[scale=0.2245]{Figures/validation/q2_compTorque.png}\hspace{-0.6cm}
		\includegraphics[scale=0.2245]{Figures/validation/q2_comp4.png}\\\vspace{-0.1cm}
		\caption{ Comparison between the real and simulated pendulum angular position (represented as tubes) using DynIbex library at different initial conditions: (a) $x_0=[0,0,-100^\circ,0]^T$ and $u=0Nm$; (b) $x_0=[0,0,50^\circ,0]^T$ and $u=0Nm$; (c)  $x_0=[0,0,-90^\circ,0]^T$ and $u=0.17Nm$ and (d) $x_0=[0,0,-45^\circ,0]^T$ and $u=-0.1Nm$.}
		\label{Experiments_valid1}
	\end{center}
\end{figure}


% For each case, same conditions are maintained for the real system and the ODE solver to get compared results. 
\vspace{-0.5cm}
In order to assess the potential  capabilities of  the predicted model, the  associated coverage rates can be computed.  We evaluate the number of cases where the measurements data are included in the simulated tubes, and then the percentage of inclusion is calculated for each tested scenarios. In other words, it is calculated by checking at each time-step whether the time variable related to the measurement is included in the simulated time interval $t_i \in [\mathbf{t_i}]$. If so, we inspect if the corresponding measurement $q_i$ is included in the simulated tube $[\mathbf{q_i}]$ given by DynIbex library. Then, the coverage rate can be approximated to $\frac{N_q}{N_t} \times 100$, where $N_t$ is  the number of times that each sampled time is incorporated in the simulated time interval, and $N_q$ is the one related to the second condition regarding $q_i$. The computed values are recapped in  Table \ref{coverge_rate}. By the way, this criterion is  draconian to test which explain the small computed percentages. This is due to the timing lag   between the measurement and estimated model. Even this issue, the coverage ratios obtained using IA parameters indicates a good  match with the measurements process. It also confirms that the model is estimated with high precision when the uncertainties are accounted.

\begin{center}
	\begin{table}[h!]
		\vspace{-0.3cm}
		\caption{Degree of coverage between the model and physical reality.}
		\centering
		\begin{tabular}{l c c c c}
			Scenario &  (a) & (b) &  (c) &  (d) \\
			\hline
			With IA parameters & $46\%$ & $61\%$ & $38\%$ & $25\%$ \\ 
			With LSMI parameters & $41\%$ & $34\%$ & $20\%$ & $21\%$ \\  \hline
		\end{tabular}
		\label{coverge_rate}
	\end{table}
\end{center}

\vspace{-0.25cm}
To sum up, the experimental results provided by the guaranteed identification of viscous friction parameters are more relevant than the ones obtained by the classical LSMI. Indeed, the estimation approach based on IA and set inversion techniques allows us  to successfully identify these coefficients as intervals taking into account all the errors and uncertainties related to the sensors and the system modeling. Despite of its proven effectiveness to more closely approximate the model output to the real output, it  still has some drawbacks. The main ones are related globally to the computation time, which depends mainly on the number of parameters and the initial domain. This issue can be enhanced with a contraction approach \cite{c22}.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{CONCLUSION AND FUTURE WORK}
\label{Sec4}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This paper describes a guaranteed identification of viscous friction parameters for our new experimental device (nonlinear inverted pendulum). This method is compared to the traditional approach least square method identification  (LSMI), which is principally used to approximate the system's inertial and geometric variables. From a practical perspective, this guaranteed approach  leads to consider all the system uncertainties coming, mainly, from its design and modeling, including also the collected data errors. In order to deal with these uncertainties, we exploit in this paper the interval analysis tools (IA)  and set inversion technique (SIVIA algorithm) to determine the acceptable set of viscous friction coefficients. These latter are adjusted in a feasible way  such that the theoretical value of the
model is consistent with the collected measurement data and their uncertainties. Finally, we show  experimental results obtained by implementing both methods  in the real inverted pendulum, as well as simulation results under DynIbex library, which is applied  to simulate the high-order ordinary differential equations (ODEs) detailing the behavior of our system, this constitutes also one of the main contributions of this paper. The results prove that the capabilities of the guaranteed identification are much more interesting and promising than the deterministic LSMI to approximate the mathematical model to the physical reality.



The  next  work  will  focus  on  the design of a new guaranteed and robust controller to stabilize the pendulum arm in the upright position. Relying on the dynamic model, a constrained nonlinear model predictive controller (NMPC) will be synthesized and applied in a guaranteed way  so as to anticipate future changes in set-points and handle all the constraints, which are critical and necessary for the safety and stability of the system. This work will be principally  inspired by our previous works presented in \cite{c2}\cite{c3}.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{APPENDIX}
In this appendix, we remind some  mathematical notions related to IA, used particularly to solve  linear and/or nonlinear problems in a guaranteed way. For further details on this subject, we encourage  the readers to refer to \cite{c18} and \cite{c19}. 
\subsection{Interval Analysis : the main notations and definitions}
\begin{itemize}
	\item A scalar interval of $\mathbb{R}$ is denoted as $[x]=\left[\underline{x},\overline{x}\right]$ where $\underline{x}$ is the lower bound and $\overline{x}$ is the upper bound. Any interval  of $\mathbb{R}$ is finite, closed and connected subsets.
	\item An interval vector, or a box $\left[\mathbf{x}\right]$ of $\mathbb{R}^n$  is a Cartesian
	product of $n$ intervals, i.e. $ \left[\mathbf{x}\right] =[x_1]\times \ldots \times [x_n]$. The set of all boxes of $\mathbb{R}^n$ is denoted by $\mathbb{IR}^n$.
	\item The width $w(\left[\mathbf{x}\right])$ of a box $\left[\mathbf{x}\right]$ is the length of its largest side.
	\item The interval function $[f]$ from $\mathbb{IR}^n$ to $\mathbb{IR}^m$, is an inclusion function of $f$ if $\forall[\mathbf{x}] \in \mathbb{I} \mathbb{R}^{n}, \quad[f]([\mathbf{x}]) \supseteq\{f(x), x \in[\mathbf{x}]\}$.
	\item A subpaving of $\mathbb{R}^{n}$ is a list of non-overlapping boxes of $\mathbb{R}^{n}$, enabling to bracket any compact set between a list of inner and outer subpavings.
\end{itemize}




\subsection{Set Inversion and SIVIA algorithm}
If $f: \mathbb{R}^{n} \rightarrow \mathbb{R}^{m} \text { and } {Y} \subset \mathbb{R}^{m}$. A set inversion problem is described by the reciprocal image of a set $Y$ by the function $f$ that can be written as, 
\begin{equation}
	\mathbb{P}=\left\{{p} \in \mathbb{R}^{n} \mid {f}({p}) \in {Y}\right\}={f}^{-1}({Y})
\end{equation}

The set $\mathbb{P}$ should be computed in a guaranteed way relying on the IA tools. 


\begin{figure}[!h]
	\vspace{-0.3cm}
	\begin{center}
		\includegraphics[scale=0.7]{Figures/SIVIA.pdf}\hspace{-0.3cm}
			\vspace{-0.4cm}
		\caption{Illustration of the  principal of the set inversion problem - Inclusion tests: infeasible (blue), undetermined (green), feasible (red).}
		\label{Siv} 
	\end{center}
\end{figure}

	\vspace{-0.6cm}
For this purpose, SIVIA algorithm is a proficient set-inversion problem solver via IA. It allows us to get the  set $\mathbb{P}$ in an effective and a guaranteed way, through successive and recursive bisections of the prior interest domain. To obtain the feasible set $\mathbb{P}$ three kinds of union boxes can be distinguished, 
\begin{enumerate}
	\item the feasible subpavings that  belong to the set $\mathbb{P}$ and satisfy this implication test $[{f}]([{\mathbf{p}}]) \subset {Y} \quad \Rightarrow \quad[{\mathbf{p}}] \subset \mathbb{P}$ (denoted $[{\mathbf{p_i}}]$ in Fig. \ref{Siv}).  
	\item the infeasible subpavings that make the empty intersection with $\mathbb{P}$, i.e. $[{f}]([{\mathbf{p}}]) \cap {Y}=\emptyset \Rightarrow[{\mathbf{p}}] \cap \mathbb{P}=\emptyset$ (denoted $[{\mathbf{p_o}}]$ in Fig. \ref{Siv}).
	\item the undetermined subpavings (or penumbra) for which nothing can be decided, and will be bisected, except if their width is too small (denoted $[{\mathbf{p_u}}]$ in Fig. \ref{Siv}).
\end{enumerate}


%\subsection{Contractor and Separators}
\vspace{-0.15cm}
\section*{ACKNOWLEDGMENT}

This  research  work   is funded   by the French LABEX Digicosme project ``PREGARI" (CRa 19006).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\addtolength{\textheight}{-20cm}  

\begin{thebibliography}{99}

\bibitem{c2x2} F. Al-Bender and J. Swevers. Characterization of friction force dynamics. In: IEEE Control Systems Magazine, 28(6), pp. 64-81, 2008.
\bibitem{c20} H. Olsson,  K. J. Astrom, C. C. De Wit,  M. Gafvert, and P. Lischinsky. Friction models and friction compensation. In: European Journal of Control, 4(3), pp. 176-195, 1998.
\bibitem{c23} K. Yoshida. Swing-up control of an inverted pendulum by energy-based methods. In Proceedings of the American Control Conference (ACC) (Cat. No. 99CH36251), Vol. 6, pp. 4045-4047, IEEE, 1999.
\bibitem{c10}  G. Abba, and P. Sardain. Modeling of frictions in the transmission elements of a robot axis for its identification. In IFAC Proceedings Volumes, 38(1), pp. 7-12, 2005.


\bibitem{c11}  I. C. Bogdan and G. Abba. Identification of the servomechanism used for micro-displacement. In IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 1986-1991, IEEE, 2009.
\bibitem{c100}R. H. Hensen,   M. J. van de Molengraft and M. Steinbuch. Frequency domain identification of dynamic friction model parameters. In IEEE Transactions on Control Systems Technology, pp. 191-196, 2002.
\bibitem{c7} L. Jaulin  and E. Walter. Set inversion via interval analysis for nonlinear bounded-error estimation. In Automatica, 29(4), pp.1053-1064, 1993.

\newpage
\bibitem{c777} D. Meizel,  O. Lévêque,  L. Jaulin,  and E. Walter. Initial localization by set inversion. In IEEE transactions on robotics and Automation, 18(6), pp. 966-971, 2002.
\bibitem{c15} I. Braems,  F. Berthier,  L. Jaulin,  M. Kieffer and E. Walter. Guaranteed estimation of electrochemical parameters by set inversion using interval analysis. In Journal of Electroanalytical Chemistry, 495(1), pp. 1-9, 2000.


\bibitem{c155} T. Kernicky,  M. Whelan  and E. Al-Shaer. Dynamic identification of axial force and boundary restraints in tie rods and cables with uncertainty quantification using Set Inversion Via Interval Analysis. In: Journal of Sound and Vibration, 423, pp. 401-420, 2018.	
\bibitem{c14} M. S. Madi,  K. Khayati and P. Bigras. Parameter estimation for the LuGre friction model using interval analysis and set inversion. In 2004 IEEE International Conference on Systems, Man and Cybernetics, IEEE Cat. No. 04CH37583, Vol. 1, pp. 428-433, 2004.


\bibitem{l1}  R. Paulen, M. E. Villanueva and B. Chachuat. Guaranteed parameter estimation of non-linear dynamic systems using high-order bounding techniques with domain and CPU-time reduction strategies. IMA - Mathematical Control and Information, 33(3), pp. 563-587, 2016.


\bibitem{l2} B. Chachuat, B. Houska,  R. Paulen,   N. Peric,  J. Rajyaguru and M. E.  Villanueva. Set-theoretic approaches in analysis, estimation and control of nonlinear systems. IFAC-PapersOnLine,  pp. 981-995, 2015.

\bibitem{c2} J. Alexandre Dit Sandretto,  Reliable NonLinear Model-Predictive Control via Validated Simulation. In  American Control Conference (ACC), IEEE, pp. 609-614, 2018. 

\bibitem{c3} M. Fnadi, W. Du, F. Plumet, F. Benamar. Model Predictive Control based Dynamic Path Tracking of a Four-Wheel Steering Mobile Robot. In 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 4518-4523, IEEE, 2019.


\bibitem{c4} M. Gafvert. Modelling the furuta pendulum. Department of Automatic Control, Lund Institute of Technology (LTH), 2006.


\bibitem{c12} J. Alexandre Dit Sandretto and A. Chapoutot. Validated explicit and implicit Runge-Kutta methods, 2016. \textit{https://perso.ensta-paris.fr/~chapoutot/dynibex/}.

\bibitem{c22} J. Alexandre Dit Sandretto, G. Trombettoni, D. Daney and G. Chabert. Certified calibration of a cable-driven robot using interval contractor programming. In Computational Kinematics,  Springer, pp. 209-217, 2014.


\bibitem{c18} R. E. Moore. Interval analysis. Englewood Cliffs: Prentice-Hall, Vol. 4, 1966.

\bibitem{c19} L. Jaulin, M. Kieffer,  O. Didrit,  and  E. Walter. Interval analysis. In: Applied interval analysis, pp. 11-43. Springer, London, 2001.



\end{thebibliography}




\end{document}