Optimal ARMAX Model Order Identification of Dynamic Systems

I. INTRODUCTION

Operational Transfer Functions (TFs) or Frequency Response Functions (FRFs) of mechanical structures play a vital role in understanding the dynamic characteristics of the systems and in solving general vibration problems during the operational process. They constitute an effective tool aiding the extraction of modal parameters. Estimating the transfer functions of a mechanical system has thus become an important task in many engineering applications. Different representations of transfer functions are crucial in the description and analysis of system properties. In industrial applications, a measurement of the transfer functions defining the structure properties in the frequency domain can be implemented using vibration instrumentations. Different methodologies are proposed in the literature with the aim of estimating operational transfer functions, with the most common applying the Fourier analysis. The Empirical Transfer Function Estimate (ETFE) is a natural nonparametric method that identifies transfer functions by taking the ratios of the Fourier transform of the outputs to those of the inputs [1]. However, this method requires more data points and raw ETFE estimates are generally not accurate enough. With these estimates, the variance does not decrease as the number of data points increases because they contain no information compression feature. Researchers have conducted various experimental studies on structural dynamics under operational conditions. In [2], the FRFs of a flexible joint industrial manipulator with a serial kinematic were identified based on a non-parametric closed loop. However, due to the nonlinearities of the robot, the method faced a challenge in eliminating disturbances in the estimated FRFs. Operational Modal Analysis (OMA) is another approach for identifying the modal properties of the structure using vibration data obtained under operating conditions. Yili Peng et al. [3] identified in-process FRFs based on the OMA and Experimental Modal Analysis (EMA), which uses the natural frequencies and damping ratios to build FRFs under operating conditions. A simulation of a three-degree-of-freedom-mass-spring-damper system and experiments on a machine tool are adopted to verify the proposed method. Similarly, Zaghbani et al. [4] used OMA in the identification of the dynamics of a milling machine under a cutting process work. At the same time, another method was presented in [5] to generate FRFs from identified poles and zeros in the low-frequency domain. Recently, Coppotelli et al. [6] proposed an approach for estimating FRFs from operational data by changing different mass and stiffness distributions. This method also allows evaluating the modal parameters of the structure via operational modal testing. Conversely, Özşahin et al. [7] introduced a new technique to calculate the variation in tool point FRFs under different working conditions by using an inverse analysis of self-excited chatter vibration. In their method, chatter frequencies were experimentally determined and applied to estimate tool point FRFs on 5-axis milling machine via the relation between the measured force inputs and acceleration outputs. However, the tool point FRFs are not well estimated at high spindle speeds due to the presence at those speeds of a low signal-to-noise ratio and the bandwidth limitation of the dynamometer. Another in-process FRF identification approach of the spindle structure was presented in [8]. In that case, the tooltip FRFs were identified under operational conditions based on an inverse solution of critical stability limits. The method is helpful for predicting the stability of the tool holders when the direct measurement of the tool point FRFs is uncompromised. Parametric estimation methods constitute another system identification class. In these methods, it is suggested to use time series modeling for the mathematical description of the transfer functions. It combines the advantage and information obtained from both measurements and theoretical modeling. Depending on the availability of the measurement signals, the Auto Regressive model (AR) [9, 10] or the Auto Regressive Moving Average model (ARMA) [11] can be used if only the output is available. In [12,

13], a modal analysis was conducted in different industrial structures based on three Auto-Regressive Moving Average methods, namely, the recursive least-square, output error, and corrected covariance matrix methods to determine the optimal model order. Conversely, in the case of measurable or identifiable excitation forces, the Auto-Regressive eXogenous (ARX) model [14] can exploit, by assuming that the model's errors and disturbances are white noise. However, because of the unavoidable noise contaminated in the measured signals, the quality of estimated FRFs can be adversely affected by noise originating from the test environment. When the system operates in an industrial condition with a lot of disturbance, identifying the transfer functions of a complex structure may become difficult.

In this paper, we present an original method designed for automatically extracting the modal parameters from identified transfer functions based on the concept of the optimal ARMAX model. Particular attention is paid to selecting optimal model orders, which can closely reflect the dynamic system. The work contributes to the determination of a model order based on the estimated transfer functions, by using the framework of the ARMAX model. The proposed method is experimentally applied to a robot during its grinding operation, and the results are compared to those of the original ARX model. The measured grinding forces may be considered the exogenous inputs excitation, and the disturbances of the system are taken into account by adding the Moving Average part into the model. The estimated orders are verified based on the most common selection criteria, such as the Akaike Information Criteria (AIC) [15], the Bayesian Information Criteria (BIC) [16], and the Noise Order Factor (NOF) [17]. In this study, the ARMAX model is expressed in a convenient way for computation at the low orders, which gives a more parsimonious representation and helps improve the modeling performance, with less computational complexity. We have organised the rest of this paper in the following way. The motivation for the research is established in Section 2 through a detailed description of the time series modeling, with a focus on the ARX and ARMAX models. Section 3 proposes an original method to determine an optimal model order of the mechanical system. Experiments are then conducted on the flexible manipulator, SCOMPI, under grinding operation to validate the proposed methods in Section 4, followed by the identification procedure and the results. We end by drawing several conclusions from this research.

II. TIME SERIES MODELING

System identification is the art of modeling a dynamical system from raw time series data. We consider the problem of estimating a dynamic system model based on the measurement of an N points input-output data, which will be pre-classified into input $\mathbf{u}(t)\in \mathbb{R},t=1,...,T$ and output $y(t)\in \mathbb{R},t=1,...,T:$

Various representations of linear time series such as Auto-Regressive (AR), Auto-Regressive Moving Average (ARMA), Auto-Regressive eXogenous (ARX), and Auto-Regressive Moving Average eXogenous (ARMAX) can be employed to extract dynamic parameters $[18]$. Since there are various time series data types, we should choose an appropriate model. In general, such models are based on an Auto-Regressive (AR) part or output, an eXogenous (X) part or input, or a Moving Average (MA) part or error term, depending on the situation. The AR model is the simplest time series representation, which linearly depends on output data (the vibration responses). In the availability of both input (the measurable and known excitation force) and output data, the ARX model is usable. It is possible to combine these models with the MA term and produce the ARMA representation for the output-only cases and the ARMAX model for the input-output conditions. Once the most appropriate modal structure is selected, we can apply a model to the measurement by minimizing certain criteria:

where $\theta$ is the unknown parameter vector of the parametric model structure.

In automatic control applications, given the current state and input signal, the model can be applied to predict the output of the system by choosing a cost function in the form.

where $\mathbf{L}(q)$ represents a filter that removes unwanted properties in the measurement data, and $l(.)$ is a convex function.

A common representation of the Linear Time-Invariant (LTI) system can be expressed in the from of the linear transfer function model:

where q is the forward shift operator, that is, $q^{-k}.\mathbf{y}(t)=\mathbf{y}(t-k)$. Here, $\mathbf{y}(t)$ is a $n_y$ dimensional vector of output, $\mathbf{u}(t)$ is a $n_u$ dimensional vector of input, and $\mathbf{w}(t)$ is the disturbance sequence with an appropriate dimension and assumed to be an independent and identically distributed stochastic process, respectively. Furthermore, the transfer functions $\mathbf{G}(q,\theta )$ and $\mathbf{H}(q,\theta )$ are rational functions in the backward shift operator q, and the coefficients are given by the elements of the parameter vector $\theta$.

The predictor associated with the output is given by [18]

This model structure is quite general, but we can develop some special cases. A simple case is the ARX model structure, which is:

that can be rewritten in a more general polynomial form as:

where

are autoregressive and exogenous matrix parameters, with I denoting the identity matrix. $n_a$, $n_b$, and $n_k$ are the orders of the ARX model, $n_a$ is equal to the number of poles and $n_b$ is the number of zeros, while $n_k$ is the pure time delay in the system. Since $\mathbf{G}(q,\theta )=\mathbf{B}(q)/\mathbf{A}(q)$ and $\mathbf{H}(q,\theta )=1/\mathbf{A}(q)$, the predictor of the output can be written as:

where

When the noise is assumed to be a white Gaussian process with zero means, which is uncorrelated with the regressors, the model parameters $\theta$ are estimated via the Least-Square (LS) estimator [19]:

The variance $\mathbf{y}(t)-\mathit{\phi }(t{)}^T\mathit{\theta }$ represents the remaining un-modelled behavior of the data.

The corresponding transfer function is $\mathbf{G}(q,\theta )=\mathbf{B}(q)/\mathbf{A}(q)$. (15-a)

However, in this case, only the deterministic part of equation (1) is estimated by considering no noise $\mathbf{H}(q,\theta )$. If a noise model is sought, additional steps are needed, assuming that the noise is described by a Moving Average process $\mathbf{C}(q)$, which results in an ARMAX structure:

and its polynomial form

where

is the polynomial of order $n_c$ which represents for Moving Average term.

(18) In contrast with the simpler ARX model, this presentation form offers a noisy transfer function $\mathbf{H}(q,\theta )=\mathbf{C}(q)/\mathbf{A}(q)$, which allows representing different types of noise characteristics through a proper choice of the MA polynomial term. In engineering applications, it is unavoidable when environmental noise contaminated in measured data. Therefore, a parametric system identification algorithm should be adopted to identify the modal parameters from the noisy data. Under this condition, the ARMAX model with real-time identification proves to be efficient. Figure 1 presents the structures of both ARX and ARMAX models.

If only the ARX model is use, the noise model is described as $\mathbf{H}(q)=1/\mathbf{A}(q)$, where $\mathbf{A}(q)$ is also used as the denominator of $\mathbf{G}(q)$ describing the dynamic model. This model implies that the polynomial form must be an average estimate of the poles of $\mathbf{G}(q)$ and $\mathbf{H}(q)$. However, the noise model can be better estimated using the ARMAX model, as described by the numerator polynomial term $\mathbf{C}(q)$. When parametric methods use for model estimation, the transfer functions to be estimated are defined as a function of a parameter vector ${\theta }_N\triangleq (a_1\dots a_{na}{b}_1\dots b_{nb}{c}_1\dots c_{nc}{)}^T$. For the identification of the value ${\hat{\mathit{\theta }}}_N$, $\overline{\mathbf{G}}(q,{\overline{\theta }}_N)$ and $\hat{\mathbf{H}}(q,{\hat{\theta }}_N)$ are closest to $\mathbf{G}(q,{\theta }_N)$, and $\mathbf{H}(q,{\theta }_N)$. In [20], the authors showed that under reasonable conditions,

where and $\mathit{\epsilon }(t,\theta )$ is prediction error and defined in equation (4). According to [21, 22], with this definition of ${\theta }_N^{\ast }$, it is possible to split the total estimation error between the true frequency response function $\mathbf{G}(q)$ and the estimated one $\mathbf{G}(q,\widehat{{\theta }_N})$ into two parts as follows:

III. A MODEL ORDER DETERMINATION APPROACH

(21) As we can see in equation (21), the errors in the estimated FRF have two components. The first contribution $[\mathbf{G}(q)-\mathbf{G}(q,{\theta }_N^{\ast })]$ is the bias contribution, a deterministic quantity due to the modeling error. If the selected model has a lower or a higher complexity than the true system, this bias error will be present at some frequencies. Choosing the model order should be flexible enough to allow a good fit to the measurement data but adequately constrained so that noise does not invoke unsuitable models. Consequently, selecting an optimal model order respecting this compromise is an important issue in system identification. This concern is analyzed in the following section to select the optimal model orders of an ARMAX model. The second contribution $[\mathbf{G}(q,{\theta }_N^{\ast })-\mathbf{G}(q,{\theta }_N^{-})]$ represents the noise or variance errors, which are due to noise in the measured input and output data. It is a random variable that will disappear if there is no noise or if the number of data tends to infinity.

When time series modeling ARX or ARMAX models employed, the performances may be affected by selecting the model order. Choosing a sufficient and correct model order has always been a challenging issue. Once the model orders are properly selected, the models successfully represent the underlying phenomenon with the lowest complexity. This method aims not only to find a model capable of describing a specific set of data but is also helpful for the validation of the inference procedures. Consequently, there is a need to develop a reliable method to identify the orders of AR, MA, and eXogenous polynomials. In general, most of the mechanical structures are operated in the low-frequency range with limited bandwidth. Having a model with orders that are too high may lead to an overfitting problem as it includes too much irrelevant oscillation information and generates high computational costs, and a model with orders that are too low will not be solid enough to capture the underlying physical system dynamics. For its part, a model with an appropriate order can precisely describe the dynamic characteristics of the system. Because of its important role in system identification, model order determination has attracted much attention in the literature, with researchers proposing different criteria for order determination [23]. The final prediction error (FPE) criterion was originally proposed by Akaike [24] for determining the AR order and was extended to the ARMA model by [25]. After adding the inflating effects of estimated coefficients, the optimal order was chosen by minimizin the one-ste-ahead mean square forecast error. A method based on the eigenvalues of a modified covariance matrix, which is robust to noise levels, was proposed in [12, 13] for determining the model order. Some other concept as information theory like Akaike's Information Criterion (AIC) [15], Bayesian Information Criterion (BIC) [16] or Minimum Description Length (MDL) [26] were developed to produce an estimate model order. Among these techniques, the AIC is a heuristic approach, which has attracted much attention in the literature. His technique penalizes the likelihood of the number of parameters in the model by attempting to choose the most suitable model order.

Considering an N-dimensional time-series data, the AIC is given by equation (22):

where N denotes the number of data points, $(z)$ is a dimension associated with the vector of unknown parameters to be estimated and

$\hat{\mathrm{\Sigma }}$ is the covariance matrix of the innovation sequence associated with the estimated coefficients, the w[t] is innovation square, or the model error.

When the AIC value is at a minimum, we obtained the most suitable order. A minimum AIC is theoretically situated at a sufficient value of $(z)$ that best represents the dimension of the unknown parameters. For an ARMAX $(n_a,n_b,n_c,n_k)(p)$ model, $(z)$ would typically be equal to $(n_a+n_b+n_c+n_k)$, with $n_a,n_b,n_c$ and $n_k$ orders for its AR, MA, eXogenous components and time delay, respectively, while p represents the number of orthonormal functions by which each of these components is multiplied. Here, it should be noted that although there would be many possible combinations of $n_a,n_b,n_c$ and $n_k$ that can produce the same adequate value of $(z)$, only the right combination would yield the smallest AIC value. The z value may be defined as:

The AIC corresponding to an ARMAX $(n_a,n_b,n_c,n_k)(p)$ model is written here as:

where z is the number of scalar parameters in the ARMAX model.

However, in many cases, AIC does not give an optimal order. [27] showed empirical evidence that AIC tends to pick models which are over-parameterized. The BIC overcomes this shortcoming by including an additional term that

(25) In this paper, it is written as:

penalizes the model complexity and enhances the procedure, which are based on the same concepts governing the AIC but are better suited for large data sequences. The BIC criterion has the following general form:

These criteria rely on the evolution of the error covariance, which monotonically diminishes concerning the model order. It asymptotically chooses the correct order model if the underlying multiple time series has high dimensions but tends to overestimate the model order as the data length increases. Thus, selected model orders can be greater than the optimal model orders.

However, attention has recently shifted to the equally important problem of bias resulting from

(27) under-modeling. Wahlberg and Ljung [28] have conducted excessive research on the distribution of bias and variance in the estimated transfer function.

In this paper, we present an improved approach to determining time series model orders based on means square errors of the estimated transfer function. This approach is different from traditional criteria such as AIC and BIC. We transformed an averaged frequency means square error on the estimated transfer function into a mean square output prediction error criterion.

The proposed method allows extracting the modal residuals with a sufficient order that guarantees the extraction of uncorrelated residual samples and the avoidance of an overfitting problem. The selection of the optimal order is based on a minimal variance of the total mean square error (1, 2)

\varepsilon\left{\left|$\mathbf{G}$\left(q\right)-$\mathbf{G}$\left(q,${\hat{\theta }}_N$\right)\right|^{2}

$\mathbf{P}$ \triangleq $\frac{1}{2\pi }$ \int_{-\pi}^{\pi} \hat{\mathbf{E}}{p{o}}(\omega) ${\mathbf{D}}_u$(\omega) d\omega

The optimal order obtained when:

The aim of this technique is converting the bias error into a random variable by ascribing a prior distribution to it. Consequently, we obtain an estimate for the average characteristics of the total errors. We assume that the transfer function represented by ${\mathbf{G}}_T(q)$ is a stochastic process indexed by the variable $\omega$, and given some value of ${\theta }_0$, it can be decomposed as the sum of a parameterized ${\mathbf{G}}_T(q,{\theta }_0)$ plus the residual ${\mathbf{G}}_{\mathrm{\Delta }}(q)$:

where ${\mathbf{G}}_{\mathrm{\Delta }}(q)$ is a zero-mean stochastic process:

Each system will provide one realization and analogous to the embedding of the single noise realization in a stochastic process. In the framework of an ARMAX model, for ease of implementation, we restrict our consideration to the linear model structures of the Single Input Single Output (SISO) case. The question of computing asymptotic variance expressions in the presence of under-modeling for transfer function prediction errors in frequency domain identification is addressed in this paper. The discrete linear transfer function of the ARMAX model can be expressed in the following form:

The true transfer function can be expressed by the relationship between the measured input and output of the system:

The ${\mathbf{E}}_{p_o}(\omega )$ means square error between the true and estimated transfer function in (28) is calculated by the following equations:

This index is an effective criterion since it includes the stochastic participation in the denominator by adding an input power spectral density ${\mathbf{D}}_u(\omega )$.

We converted a average frequency means square error criterion on the estimated transfer functions into a mean square output prediction error criterion. Due to the presence of measurement errors, the proposed method is adopted to determine the optimal model orders and proved its robustness to parameter uncertainties. This modification bridged the connection between classical variable selection criteria such as AIC and BIC and the non-concave penalized likelihood methodology that allows greater flexibility in choosing the desired models. The goodness of its performance will be assessed in the next section through a comparison with traditional validation methods. The optimal model will be tested by the quality of the residuals, the histogram, and autocorrelation functions.

IV. INDUSTRIAL APPLICATION ON A SCOMPI ROBOT

4.1 Description of the test structure

Flexible manipulators are employed in the maintenance of large hydro electronic equipment as they represent an effective solution for repair jobs $[29]$. Despite their attractive properties, controlling lightweight robot manipulators is still a challenging task. Their flexibility induces structure vibration that deteriorates the trajectory tracking accuracy and may lead to instability issues. In this part, the present method is adopted to identify the dynamic behavior of a light, portable, track-based multi-process manipulator named SCOMPI (Super-COMPact Ireq) $[30,31]$. The dynamical transfer functions and modal parameters must be monitored during the grinding process to control vibrations. Figure 2 presents the working envelop and the construction of the robot, with its links and joints.

4.2 Experimental procedure

The tool contact point Frequency Response Functions (FRFs) and their modal parameters represent the key for solving the dynamics analysis. In this paper, the actual grinding forces are used as the input excitations in estimating the FRFs. The grinding test is performed at the 3225 rpm rotating speed with an average 0.08 mm axial depth of grinding cut. The grinding work is realized on a hard steel EN31-64HNC workpiece with the dimensions of 150 x 7 x 48 mm. During the cutting operation, cutting forces are measured with a type CH8408 3-axis Kistler dynamometer, which is directly attached under the workpiece to record the normal force direction ( $F_x$ ), the tangential force direction ( $F_y$ ) and the axial force direction ( $F_z$ ). At the same time, three

PCB-352C34 piezoelectric sensors with a sensitivity of 5.29 mV/G are used to measure accelerations at the robot's end-effector ( $S_1-S_3$ ) to capture the vibration in three directions. The measurement is conducted through the LMS data acquisition system for 10 s at a 512 Hz sampling frequency. Figure 3 shows the actual experimental setup of the SCOMPI robot under a real grinding operation. An LMS test lab system was used for acquiring the data in real-time during experiments, as well as for calculating averaged FRFs, which were used to validate the estimated results of the proposed approach.

The structure is excited through a grinding operation, during which all the effects related to the rotating tool or the grinding process are considered. FRFs are obtained during the cutting operation through the relations of the cutting forces and the vibration responses. The measured cutting forces are taken as the excitation sources of the system. In this paper, the operational FRFs are determined from the relation between the cutting forces and responses through the time series ARMAX modeling. The detailed experimental description is provided in Table 1.

{"image_source":{"path":"images/1224bcb25343ba3eb33066d8cc1267e3c9d66f1bd6592013998e754100a04a92.jpg"},"content":"","chart_caption":[{"type":"text","content":"(b)"},{"type":"text","content":"(c)"},{"type":"text","content":"Fig. 4: Measured data on SCOMPI robot"}],"chart_footnote":[]}

V. IDENTIFICATION PROCEDURE

5.1 Model orders estimation

The parametric identification of the structural dynamics is based on force and response signals with a 10s sample length. The modeling strategy consists of the successive fitting of the ARMAX $(n_a,n_b,n_c,n_k)$ model. First, the model orders are selected based on the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). Figure 5 plots result from the AIC and BIC techniques respectively obtained by directly fitting the ARMAX model of increasing orders p = 1 - 60 to the different sets of experimental data.

{"image_source":{"path":"images/460269222dd2a4edbe38ec29248f569f242563e5184a8abdf3c069223a393eae.jpg"},"content":"","chart_caption":[],"chart_footnote":[]} {"image_source":{"path":"images/0c445fd2d4556a6e934b50888d4379f08b80693950021b0cd705687bc539cc11.jpg"},"content":"","chart_caption":[],"chart_footnote":[]}

{"image_source":{"path":"images/d4cdb76f501e83ea0b2245cb0151e381ed2351c01a9c4ac2feed31a8453a1b11.jpg"},"content":"","chart_caption":[{"type":"text","content":"Fig. 5: Model order selection based on AIC and BIC criteria for different sets of data: (a) X data, (b) Y data, (c) Z data"},{"type":"text","content":"Optimal ARMAX Model Order Identification of Dynamic Systems"}],"chart_footnote":[]} By dividing the experiment into different test directions based on the operation of the robot, we can decompose and characterize the dynamical behavior of the system for each direction. From Figure 5, both the AIC and BIC methods decrease with the model order, and a minimum may be assumed close to 47. However, the main limitation of using such techniques is that they may suggest different model orders and not determine the optimal orders. Thus, these values must be judged carefully.

{"image_source":{"path":"images/3a4f7d042378c867ed1d31c9a620da00d7f4e22da20e5fdfa9b523a4792152c2.jpg"},"content":"","chart_caption":[{"type":"text","content":"(b)"}],"chart_footnote":[]} By applying the proposed method in selecting an optimal model order, the orders for which all curves lead to a stable value are identified. Figure 6 illustrates an optimal order, defined as the smallest order value. The point of convergence starts at around order of 10 in all three directions.

{"image_source":{"path":"images/77fa5dd35202760338168c30fcce22379a61fca99c6264460c53165e62c03298.jpg"},"content":"","chart_caption":[],"chart_footnote":[]} {"image_source":{"path":"images/f77e01b23f9be33ef9a890b1503c9370fe20ddd4930b37dbb6868a55ccdcfeb5.jpg"},"content":"","chart_caption":[{"type":"text","content":"(c)"},{"type":"text","content":"Fig. 6: Model orders selection based on the proposed approach for different sets of data: (a) X data, (b) Y data, (c) Z data"}],"chart_footnote":[]}

To compare the efficiency of the present method, Vu et al. [17] proposed a technique for determining an efficient model order $p_{eff}$ based on the analysis of the Noise-to-Signal Ratio (NSR).

A Noise-ratio Order Factor (NOF) is calculated as a variation of the NSR between two successive orders:

The estimated $NSR$ is given in equation (37) based on the trace norm part of the error covariance matrices $\hat{M}$ and the estimated deterministic $\hat{K}$.

The NOF is a representative factor for the convergence of the NSR, which changes dramatically at low orders and converges at high orders. Since this criterion is positive and close to zero, the convergence may be assumed close to 10, considered as an efficient model order (Figure 7).

{"image_source":{"path":"images/3fe7f94520f2c4be312649e1a7a444311b128ff504476442dafc6b6a06ec057f.jpg"},"content":"","chart_caption":[{"type":"text","content":"(c)"}],"chart_footnote":[]} {"image_source":{"path":"images/64ac9d6f22660171cb2ade24582f557e7cf60db44cdecca403ca8397d77ad81e.jpg"},"content":"","chart_caption":[{"type":"text","content":"(b)"}],"chart_footnote":[]}

{"image_source":{"path":"images/8679ef591171c2376bbdac4e3cf6371ba008f52c7f6584c06b14244788ccc9da.jpg"},"content":"","chart_caption":[{"type":"text","content":"Fig. 7: NOF evolution and efficient model orders selection for different sets of data: (a) X data, (b) Y data, (c) Z data"}],"chart_footnote":[]}

Theoretically, the modeling of a complex structure like the SCOMPI robot should result in a high-order model. Based on the AIC and BIC criteria, a 47-order model should be selected as a suitable choice. However, a lower order could be chosen with all essential characteristics of the real system preserved. By comparison, the Noise-to-Signal Ratio and the proposed method are selected with a model order of around 10. However, in the case of the complex ARMAX model, it is characterized by three different orders, the model order estimation is not straightforward. In experimental modal analysis, the orders $n_a$, $n_b$, and $n_c$ depend on the model parameterization. According to [32, 33], the choice of $n_b$ is a function of the type of response measurements used and the inter-sample behavior of the data. Moore et al. [34] suggest an ARMAX ( $p,p,p$ ) model in which $p=n_a=n_b=n_c$, for the case of a vibration acceleration measurement, or ARMAX( $p$, $p-1$, $p$ ), for the case of a vibration displacement or velocity measurement with an appropriate time delay. As can be noted, since structural vibrations are usually measured in terms of acceleration rather than displacement and velocity in actual experiments, we would choose $p_{optimal}=n_a=n_b=10$.

The Moving Average order is initially set equal to the Auto Regressive part since the resulting noise model has the flexibility of representing several stochastic processes, including white noise. There is some experimental evidence in the structural systems $[34,35]$, which indicates that for low noise levels, the required MA order is often smaller or equal to the AR one. The order $n_c$ of the MA matrix and time delay order $n_k$ are dependent on the noise present in the system, and generally, no information on the nature of this disturbance is available. Therefore, the extracted required value of the MA order, as well as an appropriate time delay order, will be carefully examined in this study. This can be done by initially setting $n_a=n_b=10$, and then selecting the best model by using a proposed criterion to test the effectiveness of changing $n_c$ and $n_k$ from the set of $\{1,2,...,20\}$, resulting in the estimated results presented in Figure 8.

{"image_source":{"path":"images/f9c2f779aa430f7b57483d451e169dc83cca446951502b9cd16d9f543b2a416e.jpg"},"content":"","chart_caption":[],"chart_footnote":[]} {"image_source":{"path":"images/71591a95f2eb05e0f3851abb104d1b384f9b9b9a9e2065cdd73173320dac875d.jpg"},"content":"","chart_caption":[{"type":"text","content":"Fig. 8: Moving average order (a) and time delay order (b) selection"}],"chart_footnote":[]}

Based on these results, the set of orders $(10,10,8,2)$ is chosen as consisting of the smallest order values, which can be applied to fit the data with a negligible discrepancy and can be effectively utilized for modal analysis. However, a further step needs to be validated to assess the adequacy of the estimated model. Several diagnostic checks can be used to decide whether the ARMAX model is adequate based on the residuals, which characterized by an uncorrelated sequence. Figures 9 and 10 display the histogram, {"image_source":{"path":"images/c3c75fef25197caa08e84a44ce8801bf01669f3cd7081301c128cef2b165221f.jpg"},"content":"","chart_caption":[],"chart_footnote":[]}

Quantile-Quantile (QQ) plot, Autocorrelation Function (ACF), and Partial Autocorrelation Function (PACF) of the residuals of the ARMAX model. The histogram is unit-modal and symmetric around zero. From the QQ-plot, the residual approximately fit a straight line, and this can be assumed normal. The values of both ACF and PACF are located roughly between the upper and lower bounds of a confidence interval. Therefore, the residual can be assumed independent, and the selected orders can be used.

{"image_source":{"path":"images/799d03df44f245b5a081b4629912fb58e57d56e8d0122505cd71c5bcc9c990af.jpg"},"content":"","chart_caption":[],"chart_footnote":[]} {"image_source":{"path":"images/b7cf49efb81781d8264d4ec97bf9117f02860c1330cdd6738b15c4f12683e4c8.jpg"},"content":"","chart_caption":[{"type":"text","content":"Fig. 9: Residual errors histogram and normal probability plot of an estimated ARMAX (10,10,8,2) model"}],"chart_footnote":[]}

{"image_source":{"path":"images/d72b8e4de1cb6a7fd5a951594b538f2de2ab0457c191ad5aed4a53470bc17ddf.jpg"},"content":"","chart_caption":[{"type":"text","content":"Fig. 10: Autocorrelation function (ACF) and Partial Autocorrelation Function (PACF) of the residual ARMAX (10,10,8,2)"}],"chart_footnote":[{"type":"text","content":"Optimal ARMAX Model Order Identification of Dynamic Systems"}]}

The above identification procedure leads to an optimal ARMAX (10,10,8,2) model, which is selected as an adequate model for structural analysis, and for the extraction of modal parameters. The next section evaluates the performance and effectiveness of this dominant reduced model, in which all the essential characteristics of the real system can be retained.

5.2 Frequency Response Functions identification

Based on the discussion above, once the excitation and vibration measurement data have been selected with the appropriate orders, an ARMAX model needs to be estimated within the model structure. In this section, the two parametric transfer function estimators for assessing the dynamic flexible manipulator are analyzed. The FRFs estimation is performed using two parametric models, ARX ( $n_a$, $n_b$, $n_k$ ) and ARMAX ( $n_a$, $n_b$, $n_c$, $n_k$ ), based on data records from the structure. A scalar Single Input – Single Output (SISO), ARX and ARMAX models are used at a time. The input signals to each model are the force signals, and the output signals are measured accelerations, respectively. The assessments of an ARX model, and an ARMAX model, are undertaken based on an experimental SCOMPI structure during the grinding operation. Figures 11-19 demonstrate the estimated Frequency Response Functions obtained by the ARX and ARMAX models, which are compared to those measured by the LMS system.