Identification of Linear Time-Invariant Systems: A Least Squares of Orthogonal Distances Approach

Abstract

This work describes the parameter identification of servo systems using the least squares of orthogonal distances method. The parameter identification problem was reconsidered as data fitting to a plane, which in turn corresponds to a nonlinear minimization problem. Three models of a servo system, having one, two, and three parameters, were experimentally identified using both the classic least squares and the least squares of orthogonal distances. The models with two and three parameters were identified through numerical routines. The servo system model with a single parameter only considered the input gain. In this particular case, the analytical conditions for finding the critical points and for determining the existence of a minimum were presented, and the estimate of the input gain was obtained by solving a simple quadratic equation whose coefficients depended on measured data. The results showed that as opposed to the least squares method, the least squares of orthogonal distances method experimentally produced consistent estimates without regard for the classic persistency-of-excitation condition. Moreover, the parameter estimates of the least squares of orthogonal distances method produced the best tracking performance when they were used to compute a trajectory-tracking controller.

1. Introduction

The aims of a mathematical model are to make forecasts and deductions about the behavior of the system that it is attempting to describe, and then, it must accurately reproduce its measured variables. However, since the mathematical model is developed based on empirical knowledge and physical laws, most of the parameters used to build the model may be unknown. Consequently, it is necessary to determine the values of the parameters to produce the best fit between the signals produced by the model and the measured data. This procedure is called parameter identification or parameter estimation, and it corresponds to a linear estimation if a model is a linear combination of linearly independent functions; otherwise, the estimation is nonlinear in nature [1,2,3,4,5,6].

A widely used approach to solve parameter identification problems and curve fitting is the least squares (LS) method. It was attributed to Gauss, but the first published work on this subject was by Legendre at the beginning of the 19th century [7]. The seminal papers of R. J. Adcock in 1877 [8] and 1878 [9] proposed a variant of the LS method, namely the least squares of the orthogonal distances (LSOD).

In order to appreciate the differences between the LS and the LSOD algorithms, consider Figure 1, where it is desired to fit a set of experimental data to a straight line. The standard LS algorithm minimized the sum of the squares of the vertical distances between the experimental data and the straight line, and consequently, it only considered errors in the y variable. In contrast, the LSOD algorithm minimizes the sum of the squares of the orthogonal distances between the experimental data and the straight line; hence, it accounted for errors in both the x and y variables.

The LSOD method has been mainly used for fitting data to geometric objects, including lines, planes, spheres, ellipses, and hyperbolas [7,10,11,12]. It has also been applied in several problems, including parameter identification of tire models [13] and the modeling of sorption data [14]. It has been reported that LSOD method outperformed the classic LS method when they were used for model-fitting of experimental data [15]. However, parameter identification of dynamical systems is a classic subject, as it has been covered in many references and textbooks, and the LS algorithm has probably been the most widely used method to this end [2,3,4,5,6,16,17,18,19].

A practical case of interest has been the parameter identification in servo systems, which are used in many applications, including motion control, robotics, solar-tracking systems, and prosthetics. However, advanced control algorithms for servo systems require previous knowledge of their parameters. This makes it necessary to use the most suitable method for parameter identification. The parameter identification of servo systems may be solved through different techniques, for example, employing the classical least squares (LS) method [16,20,21,22,23]. Other techniques include Kalman filtering [24], gradient algorithms [25], particle swarm optimization (PSO) [26], total least squares (TLS), and the LSOD method [27].

This work proposes a parameter identification procedure for servo systems using the LSOD method, which recasts the parameter identification problem as data fitting to a plane. Then, the parameter identification could be accomplished by solving a nonlinear minimization problem for which the cost function and the conditions for finding a critical point were described. Then, finding the critical points was equivalent to obtaining the parameter estimates. Three models of the servo systems, having one, two, and three parameters, were considered for parameter identification. The two and three parameter models were solved using numerical routines, whereas for the model with one parameter, which corresponded to the input gain, the parameter identification problem was equivalent to the data fitting to a line, and analytical conditions for obtaining a critical point and to determine if it corresponded to a minimum were derived. Moreover, the estimate of the input gain was obtained by solving a simple quadratic equation whose coefficients depended on measured data. To the best of the authors’ knowledge, this last result had no equivalent when using the LS algorithm.

The key points of this work were the following:

• We presented a parameter identification scheme using the LSOD method for servo systems.

• It was easily extended for the parameter identification of n-order LTI systems subjected to constant disturbances, assuming that only input and output data were available from measurements.

• A cost function depending on the orthogonal distances between the data and a plane was derived, where the latter was built from a linear parametrization.

• Conditions for obtaining a critical point, i.e., the parameter estimates, as well as the numerical routine to solve the nonlinear minimization problem were described.

• For the model of a servo system containing only an unknown parameter, the input gain, the simple analytical conditions for computing the critical point and to establish if the point produced a minimum of the cost function were derived. The LSOD method provides the possibility of using smooth excitation signals with low harmonic content, which are more suitable for applications where aggressive excitation signals cannot be used.

• The proposed parameter identification methodology was experimentally compared with the LS method. Three servo system models were estimated using excitation signals with different frequency spectra, and the parameter estimates obtained with both methods were used to design a trajectory-tracking controller.

The work is structured as follows: Section 2 describes the linear regression for perturbed LTI systems, assuming that only input and output data are available. The LS method is briefly described in Section 3. Section 4 is devoted to the LSOD method. The three models of a servo system are presented in Section 5. The parameter estimates obtained through the LS and LSOD methods using experimental data and the design and testing of three trajectory-tracking controllers are presented, respectively, in Section 6 and Section 7. The main conclusions of this work are described in Section 8.

2. Linear Regression for Perturbed LTI Systems

The general problem of parameter identification was defined as follows. Given the structure of a model and a set of measured data points, the problem consisted of estimating the unknown model parameters, so the output of the model computed using the estimates would match the data in an optimal manner [2].

Let the mathematical model of a linear-time invariant (LTI) system be as follows:

(1)$\begin{array}{ccc}\hfill A\left(D\right)y& =& B\left(D\right)u+d;D^i=\frac{d^i}{d{t}^i}\hfill \\ \hfill A\left(D\right)& =& D^n+{\alpha }_1{D}^{n-1}+\cdots +{\alpha }_n;n\ge 1\hfill \\ \hfill B\left(D\right)& =& {\beta }_0{D}^m+{\beta }_1{D}^{m-1}+\cdots +{\beta }_{m}n\ge m\ge 0\hfill \end{array}$

where it is subjected to a constant disturbance d, with input u, and output y. The aim was to estimate the constant unknown parameters ${\alpha }_1,\dots ,{\alpha }_n$, ${\beta }_0,{\beta }_1,\dots ,{\beta }_m$ and the disturbance d. A linear regression for parameter identification purposes was obtained by rewriting (1) as follows:

(2)$\begin{array}{ccc}\hfill y^{\left(n\right)}& =& -{\alpha }_1{y}^{(n-1)}-\cdots -{\alpha }_{n}y+{\beta }_0{u}^{\left(m\right)}+{\beta }_1{u}^{(m-1)}+\dots +{\beta }_{m}u+d,\hfill \end{array}$

(3)$\begin{array}{ccc}& =& {\theta }^{\top }\mathrm{\Phi }, \hfill \end{array}$

where

$\begin{array}{ccc}\hfill \theta & =& {\left[{\alpha }_1,\dots ,{\alpha }_n,{\beta }_0,{\beta }_1,\dots ,{\beta }_m,d\right]}^{\top },\hfill \\ \hfill \mathrm{\Phi }& =& {\left[-y^{(n-1)},\dots ,-y,u^{\left(m\right)},u^{(m-1)},\dots ,u,1\right]}^{\top }.\hfill \end{array}$

Note that the regression vector $\mathrm{\Phi }$ depends on the system input and output and their time derivatives. In practice, these derivatives are seldom available, thus precluding the use of (2).

If only the input u and the output y are available from measurements, it is current practice to use linear filters to obtain a linear regression for the unknown parameters but with a regression vector, depending only on input–output measurements [17,18]. To this end, we defined the transfer function of the following asymptotically stable filter:

(4)$F\left(s\right)=\frac{f_n}{s^n+f_1{s}^{n-1}+\dots +f_n},$

and the filtered variables $u_f$ and $y_f$,

(5)$\begin{array}{ccc}\hfill \mathcal{L}\left\{u_f\right\}& =& F\left(s\right)\mathcal{L}\left\{u\right\},\hfill \end{array}$

(6)$\begin{array}{ccc}\hfill \mathcal{L}\left\{y_f\right\}& =& F\left(s\right)\mathcal{L}\left\{y\right\},\hfill \end{array}$

where $\mathcal{L}\{·\}$ stands for the Laplace transform and s is a complex variable.

The expressions (5) and (6) found the next state space representations:

(7)$\begin{array}{ccc}\hfill {\dot{z}}_1& =& A_f{z}_1+b_{f}u,\hfill \\ \hfill u_f& =& c_f^{\top }{z}_1,\hfill \end{array}$

(8)$\begin{array}{ccc}\hfill {\dot{z}}_2& =& A_f{z}_2+b_{f}y,\hfill \\ \hfill y_f& =& c_f^{\top }{z}_2,\hfill \end{array}$

where

(9)$\begin{array}{ccc}\hfill A_f& =& \left[\begin{array}{ccccc}-f_1& 1& 0& \cdots & 0\\ ⋮& & \ddots & & \\ -f_{n-1}& 0& 0& \cdots & 1\\ -f_n& 0& 0& \cdots & 0\end{array}\right];b_f=\left[\begin{array}{c}0\\ ⋮\\ 0\\ f_n\end{array}\right];c_f=\left[\begin{array}{c}1\\ ⋮\\ 0\\ 0\end{array}\right],\hfill \\ \hfill F\left(s\right)& =& c_f^{\top }{(sI-A_f)}^{-1}{b}_f,\hfill \end{array}$

with $f_i>0,i=1,\dots ,n$ and $det\left(A_f\right)={(-1)}^{n+2}{f}_n$. Hence, $A_f$ is invertible and Hurwitz-stable.

However, define

(10)$z_3=A\left(D\right)z_2-B\left(D\right)z_1,$

taking the time-derivative of (10) and substituting (7) and (8) produce

(11)$\begin{array}{ccc}\hfill {\dot{z}}_3& =& A\left(D\right){\dot{z}}_2-B\left(D\right){\dot{z}}_1,\hfill \\ & =& A\left(D\right)(A_f{z}_2+b_{f}y)-B\left(D\right)(A_f{z}_1+b_{f}u),\hfill \\ & =& A_f\left[A\left(D\right)z_2-B\left(D\right)z_1\right]+b_f\left[A\left(D\right)y-B\left(D\right)u\right],\hfill \end{array}$

which results in

(12)${\dot{z}}_3=A_f{z}_3+b_{f}d,$

by substituting Equations (1) and (10). The solution of Equation (12) for a constant disturbance d was shown in [28].

(13)$z_3=e^{A_{f}t}{z}_3\left(0\right)+A_f^{-1}(e^{A_{f}t}-I)b_{f}d.$

Substituting Equation (10), multiplying both sides of Equation (13) by $c_f^{\top }$, substituting $u_f$ and $y_f$ in Equation (7) and Equation (8), respectively, and, assuming zero initial conditions, finally yielded

(14)$\begin{array}{c}\hfill A\left(D\right)y_f=B\left(D\right)u_f+d_f,\end{array}$

with $d_f=c_f^{\top }{A}_f^{-1}(e^{A_{f}t}-I)b_{f}d$.

The term $d_f$ was decomposed as $d_f=d_{ft}+d_{fs}$ with $d_{ft}=c_f^{\top }{A}_f^{-1}{e}^{A_{f}t}{b}_{f}d$ and $d_{fs}=-c_f^{\top }{A}_f^{-1}{b}_{f}d$. The term $d_{ft}$ exponentially converged to zero because $A_f$ was Hurwitz-stable. However, $d_{fs}$ corresponded to the steady-state value of $d_f$, which was computed using Equations (4) and (9), and the final value theorem in [28], i.e.,

$\begin{array}{ccc}\hfill d_{fs}& =& \underset{s\to 0}{lim}sF\left(s\right)\left(\frac{d}{s}\right)=d,\hfill \\ & =& \underset{s\to 0}{lim}s{c}_f^{\top }{(sI-A_f)}^{-1}{b}_f\left(\frac{d}{s}\right)=-c_f^{\top }{A}_f^{-1}{b}_{f}d.\hfill \end{array}$

The results above showed that $-c_f^{\top }{A}_f^{-1}{b}_f=1$ and $d_{fs}=d$, which allowed us to express Equation (14) as follows:

(15)$A\left(D\right)y_f=B\left(D\right)u_f+d+d_{ft}.$

A linear regression followed from Equation (15)

(16)$\begin{array}{ccc}\hfill y_f^{\left(n\right)}& =& -{\alpha }_1{y}_f^{(n-1)}-\cdots -{\alpha }_n{y}_f+{\beta }_0{u}_f^{\left(m\right)}+{\beta }_1{u}_f^{(m-1)}+\cdots +{\beta }_m{u}_f+d+d_{ft},\hfill \end{array}$

(17)$\begin{array}{ccc}& =& {\theta }^{\top }{\mathrm{\Phi }}_f+d_{ft}, \hfill \end{array}$

where

$\begin{array}{ccc}\hfill \theta & =& {\left[{\alpha }_1,\dots ,{\alpha }_n,{\beta }_0,{\beta }_1,\dots ,{\beta }_m,d\right]}^{\top },\hfill \\ \hfill {\mathrm{\Phi }}_f& =& {\left[-y_f^{(n-1)},\dots ,-y_f,u_f^{\left(m\right)},u_f^{(m-1)},\dots ,u_f,1\right]}^{\top }.\hfill \end{array}$

Since the term $d_{ft}$ decayed exponentially to zero, it was not considered in subsequent developments. Equations (3) and (17) have similar structures, since vector ${\mathrm{\Phi }}_f$ only depends on $u_f$ and $y_f$, all entries as well as $y_f^{\left(n\right)}$ are obtained through the next filters due to having the same characteristic polynomial as Equation (4).

(18)$\begin{array}{ccc}\hfill \frac{\mathcal{L}\left\{u_f^{\left(i\right)}\right\}}{\mathcal{L}\left\{u\right\}}& =& \frac{f_n{s}^i}{s^n+f_1{s}^{n-1}+\dots +f_n};i=0,\dots ,m,\hfill \end{array}$

(19)$\begin{array}{ccc}\hfill \frac{\mathcal{L}\left\{y_f^{\left(i\right)}\right\}}{\mathcal{L}\left\{y\right\}}& =& \frac{f_n{s}^i}{s^n+f_1{s}^{n-1}+\dots +f_n};i=0,\dots ,n.\hfill \end{array}$

3. Parameter Identification Using the LS Method

Consider the model (16) with the following unknown parameters: ${\alpha }_1,\dots ,{\alpha }_n,$${\beta }_0,{\beta }_1,\dots ,{\beta }_m,d$. If the $\rho$ measurements of the variables $y_{f_j}^{\left(n\right)},y_{f_j}^{(n-1)},\dots ,{\dot{y}}_{f_j},y_{f_j},u_{f_j}^{\left(m\right)},u_{f_j}^{(m-1)},$$\dots ,{\dot{u}}_{f_j},u_{f_j}$ with $j=1,\dots ,\rho$, are available, then, the equation to estimate the unknown parameters is the follows:

(20)$\overline{\mathbf{r}}=\mathbf{y}-\mathbf{Ap},$

where

$\begin{array}{c}\hfill \overline{\mathbf{r}}=\left[\begin{array}{c}{\overline{r}}_1\\ {\overline{r}}_2\\ ⋮\\ {\overline{r}}_{\rho }\end{array}\right],\mathbf{y}=\left[\begin{array}{c}{y}_{f_1}^{\left(n\right)}\\ y_{f_2}^{\left(n\right)}\\ ⋮\\ y_{f_{\rho }}^{\left(n\right)}\end{array}\right],\mathbf{A}=\left[\begin{array}{ccccccc}-y_{f_1}^{(n-1)}& \cdots & -y_{f_1}& u_{f_1}^{\left(m\right)}& \cdots & u_{f_1}& 1\\ -y_{f_2}^{(n-1)}& \cdots & -y_{f_2}& u_{f_2}^{\left(m\right)}& \cdots & u_{f_2}& 1\\ ⋮& & ⋮& ⋮& & ⋮& ⋮\\ -y_{f_{\rho }}^{(n-1)}& \cdots & -y_{f_{\rho }}& u_{f_{\rho }}^{\left(m\right)}& \cdots & u_{f_{\rho }}& 1\end{array}\right],\end{array}$

and

$\begin{array}{c}\hfill \mathbf{p}={\left[\begin{array}{ccccccc}{\alpha }_1& \cdots & {\alpha }_n& {\beta }_0& \cdots & {\beta }_m& d\end{array}\right]}^{\top }.\end{array}$

In the standard LS method, the goal is to compute $\mathbf{p}\in {\mathbf{R}}^{\kappa }$, $\kappa =n+m+2$, such that

(21)$J_{LS}=\underset{\mathbf{p}\in {\mathbf{R}}^{\kappa }}{min}{\parallel \mathbf{y}-\mathbf{Ap}\parallel }^2.$

where the operator $\parallel ·\parallel$ stands for the Euclidean norm. The value $\mathbf{p}$ minimizing $J_{LS}$ was provided by [19].

(22)$\begin{array}{ccc}\hfill \mathbf{p}& =& {\mathbf{A}}^{†}\mathbf{y}\hfill \end{array}$

(23)$\begin{array}{ccc}\hfill {\mathbf{A}}^{†}& =& {\left({\mathbf{A}}^{\top }\mathbf{A}\right)}^{-1}{\mathbf{A}}^{\top }=V\left[\begin{array}{cc}{\Sigma }_q^{-1}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]U^{\top }\hfill \end{array}$

The term ${\mathbf{A}}^{†}$ stood for the Moore–Penrose pseudo-inverse of $\mathbf{A}$, which was defined by the Singular value decomposition (SVD), whose properties are described in the next theorem.

Theorem 1

([29]). Let any matrix $\mathbf{A}\in {\mathbb{R}}^{m\times n}$. The SVD of $\mathbf{A}$ is denoted by the following:

$\begin{array}{c}\hfill \mathbf{A}=U\left[\begin{array}{cc}{\Sigma }_q& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]V^{\top },\end{array}$

where $U\in {\mathbb{R}}^{m\times m}$ and $V\in {\mathbb{R}}^{n\times n}$ are orthogonal matrices, and ${\Sigma }_q=\mathrm{diag}({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_q)$. The numbers ${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_q>0$ are called the singular values of $\mathbf{A}$.

4. Parameter Identification Using the LSOD Method

In order to apply the LSOD method, the model (16) has been rewritten, as follows:

(24)$\begin{array}{c}\hfill x_1+{\alpha }_1{x}_2+\cdots +{\alpha }_n{x}_k+{\beta }_0{x}_{k+1}+{\beta }_1{x}_{k+2}+\cdots +{\beta }_m{x}_{\eta }-d=0,\end{array}$

with $x_1=y_f^{\left(n\right)},x_2=y_f^{(n-1)},\dots ,x_k=y_f,x_{k+1}=-u_f^{\left(m\right)},x_{k+2}=-u_f^{(m-1)},\dots ,x_{\eta }=-u_f$ and unknown parameters ${\alpha }_1,\dots ,{\alpha }_n$, ${\beta }_0,{\beta }_1,\dots ,{\beta }_m,d$, which corresponds to the equation describing a hyperplane.

Note that the LSOD method minimizes the sum of the squares of the orthogonal distances from the experimental data points to an approximation function. Then, the parameter identification of an LTI system, parameterized according to Equation (24), and the application of the LSOD method is equivalent to solve a hyperplane-fitting problem. To illustrate this point, Figure 2 shows the orthogonal and vertical distances used by the LSOD and LS methods, respectively, corresponding to the plane fitting in ${\mathbb{R}}^3$.

4.1. The LSOD Method for Plane Fitting

Let the plane ${\pi }_{\eta }\in {\mathbb{R}}^{\eta }$ be defined as follows:

(25)$\begin{array}{c}\hfill {\pi }_{\eta }:a_1{x}_1+a_2{x}_2+a_3{x}_3+\cdots +a_{\eta }{x}_{\eta }+a_0=0,\end{array}$

with the normal vector $\mathbf{N}={\left[\begin{array}{ccccc}{a}_1& a_2& a_3& \cdots & a_{\eta }\end{array}\right]}^{\top }$, and let the experimental data vector $Q_i={\left(x_{1i},x_{2i},\dots ,x_{\eta i}\right)}^{\top }$ with $i=1,2,\dots ,\rho$. Then, the orthogonal distance $\parallel {\mathbf{d}}_i\parallel$ between $Q_i$ and a point on the plane ${\pi }_{\eta }$ could be defined as follows:

Definition 1.

The orthogonal distance $\parallel {\mathbf{d}}_i\parallel$ from the point $Q_i$ to the plane ${\pi }_{\eta }$ is the distance from $Q_i$ to the intersection point $P_i$ between the plane ${\pi }_{\eta }$ and the straight line ${\mathcal{L}}_i$, which is defined along $\mathbf{N}$ and passes through $Q_i$.

In the above definition, $P_i$ is a vector on the plane ${\pi }_{\eta }$. Therefore, the square of the orthogonal distance between $Q_i$ and $P_i$ was the following:

(26)$\begin{array}{c}\hfill \parallel {\mathbf{d}}_i{\parallel }^2={\parallel Q_i-P_i\parallel }^2\end{array}$

and the straight line ${\mathcal{L}}_i$ was defined as the following:

(27)$\begin{array}{c}\hfill {\mathcal{L}}_i=\left\{Q_i+\zeta \mathbf{N}:\zeta \in \mathbb{R}\right\}.\end{array}$

The intersection between ${\pi }_{\eta }$ and ${\mathcal{L}}_i$ was expressed as follows:

$\begin{array}{c}\hfill {\pi }_{\eta }\cap {\mathcal{L}}_i=\left\{(x_{1i}+a_1\zeta ,\dots ,x_{\eta i}+a_{\eta }\zeta )\mid a_1(x_{1i}+a_1\zeta )+\cdots +a_{\eta }(x_{\eta i}+a_{\eta }\zeta )+a_0=0\right\},\end{array}$

and the intersection point was defined by $P_i={(x_{1i}+a_1\zeta ,x_{2i}+a_2\zeta ,\dots ,x_{\eta i}+a_{\eta }\zeta )}^{\top }$ with

$\begin{array}{c}\hfill \zeta =-\frac{a_1{x}_{1i}+a_2{x}_{2i}+\cdots +a_{\eta }{x}_{\eta i}+a_0}{a_1^2+a_2^2+a_3^2+\cdots +a_{\eta }^2}.\end{array}$

Therefore, substituting $Q_i$ and $P_i$ into Equation (26) yielded the following:

$\begin{array}{c}\hfill \parallel {\mathbf{d}}_i{\parallel }^2=\frac{{(a_1{x}_{1i}+a_2{x}_{2i}+\cdots +a_{\eta }{x}_{\eta i}+a_0)}^2}{a_1^2+a_2^2+a_3^2+\cdots +a_{\eta }^2}=\frac{r_i^2}{a_1^2+a_2^2+a_3^2+\cdots +a_{\eta }^2},\end{array}$

where $r_i=a_1{x}_{1i}+a_2{x}_{2i}+\cdots +a_{\eta }{x}_{\eta i}+a_0$.

Now, consider the vector of orthogonal distances $\mathbf{d}={\left[\begin{array}{cccc}\parallel {\mathbf{d}}_1\parallel & \parallel {\mathbf{d}}_2\parallel & \dots & \parallel {\mathbf{d}}_{\rho }\parallel \end{array}\right]}^{\top }$. Then, the cost function to be minimized using the LSOD method would be the following:

(28)$\begin{array}{c}\hfill J_{LSOD}={\parallel \mathbf{d}\parallel }^2=\frac{1}{a_1^2+a_2^2+a_3^2+\cdots +a_{\eta }^2}\sum _{i=1}^{\rho }{r}_i^2.\end{array}$

Therefore, to estimate the parameters $a_0,a_1,\dots ,a_{\eta }$ of the plane (25) ${\pi }_{\eta }\in {\mathbb{R}}^{\eta }$ using the LSOD method, the function (28) must be minimized. Hence, setting the first derivatives of (28) to zero, with respect to $a_0,a_1,\dots ,a_{\eta }$, yielded:

(29)$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial J_{LSOD}}{\partial a_1}}& =& {\displaystyle \frac{2}{a_1^2+a_2^2+a_3^2+\cdots +a_{\eta }^2}}\left[\sum _{i=1}^{\rho }{r}_i{x}_{1i}-{\displaystyle \frac{a_1}{a_1^2+a_2^2+a_3^2+\cdots +a_{\eta }^2}}\sum _{i=1}^{\rho }{r}_i^2\right]=0\hfill \\ & ⋮& \hfill \\ \hfill {\displaystyle \frac{\partial J_{LSOD}}{\partial a_{\eta }}}& =& {\displaystyle \frac{2}{a_1^2+a_2^2+a_3^2+\cdots +a_{\eta }^2}}\left[\sum _{i=1}^{\rho }{r}_i{x}_{\eta i}-{\displaystyle \frac{a_{\eta }}{a_1^2+a_2^2+a_3^2+\cdots +a_{\eta }^2}}\sum _{i=1}^{\rho }{r}_i^2\right]=0\hfill \\ \hfill {\displaystyle \frac{\partial J_{LSOD}}{\partial a_0}}& =& {\displaystyle \frac{2}{a_1^2+a_2^2+a_3^2+\cdots +a_{\eta }^2}}\sum _{i=1}^{\rho }{r}_i=0.\hfill \end{array}$

Since the expression (29) become to a set of nonlinear algebraic equations, they must be solved for $a_0,a_1,\dots ,a_{\eta }$ in order to find a critical point and, then, to obtain a minimum for the cost function (28). From the above, it was clear that the parameter identification problem solved through the LSOD method turns to a nonlinear minimization problem. Regarding this issue, there are several methods to solve this problem, including the Newton method [2,29,30]. The MATLAB routine fminunc is based on the Newton method [31], which was used in this work to solve the set (29).

Finally, the cost function to be minimized in order to estimate the unknown parameters of an LTI system, parameterized according to (24) by the LSOD method, was provided by the following:

(30)$\begin{array}{c}\hfill J_{LSOD}={\parallel \mathbf{d}\parallel }^2=\frac{1}{1+{\alpha }_1^2+\cdots +{\alpha }_n^2+{\beta }_0^2+{\beta }_1^2+\cdots +{\beta }_m^2}\sum _{i=1}^{\rho }{r}_i^2,\end{array}$

with $r_i=x_{1i}+{\alpha }_1{x}_{2i}+\cdots +{\alpha }_n{x}_{ki}+{\beta }_0{x}_{k+1,i}+{\beta }_1{x}_{k+2,i}+\cdots +{\beta }_m{x}_{\eta i}-d$.

4.2. The LSOD Method for a Straight Line

A particular case of the LTI system model (1) was the following:

(31)$\begin{array}{c}\hfill \ddot{y}={\beta }_{0}u+d.\end{array}$

Its corresponding filtered version was the following:

(32)$\begin{array}{c}\hfill {\ddot{y}}_f={\beta }_0{u}_f+d.\end{array}$

The goal was to estimate ${\beta }_0$ using the LSOD method according to ${\ddot{y}}_f$ and $u_f$ measurements. To obtain a suitable analytical solution, instead of resorting to numerical methods, the parameter estimation developed in the next paragraphs only considered the estimation of ${\beta }_0$ and did not consider the disturbance d. Therefore, the problem was reduced to the data fitting to a straight line, passing through the origin, as defined by the following:

(33)$\begin{array}{c}\hfill \ell =\left\{(x_1,x_2)\in {\mathbb{R}}^2\mid x_1+{\beta }_0{x}_2=0\right\},\end{array}$

with $x_1={\ddot{y}}_f$ and $x_2=-u_f$. Hence, applying the LSOD method to estimate the unknown parameters of system (31) yielded the cost function:

(34)$\begin{array}{ccc}\hfill J_{LSOD}& =& {\parallel \mathbf{d}\parallel }^2=\frac{1}{1+{\beta }_0^2}\sum _{i=1}^{\rho }{r}_i^2,\hfill \end{array}$

(35)$\begin{array}{ccc}\hfill r_i& =& x_{1i}+{\beta }_0{x}_{2i},\hfill \end{array}$

which was a special case of (30), and the estimation of the input gain ${\beta }_0$ was also a nonlinear problem. Nevertheless, its simplicity allowed the first and second derivatives of (34) to be explicitly found, with respect to ${\beta }_0$.

The first derivative of Equation (34), with respect to ${\beta }_0$, was given by the following:

(36)$\begin{array}{c}\hfill J_{LSOD}^{{}^{\prime }}\left({\beta }_0\right):=\frac{\partial J_{LSOD}}{\partial {\beta }_0}=-\frac{2}{1+{\beta }_0^2}\left[-\sum _{i=1}^{\rho }{r}_i{x}_{2i}+\frac{{\beta }_0}{1+{\beta }_0^2}\sum _{i=1}^{\rho }{r}_i^2\right].\end{array}$

Equating (36) to zero and solving for ${\beta }_0$ led to a second-degree equation:

(37)$\begin{array}{c}\hfill {\mu }_1{\beta }_0^2+{\mu }_2{\beta }_0+{\mu }_3=0,\end{array}$

with

$\begin{array}{ccc}\hfill {\mu }_1& =& \sum _{i=1}^{\rho }{x}_{1i}{x}_{2i}\hfill \\ & =& {\mathbf{x}}_1^{\top }{\mathbf{x}}_2,\hfill \\ \hfill \\ \hfill {\mu }_2& =& \sum _{i=1}^{\rho }{x}_{1i}^2-\sum _{i=1}^{\rho }{x}_{2i}^2\hfill \\ & =& \parallel {\mathbf{x}}_1{\parallel }^2-{\parallel {\mathbf{x}}_2\parallel }^2,\hfill \\ \hfill \\ \hfill {\mu }_3& =& -\sum _{i=1}^{\rho }{x}_{1i}{x}_{2i}=-{\mu }_1,\hfill \end{array}$

where ${\mathbf{x}}_1={[x_{11},\dots ,x_{1\rho }]}^{\top }$ and ${\mathbf{x}}_2={[x_{21},\dots ,x_{2\rho }]}^{\top }$ are the measurement vectors.

The solutions of Equation (37), which corresponds to the critical points, were provided by the following:

(38)$\begin{array}{ccc}\hfill {\beta }_{0_{1,2}}^{\ast }& =& \frac{-{\mu }_2\pm \sqrt{{\mu }_2^2-4{\mu }_1{\mu }_3}}{2{\mu }_1},\hfill \\ & =& -\frac{{\mu }_2}{2{\mu }_1}\pm \sqrt{\frac{{\mu }_2^2}{4{\mu }_1^2}+1},\hfill \\ & =& -\frac{\parallel {\mathbf{x}}_1{\parallel }^2-{\parallel {\mathbf{x}}_2\parallel }^2}{2\left({\mathbf{x}}_1^{\top }{\mathbf{x}}_2\right)}\pm \sqrt{\frac{(\parallel {\mathbf{x}}_1{\parallel }^2-\parallel {\mathbf{x}}_2{{\parallel }^2)}^2}{4{\left({\mathbf{x}}_1^{\top }{\mathbf{x}}_2\right)}^2}+1},\hfill \end{array}$

with the condition

(39)${\mathbf{x}}_1^{\top }{\mathbf{x}}_2\ne 0.$

Then, the solutions ${\beta }_{0_{1,2}}^{\ast }$ are always real, and they exist as long as the measurement vectors ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$ are not orthogonal. To verify that one of these solutions corresponds to a minimum in (34), the second derivative test was employed.

Definition 2.

Suppose $J_{LSOD}^{{}^{″}}$ is continuous near the critical point ${\beta }_0^{\ast }$

• If $J_{LSOD}^{{}^{\prime }}\left({\beta }_0^{\ast }\right)=0$ and $J_{LSOD}^{{}^{″}}\left({\beta }_0^{\ast }\right)>0$, then $J_{LSOD}$ has a local minimum at ${\beta }_0^{\ast }$.

• If $J_{LSOD}^{{}^{\prime }}\left({\beta }_0^{\ast }\right)=0$ and $J_{LSOD}^{{}^{″}}\left({\beta }_0^{\ast }\right)<0$, then $J_{LSOD}$ has a local maximum at ${\beta }_0^{\ast }$.

The second derivative of (34), with respect to ${\beta }_0$, was expressed, as follows:

(40)$\begin{array}{c}\hfill J_{LSOD}^{{}^{″}}\left({\beta }_0^{\ast }\right):=\frac{{\partial }^2{J}_{LSOD}}{\partial {\beta }_0^2}=\frac{2}{1+{\beta }_0^2}\left[\sum _{i=1}^{\rho }{x}_{2i}^2-\frac{4{\beta }_0}{1+{\beta }_0^2}\sum _{i=1}^{\rho }{r}_i{x}_{2i}-\frac{1-3{\beta }_0^2}{{(1+{\beta }_0^2)}^2}\sum _{i=1}^{\rho }{r}_i^2\right].\end{array}$

Equating the right-hand side of (36) to zero produced the following:

(41)$\begin{array}{c}\hfill \sum _{i=1}^{\rho }{r}_i{x}_{2i}=\frac{{\beta }_0}{1+{\beta }_0^2}\sum _{i=1}^{\rho }{r}_i^2.\end{array}$

Substituting (41) into (40) and using (34) produced:

(42)$\begin{array}{ccc}\hfill J_{LSOD}^{{}^{″}}& =& \frac{2}{1+{\beta }_0^2}\left[\sum _{i=1}^{\rho }{x}_{2i}^2-\frac{1}{1+{\beta }_0^2}\sum _{i=1}^{\rho }{r}_i^2\right]\hfill \\ & =& \frac{2}{1+{\beta }_0^2}\left[\parallel {\mathbf{x}}_2^2\parallel -J_{LSOD}\right].\hfill \end{array}$

Therefore, a critical point ${\beta }_0^{\ast }$ of (34) corresponds to a minimum if

(43)$\parallel {\mathbf{x}}_2^2\parallel -J_{LSOD}\left({\beta }_0^{\ast }\right)>0$

where $J_{LSOD}\left({\beta }_0^{\ast }\right)$ denotes the cost function evaluated at ${\beta }_0^{\ast }$.

5. Servo System Mathematical Models and Parameter Identification Procedure

Servo systems are widely used in CNC machines, industrial processes, assembly robots, electric elevators, electric and hybrid vehicles, trains, and many other areas. Although there has not been a general mathematical model for servo systems, most share similar structures, and certain parameters, such as the input gain due to the power amplifier driving the motor, have more influence in the servo system performance. Moreover, if the requirements of its performance increase, so does the complexity of the mathematical model, which then requires more elaborate control algorithms.

Several key parameters had to be considered when designing servo control systems, including friction phenomena, internal and external disturbances, load changes, un-modeled dynamics, and measurement noise [32,33]. In this work, three servo system models were presented for parameter identification purposes. The first model considered the servo system input gain, a friction term, and a constant disturbance, as the parameters to be estimated. The second model was built by ignoring the disturbance in the first model, and the third model only considered the estimation of the input gain. Moreover, the relevance of the parameter estimates was also assessed by designing feedback control laws for trajectory tracking.

Thus, having several models allows a practitioner to choose one of them for a specific task. For example, some control algorithms, including active disturbance rejection control (ADRC) [34,35], only need the input gain for their implementation, and then, it is only necessary to identify this term, and the model considering only the input gain would be more than adequate for this task.

5.1. First Mathematical Model

Figure 3 depicts a block diagram of a servo system model, which consists of a direct current (DC) motor driven by a power amplifier, which works in current mode though a proportional current controller. The variable y corresponds to the servo angular position. Signal u is the input control signal, $K_A$ is the power amplifier gain, $K_p$ is the proportional gain of the current controller, the constant $K_{cl}$ is the current loop feedback gain, and the constant $K_{H-B}$ is the gain of the H-Bridge in the power amplifier. The variable $\nu$ is the motor input voltage, the constants $R_A$, $L_A$, $K_{\tau }$, and $K_{EMF}$ are respectively, the armature resistance, the armature inductance, the torque constant, and the back electromotive force constant. The parameters J and B are the lumped DC motor and load inertia and viscous friction coefficient respectively. The term ${\tau }_d$ corresponds to constant disturbances that include parasitic constant voltages produced inside the power amplifier. The task of the current loop is to keep the armature current $i_a$ proportional to the control signal u, and to speed up the electrical dynamics of the servo system. Due to this feature, the electrical dynamics could have been overlooked at low frequencies.

Thus, by ignoring the electrical dynamics of the servo system and assuming that the external disturbances were constant, we obtained the following model:

(44)$\ddot{y}+{\alpha }_1\dot{y}={\beta }_{0}u+d$

where

$\begin{array}{cc}\hfill {\alpha }_1={\displaystyle \frac{B}{J}}+{\displaystyle \frac{K_{\tau }{K}_{EMF}}{J(R_a+K_p{K}_{H-B}{K}_{cl})}},& {\beta }_0={\displaystyle \frac{K_A{K}_p{K}_{H-B}{K}_{\tau }}{J(R_a+K_p{K}_{H-B}{K}_{cl})}},d={\displaystyle \frac{R_a{\tau }_d}{J(R_a+K_p{K}_{H-B}{K}_{cl})}}.\end{array}$

Then, ${\alpha }_1$ is related to the damping due to the mechanical friction and the back electromotive force, ${\beta }_0$ is the input gain, and d is a function of the constant disturbance.

According to (17), the corresponding filtered version of (44) was the following:

(45)${\ddot{y}}_f+{\alpha }_1{\dot{y}}_f={\beta }_0{u}_f+d$

and the exponential decaying term $d_{ft}$, due to the filtering, was neglected.

The approximation plane for (45) was obtained by defining the change of variables: $(x_1,x_2,x_3)=({\ddot{y}}_f,-u_f,{\dot{y}}_f)$

(46)$\begin{array}{c}\hfill {\pi }_3=\left\{(x_1,x_2,x_3)\in {\mathbb{R}}^3\mid x_1+{\beta }_0{x}_2+{\alpha }_1{x}_3-d=0\right\},\end{array}$

with normal vector $\mathbf{N}={\left[\begin{array}{ccc}1& {\beta }_0& {\alpha }_1\end{array}\right]}^{\top }$. Applying the LSOD method to the plane (46), the cost function to be minimized to estimate ${\beta }_0,{\alpha }_1$ and d was determined by the following:

(47)$\begin{array}{c}\hfill J_{LSOD}={\parallel \mathbf{d}\parallel }^2=\frac{1}{1+{\beta }_0^2+{\alpha }_1^2}\sum _{i=1}^{\rho }{r}_i^2,\end{array}$

The residuals corresponded to $r_i=x_{1i}+{\beta }_0{x}_{2i}+{\alpha }_1{x}_{3i}-d$.

The LS method, as described in Section 3, was applied to fit the plane (46), considering $\rho$ measurements, which resulted in the following:

(48)$\overline{\mathbf{r}}=\left[\begin{array}{c}{\overline{r}}_1\\ ⋮\\ {\overline{r}}_{\rho }\end{array}\right]\in {\mathbb{R}}^{\rho \times 1},\mathbf{y}=\left[\begin{array}{c}{x}_{11}\\ ⋮\\ x_{1\rho }\end{array}\right]\in {\mathbb{R}}^{\rho \times 1},{\mathbf{A}}_3=\left[\begin{array}{ccc}-x_{31}& -x_{21}& 1\\ ⋮& ⋮& ⋮\\ -x_{3\rho }& -x_{2\rho }& 1\end{array}\right]\in {\mathbb{R}}^{\rho \times 3},{\mathbf{p}}_3=\left[\begin{array}{c}{\alpha }_1\\ {\beta }_0\\ d\end{array}\right]\in {\mathbb{R}}^{3\times 1}.$

5.2. Second Mathematical Model

The second mathematical model to be identified was obtained by neglecting the disturbance term d in (44), i.e.,

(49)$\ddot{y}+{\alpha }_1\dot{y}={\beta }_{0}u,$

and the corresponding filtered version was provided by

(50)${\ddot{y}}_f+{\alpha }_1{\dot{y}}_f={\beta }_0{u}_f.$

The approximation plane with the change of variable $(x_1,x_2,x_3)=({\ddot{y}}_f,-u_f,{\dot{y}}_f)$ to be identified was the following:

(51)$\begin{array}{c}\hfill {\pi }_2=\left\{(x_1,x_2,x_3)\in {\mathbb{R}}^3\mid x_1+{\beta }_0{x}_2+{\alpha }_1{x}_3=0\right\},\end{array}$

with the normal vector $\mathbf{N}={\left[\begin{array}{ccc}1& {\beta }_0& {\alpha }_1\end{array}\right]}^{\top }$.

The LSOD method applied to the plane (51) to estimate ${\beta }_0$ and ${\alpha }_1$ resulted in the cost function:

(52)$\begin{array}{c}\hfill J_{LSOD}={\parallel \mathbf{d}\parallel }^2=\frac{1}{1+{\beta }_0^2+{\alpha }_1^2}\sum _{i=1}^{\rho }{r}_i^2,\end{array}$

and the residuals were $r_i=x_{1i}+{\beta }_0{x}_{2i}+{\alpha }_1{x}_{3i}$.

Now, considering $\rho$ measurements, we applied the LS method to (51) yielded the following:

(53)$\overline{\mathbf{r}}=\left[\begin{array}{c}{\overline{r}}_1\\ ⋮\\ {\overline{r}}_{\rho }\end{array}\right]\in {\mathbb{R}}^{\rho \times 1},\mathbf{y}=\left[\begin{array}{c}{x}_{11}\\ ⋮\\ x_{1\rho }\end{array}\right]\in {\mathbb{R}}^{\rho \times 1},{\mathbf{A}}_2=\left[\begin{array}{cc}-x_{31}& -x_{21}\\ ⋮& ⋮\\ -x_{3\rho }& -x_{2\rho }\end{array}\right]\in {\mathbb{R}}^{\rho \times 2},{\mathbf{p}}_2=\left[\begin{array}{c}{\alpha }_1\\ {\beta }_0\end{array}\right]\in {\mathbb{R}}^{2\times 1}.$

5.3. Third Mathematical Model

Finally, if only the input gain of the servo system needs to be estimated, the ${\alpha }_1\dot{y}$ and d terms are neglected in (44), thus yielding:

(54)$\ddot{y}={\beta }_{0}u,$

The filtered version of the above model is:

(55)${\ddot{y}}_f={\beta }_0{u}_f.$

Since the mathematical model (55) was a double-integrator system and considering the change of variable $(x_1,x_2)=({\ddot{y}}_f,-u_f)$, the straight line to be identified was the following:

(56)$\begin{array}{c}\hfill \ell =\left\{(x_1,x_2)\in {\mathbb{R}}^2\mid x_1+{\beta }_0{x}_2=0\right\}\end{array}$

and the function to be minimized to estimate ${\beta }_0$ of the model (56) by the LSOD method was (34).

The LS method applied to fit the line (56) resulted in the following terms, considering $\rho$ measurements:

(57)$\overline{\mathbf{r}}=\left[\begin{array}{c}{\overline{r}}_1\\ ⋮\\ {\overline{r}}_{\rho }\end{array}\right]\in {\mathbb{R}}^{\rho \times 1},\mathbf{y}=\left[\begin{array}{c}{x}_{11}\\ ⋮\\ x_{1\rho }\end{array}\right]\in {\mathbb{R}}^{\rho \times 1},{\mathbf{A}}_1=\left[\begin{array}{c}-x_{21}\\ ⋮\\ -x_{2\rho }\end{array}\right]\in {\mathbb{R}}^{\rho \times 1},{\mathbf{p}}_1={\beta }_0\in \mathbb{R}.$

As was previously mentioned, parameter identification using the LSOD method corresponds to a nonlinear minimization problem. Then, the minimization of the cost functions (47) and (52) was carried out by the fminunc MATLAB routine, whereas the cost function (34) was minimized using the LSOD analytical method, as described in Section 4.2. However, the parameter identification of the planes (46) and (51), and the line (56), was also carried out by the LS method through (22), together with (48), (53), and (57).

6. Servo System Parameter Identification Using Experimental Measurements

The experimental environment used for data acquisition and real-time experiments can be see in Appendix A. Since the models (44), (49), and (54) had similar information about the behavior of the servo system, the aim of the experiments was twofold: first, to estimate the parameters of these models, and second, to evaluate the identified models by designing and evaluating a trajectory-tracking controller.

However, it is well known that for parameter identification using the LS algorithm, a persistency-of-excitation (PE) condition has to be met in order to ensure the invertibility of ${\mathbf{A}}^{\top }\mathbf{A}$ in (23) [4,5,6]. To fulfill this condition, a test signal with a wide frequency spectrum is typically used. Examples of these signals include filtered white noise and pseudo-random binary signals. For the LSOD algorithm, and to the best of the authors’ knowledge, there is no previous work on how to choose an excitation test signal.

This work considered two test signals, namely filtered white noise and a sine wave. This choice allowed us to study the behavior of the LS and LSOD algorithms on signals with radically different spectra.

In order to perform parameter identification, a proportional derivative (PD) feedback control law was proposed to stabilize the servo system:

(58)$u=K_{p}e-K_d{\dot{y}}_e$

where e denotes the position error defined as $e=y_r-y$, $y_r$ is the reference, the position y is obtained by an optical encoder, and ${\dot{y}}_e$ is an estimate obtained from the motor position measurements through the next filter

(59)$\frac{\mathcal{L}\left\{{\dot{y}}_e\right\}}{\mathcal{L}\left\{y\right\}}=\left[\frac{300s}{s+300}\right]\left[\frac{300}{s+300}\right],$

where $\mathcal{L}$$\{·\}$ corresponds to the Laplace operator. The PD controller gains were adjusted empirically until obtaining an acceptable response from position, and then the proportional and derivative gains were set to $K_p=2$ and $K_d=0.09$, respectively.

The measurements of the servo system were processed using the filters (18) and (19) to obtain estimates of the angular velocity and acceleration, respectively, which were not available from measurements (see Figure 4). The poles of the filters must be real repeated value to avoid overshoots and they must have fast dynamics. For this study, the filters were defined as follows:

$\begin{array}{c}{y}_j\Rightarrow {\displaystyle \frac{400s}{s^2+40s+400}}\Rightarrow {\dot{y}}_{jf}\hfill \\ \\ y_j\Rightarrow {\displaystyle \frac{400s^2}{s^2+40s+400}}\Rightarrow {\ddot{y}}_{jf}\hfill \\ \\ u_j\Rightarrow {\displaystyle \frac{400}{s^2+40s+400}}\Rightarrow u_{jf}\hfill \end{array}$

In order to apply the LS and the LSOD methods, the values of $u_{jf}$, ${\dot{y}}_{jf}$, ${\ddot{y}}_{jf},j=1,\cdots ,\rho$ were acquired at the time instants $t_1,\cdots ,t_m$ (see Figure 4). For the LS method, this allowed us to build ${\mathbf{A}}_1$, ${\mathbf{A}}_2$, ${\mathbf{A}}_3$ and $\mathbf{y}$ in (48), (53), and (57).

Experimental Identification Results

Table 1, Table 2 and Table 3 show the results of the parameter identification of the models (45), (50) and (55), where ${\beta }_0^{\ast }$, ${\alpha }_1^{\ast }$, and $d^{\ast }$ are the estimates of ${\beta }_0$, ${\alpha }_1$ and d, respectively. The white noise power and the sample time in the Simulink band-limited white noise block were $0.01$ and $0.01$, respectively. The signal generator block producing a sine-wave signal was used with an amplitude of 0.5 s and a frequency of 0.25 Hertz. A zero-order hold (ZOH), with a sample time of $0.05$, was employed to obtain the samples and the experiments in each case were 10 s in duration, which provided $\rho =200$ samples.

The parameter estimates produced by the LS and LSOD methods using the filtered white-noise signal were quite similar for the thee models. A well-known aspect of the PE condition in gradient and LS algorithms for adaptive control and parameter identification was that it guaranteed estimates that were close to those obtained in the absence of disturbances [17,18]. Therefore, the results obtained by the LS method for the three models were not surprising.

However, in the case of the sine-wave signal, the results varied. It is worth noting that the estimate ${\beta }_0^{\ast }$ produced by the LSOD algorithm remained close to a value of 40, regardless of the excitation signal and the identified model, which contrasted with the results obtained using the LS method. In addition, as shown in Table 1 and Table 2, there was also an increase in the estimation of parameter ${\alpha }_1^{\ast }$ when the sine wave and the LSOD method were used in comparison with the white-noise results.

However, the LS method produced very poor estimates of ${\beta }_0^{\ast }$ for the three models when the sine-wave signal was employed as excitation. A possible explanation of these outcomes was that a sine wave has two spectral lines, which is theoretically enough to fulfill the persistent excitation (PE) condition in order to identify a single parameter [17]. However, the presence of disturbances and un-modeled terms introduced a significant bias in the estimate. However, the LS algorithm also produced poor estimates for the models with two and three parameters. In this regard, the LS method is an ill-conditioned problem [29], and it requires an input signal that satisfies the PE condition in order to obtain consistent estimates [4].

In addition, for the estimation of ${\beta }_0^{\ast }$ in model (54) using the LSOD method, the solution was the result of the quadratic Equation (37), and only the orthogonality condition (39) had to be met to obtain an estimate, see Table 4. This condition did not appear to be related to the frequency content of the excitation signal, thus indicating the possibility of using the LSOD method with smooth excitation signals in applications where more aggressive signals with high-frequency content cannot be employed for parameter identification purposes.

7. Real-Time Position Control and Servo System Model Validation

The estimated models (46), (51), and (56) were evaluated by using their parameter estimates for designing trajectory tracking feedback control laws. To evaluate the closed-loop performance of the integral of the squared errors ($ISE$), the integral of the absolute value of the error ($IAE$), the integral of the absolute value of the control ($IAC$), and the integral of the absolute value of the control variation ($IACV$) indexes were employed. The $IACV$ evaluates the chattering level and high-frequency content of the control signal. These performance indexes [36] are defined as follows:

$\begin{array}{c}{\displaystyle ISE={\int }_{{\tau }_1}^{{\tau }_2}{e}^2\left(t\right)dt,IAE={\int }_{{\tau }_1}^{{\tau }_2}\left|e\left(t\right)\right|dt,}\\ \\ {\displaystyle IACV={\int }_{{\tau }_1}^{{\tau }_2}\left|{\displaystyle \frac{du\left(t\right)}{dt}}\right|dt,IAC={\int }_{{\tau }_1}^{{\tau }_2}\left|u\left(t\right)\right|dt.}\end{array}$

The time window for evaluating the above indexes was defined by ${\tau }_1=0$ s and ${\tau }_2=20$ s.

The time-varying reference $y_r$ and its time derivatives used in all the experiments were the following:

$\begin{array}{ccc}\hfill y_r& =& y_1\hfill \\ \hfill {\dot{y}}_r& =& {\dot{y}}_1\hfill \\ \hfill {\ddot{y}}_r& =& a\omega \pi {\dot{y}}_2\hfill \end{array}$

which corresponds to the solution of the Duffing oscillator equation

$\begin{array}{ccc}\hfill {\dot{y}}_1& =& a\omega \pi y_2,\hfill \\ \hfill {\dot{y}}_2& =& a\omega \pi \left[-0.25y_2+y_1-1.05y_1^3+0.3sin\left(\omega \pi t\right)\right],\hfill \end{array}$

with $a=0.5$, $\omega =0.55$. The gains in all the cases described in the next subsections were tuned as $K_p=225$ and $K_d=21$.

7.1. Servo System Model with Three Parameter Estimates

Consider the servo system model (44), $\ddot{y}=-{\alpha }_1\dot{y}+{\beta }_0{u}_3+d$, and its parameter estimates depicted in Table 1. The control law used for the trajectory-tracking task was the following:

(60)$u_3={\displaystyle \frac{1}{{\beta }_0^{\ast }}}\left[{\ddot{y}}_r+K_p(y_r-y)+K_d\left({\dot{y}}_r-{\dot{y}}_e\right)+{\alpha }_1^{\ast }{\dot{y}}_e-d^{\ast }\right],$

Table 5 continues the performance indexes of the experiments using the control law (60). The closed-loop responses, the position errors, and the control signals are depicted in Figure 5, Figure 6 and Figure 7, respectively. The experimental results in Table 5 showed that the smallest values of the $ISE$ and $IAE$ performance indexes were obtained with the parameters produced by the LSOD method.

We noted that the control law (60), computed with the parameter estimates from the LS method under sine-wave excitation, produced the highest value of the $IACV$ index. This was due to the poor estimate of the input gain ${\beta }_0$, which was the lowest in Table 1. Moreover, the values of the $IAC$ index were similar for the two identification methods.

7.2. Servo System Model with Two Parameter Estimates

The servo system model (49), $\ddot{y}=-{\alpha }_1\dot{y}+{\beta }_0{u}_2$ and the parameters in Table 2 were considered in this study. The control law used for the trajectory-tracking task was provided by the following:

(61)$u_2={\displaystyle \frac{1}{{\beta }_0^{\ast }}}\left[{\ddot{y}}_r+K_p(y_r-y)+K_d\left({\dot{y}}_r-{\dot{y}}_e\right)+{\alpha }_1^{\ast }{\dot{y}}_e\right]$

The closed-loop responses, the position errors, and the control signals are shown in Figure 8, Figure 9 and Figure 10, respectively. The experimental results from Table 6 indicated that the control law, computed using the parameter estimates, generated by the filtered white-noise signal and provided by the LSOD method produced the lowest values of the $ISE$ and $IAE$ indexes. However, when the parameters obtained with the sine-wave signal were employed in the control law, the $ISE$ index, related to the LS algorithm, was the lowest but at the expense of an increased value of the $IACV$ index, i.e., with a control signal having significant high-frequency components. Table 2 shows that the value ${\beta }_0^{\ast }$ was too small, and then its inverse was large, thus producing a larger control signal. Further, the LSOD method produced acceptable performance with lower values of control signal activity, according to the $IACV$ index. Finally, we noted that all the index values related to the LSOD method did not exhibit significant differences for the two excitation signals.

7.3. Servo System Model with One Estimated Parameter

The next control law is designed for the servo system model (54), $\ddot{y}={\beta }_0{u}_1$, with the parameter values displayed in Table 3

(62)$u_1={\displaystyle \frac{1}{{\beta }_0^{\ast }}}\left[{\ddot{y}}_r+K_p(y_r-y)+K_d\left({\dot{y}}_r-{\dot{y}}_e\right)\right]$

Figure 11, Figure 12 and Figure 13 depict the closed-loop responses, the position errors, and control signals, respectively.

According to Table 7, the performance indexes corresponding to the white-noise signal were similar for both the LS and LSOD methods. We noted that the $ISE$ and $IAE$ indexes, obtained with parameter estimates using the LS method with a sine-wave signal, were the smallest, but the values of the $IACV$ and $IAC$ indexes were the highest. For model (49), a very low value of the estimate ${\beta }_0^{\ast }$ had a larger inverse, thus producing a large control activity. As a consequence, the control displayed large oscillations, as it could be verified in Figure 13. These signals could damage the servo system and/or the load it drives. We noted the LSOD method produced essentially the same performance, regardless of the excitation signal used for obtaining the parameter estimates.

7.4. Analysis of the Results Produced by the LSOD Method

Here, we analyzed the parameter estimates and the closed-loop performance obtained with the LSOD method, regarding the servo system models and the excitation signal.

Table 8 continues all the parameter estimates produced by the LSOD method. One interesting outcome was that the value of the input gain estimate ${\beta }_0^{\ast }$ remained in the range 38–43 for all the models and excitation signals. Furthermore, a consequence of the low variability of this estimate was that the closed-loop performance was similar for the same control law regardless of the excitation signal, as it could be verified in Table 9 and Table 10. These tables also showed that designing a control law using a model with more parameters improved performance. For example, the values of the $ISE$ and $IAE$ indexes considerably decreased using control law $u_3$ with three parameters, as compared to the control law $u_1$, which only considered the input gain. These results were not surprising because more information was used for designing control laws $u_2$ and $u_3$. However, the $IACV$ and $IAC$ indexes had similar values in all the cases.

We should also emphasize that for the same control law, the choice of the excitation signal for parameter estimation appeared to have no influence on the closed-loop performance. A practical consequence of this was that the LSOD method could be used in systems where applying excitation signals with high-frequency content could be impractical or dangerous. Examples include quad-rotors, mobile robots, and servo systems driving robot manipulators, where the robot joints are driven though gearboxes.

8. Conclusions

In the previous sections, the least squares of orthogonal distances (LSOD) method was evaluated for the parameter identification of servo systems. The main conclusions of this study included the following:

• The LSOD method could be a practical alternative to the classic least squares for parameter identification of linear time-invariant systems. The parameter identification problem was transformed into data fitting to a hyperplane, which was subsequently solved as a nonlinear minimization problem.

• The LSOD method required well-established numerical routines based on the Newton method.

• For the parameter identification of the servo model having only one parameter to identify, i.e., the input gain, the parameter estimate was computed analytically by solving a simple quadratic equation.

• Three models of a servo system with one, two and three parameters were estimated using a sine wave and a filtered white noise signals. Tracking control laws were designed and experimentally tested for the three identified models. The overall best tracking performance was obtained when the tracking control laws were computed using the parameter estimates produced by the LSOD method irrespective of the signal used for exciting the servo system and the number of identified parameters.

The LSOD estimates showed consistency, regardless of the excitation signal spectra, which suggested the possibility of estimating the input gain in second-order electro-mechanical systems using this method, especially when it might be impractical or even dangerous to use aggressive test signals for parameter identification purposes. Future research could include applying the LSOD method for parameter identification in nonlinear systems and determining the fewest terms required in the mathematical models, using the methodology of Burton et al. [37] and the LSOD method.

Author Contributions

Conceptualization, C.V.-J. and R.G.; methodology, L.A.C.-C. and L.L.; software, L.A.C.-C. and L.L.; validation, C.V.-J., R.G. and E.A.; formal analysis, L.A.C.-C. and L.L.; investigation, L.A.C.-C., L.L. and E.A.; writing—original draft preparation, R.G., L.A.C.-C. and L.L.; writing—review and editing, E.A., L.A.C.-C. and L.L.; supervision, C.V.-J. and R.G. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Acknowledgments

Luis Luna and Erick Asiain want to thank the support of Martha Hernandez, Edgar Ramirez and Natisma Lopez from CECyT 9-IPN. The work of Luis Cantera, Luis Luna, Erick Asiain was supported by scholarships from CONACyT-MEXICO.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Generate image description

Figure 1. Straight line fitting problem. (a) Classic least squares algorithm.; (b) Least squares of orthogonal distances algorithm.

Generate image description

Figure 2. Approximation plane.

Generate image description

Figure 3. Block diagram of a servo system model.

Generate image description

Figure 4. Servo system signal filtering block diagram.

Generate image description

Figure 5. Reference and servo system angular positions produced by control law (60) with parameter estimates [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Generate image description

Figure 6. Position errors produced by control law (60) with parameter estimates [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Generate image description

Figure 7. Control signals produced by control law (60) with parameter estimates [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Generate image description

Figure 8. Reference and servo system angular positions produced by control law (61) with parameter estimates [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Generate image description

Figure 9. Position errors produced by control law (61) with parameter estimates [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Generate image description

Figure 10. Control signals produced by control law (61) with parameter estimates [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Generate image description

Figure 11. Reference and servo system angular positions produced by control law (62) with the parameter estimate [Forumla omitted. See PDF.].

Generate image description

Figure 12. Position errors produced by control law (62) with the parameter estimate [Forumla omitted. See PDF.].

Generate image description

Figure 13. Control signals produced by control law (62) with the parameter estimate [Forumla omitted. See PDF.].

Table 1

Parameter identification of the servo system model (44).

Model Parameter	White-Noise Signal		Sine-Wave Signal
	LS	LSOD	LS	LSOD
${\beta }_0^{\ast }$	40.809663	41.402230	15.300981	38.682501
${\alpha }_1^{\ast }$	2.142594	2.174090	2.269066	5.868168
$d^{\ast }$	−1.653451	−1.678005	−0.841450	−2.008895

Table 2

Parameter identification of the servo system model (49).

Model Parameter	White-Noise Signal		Sine-Wave Signal
	LS	LSOD	LS	LSOD
${\beta }_0^{\ast }$	39.924891	41.380829	7.520022	42.460353
${\alpha }_1^{\ast }$	2.111703	2.189721	1.123041	6.704563

Table 3

Parameter identification of the servo system model (54).

Model Parameter	White-Noise Signal		Sine Wave Signal
	LS	LSOD	LS	LSOD
${\beta }_0^{\ast }$	37.679891	41.461638	3.484775	39.823089

Table 4

Orthogonality condition ${\mathbf{x}}_1^{\top }{\mathbf{x}}_2$ (39) and the value of Minimum condition $J_{LSOD}^{{}^{″}}\left({\beta }_0^{\ast }\right)$ (42) for the parameter identification of the servo system model (54) using the LSOD method.

Condition	White-Noise Signal	Sine-Wave Signal
${\mathbf{x}}_1^{\top }{\mathbf{x}}_2$	617.349730	8.607933
$J_{LSOD}^{{}^{″}}\left({\beta }_0^{\ast }\right)$	0.017312	0.000272

Table 5

Performance indexes produced by control law (60) with parameter estimates ${\beta }_0^{\ast }$, ${\alpha }_1^{\ast }$ and $d^{\ast }$.

Index	White-Noise Signal		Sine-Wave Signal
	LS	LSOD	LS	LSOD
$ISE$	18.6215	16.4898	13.8181	11.9049
$IAE$	16.0305	15.1407	13.6463	12.4864
$IACV$	9.7414	9.9285	23.1436	10.0730
$IAC$	1.7290	1.6977	1.7812	1.7345

Table 6

Performance indexes produced by control law (61) with parameter estimates ${\beta }_0^{\ast }$ and ${\alpha }_1^{\ast }$.

Index	White-Noise Signal		Sine-Wave Signal
	LS	LSOD	LS	LSOD
$ISE$	34.7949	29.2900	16.0997	24.0802
$IAE$	21.8799	20.0538	14.34865	18.57136
$IACV$	9.94869	9.8466	45.6438	9.1073
$IAC$	1.7231	1.6963	1.8044	1.6745

Table 7

Performance indexes produced by control law (62) with the parameter estimates ${\beta }_0^{\ast }$.

Index	White-Noise Signal		Sine-Wave Signal
	LS	LSOD	LS	LSOD
$ISE$	45.1281	45.3952	16.4636	47.1772
$IAE$	25.1922	25.5792	13.9755	25.8390
$IACV$	10.6522	10.3371	200.3603	10.2424
$IAC$	1.7399	1.7053	2.6660	1.7106

Table 8

Parameter estimates obtained by the LSOD method.

ParameterEstimates	White-Noise Signal			Sine-Wave Signal
	Model 3	Model 2	Model 1	Model 3	Model 2	Model 1
${\beta }_0^{\ast }$	41.461638	41.380829	41.402230	39.823089	42.460353	38.682501
${\alpha }_1^{\ast }$	–	2.189721	2.174090	–	6.704563	5.868168
$d^{\ast }$	–	–	−1.678005	–	–	−2.008895

Table 9

Performance indexes using the LSOD method and white noise as an excitation signal using the control laws (60)–(62).

Control Law	Parameter Estimates	$\mathit{ISE}$	$\mathit{IAE}$	$\mathit{IACV}$	$\mathit{IAC}$
$u_1$	${\beta }_0^{\ast }$	45.3952	25.5792	10.3371	1.7053
$u_2$	${\beta }_0^{\ast }$, ${\alpha }_1^{\ast }$	29.29	20.0538	9.8466	1.6963
$u_3$	${\beta }_0^{\ast }$, ${\alpha }_1^{\ast }$, $d^{\ast }$	16.4898	15.1407	9.9285	1.6977

Table 10

Performance indexes using the LSOD method and a sine wave as excitation signal using the control laws (60)–(62).

Control Law	ParameterEstimates	$\mathit{ISE}$	$\mathit{IAE}$	$\mathit{IACV}$	$\mathit{IAC}$
$u_1$	${\beta }_0^{\ast }$	47.1772	25.8390	10.2424	1.7106
$u_2$	${\beta }_0^{\ast }$, ${\alpha }_1^{\ast }$	24.0802	18.57136	9.1073	1.6745
$u_3$	${\beta }_0^{\ast }$, ${\alpha }_1^{\ast }$, $d^{\ast }$	11.9049	12.4864	10.0730	1.7345

Appendix A. Servo System Experimental Setup

The servo system setup used in the experiments is shown in Figure A1. It consisted of the following components:

• A Quanser QPID data acquisition card.

• A DC tachogenerator plus an optical encoder from Servotek model 1024PTSA-7388F-1, which provided velocity and position measurements.

• Two mechanically coupled brushed DC motors from Clifton Precision model JDTH-2250-BQ-IC. One was used for parameter identification, while the other was used as a load.

• A pair of amplifiers, model 30A20AC PWM, from Advanced Motion Controls.

• The control algorithms were implemented using the MATLAB/Simulink programming platform under the Quanser QUARC real-time control environment and using a sampling period of 1 ms and the Euler-ode1 integration method.

Generate image description

Figure A1. Experimental Setup.

References

1. Bard, Y. Nonlinear Parameter Estimation; Academic Press: New York, NY, USA, 1974.

2. Englezos, P.; Kalogerakis, N. Applied Parameter Estimation for Chemical Engineers; CRC Press: Boca Raton, FL, USA, 2000.

3. Sorenson, H.W. Parameter Estimation: Principles and Problems; M. Dekker: New York, NY, USA, 1980; Volume 9.

4. Åström, K.J.; Eykhoff, P. System identification—A survey. Automatica; 1971; 7, pp. 123-162. [DOI: https://dx.doi.org/10.1016/0005-1098(71)90059-8]

5. Ljung, L. System Identification: Theory for User; Prentice Hall: Upper Saddle River, NJ, USA, 1999.

6. Isermann, R.; Münchhof, M. Identification of Dynamic Systems: An Introduction with Applications; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2010.

7. Petras, I.; Podlubny, I. Least Squares or Least Circles? A comparison of classical regression and orthogonal regression. Chance; 2010; 23, pp. 38-42. [DOI: https://dx.doi.org/10.1080/09332480.2010.10739804]

8. Adcock, R.J. Note on the method of least squares. Analyst; 1877; 4, pp. 183-184. [DOI: https://dx.doi.org/10.2307/2635777]

9. Adcock, R.J. A problem in least squares. Analyst; 1878; 5, pp. 53-54. [DOI: https://dx.doi.org/10.2307/2635758]

10. Ahn, S.J. Least Squares Orthogonal Distance Fitting of Curves and Surfaces in Space; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2004; Volume 3151.

11. Ahn, S.J.; Rauh, W.; Warnecke, H.J. Least-squares orthogonal distances fitting of circle, sphere, ellipse, hyperbola, and parabola. Pattern Recognit.; 2001; 34, pp. 2283-2303. [DOI: https://dx.doi.org/10.1016/S0031-3203(00)00152-7]

12. Penner, A. Fitting Splines to a Parametric Function; Springer: Berlin/Heidelberg, Germany, 2019.

13. López, A.; Olazagoitia, J.L.; Marzal, F.; Rubio, M.R. Optimal parameter estimation in semi-empirical tire models. Proc. Inst. Mech. Eng. Part D J. Automob. Eng.; 2019; 233, pp. 73-87. [DOI: https://dx.doi.org/10.1177/0954407018779851]

14. Poch, J.; Villaescusa, I. Orthogonal distance regression: A good alternative to least squares for modeling sorption data. J. Chem. Eng. Data; 2012; 57, pp. 490-499. [DOI: https://dx.doi.org/10.1021/je201070u]

15. Keles, T. Comparison of Classical Least Squares and Orthogonal Regression in measurement error models. Int. Online J. Educ. Sci.; 2018; 10, pp. 200-214.

16. Åström, K.J. Lectures on the Identification Problem: The Least Squares Method; Technical Report Division of Automatic Control, Lund Institute of Technology: Lund, Sweden, 1968.

17. Sastry, S.; Bodson, M. Adaptive Control: Stability, Convergence, and Robustness; Prentice-Hall, Inc.: Hoboken, NJ, USA, 1989.

18. Ioannou, P.A.; Sun, J. Robust Adaptive Control; Prentice-Hall, Inc.: Hoboken, NJ, USA, 1995.

19. Keesman, K.J. System Identification: An Introduction; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2011.

20. Garrido, R.; Miranda, R. DC servomechanism parameter identification: A closed loop input error approach. ISA Trans.; 2012; 51, pp. 42-49. [DOI: https://dx.doi.org/10.1016/j.isatra.2011.07.003] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21868009]

21. Abdulle, A.; Wanner, G. 200 years of least squares method. Elem. Der Math.; 2002; 57, pp. 45-60. [DOI: https://dx.doi.org/10.1007/PL00000559]

22. Bjorck, A. Numerical Methods for Least Squares Problems; SIAM: Philadelphia, PA, USA, 1996; Volume 51.

23. Lawson, C.L.; Hanson, R.J. Solving Least Squares Problems; SIAM: Philadelphia, PA, USA, 1995.

24. Terzic, B.; Jadric, M. Design and implementation of the extended Kalman filter for the speed and rotor position estimation of brushless DC motor. IEEE Trans. Ind. Electron.; 2001; 48, pp. 1065-1073. [DOI: https://dx.doi.org/10.1109/41.969385]

25. Andreescu, G.D.; Rabinovici, R. Torque-speed adaptive observer and inertia identification without current transducers for control of electric drives. Proc. Rom. Acad. Ser. A; 2003; 4, 8.

26. Cortez-Vega, R.; Maldonado, J.; Garrido, R. Parameter Identification using PSO under measurement noise conditions. Proceedings of the 2019 6th International Conference on Control, Decision and Information Technologies (CoDIT); Paris, France, 23–26 April 2019; pp. 103-108.

27. Cantera, L.A.C.; Luna, L.; Vargas-Jarillo, C.; Garrido, R. Parameter estimation of a linear ultrasonic motor using the least squares of orthogonal distances algorithm. Proceedings of the 2019 16th International Conference on Electrical Engineering, Computing Science and Automatic Control (CCE); Mexico City, Mexico, 9–11 November 2019; pp. 1-6.

28. Ogata, K. Modern Control Engineering; Prentice & Hall Inc.: Upper Saddle River, NJ, USA, 2002.

29. Atkinson, K.E. An Introduction to Numerical Analysis; John Wiley & Sons: New York, NY, USA, 2008.

30. Dennis, J.E., Jr.; Schnabel, R.B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; SIAM: Hoboken, NJ, USA, 1996.

31. MathWorks. Optimization Toolbox User’s Guide; MathWorks: Natick, MA, USA, 2020.

32. Firoozian, R. Servo Motors and Industrial Control Theory; Springer: Berlin/Heidelberg, Germany, 2014.

33. Nakamura, M.; Goto, S.; Kyura, N. Mechatronic Servo System Control: Problems in Industries and Their Theoretical Solutions; Springer: Berlin/Heidelberg, Germany, 2004; Volume 300.

34. Han, J. From PID to active disturbance rejection control. IEEE Trans. Ind. Electron.; 2009; 56, pp. 900-906. [DOI: https://dx.doi.org/10.1109/TIE.2008.2011621]

35. Sariyildiz, E.; Ohnishi, K. Stability and robustness of disturbance-observer-based motion control systems. IEEE Trans. Ind. Electron.; 2014; 62, pp. 414-422. [DOI: https://dx.doi.org/10.1109/TIE.2014.2327009]

36. Dorf, R.C.; Bishop, R.H. Modern Control Systems; Pearson: New York, NY, USA, 1998.

37. Brunton, S.L.; Proctor, J.L.; Kutz, J.N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl. Acad. Sci. USA; 2016; 113, pp. 3932-3937. [DOI: https://dx.doi.org/10.1073/pnas.1517384113] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27035946]