# Model identification of a cable-based mechanical transmission for robotics using a Best Linear Approximation approach

Bassem Boukhebouz, Guillaume Mercère, Mathieu Grossard, Édouard Laroche

► To cite this version:

Bassem Boukhebouz, Guillaume Mercère, Mathieu Grossard, Édouard Laroche. Model identification of a cable-based mechanical transmission for robotics using a Best Linear Approximation approach. Control Engineering Practice, 2025, 164, pp.106457. 10.1016/j.conengprac.2025.106457. hal-05112170

HAL Id: hal-05112170

https://hal.science/hal-05112170v1

Submitted on 4 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Model identification of a cable-based mechanical transmission for robotics using a Best Linear Approximation approach ⋆

Bassem Boukhebouz a,b, Guillaume Mercère a, c, Mathieu Grossard c, Édouard Laroche b, d, ∗

ARTICLE INFO

Keywords:

Data-based modeling

Model learning

Frequency domain

Multisine excitation

Best Linear Approximation (BLA)

Flexible transmission

A B S T R A C T

This article presents an experimental approach for frequency-based identification and excitation signal design to be used for identifying the mechanical motor-to-joint cable-based transmission of a robotic actuation system. By using the Best Linear Approximation (BLA) theory, frequency responses are estimated and nonlinear distortions are mitigated. The identification process is performed in closed-loop to ensure system stability, with multisine excitation signals designed to respect system constraints. The paper presents a technique for enhancing the accuracy of estimated frequency responses when dealing with signal saturation in a closed-loop system. This approach involves shaping the multisine excitation to enhance the Crest Factor of the control signal. The study includes experimental results from non-parametric frequency identification conducted on a test bench, and a parametric model is derived using a non-convex optimization method.

1. Introduction

1.1. Background of the study

Manipulator robots have transformed various industries, including manufacturing, healthcare, automotive, and aerospace, by providing efficient solutions that simplify and enhance human tasks. These robots are capable of handling both simple and complex operations with key attributes like precision, speed, and repetitive accuracy. They incorporate various components such as sensors, actuators, mechanical parts, and control systems. However, despite their significant advantages, some tasks remain difficult to fully automate due to their complexity or limited cost-effectiveness (Abroug & Laroche, 2015). As a result, human intervention, expertise, and dexterity are still crucial in certain situations. This highlights the need for further advancements and exploration of new technologies to address these limitations and push the boundaries of automation in manipulator robotics.

In recent years, a major shift has occurred in the field, focusing on combining the strength and precision of machines with human intelligence and dexterity. This shift has led to the development of a new generation of robots known as cobots (Goodrich & Schultz, 2007), designed specifically to work alongside human operators (Colgate & Hogan, 1988). A key feature of cobots is their ability to physically interact with humans through dedicated interaction interfaces (Goodrich &

Schultz, 2007). This interaction, often called robotic co-manipulation, involves both the operator and robot jointly manipulating objects.

Unlike traditional industrial robots that prioritize speed and preci-
sion, cobots emphasize assisting operators according to specific criteria.
For instance, some cobots feature intrinsically reversible actuation
systems, which enable external contact forces to be transmitted from
the motor to the joints. An example studied in this work is the screw-
cable actuation system, which offers low inertia and friction, ensuring
high sensitivity to external forces. This type of system has been applied
successfully in various robotics fields, including teleoperation (Gosselin
& Grossard, 2023), multifingered grippers (Martin & Grossard, 2014),
and collaborative robot arms (Makarov, Grossard, Rodriguez-Ayerbe, &
Dumur, 2016). Consisting of an electric motor and a cable-based trans-
mission mechanism, these actuation systems provide high efficiency
and are widely used in collaborative robotics.

1.2. Problem formulation

The presence of cables within screw-cable actuation system intro-
duces mechanical flexibilities that can affect both static and dynamic
positioning precision (Fichera & Grossard, 2017). Additionally, due to
their specific design, nonlinear phenomena manifest in their dynamics,
such as piece-wise multizone operation principle (Fichera & Grossard,

<footnote>ª LIAS Laboratory, University of Poitiers, France</footnote>

<footnote>ᵇICube Laboratory, University of Strasbourg, France</footnote>

<footnote>ᶜInteractive Robotics Laboratory, CEA LIST, Palaiseau, France</footnote>

<footnote>⋆ This work has been supported by a 3-year PhD funding by Commissariat à l'énergie atomique et aux énergies alternatives.</footnote>

<footnote>* Corresponding author.

*E-mail addresses:* bassem.boukhebouz@expleogroup.com (B. Boukhebouz), laroche@unistra.fr (É. Laroche).</footnote>
2017) and Coulomb friction (Garrec, P.Friconneau, Measson, & Perrot, 2008). When dealing with such systems, especially in the presence of operators in proximity or on direct contact, controller synthesis must account for the human-robot coupling and cable-induced flexibilities. This requires the use of advanced control synthesis methods (Abroug & Laroche, 2015), commonly referred to as *model-based* approaches, which require the development of accurate dynamic models for these systems.

Analytical modeling of robot dynamics, as the first step of model-based controller design, when developed with a deep understanding of the robot design, can effectively describe the dynamics of a robotic system. However, the complexity of these systems, coupled with the presence of flexible modes within the bandwidth of interest, makes analytical modeling a highly specialized task that demands extensive expertise and system-specific knowledge. Another approach to modeling robotic systems involves the use of input-output measurements to construct mathematical models (Ljung, 1986; Pintelon & Schoukens, 2012; Söderström & Stoica, 1988). This technique, known as *system identification* and also called *model learning* or *data-based modeling*, focuses on deriving mathematical models based on the input and output data of a process. By leveraging system identification, it becomes possible to predict future system states, perform simulations, and synthesize controllers.

The identification of robotic systems depends mainly on the instrumentation of the system, which makes each problem of identification in robotics specific. The robotic system considered in this work is at the same time nonlinear, unstable and resonant. The identification of such a system is thus a challenging task, as well as to find appropriate excitations and efficient identification methods. In this work, it was chosen to work in the frequency domain by estimating frequency responses. This choice benefits from the advantages that non-parametric frequency responses offer in an identification process, namely, the ability to: (i) provide knowledge about the input-output behavior of the system without requiring strong prior knowledge about the system to be identified; (ii) capture resonances that occur in the frequency band of interest; (iii) serve as frequency samples for estimating parametric models.

For position-based control, the identification of such actuated mechanical transmissions is performed in closed loop to ensure their stability. In a closed-loop identification process, one of the challenges is to be able to excite the system while respecting its constraints, such as actuators saturations. These saturations limit the amplitude of the excitation signal. Therefore, on the one hand, high amplitudes of excitation must be avoided to prevent signal saturation. On the other hand, properly exciting the entire dynamic of the robot to estimate its physical parameters requires an excitation signal with a good signal-to-noise ratio (SNR) in the frequency band of interest. Often, the excitation of industrial robots is injected as a disturbance at the control level, such as motor torque. However, in the present case, due to the particularity of the experimental robotic system, the excitation is applied as a motor position reference. Therefore, it is important to design a reference excitation signal that takes into account the saturations on the control and allows for rich measurements to enhance identification.

### 1.3. Objectives

This paper proposes an experimental data-driven approach to identify such robotic mechanical transmission, using frequency-based techniques and reference excitation design. Certain systems exhibit significant nonlinearities that can affect the properties of the identified model, like passivity<fnref n="1" /> (Abroug & Laroche, 2015). Therefore, it is crucial

to analyze and mitigate the impact of nonlinearities to uncover the underlying linear behavior of the system.

To achieve the objective of minimizing the impact of nonlinear distortion, an advanced approach based on Best Linear Approximation (BLA) theory (Pintelon & Schoukens, 2012) was used to estimate frequency responses. Thus, to perform a frequency non-parametric analysis of the experimental system, based on the BLA theory, a closed-loop excitation design approach has been proposed takes into account the constraints on the signal amplitude.

The paper is organized as follows:

* Section 2 introduces the theory of best linear approximation for systems with nonlinearities and its use for frequency-response estimation.

* Section 3 focuses on the design of excitation signals for non-parametric closed-loop model identification. The impact of phase selection on the temporal shape of excitation is discussed and the problem of closed-loop excitation design considering actuation system saturations is introduced.

* In Section 4, the experimental results obtained from non-parametric frequency identification on a test bench available at CEA-LIST are presented. We start by explaining the operation of the experimental system, followed by the hardware and software structure of the control strategy used in the test bench. Finally, the experimental results obtained using the identification approach and proposed excitation design are evaluated. The chapter concludes with a validation of the proposed approach through parametric identification based on the $H_\infty$ norm.

* In Section 5, the non-parametric results are used to tune the parameters of a physical model using an identification approach based on the $H_\infty$ norm.

* Section 6 concludes the paper and provides future perspectives.

## 2. Estimation of the best linear approximation

### 2.1. Best linear approximation

The presence of nonlinearities in robotic systems arises from several factors, including: (i) coupling between rigid segments, (ii) motor or joint friction (Hamon, Gautier, & Garrec, 2011), and (iii) transmission stiffness (Fichera & Grossard, 2017). Identifying the dynamics of such nonlinear systems poses challenges in selecting an appropriate model structure and estimating its parameters, making the task more difficult compared to linear systems (Schoukens & Ljung, 2019).

To describe the dynamics of nonlinear systems, linear models, around operating points, are often sought (Peumans, 2020; Schoukens & Ljung, 2019). This approach, known as *linear approximation*, reduces the complexity associated with identifying nonlinear models. Linear approximation can be achieved in two ways:

* by selecting the amplitude (or the root mean square (RMS) value) of the excitation signal $u(t)$ in such a way that the nonlinearities of the system remain inactive. In this case, the resulting linear model is unique and independent of the excitation amplitude (see Fig. 1),

* by not restricting the choice of excitation amplitude, thus covering a wide range of operating points and extracting more information. This approach allows the nonlinearities to become active, introducing nonlinear distortions<fnref n="2" /> in the measurements (as illustrated in Fig. 1). Unlike the first case, the obtained linear models in this approach are not unique and depend on the excitation applied to the system.

<footnote n="1">¹ Passivity is the property of a system not producing energy and not providing energy greater than what it has received from an external energy source (Wyatt, Chua, Gannett, Goknar, & Green, 1981).</footnote>

<footnote n="2">² Nonlinear distortions refer to the deviation between the output of the nonlinear system and that of its approximated linear model in a noiseless context (Pintelon & Schoukens, 2012).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Impact of excitation on the output in case of a saturation system.</figcaption></figure>

This work falls into the second situation, where the nonlinearities of the system lead to the appearance of nonlinear distortions in the measurements whatever the choice of the amplitude of the excitation (Boukhebouz, Mercère, Grossard, Lamy, & Laroche, 2020). Thus, to find a linear approximation of the nonlinear system, a dedicated approach is sought to extract the underlying linear model of the nonlinear system. Among the possible approaches, the Best Linear Approximation (Enqvist, 2005; Ljung, 1986; Mäkilä & Partington, 2004; Pintelon & Schoukens, 2012) is often used to linearize the behavior of the nonlinear system. This approximation, either in the time domain or in the frequency domain, aims at minimizing, in the least squares sense, the error between the measured output of a nonlinear system and the output of a linear model (Bendat & Piersol, 1980) as follows

$$g_{BLA} = \arg\min_{u} \mathbb{E}[|y(t) - g(t) * u(t)|^2], \quad (1)$$

where $u(t)$ is the input of the nonlinear system, $y(t)$ is its output and $g(t)$ is its impulse response and * stands for the convolution product.

The literature (Enqvist, 2005; Pintelon & Schoukens, 2012) emphasizes that the best linear approximation relies on the characteristics of the excitation signal $u(t)$. This dependency is captured through the standard mathematical expectation operator $\mathbb{E}$ in Eq. (1). In simpler terms, the choice and properties of the excitation signal directly influence the obtained linear approximation of the system. Different excitation signals may lead to different linear approximations, as their properties (including distribution) impact the behavior of the nonlinear system.

The equivalent of Eq. (1) in the frequency domain is

$$G_{BLA}(j\omega) = \arg\min_{U} \mathbb{E}[|Y(j\omega) - G(j\omega)U(j\omega)|^2], \quad (2)$$

with $j$ the complex number, $\omega$ the frequencies contained on the excitation $u(t)$, $U(j\omega)$ and $Y(j\omega)$ the Fourier transforms of the input and output signals, respectively, and $G(j\omega)$ the frequency response of the system. As shown in Bendat and Piersol (1980), solving Equation (2) provides the frequency expression of the best linear approximation

$$G_{BLA}(j\omega) = \frac{\mathbb{E}[Y(j\omega)U^H(j\omega)]}{\mathbb{E}[|U(j\omega)|^2]} = \frac{S_{yu}(j\omega)}{S_{iuu}(j\omega)}, \quad (3)$$

where $S_{iuu}(j\omega)$ is the auto-spectral density of $u(t)$ and $S_{yu}(j\omega)$ is the inter-spectral power density between $u(t)$ and $y(t)$. The operator $(\cdot)^H$ represents the conjugate transpose.

For a periodic nonlinear system<fnref n="3" /> that can be represented by a Volterra series, the output spectrum $Y(j\omega)$ (as shown in Fig. 2) can be decomposed into three distinct components:

<figure><img src="image_3.png" /><figcaption>Fig. 2. Approximation of periodic nonlinear system by a linear model in presence of noise and nonlinear distortions.</figcaption></figure>

* $Y_{lin}(j\omega)$ represents the output of the underlying linear model. This component captures the behavior that can be attributed to linear dynamics in the system,

* $Y_s(j\omega)$ represents the contributions of nonlinear distortions,

* $V(j\omega)$ contains the contributions from noise sources present in the system. This component represents the effects of noise on the system's output spectrum.

By decomposing the output spectrum into these components, it becomes possible to analyze and understand the specific contributions of linear dynamics, nonlinear distortions, and noise, in the periodic nonlinear system. This decomposition provides insights into the different factors influencing the system's behavior and helps in the identification and characterization of the system's underlying dynamics.

Thus, considering a nonlinear system excited by a Gaussian excitation,<fnref n="4" /> denoted as $u(t)$, it can be established the following relation for the output spectrum $Y(j\omega)$:

$$Y(j\omega) = Y_{lin}(j\omega) + Y_s(j\omega) + V(j\omega). \quad (4)$$

In the given context, $\frac{Y_{lin}(j\omega)}{U(j\omega)}$ represents the best linear approximation, denoted as $G_{BLA}(j\omega)$. The noises, $V(j\omega)$, are aperiodic and independent of the input signal $u(t)$. They are specifically measurement noise and are often modeled as additive, stationary, and zero mean Gaussian white or colored noise (Pintelon & Schoukens, 2012). On the other hand, the nonlinear distortions, $Y_s(j\omega)$, depend on the specific input signal $u(t)$ used. Pintelon and Schoukens (2012) have extensively studied the characteristics of these nonlinear distortions and have identified several properties. The most significant are outlined below:

* the nonlinear distortions have zero mean and are uncorrelated with $u(t)$,

* for a given periodic system, the nonlinear distortions are also periodic with the same periodicity as $u(t)$,

* the nonlinear distortions follow, asymptotically, a normal distribution.

Various approaches have been proposed in the literature for estimating the optimal linear approximation $G_{BLA}(j\omega)$. One notable approach, introduced by Pintelon and Schoukens (2012), is the robust method. This approach leverages Gaussian excitation by performing multiple realizations.<fnref n="5" /> The motivation behind this approach stems from the periodic nature of nonlinear distortions. Since nonlinear distortions exhibit periodicity and a zero mean, conducting several experiments with different realizations allows for averaging over the realizations and mitigating the impact of nonlinearities (Pintelon, Rolain, Vandersteen, & Schoukens, 2004). Additionally, as noises (aperiodic) and nonlinear distortions (periodic) inherently possess distinct characteristics, it becomes possible to independently reduce their impact using periodic signals as recommended by Pintelon et al. (2004).

<footnote n="3">³ A periodic nonlinear system refers to a system in which the noise-free output maintains the same periodicity as its input (Pintelon & Schoukens, 2012).</footnote>

<footnote n="4">⁴ A Gaussian excitation refers to an excitation signal characterized by a Gaussian distribution (Pintelon & Schoukens, 2012).</footnote>

<footnote n="5">⁵ A realization refers to the random generation of an excitation signal following a specific probability distribution.</footnote>
Several types of candidate signals are used in an identification process. These include pseudo-random binary signals, chirp signals, multisine signals, and Gaussian noise. In this work, the multisine signal has been chosen, due to the advantages it offers over other signals, particularly in detecting, classifying, and quantifying nonlinearities (Pintelon & Schoukens, 2012) (note that detection and classification are not discussed in this article), as well as in reducing their impact on the estimated frequency response, as explained by Pintelon et al. (2004).

Through meticulous design of the multisine excitation properties, including frequencies, amplitudes, and phases of the sinusoidal components, it becomes feasible to discern and analyze the periodic contributions attributed to the nonlinear distortions, as will be elaborated in the subsequent sections.

## 2.2. Random-phase multisine signal

A multisine signal is generated by combining multiple sinusoidal signals that are phase-shifted with respect to each other. When considering a fixed number F of harmonics, the multisine signal can be expressed as follows:

$$u(t) = \sum_{k=1}^{F} A_k \cos(2\pi f_k t + \phi_k), \quad (5)$$

In this equation, $A_k$ represents the amplitudes of each cosine component, $f_k$ corresponds to integer multiples of the fundamental frequency $f_0$ (i.e., $f_k = k f_0$), and $\phi_k$ denotes the phase of each component.

The amplitudes, the frequencies and the phases of a multisine signal play crucial roles in determining its characteristics and properties, where:

* the amplitudes of the sinusoidal components directly influence the power spectral density (PSD) of the multisine signal (Pintelon & Schoukens, 2012). The levels of amplitudes determine the energy distribution across different frequencies. Controlling the amplitudes allows to shape the frequency response and allocate more energy to specific frequency bands of interest,

* the frequencies of the sinusoidal components define the spectral content of the multisine signal. The choice of frequencies determines which frequency bands are excited and analyzed in the system. A denser frequency grid provides higher resolution, allowing for more detailed frequency domain analysis,

* the phases of the sinusoidal components significantly influence the distribution of the multisine signal. Varying phase configurations can result in distinct distribution properties, including symmetry and skewness. In essence, the phases play a crucial role in shaping the signal's statistical characteristics and temporal behavior (Pintelon & Schoukens, 2012).

The choice of phases $\phi_k$ holds crucial importance based on the aforementioned properties of the multisine signal $u(t)$. In the context of the best linear approximation theory, the use of a signal with a Gaussian distribution is fundamental (Pintelon & Schoukens, 2012). One common approach involves randomly assigning the phases $\phi_k$ from either a normal or uniform distribution over the interval $[0; 2\pi]$. It is worth noting that both uniform and normal distributions result in a Gaussian distribution of the signal $u(t)$. Consequently, this type of signal is referred to as a random-phase multisine signal.

## 2.3. Estimation of the best linear approximation via a robust approach

The robust approach relies on the use of multisine signals with random phases. By incorporating this approach, Eq. (3) can be computed through the consideration of M independent realizations of the multisine signal. The random variation in the phases $\phi_k$ introduces diverse signals $u(t)$ and, consequently, diverse responses. By exciting the system with M independent realizations of the multisine signal, the

impact of both nonlinear distortions and noise contributions can be effectively reduced through averaging across the M realizations. This averaging process allows for the reduction of the influence for both $Y_s(j\omega)$ and $V(j\omega)$. In essence, the robust approach facilitates quantitative and qualitative analysis of the levels of nonlinear distortions and noise, enabling a comprehensive assessment (Pintelon et al., 2004).

When considering a multisine signal with random phase, Eq. (3) can be modified to:

$$G_{BLA}(j\omega) = \frac{\mathbb{E}[Y(j\omega)U^H(j\omega)]}{\mathbb{E}[|U(j\omega)|^2]} = \mathbb{E}\left[\frac{Y(j\omega)}{U(j\omega)}\right]. \quad (6)$$

Here, Eq. (6) is evaluated specifically for certain frequencies $\omega_k = 2\pi k/(NT_s)$, which correspond to rank k in the spectrum of the multisine signal. The variables N and $T_s$ represent the number of samples and the sampling time, respectively. By measuring the pairs $\{U(k), Y(k)\}$ over P periods and M realizations, as presented in Fig. 3, the best linear approximation of the nonlinear system's frequency response can be determined using the following equation:

$$G_{BLA}(j\omega_k) = \frac{1}{M} \sum_{m=1}^{M} G^{[m]}(j\omega_k), \quad (7)$$

where $G^{[m]}(j\omega_k)$ represents the frequency estimate for a single realization estimated using the previous relationship (6).

The contributions of noise are influenced by both the number of realizations (M) and the number of measured periods (P) (Pintelon et al., 2004). However, the impact of nonlinear distortions is solely dependent on the number of realizations (M) (Pintelon & Schoukens, 2012). In other words, by repeating the identification experiment M times with different random $\phi_k$ phases, both $\hat{\sigma}_{DN}^2$ (the variance of nonlinear distortions) and $\hat{\sigma}_B^2$ (the variance of noise) can be reduced. It is worth noting that a minimum of two periods P must be measured, and at least two realizations M need to be performed to calculate $\hat{\sigma}_{DN}^2(j\omega_k)$ and $\hat{\sigma}_B^2(j\omega_k)$.

## 2.4. Best linear approximation in closed-loop setup

Closed-loop identification offers several advantages, including ensuring the stability of the system during the identification process and enabling identification around specific operating points. However, in this configuration, the control signal $u(t)$ becomes correlated with the measurement noise $v(t)$, leading to biased estimates. To address this issue, it is crucial to use an exogenous signal (uncorrelated with noise) to avoid potential estimation biases. In the current case, due to the nature of the experimental setup under study, the exogenous excitation signal is applied as the reference signal $r(t)$. To formulate the identification problem, the general scheme of a system operating in a closed loop has been considered (see Fig. 4).

In the case of nonlinear systems, apart from noise, nonlinear distortions affecting the output $y(t)$ are also fed back into the control loop. Consequently, the input $u(t)$ is influenced by the nonlinear distortions of the system. To address this issue, several closed-loop identification techniques have been proposed, as documented by Friis, Tarpø, Katsanos, and Brincker (2021), Gustavsson, Ljung, and Söderström (1977), Van den Hof (1998), Ljung (1986), Pintelon, Louarroudi, and Lataire (2013), Pintelon and Schoukens (2013, 2013), Pintelon, Schoukens, and Lataire (2020), Söderström and Stoica (1988), Wellstead (1977). Among these techniques, the indirect method proposed by Gustavsson et al. (1977) and employed later by Pintelon and Schoukens (2013) in the BLA context, offers a solution for obtaining the best non-parametric linear approximation of a nonlinear system operating in a closed loop. The indirect method involves estimating two independent transfers: the transfer between $r(t)$ and $u(t)$ and the transfer between $r(t)$ and $y(t)$, enabling the identification of the optimal linear approximation as demonstrated by Pintelon and Schoukens (2013). The equation for
<figure><img src="image_4.png" /><figcaption>Fig. 3. Robust approach to estimate the best linear approximation (Pintelon & Schoukens, 2012).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. General scheme of a closed-loop system.</figcaption></figure>

the best linear closed-loop approximation, as shown in Eq. (3), can be expressed as follows:

$$G_{BLA}(j\omega) = \frac{S_{yr}(j\omega)/S_{rr}(j\omega)}{S_{ur}(j\omega)/S_{rr}(j\omega)} = \frac{G_{r\to y}(j\omega)}{G_{r\to u}(j\omega)}. \qquad (8)$$

Here, $r(t)$ represents a multisine signal that is decorrelated from both the noise $v(t)$ and nonlinear distortions. Eq. (8) indicates that an unbiased estimate of $G_{BLA}(j\omega)$ can be obtained by calculating the ratio between the best linear approximations $G_{r\to y}(j\omega)$ and $G_{r\to u}(j\omega)$. It is important to note that the robust approach described in Section 2.3 can be employed to calculate $G_{r\to y}(j\omega)$ and $G_{r\to u}(j\omega)$, thereby reducing the impact of noise and nonlinear distortions during the estimation of the closed-loop frequency response $G_{BLA}(j\omega)$.

### 3. Design of multisine excitations for closed-loop identification

Pintelon and Schoukens (2012) suggest choosing random phases in the interval $[0; 2\pi]$ for estimating the best linear approximation of the frequency response. However, random phases can lead to a multisine signal with potentially high amplitudes, which can generate a control signal that exceeds the saturation level, thus resulting in nonlinear distortions. This is often the case in robotic systems where saturations are active in internal control loops to protect actuators from excessive current or speed. Therefore, in such cases, the objective of generating excitation is to find an admissible excitation that maximizes the quality of non-parametric frequency response estimation while respecting the imposed constraints.

A simple solution to respect the system constraints and prevent excitation from saturating is to normalize it. This type of normalization involves scaling all amplitude values so that the excitation signal is

lower than or equal to the saturation value. Unfortunately, such a solution degrades the crest factor (CF) value, which is a metric that quantifies the power richness of a signal often used in constrained non-parametric identification (Boyd, 1986; Schoukens, Guillaume, & Pintelon, 1991).

For a signal $u(t)$ defined over the time interval $[0, T]$, its crest factor is defined as follows:

$$CF = \frac{\max_{t \in [0, T]} |u(t)|}{\sqrt{\frac{1}{T} \int_{0}^{T} |u(t)|^2 dt}}. \qquad (9)$$

The ratio in Eq. (9) quantifies the amount of energy that is injected into the system, for a given amplitude: the square of the CF is inversely proportional to the energy of the signal. Therefore, a lower crest factor indicates a higher energy injection into the system. In order to obtain a powerful signal for estimating the frequency response of a system while respecting its limits, especially in the case of multisine signals, achieving a minimal crest factor through appropriate phase selection is often necessary (Guillaume, Schoukens, Pintelon, & Kollar, 1991; Van Den Bos, 1987). Several research studies (Guillaume et al., 1991; Ojarand & Min, 2017; Schoukens et al., 1991; Van Den Bos, 1987; Van der Ouderaa, Schoukens, & Renneboog, 1988; Yang et al., 2015) have demonstrated the possibility of reducing the crest factor of a multisine signal using iterative numerical methods, albeit with increased computation time. Among these methods, the approach proposed by Van der Ouderaa et al. (1988) consistently produces improved results regardless of the initialization of the multisine signal's parameters. Van der Ouderaa's method involves adjusting the phases of the multisine signal using a heuristic approach inspired by the Gerchberg-Saxton algorithm (Gerchberg, 1972). This iterative method performs clipping of the multisine signal in the time domain. For more detailed information on the clipping technique, readers can refer to Boukhebouz, Mercère, Grossard, and Laroche (2021).

The optimization of the crest factor mentioned above has demonstrated great effectiveness in obtaining a power-rich multisine excitation that respects the limitations imposed in an open-loop scheme. However, for many systems, the identification process is performed in closed-loop due to the following reasons:

* the instability of the open-loop system.

* the inability to isolate the system from its controller,
<figure><img src="image_6.png" /><figcaption>Fig. 5. Schematic diagram of the crest factor optimization in closed loop setup (Boukhebouz et al., 2021).</figcaption></figure>

Moreover, in the case of the considered robotic system, the multi-sine excitation signal is not directly applied as the control signal $u(t)$ but is injected at the reference level $r(t)$. It is therefore necessary to ensure that the control signal $u(t)$ does not exceed the saturation thresholds without deteriorating the crest factor.

Addressing such a challenge in non-parametric model estimation has been studied by Boukhebouz et al. (2021). The proposed design approach for multisine excitation in closed-loop estimation iteratively adjusts the multisine excitation at the reference level $r(t)$ so that the control signal $u(t)$ converges to a desired signal $u_{des}(t)$ satisfying:

$$ \max_{t} |u_{des}(j\omega_k)| \le u_{max}, \qquad (10) $$

where $u_{max}$ represents the saturation threshold.

The approach proposed by Boukhebouz et al. (2021) can be summarized in the two following steps:

1. design a desired signal $u_{des}(t)$ with an optimized crest factor that satisfies Eq. (10),

2. shape the reference signal $r(t)$ so that the control signal $u(t)$ reaches $u_{des}(t)$ in closed-loop using the following equation:

$$ R_{des}(j\omega) = G_{r \to u}^{-1}(j\omega) U_{des}(j\omega). \qquad (11) $$

In Eq. (11), $R_{des}(j\omega)$ represents the desired reference spectrum, $U_{des}(j\omega)$ represents the desired control signal spectrum, and $G_{r \to u}(j\omega)$ denotes the frequency response between $r(t)$ and $u(t)$.

These two steps are performed offline. For the first point, given the system-imposed thresholds and the frequency band of interest, the signal $u_{des}(t)$ is defined by selecting the $A_k$, $\phi_k$, and $f_k$ parameters for $k = 1, 2, \dots, K$ using the approach proposed by Van der Ouderaa et al. (1988) to optimize the crest factor. Once the crest factor optimization is achieved, we move on to the second step, which is an iterative process depending on the closed-loop setup, including knowledge of the controller structure and the signals to be measured. The schematic diagram of Fig. 5 illustrates the overall approach.

## 4. Experimental validation of the identification approach

This section reports the implementation of the proposed approach on an experimental system. The objective is to experimentally identify a non-parametric frequency model of a robotic system. In order to achieve this objective, we begin by providing a brief description of the functioning principle of the actuation system and its control architecture in the first part. In the second part, we present the identification results obtained using the excitation generation method proposed in Section 3, along with the estimation approach based on the best linear approximation. The final part of the study is dedicated to validating the estimated frequency models by identifying grey-box parametric models.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Experimental setup (CEA-LIST).</figcaption></figure>

### 4.1. Experimental setup description

The experimental setup, depicted in Fig. 6, comprises a single-axis robotic arm, which incorporates the following components:

* a low-inertia synchronous motor ($J_m = 10^{-4}$ kg.m$^2$) equipped with an 11-bit incremental encoder to precisely measure the motor position $q_{mot}(t)$,

* a specially designed mechanical transmission, known as the screw-cable mechanism (Garrec, 2008),

* a robotic arm mounted at the joint pulley, with an inertia of $J_a = 2.23$ kg.m$^2$, calculated based on CAD software,

* two mechanical stoppers positioned at the top and bottom to restrict the arm's range of motion,

* an electronic board for controlling the arm's movement.

The hardware architecture of the experimental setup consists of the following components:

* a workstation responsible for generating the reference signals $r(t)$,

* a dedicated real-time operating system, which handles the calculation of digital control signals.

The mechanical transmission utilized in this experimental setup is designed to adapt the motor power to the joint's requirements in terms of available torque $\tau_m(t)$ and speed $\dot{q}(t)$, employing a specific operating mode. It consists of the following elements:
<figure><img src="image_8.png" /><figcaption>Fig. 7. Schematic diagram of the robotic arm control.</figcaption></figure>

* an upstream pulley,

* a joint pulley (downstream pulley) equipped with a 23-bit encoder for measuring the joint position $q(t)$,

* a toothed belt that connects the motor shaft to the input pulley,

* cables that establish a connection between the input and output pulleys,

* allowable joint positions of the robotic arm fall within the interval of [0; 2.3] rad due to the presence of mechanical stoppers.

For the robotic arm under consideration, the transfer functions respectively between: (i) the motor torque $\tau_m(t)$ and the motor position $y_1(t) = q_{mot}(t)$, and (ii) the motor torque $\tau_m(t)$ and the joint position $y_2(t) = q(t)$, are inherently unstable in open loop due to the presence of an integration. To stabilize $y_1(t)$ and $y_2(t)$, a cascade servo-control approach is employed, as illustrated in Fig. 7. The controller structure can be divided into two parts:

* an internal loop featuring a PID controller in charge of regulating the motor torque $\tau_m(t)$,

* an external position loop implemented, where the positions are controlled by another PID controller.

The motor position is measured using an encoder and transmitted to the workstation with a sampling period of $T_s = 1$ ms.

In terms of control, the key focus lies on studying the transfer functions between the motor torque $\tau_m(t)$ and four variables: the motor position $q_{mot}(t)$, motor velocity $\dot{q}_{mot}(t)$, joint position $q(t)$, and joint velocity $\dot{q}(t)$. Velocities are commonly reconstructed using a suitable high-resolution numerical differentiation method.

To comprehensively analyze the behavior of the robotic axis, the dynamic contributions are divided into two parts: upstream and downstream of the flexible transmission. For this purpose, two transfer functions of interest are defined for position control. Firstly, the transfer function $G_{\tau_m \to q_{mot}}(s)$ is referred to as the collocated transfer function.<fnref n="6" /> This transfer function characterizes the behavior of the system upstream. Secondly, the transfer function $G_{\tau_m \to q}(s)$ is known as the non-collocated transfer function.<fnref n="7" /> This transfer function captures the dynamics of the system downstream.

The considered robotic arm is subjected to several sources of nonlinearities, including dry friction (Hamon et al., 2011), gravity, cable stiffness (Fichera & Grossard, 2017), and actuators saturations. These nonlinearities play a significant role in the selection of excitation signals and can introduce distortions in the frequency responses, thereby impacting the estimation process. For instance, friction can obscure flexible modes, while variable cable stiffness can cause shifts in the frequency responses (Boukhebouz et al., 2020).

In this work, the only nonlinear compensation activated in the control algorithm is gravity. This deliberate choice enables to identify and assess the impact of the nonlinearities on the estimation of frequency responses. By examining the system's behavior in the presence of these nonlinear effects, we aim to gain insights into their influence on the overall dynamics.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Arm configured at position $q = 1.23$ rad.</figcaption></figure>

## 4.2. Estimation of the best linear approximation of the collocated transfer function

For the frequency study, the focus lies in estimating the open-loop transfer function from $\tau_m$ to $q_{mot}$ that is critical for ensuring the stability of the robotic systems (De Luca & Book, 2008). Through prior analysis of the setup, it has been determined that the frequency responses are independent from the joint position. Based on this insight, we have deliberately chosen to stabilize the arm at a horizontal position (as shown in Fig. 8), where $q = 1.23$ rad, using the position feedback loop. By operating the arm around this specific position, we aim at estimating the frequency responses in the proximity of this operating point, which facilitates the analysis of the system dynamics.

The excitation signal $r(t)$ used in the present study is a multisine signal, carefully designed to maximize information content in the data set. The phases of the harmonics are chosen based on the procedure detailed in Section 3, adhering to motor torque saturation constraints with a maximum allowable torque of $u_{max} = 1.5$ N.m. The multisine excitation spans frequencies from 1 Hz to 60 Hz, selected to capture significant system dynamics (potential resonance frequencies). A frequency resolution of $f_{res} = 0.1$ Hz was chosen to balance FRF accuracy with computational efficiency, and a sampling period of $T_s = 1$ ms is used, which adheres to the Nyquist theorem and prevents aliasing.

Initially, the system is excited using a low-amplitude reference signal $A_k = 0.004$ rad to prevent motor torque saturation. Six periods of the signals $r(t)$, $u(t)$, and $q_{mot}(t)$ are then collected, but only the last four periods are retained for estimating the frequency responses (i.e., $P = 4$). This selection eliminates transient effects. By calculating the spectra $R(j\omega_k)$ and $U(j\omega_k)$, the frequency response $G_{r \to u}(j\omega_k)$ is estimated, to be used for the closed-loop excitation generation.

In the initial experiment, a crest factor value of $u(t)$ equal to 3.03 has been obtained. This value is relatively high and could be reduced to improve the quality of the multisine signal. As explained by Ojarand and Min (2017), a lower crest factor generally signifies a better distribution of signal power, which enhances the signal's effectiveness. To improve the quality of the signal $u(t)$, we have developed the following approach:

* use the initial estimate of $G_{r \to u}(j\omega_k)$ obtained from the previous experiment,

<footnote n="6">⁶ A collocated transfer occurs when the sensor and actuator are physically placed in the same location, enabling control at the measurement point.</footnote>

<footnote n="7">⁷ When the sensor and actuator are placed in different positions, it is referred to as a non-collocated transfer.</footnote>
Desired control signal with CF=1.54

<figure><img src="image_10.png" /><figcaption>Fig. 9. Time evolution of $u_{des}(t)$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Crest factor variation of $u(t)$.</figcaption></figure>

* design a desired input signal $u_{des}(t)$ (refer to Fig. 9) using an optimization process outlined in Section 3. This optimization aims to maximize the crest factor of $u_{des}(t)$,

The goal of this iterative procedure is to find a suitable signal that results in a control signal $u(t)$ that is as close as possible to the desired signal $u_{des}(t)$. This iterative process allows to improve the resemblance between the actual control signal and the desired one, enhancing the quality of $u(t)$ for subsequent analysis.

After conducting twelve iterations of the excitation generation algorithm, the following results are observed:

* Fig. 10 represents the variations in the crest factor of the control signal $u(t)$. The best value for the crest factor is achieved at the third iteration, where CF = 2.16. Initially starting at CF = 3.03, the crest factor of $u(t)$ decreases to 2.16, which represents a 29% reduction, after three iterations. It then increases slightly in the fourth iteration and eventually stabilizes around an average value of 2.5 after twelve iterations,

* Fig. 11 presents the frequency response $G_{r→u}(j\omega_k)$ for each iteration. The initial frequency response (depicted by blue triangles) exhibits some low-frequency outliers and fails to accurately capture the flexible modes at frequencies of 4 Hz and 10 Hz. However, starting from the second iteration, the frequency responses become nearly identical and enable more precise detection of the flexible modes compared to the first iteration.

These results demonstrate the effectiveness of the iterative excitation generation process in improving the quality of the control signal

Frequency response $G_{r→u}(j\omega)$ for 12 iterations of the algorithm

<figure><img src="image_12.png" /><figcaption>Fig. 11. Frequency response $G_{r→u}(j\omega_k)$ obtained for 12 iterations.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Spectra of $r(t)$ and $u(t)$ of the first iteration compared with the spectra of the third iteration.</figcaption></figure>

$u(t)$ and refining the estimation of the frequency response $G_{r→u}(j\omega_k)$, particularly in accurately capturing the system's flexible modes.

After evaluating the results, it can be observed that the best excitation signal is obtained at the third iteration, with a crest factor of 2.16. To assess the impact of the excitation signal tuning algorithm, we compare the reference spectrum $R(j\omega_k)$ and the input spectrum $U(j\omega_k)$ for the first and third iterations in Fig. 12.

In Fig. 12a, one can observe that the excitation signal tuning algorithm has modified the initial reference spectrum (represented by blue circles) from a flat spectrum to a new spectrum represented by black circles. This adjustment results in a flatter spectrum for the $u(t)$ signal. Additionally, Fig. 12b illustrates that the amplitude spectrum of the $u(t)$ signal has increased, particularly at low frequencies, where the spectrum has nearly doubled. These results highlight the effectiveness of the proposed approach, as a power-rich closed-loop $u(t)$ signal that avoids exceeding saturation thresholds has been successfully generated. The modification of the excitation signal shape has resulted in a more balanced and enhanced input spectrum.

In order to estimate the BLA of $G_{r→u}(j\omega_k)$, M desired input signals $u_{des}(t)$, with optimized crest factors, have been generated. Using the estimate of $G_{r→u}(j\omega_k)$ obtained at the third iteration, M reference signals $r(t)$ are calculated, based on the following equation:

$$R(j\omega_k) = \hat{G}_{r→u}^{-1}(j\omega_k) U_{des}(j\omega_k), \qquad (12)$$

where $R(j\omega_k)$ is the reference spectrum, $\hat{G}_{r→u}(j\omega_k)$ is the estimated frequency response from the third iteration, and $U_{des}(j\omega_k)$ is the desired input spectrum. By inverting the estimated frequency response and multiplying it with the desired input spectrum, the reference spectrum $R(j\omega_k)$ is obtained for each frequency $\omega_k$.
<figure><img src="image_14.png" /><figcaption>Fig. 13. Non-parametric analysis of nonlinear distortions and noise.</figcaption></figure>

The results for **P** = 4 and **M** = 2 are presented in Fig. 13a. The estimated frequency response $G_{\tau_m \to q_{mot}}(j\omega_k)$ is plotted in blue, while the standard deviation of nonlinear distortions $\sigma_{DN}(j\omega)$ is presented as red circles, and the standard deviation of noise $\sigma_B(j\omega)$ is shown as green dots.

From the plots, it can be observed that the nonlinear distortions are nearly at the same level as the frequency response at low frequencies and between 2 Hz and 10 Hz. This leads to distortions in the estimated $G_{\tau_m \to q_{mot}}(j\omega_k)$, especially at resonance and anti-resonance frequencies, as illustrated in Fig. 13a. Additionally, it can be seen on the plot that the level of nonlinear distortions is approximately 28 dB higher than the noise level. Based on this observation, it can be concluded that, for this particular system, the level of nonlinearities is more significant compared to the noise level. The presence of these nonlinearities poses challenges in accurately estimating the frequency response $G_{\tau_m \to q_{mot}}(j\omega_k)$, particularly in the presence of resonance and anti-resonance frequencies.

To mitigate the impact of nonlinearities on the estimated frequency response, a second experiment has been designed with a higher value of **M** = 40 and the same number of periods **P** = 4.

The experimental results (presented in Fig. 13b) demonstrate that increasing **M** to 40 significantly reduces the nonlinear distortions. The standard deviation of the distortions (represented by red circles in Fig. 13b) is observed to be 37 dB below the frequency response in black. This reduction is achieved by averaging over the **M** realizations, leading to a more accurate estimation. Moreover, the reduction in nonlinear effects enables better distinction between resonance and anti-resonance, as illustrated in the frequency responses of Fig. 14.

Based on these results and the analysis of standard deviations of nonlinear distortions and noise, we conclude that the frequency response in black represents an accurate linear approximation for the transfer between motor torque and motor position. The gain and phase characteristics of this frequency response align with the typical behavior of this type of actuation system, where anti-resonance is followed by resonance (De Luca & Book, 2008). In this specific case, the first resonance occurs at a frequency of 10.1 Hz, while the anti-resonance is observed at a frequency of 7.5 Hz. Additionally, a second flexible mode emerges in the frequency range of 40 to 50 Hz (see Fig. 14).

## 5. Validation of results through a parametric estimation

The frequency model can be used directly for some applications, such as passivity analysis (Abroug & Laroche, 2015) or modal analysis (Guillaume et al., 2003). For other applications, such as control, starting by identifying a non-parametric model by the proposed method, is an efficient first step towards the identification of a physical model. In this section, such an identification of a physical model, based

<figure><img src="image_15.png" /><figcaption>Fig. 14. Comparison between $G_{\tau_m \to q_{mot}}(j\omega_k)$ obtained for **M**=2 and for **M**=40.</figcaption></figure>

on the previous non-parametric model, is made using the frequency response $G_{\tau_m \to q_{mot}}(j\omega_k)$ between motor torque and motor position. The validity of the non-parametric analysis is confirmed by comparing the estimated physical parameters with prior knowledge of the setup.

### 5.1. Physical model of the experimental setup

Transfer $G_{\tau_m \to q_{mot}}(j\omega_k)$ includes a pure integration between speed and position. Considering the speed $\dot{q}_{mot}(t)$ instead of the position as the output of the model allows to reduce the order of the model. To achieve this, the frequency response $G_{\tau_m \to q_{mot}}(j\omega_k)$ has been multiplied by $j\omega$, resulting in the frequency response between motor torque and motor speed $G_V(j\omega_k) = j\omega_k \cdot G_{\tau_m \to q_{mot}}(j\omega_k)$.

A linear time-invariant (LTI) model that governs the dynamics of $G_V(s)$ has been developed, represented by the following set of equations:

$$ \begin{cases} J_a \ddot{q} + f_{vt} (\dot{q} - \dot{q}_{mot}) + k (q - q_{mot}) + \tau_{besi}(q) = 0, \\ J_m \ddot{q}_{mot} + f_{vm} \dot{q}_{mot} - f_{vt} (\dot{q} - \dot{q}_{mot}) - k (q - q_{mot}) = \tau_m, \end{cases} \quad (13) $$

based on the physical parameters $J_a$ (inertia of the arm), $J_m$ (inertia of the motor), $k$ (stiffness coefficient), $f_{vt}$ (friction coefficient of the transmission), and $f_{vm}$ (friction coefficient of the motor).

The following minimal state-space representation has been chosen: $x = (q - q_{mot} \quad \dot{q} \quad \dot{q}_{mot})^T$, the input vector $u(t) = \tau_m(t)$, and the output vector $y(t) = \dot{q}_{mot}(t)$. Assuming that $\tau_g(q)$ is perfectly compensated by the control, the system can be represented by the following state-space model:

$$ \begin{aligned} \dot{x}(t) &= A(\theta)x(t) + B(\theta)u(t), \\ y(t) &= C(\theta)x(t), \end{aligned} \quad (14) $$

where matrices $A(\theta)$, $B(\theta)$, and $C(\theta)$ are defined as:

$$ A(\theta) = \begin{bmatrix} 0 & 1 & -1 \\ \theta_1 & \theta_2 & -\theta_2 \\ \theta_3 & \theta_4 & \theta_5 \end{bmatrix}, \quad B(\theta) = \begin{bmatrix} 0 \\ 0 \\ \theta_6 \end{bmatrix}, \quad C(\theta) = \begin{bmatrix} 0 & 0 & 1 \end{bmatrix}, \quad (15) $$

with the parameter vector $\theta \in \mathbb{R}^6$ is defined as:

$$ \begin{aligned} \theta_1 &= -\frac{k}{J_a}, & \theta_2 &= -\frac{f_{vt}}{J_a}, & \theta_3 &= \frac{k}{J_m}, \\ \theta_4 &= \frac{f_{vt}}{J_m}, & \theta_5 &= -\frac{f_{vm} + f_{vt}}{J_m}, & \theta_6 &= \frac{1}{J_m}. \end{aligned} \quad (16) $$

The matrices $A$, $B$, and $C$ contain the parameter values $\theta_i$, as well as 0's and 1's that do not require estimation. This grey-box model, which is affine in its parameters, is identifiable using the method described by Glover and Willems (1974), which ensures local and global identifiability by leveraging parameter mappings from the system's transfer function. Consequently, the estimation of its parameter vector $\theta$ yields a unique solution. This grey-box model is denoted $G_{GB}(\theta, s) = C(\theta)(sI - A(\theta))^{-1}B(\theta)$ in the sequel.
## 5.2. $H_\infty$ norm approach to estimate state space model parameters

Let $F(\theta, \Omega) \in C^{n_{F_{max}} \times 1}$ be the complex vector defined as:

$$F(\theta, \Omega) = [F(\theta, j\omega_1) \quad F(\theta, j\omega_2) \quad \dots \quad F(\theta, j\omega_{n_{F_{max}}})]^T, \quad (17)$$

where $F(\theta, j\omega_k) = G_v(j\omega_k) - G_{GB}(\theta, j\omega_k)$ is the complex gain at frequency $\omega_k$ and $I$ is the identity matrix.

The problem of estimating the parameter vector $\theta$ from the frequency response $G_v(j\omega)$ can be formulated as follows: given the best linear approximation of the frequency response $G_v(j\omega_k)$, estimate the parameter vector $\theta$ of the grey-box LTI state model, defined by Eq. (15), by minimizing the following cost function:

$$V_{H_{\infty}}(\theta) = \max_{\omega_k} |F(\theta, j\omega_k)|. \quad (18)$$

This cost-function corresponds to the $H_\infty$ norm of $F(\theta, s)$ evaluated over the set of frequencies $\omega_k$, following the works by Bergamasco and Lovera (2013), Vizer and Mercere (2014), Vizer, Mercère, and Laroche (2015), Vizer, Mercère, Prot, and Laroche (2016). These techniques leverage developments in nonsmooth and nonconvex optimization (Apkarian & Noll, 2006; Apkarian, Noll, & Prot, 2008, 2009; Fuduli, Gaudioso, & Giallombardo, 2004) to enable the synthesis of controllers based on the $H_\infty$ norm. Unlike classical approaches that rely on the Kalman-Yakubovich-Popov lemma (Rantzer, 1996) to design full-order and unstructured controllers, these new techniques offer the capability to design structured and reduced-order controllers.

In the identification approach, inspired by Vizer et al. (2016), the use of a black-box model is necessary. However, since the cost function (18) is defined in terms of the frequency response $G_v(j\omega)$, an intermediate step is required to address this issue. The identification procedure is carried out in two stages:

1. estimate a stable black-box model $G_{BB}(s)$ from the frequency response $G_v(j\omega)$,

2. tune the parameters vector $\theta$ in order to minimize (18) using the proximity control algorithm, where $F(\theta, j\omega_k) = G_{BB}(j\omega_k) - G_{BG}(\theta, j\omega_k)$.

In the first stage, determining a linear black-box model from the frequency response is a relatively straightforward problem to solve. Subspace algorithms, such as the N4SID algorithm (Overschee, Moor, & Favoreel, 1997), have been successful in accurately estimating linear time-invariant state models from frequency data. The N4SID algorithm, available in Matlab's System Identification Toolbox, has been chosen to estimate black-box models.

In the second stage, the optimization of the cost function $V_{H_\infty}(\theta)$ is performed using the systune function provided by Matlab's Robust Control Toolbox (Gahinet & Apkarian, 2013). The systune function enables the tuning of transfer-function parameters to achieve desired performance specifications, such as minimizing the $H_\infty$ norm in the present case.

## 5.3. Estimation results

This section presents the results of the parameter vector estimation obtained through the optimization of the $H_\infty$ norm (using the systune function). To assess the sensitivity of the optimization algorithms to the initialization of the vector $\theta$, 1000 Monte Carlo runs are performed. For each run, the parameters $\theta_i$ are initialized by sampling from a uniform distribution over the interval $\pm 50\%$, as follows:

$$\theta_i^{\text{init}} = \theta_i^{\text{pr}} (0.5 \pm \text{rand}), \quad (19)$$

where $\theta_i^{\text{pr}}$ is chosen based on prior knowledge, as presented in Section 4.1. The function `rand` generates a random value with a uniform distribution in the interval $[0; 1]$.

As mentioned earlier, the estimation of the parameter vector $\theta$ using the `systune` function involves an intermediate step to estimate a

<table><caption>**Table 1**
Median, mean and standard deviation obtained with systune.</caption><thead><tr><th>Parameter</th><th>Median</th><th>Mean</th><th>Standard deviation</th></tr></thead><tbody><tr><td>θ̂<sub>1</sub></td><td>-2191.9</td><td>-2191.9</td><td>1.19e-07</td></tr><tr><td>θ̂<sub>2</sub></td><td>-6.476</td><td>-6.476</td><td>2.67e-09</td></tr><tr><td>θ̂<sub>3</sub></td><td>1719.4</td><td>1719.4</td><td>1.25e-07</td></tr><tr><td>θ̂<sub>4</sub></td><td>8.753</td><td>8.753</td><td>2.91e-09</td></tr><tr><td>θ̂<sub>5</sub></td><td>-11.659</td><td>-11.659</td><td>2.96e-09</td></tr><tr><td>θ̂<sub>6</sub></td><td>0.249</td><td>0.249</td><td>1.29e-11</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Histograms of normalized $\theta_{n_i}$ parameters obtained with systune.</figcaption></figure>

black-box model. In the present case, the frequency response function $G_v(j\omega)$ has been approximated with a third-order model using the n4sid function. The identified stable black-box model $G_{BN}(s)$ has the following state matrices:

$$A = \begin{bmatrix} -0.4067 & 706 & 3.519 \\ -0.9727 & -15.95 & -1689 \\ 0.00321 & 1.899 & -1.778 \end{bmatrix}, \quad B = \begin{bmatrix} -0.5417 \\ 0.01113 \\ -0.00044 \end{bmatrix}, \quad (20)$$

$$C = \begin{bmatrix} -0.4473 & 0.6371 & 0.5305 \end{bmatrix}, \quad D = 0.$$

Once the black-box model is obtained, the next step is to estimate the parameters $\theta$ of the grey-box model (15) from the black-box model (20) using the `systune` function. In 1000 Monte-Carlo runs, no outliers has been encountered, which would correspond to abnormal values of the parameters. This is a clear advantage of the $H_\infty$ norm over the more usual $H_2$, as reported by Vizer et al. (2016).

In Table 1 are presented the median value, mean value, and standard deviation of each $\theta$ estimated from the 1000 Monte Carlo runs. Additionally, Fig. 15 shows the normalized histograms, where the random variables are subtracted by their mean and divided by their standard deviation, to obtain variables with a mean of 0 and variance of 1.

From the results presented in Table 1, it can be observed that the mean is equal to the median, indicating the absence of outliers. Additionally, all the standard deviations are very small compared to the mean values, showing a low sensitivity to the initialization. The histograms exhibit Gaussian distributions, as can be seen on the quantile-quantile plot shown in Fig. 16, where the estimated values fit well to the reference, indicating a good fit to a normal distribution (Thode, 2002).

After selecting the best estimates of the parameter vector $\theta$, denoted as $\theta_{\text{systune}}$, the Bode diagrams of the identified grey-box model $G_{GB}(s)$ and the frequency response $G_v(j\omega)$ are together plotted in Fig. 17. One can see that the physical model fits well to the frequency model.

To quantify the gap between the two models, the point-wise discrepancy at each frequency $\omega$ is defined as:

$$\Delta_G(j\omega) = G_v(j\omega) - G_{GB}(j\omega).$$

and the three following metrics are evaluated:
<figure><img src="image_17.png" /><figcaption>Fig. 16. Quantile-quantile diagram for $\theta$ obtained with systune.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Comparison of frequency response with the grey-box model for the transfer between $\tau_m$ and $\dot{q}_{mot}$.</figcaption></figure>

1. The fit (in %) is an usual metric, defined as:

$$ \text{fit} = 100 \left( 1 - \frac{\| \Delta_G(j\omega) \|_2}{\| G_v(j\omega) - \bar{G}_v(j\omega) \|_2} \right) = 86.1 \, \%, \quad (21) $$

where $\bar{G}_v(j\omega)$ represents the average value of the measured response over all frequency points.

2. The Root-Mean-Square Error (RMSE) is another usual metric:

$$ \text{RMSE} = \sqrt{\frac{1}{n_{F_{\max}}} \sum_{k=1}^{n_{F_{\max}}} [\Delta_G(j\omega_k)]^2} = 0.0049. \quad (22) $$

3. The mean of the relative error, expressed in dB (MRE) is defined as:

$$ \text{MRE} = \frac{1}{n_{F_{\max}}} \sum_{k=1}^{n_{F_{\max}}} \left| \frac{\Delta_G(j\omega_k)}{\bar{G}_v(j\omega_k)} \right|_{\text{dB}} = -19.4 \, \text{dB}. \quad (23) $$

The fit and RMSE present satisfying values. But, as the amplitudes highly vary from low to high frequency, the high-amplitude data dominate in these metrics. The mean value of the relative errors (MRE), expressed in dB, enable to show that the errors remain at about -20 dB below the reference value, as can be seen in the left-hand side of Fig. 17.

Finally, the estimated values of the physical parameters are presented in Table 2. The results of parameter estimation for the grey-box model (15) closely match the a priori values of the physical parameters, particularly for the motor inertia $J_m$. These favorable results are a further indication of the effectiveness of the proposed identification methodology.

<table><caption>Table 2
Estimated physical parameters.</caption><thead><tr><td>Parameter</td><td>A priori values</td><td>Systune</td><td>Unit</td></tr></thead><tbody><tr><td>J<sub>a</sub></td><td>2.23</td><td>2.94</td><td>kg m<sup>2</sup></td></tr><tr><td>J<sub>m</sub></td><td>1e<sup>-04</sup></td><td>1.439e<sup>-04</sup></td><td>kg m<sup>2</sup></td></tr><tr><td>k</td><td>–</td><td>6994.9</td><td>N m/rad</td></tr><tr><td>f<sub>ct</sub></td><td>–</td><td>34.130</td><td>N m s/rad</td></tr><tr><td>f<sub>cm</sub></td><td>–</td><td>1.166</td><td>N m s/rad</td></tr></tbody></table>

## 6. Conclusion

The research work presented in this article focused on the frequency-based identification of cable-driven mechanical transmissions, commonly used in the design of robots characterized by flexible joints. The presence of cables in these mechanical transmissions can introduce mechanical flexibilities that affect both static and dynamic positioning accuracy. Additionally, specific design features can lead to certain nonlinear phenomena in the system dynamics, such as cable stiffness nonlinearity and friction. Thus, the first objective was to analyze and detect these nonlinearities in order to study their impact on frequency response identification. This objective was achieved by employing the theory of Best Linear Approximation (BLA) and utilizing multisine signals as excitation sources.

In the open-loop configuration, the transfer functions from the motor torque to the motor and joint positions respectively can exhibit instability, prompting the need for identification in a closed-loop setup. This approach helps avoid biases in estimating frequency responses that could arise from the correlation between measured signals (control signal and output signal), noise, and nonlinear distortions. In the closed-loop configuration, the actuator is protected by saturations on the motor torque, requiring the design of a suitable excitation that adheres to saturation thresholds, while maintaining the quality of the control signal, which is crucial for accurate identification. To achieve this, the crest factor metric has been used to quantify the control signal quality. The experimental results demonstrated that the proposed algorithm effectively satisfies system constraints while maintaining a low crest factor.

The system under investigation is influenced by two sources of nonlinearities: friction and cable stiffness. These factors introduce distortions in the frequency responses of the test bench. To mitigate these distortions, an advanced approach has been used for estimating frequency responses, based on the theory of BLA. The proposed solution consists of combining two algorithms: the first minimizes the crest factor metric of the control signal, and the second iteratively modifies the spectrum of the control signal to converge towards a desired signal by adjusting the reference spectrum.
The proposed identification methodology enabled to obtain accurate estimates of the co-located frequency response. Tests and studies on the experimental system showed that the obtained frequency response accurately represents the dynamics of the system and aligns well with analytical models dedicated to mechanical transmission axes. Furthermore, a parametric estimation based on the frequency response between the motor torque and the motor speed was performed to validate the experimental results. Using an estimation method based on the $H_\infty$ norm, the identified parametric model confirmed the relevance of the identification methodology, as the estimated physical parameters closely matched the prior knowledge of the system.

This validation on an experimental device further underscores the efficacy of the proposed experimental modeling methodology, which serves as a general identification method for this type of system and can be extended to other robotic systems.

## CRediT authorship contribution statement

**Bassem Boukhebouz:** Writing – original draft, Investigation, Validation. **Guillaume Mercère:** Writing – review & editing, Methodology, Supervision, Validation. **Mathieu Grossard:** Writing – review & editing, Investigation, Methodology, Supervision. **Édouard Laroche:** Writing – review & editing, Methodology, Supervision.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## References

Abroug, N., & Laroche, E. (2015). Transforming series elastic actuators into variable stiffness actuators thanks to structured $H_\infty$ control. In *2015 European control conference* (pp. 734–740). IEEE.

Apkarian, P., & Noll, D. (2006). Nonsmooth $H_\infty$ synthesis. *IEEE Transactions on Automatic Control*, **51**(1), 71–86.

Apkarian, P., Noll, D., & Prot, O. (2008). A trust region spectral bundle method for nonconvex eigenvalue optimization. *SIAM Journal on Optimization*, **19**, 281–306.

Apkarian, P., Noll, D., & Prot, O. (2009). A proximity control algorithm to minimize nonsmooth and nonconvex semi-infinite maximum eigenvalue functions. *Journal of Convex Analysis*, **16**, 641–666.

Bendat, J. S., & Piersol, A. G. (1980). *Engineering Applications of Correlations and Spectral Analysis*. Wiley, New York.

Bergamasco, M., & Lovera, M. (2013). State space model identification: from unstructured to structured models with an $H_\infty$ approach. *IFAC Proceedings Volumes*, **46**, 202–207.

Boukhebouz, B., Mercère, G., Grossard, M., Lamy, X., & Laroche, E. (2020). Identification of single flexible-joint robot dynamics: a nonparametric approach. In *2020 28th mediterranean conference on control and automation* (pp. 357–362).

Boukhebouz, B., Mercère, G., Grossard, M., & Laroche, E. (2021). Shaping multisine excitation for closed-loop identification of a flexible transmission. *IFAC-PapersOnLine*, **54**(7), 643–648, 19th IFAC Symposium on System Identification SYSID 2021.

Boyd, S. (1986). Multitone signals with low crest factor. *IEEE Transactions on Circuits and Systems*, **33**(10), 1018–1022.

Colgate, J. E., & Hogan, N. (1988). Robust control of dynamically interacting systems. *International Journal of Control*, **48**(1), 65–88.

De Luca, A., & Book, W. (2008). Robot with flexible elements (pp. 287–319). Springer.

Enqvist, M. (2005). *Linear Models of Nonlinear Systems*. (Ph.D. dissertation), Sweden: Dept. of Electrical Engineering, Linkoping University.

Fichera, F., & Grossard, M. (2017). On the modeling and identification of stiffness in cable-based mechanical transmissions for robot manipulators. *Mechanism and Machine Theory*, **108**, 176–190.

Friis, T., Tarpö, M., Katsanos, E. I., & Brincker, R. (2021). Best linear approximation of nonlinear and nonstationary systems using operational modal analysis. *Mechanical Systems and Signal Processing*, **152**, Article 107395.

Fuduli, A., Gaudioso, M., & Giallombardo, G. (2004). Minimizing nonconvex nonsmooth functions via cutting planes and proximity control. *SIAM Journal on Optimization*, **14**, 743–756.

Gahinet, P., & Apkarian, P. (2013). Automated tuning of gain-scheduled control systems. In *52nd IEEE conference on decision and control* (pp. 2740–2745).

Garrec, P. (2008). Screw, nut and cable transmission. U.S. Patent US7389974B2.

Garrec, P., P.Friconneau, J., Measson, Y., & Perrot, Y. (2008). ABLE, an innovative transparent exoskeleton for the upper-limb. *2008 IEEE/RSJ International Conference on Intelligent Robots and Systems*, 1483–1488.

Gerchberg, R. (1972). A practical algorithm for the determination of phase from image and diffraction plane pictures. *Optik*, **35**, 237–246.

Glover, K., & Willems, J. (1974). Parametrizations of linear dynamical systems: Canonical forms and identifiability. *IEEE Transactions on Automatic Control*, **19**(6), 640–646.

Goodrich, M., & Schultz, A. (2007). Human-robot interaction: A survey. *Foundations and Trends in Human-Computer Interaction*, **1**, 203–275.

Gosselin, F., & Grossard, M. (2023). Design and development of a dexterous teleoperation setup for nuclear waste remote manipulation. In O. Gusikhin, K. Madani, & H. Nijmeijer (Eds.), *Informatics in control, automation and robotics* (pp. 80–105). Springer International Publishing.

Guillaume, P., Schoukens, J., Pintelon, R., & Kollar, I. (1991). Crest-factor minimization using nonlinear Chebyshev approximation methods. *IEEE Transactions on Instrumentation and Measurement*, **40**(6), 982–989.

Guillaume, P., Verboven, P., Cauberghie, B., Vanlanduit, S., Parloo, E., & De Sitter, G. (2003). Frequency-domain system identification techniques for experimental and operational modal analysis. *IFAC Proceedings Volumes*, **36**(16), 1609–1614.

Gustavsson, I., Ljung, L., & Söderström, T. (1977). Identification of processes in closed loop—identifiability and accuracy aspects. *Automatica*, **13**(1), 59–75.

Hamon, P., Gautier, M., & Garrec, P. (2011). New dry friction model with load- and velocity-dependence and dynamic identification of multi-DOF robots. In *Proceedings of the IEEE international conference on robotics and automation* (pp. 1077–1084).

Van den Hof, P. (1998). Closed-loop issues in system identification. *Annual Reviews in Control*, **22**, 173–186.

Ljung, L. (1986). *System Identification: Theory for the User*. Upper Saddle River, NJ, USA: Prentice-Hall, Inc.

Makarov, M., Grossard, M., Rodriguez-Ayerbe, P., & Dumur, D. (2016). Modeling and preview $H_\infty$ control design for motion control of elastic-joint robots with uncertainties. *IEEE Transactions on Industrial Electronics*, **63**, 1–1.

Mäkilä, P., & Partington, J. (2004). Least-squares LTI approximation of nonlinear systems and quasistationarity analysis. *Automatica*, **40**(7), 1157–1169.

Martin, J., & Grossard, M. (2014). Design of a fully modular and backdrivable dexterous hand. *The International Journal of Robotics Research*, **33**(5), 783–798.

Ojarand, J., & Min, M. (2017). Recent advances in crest factor minimization of multisine. *Elektronika Ir Elektrotechnika*, **23**(2), 59–62.

Overschee, P. V., Moor, B. D., & Favoreel, W. (1997). Numerical algorithms for subspace state space system identification (n4sid). In *In applied and computational control, signals and circuits* (pp. 385–422).

Peumans, D. (2020). BLA-based design and analysis of VCO-based sigma-delta modulators (Ph.D. thesis), Vrije Universiteit Brussel, Belgium.

Pintelon, R., Louarroudi, E., & Lataire, J. (2013). Detection and quantification of the influence of time variation in closed-loop frequency-response-function measurements. *IEEE Transactions on Instrumentation and Measurement*, **62**(4), 853–863.

Pintelon, R., Rolain, Y., Vandersteen, G., & Schoukens, J. (2004). Experimental characterization of operational amplifiers: a system identification approach part II: calibration and measurements. *IEEE Transactions on Instrumentation and Measurement*, **53**(3), 863–876.

Pintelon, R., & Schoukens, J. (2012). *System identification: a frequency domain approach*. John Wiley & Sons.

Pintelon, R., & Schoukens, J. (2013). FRF measurement of nonlinear systems operating in closed loop. *IEEE Transactions on Instrumentation and Measurement*, **62**(5), 1334–1345.

Pintelon, R., Schoukens, M., & Lataire, J. (2020). Best linear approximation of nonlinear continuous-time systems subject to process noise and operating in feedback. *IEEE Transactions on Instrumentation and Measurement*, **69**(10), 8600–8612.

Rantzer, A. (1996). On the Kalman—Yakubovich—Popov lemma. *Systems & Control Letters*, **28**(1), 7–10.

Schoukens, J., Guillaume, P., & Pintelon, R. (1991). Design of multisine excitations. In *Proceedings of the international conference on control 1991: vol.1*, (pp. 638–643).

Schoukens, J., & Ljung, L. (2019). Nonlinear system identification: A user-oriented road map. *IEEE Control Systems Magazine*, **39**(6), 28–99.

Söderström, T., & Stoica, P. (Eds.), (1988). *System Identification*. Upper Saddle River, NJ, USA: Prentice-Hall, Inc.

Thode, H. C. (2002). *Statistics, textbooks and monographs, Testing For Normality*. Taylor & Francis.

Van Den Bos, A. (1987). A new method for synthesis of low-peak-factor signals. *IEEE Transactions on Acoustics, Speech and Signal Processing*, **35**(1), 120–122.

Van der Ouderaa, E., Schoukens, J., & Renneboog, J. (1988). Peak factor minimization of input and output signals of linear systems. *IEEE Transactions on Instrumentation and Measurement*, **37**(2), 207–212.

Vizer, D., & Mercere, G. (2014). $H_\infty$-Based LPV model identification from local experiments with a gap metric-based operating point selection. In *2014 European control conference* (pp. 388–393). IEEE.

Vizer, D., Mercère, G., & Laroche, E. (2015). Gray-box LPV model identification of a 2-dof surgical robotic manipulator by using an $H_\infty$-norm-based local approach. In *1st IFAC Workshop on Linear Parameter Varying Systems* LPVS 2015: vol. 48, (26), (pp. 79–84).
Vizer, D., Mercère, G., Prot, O., & Laroche, E. (2016). $H_\infty$-Norm-based optimization for the identification of gray-box LTI state-space model parameters. *Systems & Control Letters, **92***, 34–41.

Wellstead, P. (1977). Reference signals for closed-loop identification. *International Journal of Control, **26***(**6**), 945–962.

Wyatt, J., Chua, L., Gannett, J., Goknar, I., & Green, D. (1981). Energy concepts in the state-space theory of nonlinear n-ports: Part I-Passivity. *IEEE Transactions on Circuits and Systems, **28***(**1**), 48–61.

Yang, Y., Zhang, F., Tao, K., Sanchez, B., Wen, H., & Teng, Z. (2015). An improved crest factor minimization algorithm to synthesize multisines with arbitrary spectrum. *Physiological Measurement, **36***.