Commit 7b87a766773cc79429315f1add53923c818612e5

1 parent
eed078d0

### LaTex Sources Papier ACC 2021

No preview for this file type

... | ... | @@ -0,0 +1,709 @@ |

1 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

2 | +%2345678901234567890123456789012345678901234567890123456789012345678901234567890 | |

3 | +% 1 2 3 4 5 6 7 8 | |

4 | + | |

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

6 | + | |

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

8 | + | |

9 | +\IEEEoverridecommandlockouts % This command is only needed if | |

10 | + % you want to use the \thanks command | |

11 | + | |

12 | +\overrideIEEEmargins % Needed to meet printer requirements. | |

13 | + | |

14 | +%In case you encounter the following error: | |

15 | +%Error 1010 The PDF file may be corrupt (unable to open PDF file) OR | |

16 | +%Error 1000 An error occurred while parsing a contents stream. Unable to analyze the PDF file. | |

17 | +%This is a known problem with pdfLaTeX conversion filter. The file cannot be opened with acrobat reader | |

18 | +%Please use one of the alternatives below to circumvent this error by uncommenting one or the other | |

19 | +%\pdfobjcompresslevel=0 | |

20 | +%\pdfminorversion=4 | |

21 | + | |

22 | +% See the \addtolength command later in the file to balance the column lengths | |

23 | +% on the last page of the document | |

24 | + | |

25 | +% The following packages can be found on http:\\www.ctan.org | |

26 | +% The following packages can be found on http:\\www.ctan.org | |

27 | +\usepackage{graphics} % for pdf, bitmapped graphics files | |

28 | +\usepackage{epsfig} % for postscript graphics files | |

29 | +\usepackage{mathptmx} % assumes new font selection scheme installed | |

30 | +\usepackage{times} % assumes new font selection scheme installed | |

31 | +\usepackage{amsmath} % assumes amsmath package installed | |

32 | +\usepackage{amssymb} % assumes amsmath package installed | |

33 | +\usepackage{amsmath,bm} | |

34 | +\usepackage{blindtext} | |

35 | +\usepackage{lipsum} | |

36 | +\usepackage{mathtools} | |

37 | +\usepackage{relsize} | |

38 | +\usepackage{cuted} | |

39 | + | |

40 | +\usepackage{stmaryrd} | |

41 | +\usepackage[mathscr]{eucal} | |

42 | +%\usepackage{arevmath} % For math symbols | |

43 | + | |

44 | +\usepackage{graphics} % for pdf, bitmapped graphics files | |

45 | +\newcommand{\dd}{ | |

46 | + \mathop{}\mathopen{}\mathrm{d} | |

47 | +} | |

48 | +\usepackage{epsfig} % for postscript graphics files | |

49 | +\usepackage{mathptmx} % assumes new font selection scheme installed | |

50 | +\usepackage{times} % assumes new font selection scheme installed | |

51 | +\usepackage{amsmath} % assumes amsmath package installed | |

52 | +\usepackage{amssymb} % assumes amsmath package installed | |

53 | +\usepackage{amsfonts,bm} | |

54 | +\usepackage{stackengine} | |

55 | +\usepackage{stfloats} | |

56 | + | |

57 | +\usepackage{algorithm} | |

58 | +\usepackage{algorithmic} | |

59 | + | |

60 | + | |

61 | + | |

62 | +%\usepackage{graphics} % for pdf, bitmapped graphics files | |

63 | +%\usepackage{epsfig} % for postscript graphics files | |

64 | +%\usepackage{mathptmx} % assumes new font selection scheme installed | |

65 | +%\usepackage{times} % assumes new font selection scheme installed | |

66 | +%\usepackage{amsmath} % assumes amsmath package installed | |

67 | +%\usepackage{amssymb} % assumes amsmath package installed | |

68 | + | |

69 | +\title{\LARGE \bf | |

70 | +%Dynamic Identification of Nonlinear Inverted Pendulum with Quantified Uncertainties via Interval Analysis | |

71 | +Guaranteed Identification of Viscous Friction for a Nonlinear Inverted Pendulum Through Interval Analysis and Set Inversion | |

72 | +} | |

73 | + | |

74 | + | |

75 | +\author{Mohamed Fnadi$^{1}$, Julien Alexandre dit Sandretto$^{1}$, Gabriel Ballet$^{1}$, and Laurent Fribourg$^{2}$% <-this % stops a space | |

76 | +%\thanks{*This work was not supported by any organization}% <-this % stops a space | |

77 | +\thanks{$^{1}$ENSTA Paris, Institut Polytechnique de Paris, U2IS laboratory, 828 bd des mar\'echaux, 91762 | |

78 | + Palaiseau Cedex, France | |

79 | + {\tt\small \{firstname.lastname\}@ensta-paris.fr}}% | |

80 | +\thanks{$^{2}$LSV, Ecole Normale Sup\'erieure (ENS) \& CNRS, Paris-Saclay, France | |

81 | + {\tt\small fribourg@lsv.ens-cachan.fr}}% | |

82 | +} | |

83 | + | |

84 | + | |

85 | +\begin{document} | |

86 | + | |

87 | + | |

88 | + | |

89 | +\maketitle | |

90 | +\thispagestyle{empty} | |

91 | +\pagestyle{empty} | |

92 | + | |

93 | + | |

94 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

95 | +\begin{abstract} | |

96 | + | |

97 | +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. | |

98 | + | |

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

100 | + | |

101 | +\end{abstract} | |

102 | + | |

103 | + | |

104 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

105 | +\section{INTRODUCTION} | |

106 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

107 | + | |

108 | +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}. | |

109 | + | |

110 | + | |

111 | + | |

112 | + | |

113 | + | |

114 | + | |

115 | + | |

116 | +%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}. | |

117 | + | |

118 | + | |

119 | +%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}. | |

120 | + | |

121 | +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. | |

122 | + | |

123 | + | |

124 | +% 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. | |

125 | + | |

126 | + | |

127 | + 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). | |

128 | + | |

129 | +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. | |

130 | + | |

131 | +%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}. | |

132 | + | |

133 | +% while considering all the uncertainties included in the system's design and the embedded sensors (represented by intervals) | |

134 | + | |

135 | +% 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}. | |

136 | + | |

137 | + | |

138 | +% 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}. | |

139 | + | |

140 | +% which is based on the IA and set inversion tools | |

141 | +% | |

142 | +%(solutions of the set-inversion problem) from an initial research domain of parameters. | |

143 | + | |

144 | +% 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. | |

145 | +% | |

146 | +%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. | |

147 | +% | |

148 | +%by the comparisons between the nonlinear ODE solutions and their real measurements. | |

149 | +% | |

150 | +%% 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. | |

151 | +% | |

152 | +% | |

153 | +%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}. | |

154 | +% | |

155 | +%The results of this paper will serve as a basis for the achievement of these goals. | |

156 | +% | |

157 | +%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. | |

158 | + | |

159 | + | |

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

161 | + | |

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

163 | + | |

164 | +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. | |

165 | + | |

166 | + | |

167 | +%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. | |

168 | +% using a nonlinear inverted pendulum | |

169 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

170 | +\section{RECALL ON THE DYNAMIC MODELING OF AN INVERTED PENDULUM} | |

171 | +\label{Sec1} | |

172 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

173 | +\subsection{Robotic Device} | |

174 | + | |

175 | + | |

176 | +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. | |

177 | + | |

178 | +\begin{figure}[!h] | |

179 | + \vspace{-0.42cm} | |

180 | + \begin{center} | |

181 | + \includegraphics[scale=0.8]{Figures/SchemCine.pdf} | |

182 | + \vspace{-0.45cm} | |

183 | + \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).} | |

184 | + \label{Pendule} | |

185 | + \end{center} | |

186 | +\end{figure} | |

187 | + | |

188 | +\vspace{-0.71cm} | |

189 | +\subsection{Dynamic Modeling of the Inverted Pendulum} | |

190 | + As described in \cite{c4}, the equations of motion are given by, | |

191 | + \begin{equation}\label{DynEqu} | |

192 | + \Gamma={M} \ddot{{q}}+{B}({q}, \dot{{q}}) +{Q}(q), | |

193 | + \end{equation} | |

194 | + 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, | |

195 | + \begin{equation*} | |

196 | + \begin{aligned} | |

197 | + M(q) =\left[\begin{array}{ccc} \mu_1\sin(q_2)^2+\mu_2&& \mu_3 \cos(q_2) \\\\ | |

198 | + | |

199 | + \mu_3 \cos(q_2) && \mu_4\end{array}\right], | |

200 | + \end{aligned} | |

201 | + \end{equation*} | |

202 | + \begin{equation*} | |

203 | + \begin{aligned} | |

204 | +{B}({q}, \dot{{q}})= | |

205 | + \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 | |

206 | + \\\\ | |

207 | + - \mu_1\cos(q_2)\sin(q_2)\dot{q}_1^2 \end{array}\right], | |

208 | + \end{aligned} | |

209 | + \end{equation*} | |

210 | + \begin{equation*} | |

211 | + \mu_1 = l_p^2\left[ \dfrac{m_p}{4}+M\right], \hspace{0.4cm} | |

212 | + \mu_2=l_a^2[m_p+M]+\dfrac{J_m}{N_g^2}+J_a, | |

213 | + \end{equation*} | |

214 | + \begin{equation*} | |

215 | + \mu_3 = l_pl_a\left[ \dfrac{m_p}{2}+ M,\right], \hspace{0.4cm} | |

216 | + \mu_4 = l_p^2 \left[ \dfrac{m_p}{4}+M\right]+J_p,\hspace{0.15cm} | |

217 | + \end{equation*} | |

218 | + \begin{equation*} | |

219 | + \begin{aligned} | |

220 | + {Q}({q}, \dot{{q}})= | |

221 | + \left[\begin{array}{c} 0 | |

222 | + \\\\ | |

223 | + \mu_g\sin(q_2) \end{array}\right], | |

224 | + \end{aligned} | |

225 | + ~~~\mu_g=\left[\dfrac{m_p}{2}+M\right]l_pg, | |

226 | + \end{equation*} | |

227 | + 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. | |

228 | + | |

229 | + 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}$, | |

230 | + \begin{equation} | |

231 | + \begin{aligned} | |

232 | + \Gamma =\left[\begin{array}{c} \tau -f_{v_1} \dot{q}_1 \\\\ | |

233 | + -f_{v_2} \dot{q}_2\end{array}\right], | |

234 | + \end{aligned} | |

235 | + \end{equation} | |

236 | + 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, | |

237 | + \begin{equation}\label{DynEqu2} | |

238 | + y_m=f(q,\dot{q}, \ddot{q}, p), | |

239 | + \end{equation} | |

240 | + 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 | |

241 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

242 | +\section{IDENTIFICATION APPROACHES} | |

243 | +\label{Sec2} | |

244 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

245 | +The purpose of this section is to propose an identification strategy for the unknown vector ${p}$ embedding all the model parameters. | |

246 | +%$${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}$$ | |

247 | +%where $n_p$ is the number of the identified parameters. | |

248 | +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. | |

249 | + | |

250 | +%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}). | |

251 | + | |

252 | +\subsection{Torque, speed and acceleration computation } | |

253 | + | |

254 | +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). | |

255 | + | |

256 | + | |

257 | + | |

258 | + | |

259 | + | |

260 | + | |

261 | +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, | |

262 | +\begin{eqnarray}\label{torque} | |

263 | +\tau=\dfrac{N_mN_g\eta_m\eta_g(\omega_0-\omega)}{R_v} | |

264 | +\end{eqnarray} | |

265 | +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. | |

266 | + | |

267 | +\subsection{Identification with classical least square method} | |

268 | + | |

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

270 | + | |

271 | +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: | |

272 | + | |

273 | +\begin{equation} | |

274 | +\begin{aligned} | |

275 | +y_m=D_{}(q, \dot{q}, \ddot{q}) p | |

276 | +\end{aligned} | |

277 | +\end{equation} | |

278 | +where $D \in\mathbb{R}^{2\times n_p} $ is the matrix of the measured values, expressed as | |

279 | +a function of joint positions $q$, velocities $\dot{q}$ and joint acceleration vector $\ddot{q}$. It is defined as, | |

280 | +\small | |

281 | +$$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 \\\\ | |

282 | +-\cos(q_2)\sin(q_2)\dot{q}_1^2 & 0 & \ldots \end{array}\right.$$ | |

283 | +$$\left.\begin{array}{cccccc} \cos(q_2)\ddot{q}_2 - \sin(q_2)\dot{q}_2^2 & 0& 0 & \dot{{q}}_1 & 0 \\\\ | |

284 | +\cos(q_2)\ddot{q}_1 & \ddot{q}_2 & \sin(q_2) & 0 &\dot{{q}}_2 \end{array}\right],$$ | |

285 | + | |

286 | +\normalsize | |

287 | +To identify the inertial and geometric parameters, we apply the LSMI method to the augmented form relying wholly on measurement data, | |

288 | +\begin{equation} | |

289 | +\begin{aligned} | |

290 | +Y=W p+\rho | |

291 | +\end{aligned} | |

292 | +\end{equation} | |

293 | +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$. | |

294 | + | |

295 | +The optimal solution $\widehat{p}$ can be expressed as follows, | |

296 | +\begin{equation} | |

297 | +\begin{aligned} | |

298 | +\widehat{p}=\underset{p}{\operatorname{\arg\min}}\|Y-Wp\|_{_2}=\left(W^{T} W\right)^{-1} W^{T} Y=W^{+} Y | |

299 | +\end{aligned} | |

300 | +\end{equation} | |

301 | +where $\|{a}\|_{_2}$ denotes the second euclidean norm of the vector $a$, and $W^{+}$ denote the Moore-Penrose inverse matrix of $W$. | |

302 | +In practice, it is recommended to use at least $500\times n_p$ measures to identify correctly the needed parameters by LSMI. | |

303 | + | |

304 | + | |

305 | +\subsection{Guaranteed identification with uncertainty quantification} | |

306 | + | |

307 | +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). | |

308 | + | |

309 | +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. | |

310 | + | |

311 | +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, | |

312 | +%\in \mathbb{R}^{n_p} | |

313 | +\begin{eqnarray}\label{P_set}\small | |

314 | +\begin{split} | |

315 | +\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.\\ | |

316 | +&~~~~~~~~~~~~~~~~~~\left.\text { s.t. }~~ f(q(i),\dot{q}(i),\ddot{q}(i),p)~\in~ [\mathbf{y}](i) ~ \right. \Big \} | |

317 | +\end{split} | |

318 | +\end{eqnarray} | |

319 | +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$. | |

320 | + | |

321 | +Then, the feasible set of accepted parameters $\mathbb{P}$ can be derived as, | |

322 | +\begin{eqnarray}\label{P_sett} | |

323 | +\small | |

324 | +\begin{split} | |

325 | +\mathbb{P} | |

326 | +%&=\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)\\ | |

327 | +%&\left.\text { s.t. }~~ f(q(1),\dot{q}(1),\ddot{q}(1),p)-y(1)=0,~ y(1)\in [y](1)\right. ~~\&\\ | |

328 | +%%& \&\\ | |

329 | +%&~~~~~~~~~~~~~~ \vdots ~~~~~~~~~\vdots~~~~~~~~~\vdots ~~~~~~~~~\vdots~~~~~~~~~\vdots~ \\ | |

330 | +%&~~~\&~\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)\\ | |

331 | +%&\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 \},\\ | |

332 | +%\vspace{0.1cm} | |

333 | +& = \mathbb{P}_{1}\cap \mathbb{P}_{2}\ldots \cap \mathbb{P}_{n_y} | |

334 | + = \bigcap_{i=1}^{n_y} \mathbb{P}_i | |

335 | +\end{split} | |

336 | +\end{eqnarray} | |

337 | + | |

338 | +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, | |

339 | + \begin{equation}\small | |

340 | +\begin{split} | |

341 | +&~~~~~~~\underbrace{\left[\begin{array}{c} y_m(1)\\\\ \vdots \\\\ | |

342 | + y_m(n_y)\end{array}\right]}_{Y_m} | |

343 | + \in | |

344 | + \underbrace{\left[\begin{array}{c}[\mathbf{y}](1) \\\\ \vdots \\\\ \left[\mathbf{y}\right](n_y)\end{array}\right]}_{[\mathbf{Y}]},\\\\ | |

345 | +& ~\text{with} ~~y_m(i)=f(q(i),\dot{q}(i),\ddot{q}(i),p),~~\forall i\in\llbracket 1,n_y\rrbracket | |

346 | +\end{split} | |

347 | +\label{inclusion} | |

348 | +\end{equation} | |

349 | + | |

350 | +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}}$. | |

351 | +\begin{eqnarray}\label{P_settt}\small | |

352 | +\widehat{\mathbb{P}}=f^{-1}([\mathbf{Y}])\cap[\mathbf{p}] | |

353 | +\end{eqnarray} | |

354 | +% \begin{equation}\small | |

355 | +%\begin{split} | |

356 | +%\left[\begin{array}{c}f(q(1),\dot{q}(1),\ddot{q}(1),p) \\\\ \vdots \\\\ | |

357 | +%f(q(n_y),\dot{q}(n_y),\ddot{q}(n_y),p)-y(n_y)\end{array}\right] \in | |

358 | +%\left[\begin{array}{c}[y](1) \\\\ \vdots \\\\ \left[y\right](n_y)\end{array}\right], | |

359 | +%\end{split} | |

360 | +%\end{equation} | |

361 | + | |

362 | + | |

363 | + | |

364 | +%\begin{eqnarray}\label{P_settt} | |

365 | +%\begin{split} | |

366 | +%\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}}]\\ | |

367 | +%&~~~~~~\left.\text { s.t. }~~ \textbf{y}_m-\textbf{y}=0,~~ \text{with}~ \textbf{y}\in [\textbf{y}]\right. \big \}, | |

368 | +%\end{split} | |

369 | +%\end{eqnarray} | |

370 | + | |

371 | +%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. | |

372 | + | |

373 | +\begin{figure*}[b] | |

374 | + \vspace{-0.6cm} | |

375 | + \begin{minipage}[c]{.46\linewidth} | |

376 | + \begin{center} | |

377 | + \centering | |

378 | + \includegraphics[scale=0.21]{Figures/q.png}\hspace{-0.5cm} | |

379 | + \includegraphics[scale=0.21]{Figures/dq.png}\\\vspace{-0.1cm} | |

380 | + \includegraphics[scale=0.21]{Figures/ddq.png}\hspace{-0.5cm} | |

381 | + \includegraphics[scale=0.21]{Figures/tau.png}\\\vspace{-0.1cm} | |

382 | + \caption{ Uncertainty rectangles associated to the measurements of acquisition \#1: (a) angular positions; (b) Angular velocities; (c) Angular accelerations; (d) DC motor torque.} | |

383 | + \label{Experiments1} | |

384 | + \end{center} | |

385 | + \end{minipage} \hfill | |

386 | + \vspace{-0.6cm} | |

387 | + \begin{minipage}[c]{.46\linewidth} | |

388 | + \begin{center} | |

389 | + \hspace{-0.3cm} \centering | |

390 | + \includegraphics[scale=0.21]{Figures/q2.png}\hspace{-0.5cm} | |

391 | + \includegraphics[scale=0.21]{Figures/dq2.png}\\\vspace{-0.1cm} | |

392 | + \includegraphics[scale=0.21]{Figures/ddq2.png}\hspace{-0.5cm} | |

393 | + \includegraphics[scale=0.21]{Figures/tau2.png}\\\vspace{-0.1cm} | |

394 | + \caption{ Uncertainty rectangles associated to the measurements of acquisition \#2: (a) angular positions; (b) Angular velocities; (c) Angular accelerations; (d) DC motor torque.} | |

395 | + \label{Experiments2} | |

396 | + \end{center} | |

397 | + \end{minipage} | |

398 | +\end{figure*} | |

399 | + | |

400 | + | |

401 | +%\begin{algorithm} | |

402 | +% \caption{Compute the paving of the friction Parameters} | |

403 | +% \begin{algorithmic}[1] | |

404 | +% \REQUIRE $\mathbb{P}_0=[-3,3]\times[-3,3],~ \mathbb{P}=[]\times[]$, $~[q],~[\dot{q}],~[\ddot{q}],~[y_m]$%, \mathcal{P}:=\emptyset .$ | |

405 | +%% \ENSURE $y = x^n$ %\COMMENT{put some comments here} | |

406 | +% \FOR{$i = 1$ to $n_y$} | |

407 | +%% \STATE {\COMMENT{//$m\leq n^2$}} | |

408 | +% \STATE $\mathbb{P}_{[if]} \gets f_{Fwd}([q](i),[\dot{q}](i),[\ddot{q}](i),p)$ | |

409 | +% \STATE $\mathbb{P}_{[ib]} \gets f_{Bwd}([q](i),[\dot{q}](i),[\ddot{q}](i),p)$ | |

410 | +% \STATE $\mathbb{P}_{i} \gets \mathbb{P}_{[if]}\cap \mathbb{P}_{[ib]}$ | |

411 | +% \STATE Push $\mathbb{P}_i$ onto $\mathbb{P}$ | |

412 | +% \ENDFOR | |

413 | +% \STATE $\mathbb{P} \gets \mathbb{P}_{1}\cap \mathbb{P}_{2}\ldots \cap \mathbb{P}_{n_y}$ | |

414 | +% \end{algorithmic} | |

415 | +%\label{algo1} | |

416 | +%\end{algorithm} | |

417 | + | |

418 | + | |

419 | +%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. | |

420 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

421 | +\section{EXPERIMENTS} | |

422 | +\label{Sec3} | |

423 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

424 | + | |

425 | +\subsection{Identification results} | |

426 | +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). | |

427 | + | |

428 | +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}. | |

429 | + | |

430 | + | |

431 | + | |

432 | +% Many trails on the real system were performed to quantify these error boxes. | |

433 | +%In order to excite the dynamics of the inverted pendulum, | |

434 | + | |

435 | +%($[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]$). | |

436 | + | |

437 | +%Since the first acquisition | |

438 | + | |

439 | + | |

440 | +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$. | |

441 | + | |

442 | +\begin{center} | |

443 | + \begin{table}[h!] | |

444 | + \vspace{-0.05cm} | |

445 | + \caption{Identification results by LSMI method.} | |

446 | + \centering | |

447 | + \begin{tabular}{lccc} | |

448 | + Symbol & Values with \#1 & Values with \#2 & CAD Values \\ | |

449 | + \hline | |

450 | + $\mu_1[kg.m^2]$ & $1.04082e^{-3}$ & $1.06171e^{-3}$ & $1.18152e^{-3}$\\ | |

451 | + $\mu_2[kg.m^2]$ &$2.56211e^{-3}$ & $2.37309e^{-3}$ & $2.55064e^{-3}$\\ | |

452 | + $\mu_3[kg.m^2]$ & $8.21e^{-4}$ & $8.29e^{-4}$ & $7.59e^{-4}$\\ | |

453 | + $\mu_4[kg.m^2]$ & $1.31781e^{-3}$ & $1.19171e^{-3}$ & $1.18152e^{-3}$ \\ | |

454 | + $\mu_g[kg.m^2.s^{-2}]$ & $7.4175e^{-2}$ & $7.2591e^{-2}$ & $7.3575e^{-2}$\\ | |

455 | + $f_{v_1}[N.m.s]$ & $8.13e^{-2}$ & $7.02e^{-2}$ & $-$\\ | |

456 | + $f_{v_2}[N.m.s]$ & $6.42e^{-4}$ & $5.93e^{-4}$ & $-$\\ \hline | |

457 | + \end{tabular} | |

458 | + \label{Table_LSIM} | |

459 | + \end{table} | |

460 | +\end{center} | |

461 | + | |

462 | +\vspace{-0.45cm} | |

463 | +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}). | |

464 | + | |

465 | + | |

466 | + | |

467 | +\begin{center} | |

468 | + \begin{table}[h!] | |

469 | + \vspace{-0.3cm} | |

470 | + \caption{Inertial and geometric parameters with uncertainties.} | |

471 | + \centering | |

472 | + \begin{tabular}{l c c} | |

473 | + Symbol & Mean value & Uncertainty \\ | |

474 | + \hline | |

475 | + $\mu_1[kg.m^2]$ & $1.09468e^{-3}$ & $\pm 12 \%$ \\ | |

476 | + $\mu_2[kg.m^2]$ & $2.49528e^{-3}$ & $\pm 12\%$\\ | |

477 | + $\mu_3[kg.m^2]$ & $8.03e^{-4}$ & $\pm 10\%$ \\ | |

478 | + $\mu_4[kg.m^2]$ & $1.23035e^{-3}$ & $\pm 12\%$ \\ | |

479 | + $\mu_g [kg.m^2.s^{-2}]$ & $7.3447e^{-2}$ & $\pm 15\%$\\ \hline | |

480 | + \end{tabular} | |

481 | + \label{Tab1_uncetainties} | |

482 | + \end{table} | |

483 | +\end{center} | |

484 | + | |

485 | + \vspace{-0.3cm} | |

486 | +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]$. | |

487 | + | |

488 | +% These subpavings, by the way, have been easily drawn using the Visualizer for Intervals and Boxes (VIBEs) \cite{c16} | |

489 | +% (red boxes are feasible sets, yellow boxes are undetermined, and the blues ones are the infeasible sets) | |

490 | +\vspace{-0.4cm} | |

491 | +\begin{figure}[!h] | |

492 | + \centering | |

493 | + \includegraphics[scale=0.56]{Figures/Siva1.pdf}\hspace{0.cm} | |

494 | + \includegraphics[scale=0.56]{Figures/Siva.pdf}\hspace{-0.3cm} | |

495 | + \caption{ Pavings generated by SIVIA algorithm to show the feasible set $\widehat{\mathbb{P}}$ for the friction parameters. [Left] : scenario \#1. [Right] : scenario \#2.} | |

496 | + \label{Vibes} | |

497 | +\end{figure} | |

498 | + | |

499 | +\vspace{-0.55cm} | |

500 | +\subsection{Identification results validation} | |

501 | + | |

502 | +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. | |

503 | + | |

504 | + \begin{equation}\small | |

505 | + \mathrm{RMSE}=\sqrt{\dfrac{1}{n} \sum_{i=1}^{n}\left(\tau_{i}-\widehat{\tau}_{i}\right)^{2}} | |

506 | + \label{RMSE} | |

507 | + \end{equation} | |

508 | + | |

509 | + | |

510 | + | |

511 | + | |

512 | + | |

513 | + | |

514 | + | |

515 | + 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, | |

516 | +\begin{equation}\label{ODE} | |

517 | +\left\{\begin{array}{l} | |

518 | +\dot{{x}}(t)=g(t, {x}(t),u(t)) \\ | |

519 | +{x}(0) \in\left[\mathbf{{x}}_{0}\right] \subseteq {\mathbb{IR}}^{4} | |

520 | +\\ | |

521 | +{u}(0) \in\left[\mathbf{{u}}_{0}\right] \subseteq {\mathbb{IR}} | |

522 | +\end{array}\right. | |

523 | +\end{equation} | |

524 | +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. | |

525 | + | |

526 | + | |

527 | + | |

528 | +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}. | |

529 | +\vspace{-0.4cm} | |

530 | +\begin{figure}[!h] | |

531 | + \begin{center} | |

532 | + \includegraphics[scale=0.37]{Figures/validation/Torques.png}\hspace{-0.3cm} | |

533 | + \vspace{-0.3cm} | |

534 | + \caption{Comparison between measured and estimated torques via IA and LSMI methods.} | |

535 | + \label{Torques} | |

536 | + \end{center} | |

537 | +\end{figure} | |

538 | + | |

539 | +\vspace{-0.3cm} | |

540 | +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. | |

541 | +\begin{figure}[!h] | |

542 | + \vspace{-0.4cm} | |

543 | + \begin{center} | |

544 | + \centering | |

545 | + \includegraphics[scale=0.2245]{Figures/validation/q2_comp100.png}\hspace{-0.6cm} | |

546 | + \includegraphics[scale=0.2245]{Figures/validation/q2_comp.png}\\\vspace{-0.1cm} | |

547 | + \includegraphics[scale=0.2245]{Figures/validation/q2_compTorque.png}\hspace{-0.6cm} | |

548 | + \includegraphics[scale=0.2245]{Figures/validation/q2_comp4.png}\\\vspace{-0.1cm} | |

549 | + \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$.} | |

550 | + \label{Experiments_valid1} | |

551 | + \end{center} | |

552 | +\end{figure} | |

553 | + | |

554 | + | |

555 | +% For each case, same conditions are maintained for the real system and the ODE solver to get compared results. | |

556 | +\vspace{-0.5cm} | |

557 | +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. | |

558 | + | |

559 | +\begin{center} | |

560 | + \begin{table}[h!] | |

561 | + \vspace{-0.3cm} | |

562 | + \caption{Degree of coverage between the model and physical reality.} | |

563 | + \centering | |

564 | + \begin{tabular}{l c c c c} | |

565 | + Scenario & (a) & (b) & (c) & (d) \\ | |

566 | + \hline | |

567 | + With IA parameters & $46\%$ & $61\%$ & $38\%$ & $25\%$ \\ | |

568 | + With LSMI parameters & $41\%$ & $34\%$ & $20\%$ & $21\%$ \\ \hline | |

569 | + \end{tabular} | |

570 | + \label{coverge_rate} | |

571 | + \end{table} | |

572 | +\end{center} | |

573 | + | |

574 | +\vspace{-0.25cm} | |

575 | +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}. | |

576 | + | |

577 | + | |

578 | + | |

579 | + | |

580 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

581 | +\section{CONCLUSION AND FUTURE WORK} | |

582 | +\label{Sec4} | |

583 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

584 | + | |

585 | +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 | |

586 | +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. | |

587 | + | |

588 | + | |

589 | + | |

590 | +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}. | |

591 | + | |

592 | + | |

593 | + | |

594 | + | |

595 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

596 | + | |

597 | + | |

598 | + | |

599 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

600 | + | |

601 | + | |

602 | + | |

603 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

604 | +\section*{APPENDIX} | |

605 | +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}. | |

606 | +\subsection{Interval Analysis : the main notations and definitions} | |

607 | +\begin{itemize} | |

608 | + \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. | |

609 | + \item An interval vector, or a box $\left[\mathbf{x}\right]$ of $\mathbb{R}^n$ is a Cartesian | |

610 | + 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$. | |

611 | + \item The width $w(\left[\mathbf{x}\right])$ of a box $\left[\mathbf{x}\right]$ is the length of its largest side. | |

612 | + \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}]\}$. | |

613 | + \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. | |

614 | +\end{itemize} | |

615 | + | |

616 | + | |

617 | + | |

618 | + | |

619 | +\subsection{Set Inversion and SIVIA algorithm} | |

620 | +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, | |

621 | +\begin{equation} | |

622 | + \mathbb{P}=\left\{{p} \in \mathbb{R}^{n} \mid {f}({p}) \in {Y}\right\}={f}^{-1}({Y}) | |

623 | +\end{equation} | |

624 | + | |

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

626 | + | |

627 | + | |

628 | +\begin{figure}[!h] | |

629 | + \vspace{-0.3cm} | |

630 | + \begin{center} | |

631 | + \includegraphics[scale=0.7]{Figures/SIVIA.pdf}\hspace{-0.3cm} | |

632 | + \vspace{-0.4cm} | |

633 | + \caption{Illustration of the principal of the set inversion problem - Inclusion tests: infeasible (blue), undetermined (green), feasible (red).} | |

634 | + \label{Siv} | |

635 | + \end{center} | |

636 | +\end{figure} | |

637 | + | |

638 | + \vspace{-0.6cm} | |

639 | +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, | |

640 | +\begin{enumerate} | |

641 | + \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}). | |

642 | + \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}). | |

643 | + \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}). | |

644 | +\end{enumerate} | |

645 | + | |

646 | + | |

647 | +%\subsection{Contractor and Separators} | |

648 | +\vspace{-0.15cm} | |

649 | +\section*{ACKNOWLEDGMENT} | |

650 | + | |

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

652 | + | |

653 | + | |

654 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |

655 | + | |

656 | + | |

657 | +\addtolength{\textheight}{-20cm} | |

658 | + | |

659 | +\begin{thebibliography}{99} | |

660 | + | |

661 | +\bibitem{c2x2} F. Al-Bender and J. Swevers. Characterization of friction force dynamics. In: IEEE Control Systems Magazine, 28(6), pp. 64-81, 2008. | |

662 | +\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. | |

663 | +\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. | |

664 | +\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. | |

665 | + | |

666 | + | |

667 | +\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. | |

668 | +\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. | |

669 | +\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. | |

670 | + | |

671 | +\newpage | |

672 | +\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. | |

673 | +\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. | |

674 | + | |

675 | + | |

676 | +\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. | |

677 | +\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. | |

678 | + | |

679 | + | |

680 | +\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. | |

681 | + | |

682 | + | |

683 | +\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. | |

684 | + | |

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

686 | + | |

687 | +\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. | |

688 | + | |

689 | + | |

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

691 | + | |

692 | + | |

693 | +\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/}. | |

694 | + | |

695 | +\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. | |

696 | + | |

697 | + | |

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

699 | + | |

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

701 | + | |

702 | + | |

703 | + | |

704 | +\end{thebibliography} | |

705 | + | |

706 | + | |

707 | + | |

708 | + | |

709 | +\end{document} | |

... | ... |

94.8 KB

No preview for this file type

No preview for this file type

No preview for this file type

No preview for this file type

No preview for this file type