Table of Contents Author Guidelines Submit a Manuscript
International Journal of Aerospace Engineering
Volume 2018, Article ID 2615203, 9 pages
Research Article

Hankel Matrix Correlation Function-Based Subspace Identification Method for UAV Servo System

1School of Automation, Chongqing University, Chongqing, China
2Computer College, Chongqing College of Electronic Engineering, Chongqing, China

Correspondence should be addressed to Minghong She; moc.621@ehs_hm

Received 2 October 2017; Revised 24 December 2017; Accepted 14 January 2018; Published 5 April 2018

Academic Editor: Hikmat Asadov

Copyright © 2018 Minghong She and Pengju Zhao. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


For the identification problem of closed-loop subspace model, we propose a zero space projection method based on the estimation of correlation function to fill the block Hankel matrix of identification model by combining the linear algebra with geometry. By using the same projection of related data in time offset set and LQ decomposition, the multiplication operation of projection is achieved and dynamics estimation of the unknown equipment system model is obtained. Consequently, we have solved the problem of biased estimation caused when the open-loop subspace identification algorithm is applied to the closed-loop identification. A simulation example is given to show the effectiveness of the proposed approach. In final, the practicability of the identification algorithm is verified by hardware test of UAV servo system in real environment.

1. Introduction

In the past twenty years, the subspace model identification (SMI) has received great attention, not only because of its excellent convergence and simple numerical calculation, but also the suitability for the application in the estimation, prediction, and control algorithm. In the early references, most subspace identification methods have open-loop recognition characteristics. Considering stability, safety, and control-oriented identification problems, researchers have been trying to apply these subspace methods to closed-loop identification [14].

The main difficulty in closed-loop system identification is that the correlation between device input and interference leads to a bias in the parameter estimation of the system model [5]. So far, it has developed many subspace identification methods, such as the literature [611], which can obtain consistent estimates and closed-loop data. It is noticed that most subspace system identification methods are based on the input and output data in time domain, and the frequency response methods of some linear time invariant systems are often based on the signal correlation function [12, 13].

Subspace method is extended to frequency response function estimation, and continuous and discrete time models are determined by auxiliary variables [14]. Two frequency statistical properties and subspace convergence analysis methods are proposed in the literature [15]. For a linear closed-loop system, where the external input is independent of the observed noise, the cross correlation function of the output and the external input signals is equal to the cross correlation function of the input and external signals of the dynamic system [16]. By using correlation function sequence as the interface function, the important information has been carried by the correlation function in data compression which has been hidden in the sequence, and by extracting the parameter information of interface function, it provides the basis for parameter identification [17].

The above algorithms can solve the problem of correlation between input and interference of the equipment to some extent. The unbiased parameter estimation can be obtained under arbitrary noise characteristics [18]. But the subspace identification algorithm can solve the relationship between the input and interference at the same time, because the relation between the two related function sequences cannot be completely determined, which leads to the higher dimension of the input matrix of the subspace computation matrix [19].

In order to simplify the calculation and improve the calculation accuracy, we design the input deletion strategy based on the zero space projection on the basis of block Hankel matrix and related data estimation separately [20]. This paper also uses LQ framework to solve the identification process and simplifies the calculation of the algorithm [21]. In this paper, a new subspace identification algorithm based on the estimation of correlation function is proposed, and the unbiased parameter estimation of closed-loop dynamics under linear closed-loop conditions is obtained.

For the EIV (errors-in-variables) model structure, namely, the input and output, both are affected by noise pollution. Chou and Verhaegen put forward a new method of subspace identification [12]. The method eliminates noise effects by regarding the past input/output data as auxiliary variables. Gustafsson [22] changed the steps of the traditional subspace identification method and proposed a new subspace auxiliary variable method (subspace-based identification using instrumental variables (SIVs)). The algorithm presented in the literature [12] was included, and the identification precision was improved after the algorithm was modified. In terms of the algorithm itself, the input being independent of noise assumption is not involved in the two methods; thus, it seems that they can be applied for identification in a closed-loop system. However, this is not the case based on simulation examples; these two types of algorithms used in closed-loop identification do not obtain consensus estimates in some cases.

On the other hand, for closed-loop systems subject to the input and output measurement noises, [23] developed a new closed-loop subspace identification algorithm (SOPIM) by adopting the EIV model structure of SIMPCA and proposed an orthogonal projection approach to avoid identifying the parity space of the feedback controller. In the literature [24], it proposed the use of parity space and principal component analysis (SIMPCA) for EIV identification with colored input excitation can also be applied to closed-loop identification. Instead of preestimating the Markov parameters or eliminating them via noncausal projections, SIMPCA reformulates the SIM problem in parity space. In the process of algorithm implementation, SVD should be applied twice to solve the orthogonal complement space in the two methods. It is difficult to determine the dimensions of the orthogonal complement space when the system order is unknown. The contribution of this paper is to use a zero space projection method based on the estimation of correlation function to fill the block Hankel matrix of identification model by combining the linear algebra with geometry. By using the same projection of related data in time offset set and LQ decomposition, the multiplication operation of projection is achieved and dynamics estimation of the unknown equipment system model is obtained. Simulation shows that this algorithm has higher accuracy, which is another contribution of this paper.

The outline of this paper is as follows: we start in Section 2 with the statement of the problem and notation, and then a state space model based on estimation of correlation function is obtained. In Section 3, we present a block Hankel matrix, related data estimation equations, and deletion of input items based on null space projection for subspace model identification process. In Section 4, both a simulation example and a real hardware verification experiment are presented. We end this paper with our conclusions in Section 5.

2. Closed-Loop Subspace Identification

2.1. Statement of the Problem and Notation

In the UAV servo closed-loop control system, there is an unknown device, as shown in Figure 1. The device model contains deterministic part , and random parts are obtained by filtering white noise sequence with the noise filter . Therefore, the device model can be represented as

Figure 1: Schematic diagram of closed-loop system.

In the formula, and are input and output vectors, respectively, is the zero mean and covariance matrix, and is the white noise.

The system operates under the controller C and shows closed-loop characteristics, so that the closed-loop system can be stably controlled in the whole control trajectory. The resulting control output signals are as follows:

In the formula, the variance matrix is the zero mean of . The white noise sequence is filtered by . Exogenous inputs of and meet persistent excitation conditions and are independent of white noise and . The internal signals and can be represented by and , which are statistically independent quantities of arbitrary colors and known external sequences and . In Figure 1, the system scheme based on the following loop states recognition problem of this study: given the exogenous input of and , and the input and output sequences of and , the goal is to determine the unknown devices in the closed-loop system without the state space model of partial parameter description, as shown in Figure 1.

The signal and are quasi steady processes satisfying the following two conditions:

In the formula, , . E indicates expected operation. The function : is called the autocorrelation function of . Similarly, if is a quasi stationary signal, then the cross-correlation function of and is computed as

If only data is available, the estimates of the autocorrelation function and cross-correlation function can be calculated as

At the same time, it is assumed that in case, and are convergents, respectively.

2.2. State Space Model Based on Estimation of Correlation Function

The subspace identification process adopts the method of calculating the state space matrix to identify the system parameters. Therefore, the first step of the algorithm is to represent the system model into the state space model. Here, the unknown device model shown in Figure 1 is rewritten as the following state space model form:

In the formula, is the state vector of the device, and the system matrices are , , , and . If the exogenous signal or is related to the input , the output , and the noise , the cross-correlation functions , , and exist. The covariance of the states and is defined, and then the correlation function of the device can be represented as in the state space matrix.

It is assumed that the correlation functions and are known within a certain interval , so that is the estimate of the correlation function , and the definitions of and are similar. Then, the correlation function shown in (7) can be rewritten as

Since the external input instrument signals are uncorrelated with noise of and , all elements of along with will gradually tend to 0. Then, the correlation function (8) will converge to

Based on the consistency of correlation function estimation, the direct method of correlation function estimation can be obtained, instead of the original input and output data consistent with the system estimate.

3. The Correlation Function Recognition Based on Subspace Estimation

Through the correlation function estimation and zero space followed by projection, block Hankel matrix of identification framework for filling, we can obtain the scope of the extended observability matrix, which is the basic step of subspace identification methods. Then, an estimate of the dynamics of an unknown device system model can be obtained based on the same projection on the time offset set of the correlation data. Finally, the operation of projection multiplication can be realized by numerical calculation of LQ decomposition.

3.1. Block Hankel Matrix and Related Data Estimation Equations

Construct the correlation function to estimate the partitioned Hankel matrices of and , including rows and columns, as defined below:

In the formula, , , each element of and is the correlation function data, and and are user-defined indexes. According to (9), the above matrix satisfies the relation.

In the formula, vector consists of data from the state cross-correlation function.

The extended observation matrix and the lower triangular block Toeplitz matrix are defined as follows:

By applying the one-step displacement process, the displacement equation corresponding to (11) can be defined as

In the formula, the matrix can be added to the left of the with a column of zeros, as follows:

Similarly, a zero line extension is added to the bottom of the to obtain the representation of .

3.2. Deletion of Input Items Based on Null Space Projection

Although the range of the extended observable matrix and the Toeplitz matrix is contained in the of (11), the range of can be deleted by zero space projection. So the range of the extended observability matrix can be obtained. The main reason is that the null space projection; before feature extraction, the feature of zero space projection extracts the feature projection on nonzero vertex position, and the characteristics of complex single vertex extraction problem were simplified into any nonzero vertex extraction, which is a simplified algorithm of computing process. Since the orthogonal projection matrix of the row space of the partitioned Hankel matrix can be computed, the projection matrix is first defined as follows:

Theorem 1. , , and are defined by the above formula, and is defined according to (17) as follows: Then, In the formula, ,.

Proof. For , (17) is brought into (19). For , for example, , item can be obtained by extending zero column to the right side of the matrix. Similarly, can obtain by adding zero line vectors. For all , is a fixed value, and then available The required proof results can be obtained by replacing with . is not allowed here to ensure that (17) satisfies the singularity.

3.3. System Dynamics Estimation

The direct result obtained in Theorem 1 is that the same projection removes the fixed items including and from and by the right projection of the (11) and (15), and then the following equations can be obtained:

Therefore, from the (22) on the right side of the singular value decomposition, estimation is obtained by the extended observability matrix . The calculation form is

According to (22) and (23), the minimization problem of the least squares solution is as follows:

The system matrix can be calculated by using the rows of the extended observation matrix.

Based on the estimation of and , the entire problem becomes linear in the unknown and . Then, the estimated form of the output is

In the formula, the estimates of and can be obtained by solving the linear least squares method.

In the formula, where

3.4. LQ Computing Framework

According to (20), the expansion matrix, which constitutes a lower triangular matrix, that is, the matrix LQ decomposition, shows that the matrix expansion is used in the matrix of right to satisfy the characteristics. Therefore, select the LQ decomposition of matrix multiplication calculation. By constructing projection matrices and replacing (11) and (15) projections, , , and can be computed more efficiently by following LQ decomposition:

Then the output can be calculated as

Formula (33) is multiplied by to get and . The algorithm then performs and at and , respectively. In addition, the singular value decomposition is performed on to achieve the estimation of the extended observation matrix .

3.5. Subspace Model Identification Process

The proposed closed-loop subspace identification algorithm for the concrete calculation process is shown in Figure 2. First, the closed-loop prediction of the parameters of the , , and matrices is obtained based on the subspace projection method. Then, based on the matrix operation and singular value decomposition process, the values of , , , and are obtained. Finally, the system parameter matrices and are solved by the least squares method.

Figure 2: Closed-loop subspace identification process.

In Figure 2, can be obtained according to (25), and B can be calculated according to the (31-32). is the state vector matrix. , , , and are the state space matrices (for calculation process, see Section 2.3—system dynamic estimation part). and can be calculated according to (14) and (33), respectively.

4. Experimental Analysis

4.1. Model of UAV Servo Induction Motor

Several different induction UAV servo motor models have been described in the literature [6]. The model of the induction motor used in this experiment is

In the formula, is the inductance matrix, is the input matrix, and is the output matrix. And there are

In the formula, , , and . is excitation inductance, is rotor inductance, is stator inductance, is stator resistance, and is rotor resistance.

Here, the Bode amplitude diagram of the Gauss white noise of the induction motor under the Matlab environment is compared with the experimental performance. The hardware settings are the following: memory for the Kingston 8G ddr4-2400GHz processor, i7-6400HQ2.8GHz, and system simulation platform for the selection of win7 ultimate, matlab2014a. The input is generated by filtering a white noise with a fourth-order lowpass Butterworth filter that has a cutoff frequency of 0.8 times the Nyquist frequency and with the zero mean white noise signal disturbance. We select SOPIM and SIMPCA as comparison algorithms, which are shown in [23] and [24].

Experimental parameter settings are the following: , , , , and .The sampling period used is equal to . According to the experimental parameters, the following induction motor model parameter matrices can be obtained:

The number of conditions can be calculated as

It can be seen that the conditional number of the matrix is much larger than 1; therefore, the matrix is ill conditioned. In this case, the instant messaging model is chosen to illustrate the phenomenon of morbidity.

In Figures 35, three kinds of Bode amplitude tracking curves of the induction motor model are given. They are SOPIM estimation curve, SIMPCA estimation curve, and the estimation curve of our algorithm, respectively. In each figure, we perform 50 times Monte Carlo simulations of estimating the system to examine the accuracy of the algorithm in this paper.

Figure 3: SOPIM Bode amplitude tracking curve.
Figure 4: SIMPCA Bode amplitude tracking curve.
Figure 5: The algorithm of this paper Bode amplitude tracking curve.

As shown in Figures 35, the Bode amplitude tracking curve shows that the proposed algorithm outperforms the SOPIM estimation results and the SIMPCA estimation results in amplitude tracking. At the same time, SIMPCA estimation results are better than SOPIM estimation results, mainly because Column weighting for SIMPCA is introduced, while SOPIM uses orthogonal projection only. This is the main reason for zero space projection method based on the estimation of correlation function to fill the block Hankel matrix being adopted in this paper. The Bode amplitude tracking curves among the three algorithms are depicted in Figure 6, which shows the effectiveness of our method.

Figure 6: Comparison of Bode amplitude tracking curves among the three algorithms.

To confirm that the Bode amplitude tracking results can be obtained, the root locus simulations of the SOPIM algorithm, the SIMPCA algorithm, and the proposed estimation algorithm are performed. The curves obtained are shown in Figures 79.

Figure 7: SOPIM simulation and comparison of root loci.
Figure 8: SIMPCA simulation and comparison of root loci.
Figure 9: The algorithm of this paper simulation and comparison of root loci.

Figures 79 illustrate that the red “+” symbol is the real trajectory of the model, and the real trajectory is concentrated on the right side of the chart. In the comparative test of three algorithms, the SOPIM algorithm tracks the most points in the graph, which shows that the locus distribution is dispersed. There are relatively few points near the true trajectory of the “+” symbol, indicating the poor accuracy of trajectory tracking. The trajectory distribution of SIMPCA algorithm in Figure 8 is relatively concentrated. The distribution of our algorithm is the most concentrated, and the main results in the distribution of real track near the “+” sign, which indicates that this algorithm with respect to the root locus simulation has higher tracking accuracy compared with the SOPIM algorithm and SIMPCA algorithm, as shown in Figure 9.

4.2. Hardware Verification Experiment

The hardware test platform is under the condition of single-phase inverter unit and controlled rectifier of model PM201CL1A061. Platform identification algorithm is DSP27334 main control chip. The production company is TI. The simulation time connection circuit diagram is shown in Figure 10. The following parameters are P_ MW motor rated power, rated power factor cos_ , rated voltage U_ kV, rated speed of n_ r/min D_a1 = 8500 mm, the outside diameter of the stator core, stator core diameter D_i1 = 7500 mm, slot number , the stator slot width b_ mm, damping bar diameter d_ mm, damping number n_, and air cooling mode.

Figure 10: Hardware test connection diagram.

Because the model belongs to the serial working mode, the proposed algorithm can realize simultaneous identification of magnetizing inductance and rotor resistance with more relative storage. But with the development of science and technology of hard disk, the problem can be ignored. At the time , the motor speed changes from to step by step. Moreover, the rotor resistance synchronous identification error and excitation inductance identification results at this time are shown in Table 1.

Table 1: Convergence of identification results.

Table 1 shows the changes of motor speed with the step changing. The identification results for the algorithm of the motor rotor resistance identification error and the synchronous excitation inductance error change in a passage of time, because the speed of order at changes step by step. Here, we start recording the identification error from . According to the experimental results, the identification error of the algorithm decreases as time goes by, which shows that the algorithm has better convergence.

5. Concluding Remarks

This paper proposed a kind of zero space projection based on correlation function estimation of closed-loop subspace identification method. Through the correlation function estimation and zero space followed by projection, block Hankel matrix of identification framework was filled by the numerical calculation and LQ decomposition to achieve multiplication projection, get unknown device estimation of dynamic system model. The experimental results verify the effectiveness of the proposed method. The next step will focus on the closed-loop system combined with the subspace identification method for the control problem of UAV.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This study is supported by the Science and Technology Research Program of Chongqing Municipal Education Commission (Grant no. KJ1729403); Chongqing Science and Technology Communication and Popularization Project (no. cstc2017kp-ysczA0021); National Natural Science Foundation of China (no. 61403055); the Intelligent Robot Technology Research Center of Chongqing College of Electronic Engineering (no. XJPT201705); the Fundamental Research Funds for the Central Universities (Grant no. CDJXS12 17 11 01).


  1. S. Joe Qin and L. Ljung, “Closed-loop subspace identification with innovation estimation,” IFAC Proceedings, vol. 36, no. 16, pp. 861–866, 2003. View at Publisher · View at Google Scholar
  2. G. Mercère and M. Lovera, “Convergence analysis of instrumental variable recursive subspace identification algorithms,” Automatica, vol. 43, no. 8, pp. 1377–1386, 2007. View at Publisher · View at Google Scholar · View at Scopus
  3. M. Verhaegen, “Application of a subspace model identification technique to identify LTI systems operating in closed-loop,” Automatica, vol. 29, no. 4, pp. 1027–1040, 1993. View at Publisher · View at Google Scholar · View at Scopus
  4. J. Wang and S. J. Qin, “A new subspace identification approach based on principal component analysis,” Journal of Process Control, vol. 12, no. 8, pp. 841–855, 2002. View at Publisher · View at Google Scholar · View at Scopus
  5. G. van der Veen, J. W. van Wingerden, M. Bergamasco, M. Lovera, and M. Verhaegen, “Closed-loop subspace identification methods: an overview,” vol. 7, no. 10, pp. 1339–1358, 2013. View at Publisher · View at Google Scholar · View at Scopus
  6. I. Khan, D. Shan, and Q. Li, “Continuous modal parameter identification of a cable-stayed bridge based on Robustious decomposition and covariance-driven stochastic subspace identification,” Iranian Journal of Science and Technology, Transactions of Civil Engineering, vol. 40, no. 1, pp. 11–22, 2016. View at Publisher · View at Google Scholar · View at Scopus
  7. K. Jalaleddini and R. E. Kearney, “Subspace identification of SISO Hammerstein systems: application to stretch reflex identification,” IEEE Transactions on Biomedical Engineering, vol. 60, no. 10, pp. 2725–2734, 2013. View at Publisher · View at Google Scholar · View at Scopus
  8. A. Naitali and F. Giri, “Persistent excitation by deterministic signals for subspace parametric identification of MISO Hammerstein systems,” IEEE Transactions on Automatic Control, vol. 61, no. 1, pp. 258–263, 2016. View at Publisher · View at Google Scholar · View at Scopus
  9. C. S. Rao and M. Chidambaram, “Experimental application of subspace model identification of an unstable system,” International Journal of Advances in Engineering Sciences and Applied Mathematics, vol. 7, no. 1-2, pp. 70–76, 2015. View at Publisher · View at Google Scholar
  10. T. Jiang, H. Yuan, H. Jia, N. Zhou, and F. Li, “Stochastic subspace identification-based approach for tracking inter-area oscillatory modes in bulk power system utilising synchrophasor measurements,” IET Generation, Transmission & Distribution, vol. 9, no. 15, pp. 2409–2418, 2015. View at Publisher · View at Google Scholar · View at Scopus
  11. Z. Jianyuan, L. Xingfei, and T. Ziling, “The orthogonal decomposition of the recursive subspace identification method of closed loop based on,” Control and Decision, vol. 30, no. 3, pp. 490–494, 2015. View at Google Scholar
  12. C. T. Chou and M. Verhaegen, “Subspace algorithms for the identification of multivariable dynamic errors-in-variables models,” Automatica, vol. 33, no. 10, pp. 1857–1869, 1997. View at Publisher · View at Google Scholar
  13. T. Katayama and H. Tanaka, “An approach to closed-loop subspace identification by orthogonal decomposition,” Automatica, vol. 43, no. 9, pp. 1623–1630, 2007. View at Publisher · View at Google Scholar · View at Scopus
  14. L. Xinming, G. Xianwen, and L. Xinzhe, “Improved pressure identification method based on manifold subspace,” Chinese Journal of Instrumentation, vol. 35, no. 5, pp. 1079–1085, 2014. View at Google Scholar
  15. A. G. C. Junior, J. A. Riul, and P. H. M. Montenegro, “Application of the subspace identification method using the N4SID technique for a robotic manipulator,” IEEE Latin America Transactions, vol. 14, no. 4, pp. 1588–1593, 2016. View at Publisher · View at Google Scholar · View at Scopus
  16. Y. Xie, P. Liu, and G. P. Cai, “Modal parameter identification of flexible spacecraft using the covariance-driven stochastic subspace identification (SSI-COV) method,” Acta Mechanica Sinica, vol. 32, no. 4, pp. 710–719, 2016. View at Publisher · View at Google Scholar · View at Scopus
  17. R. N. Cabrera, V. M. A. Martínez, and M. A. Medina, “Analysis of subspace identification methods based on the estimation of the system matrices,” IEEE Latin America Transactions, vol. 13, no. 4, pp. 1068–1076, 2015. View at Publisher · View at Google Scholar · View at Scopus
  18. X. Wu, B. Huang, L. Wang, and J. Zhang, “GPU-based parallel design of the hyperspectral signal subspace identification by minimum error (HySime),” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 9, no. 9, pp. 4400–4406, 2016. View at Publisher · View at Google Scholar · View at Scopus
  19. S. A. Nezam Sarmadi and V. Venkatasubramanian, “Electromechanical mode estimation using recursive adaptive stochastic subspace identification,” IEEE Transactions on Power Systems, vol. 29, no. 1, pp. 349–358, 2014. View at Publisher · View at Google Scholar · View at Scopus
  20. D. Wang, F. Ding, and L. Ximei, “Least squares algorithm for an input nonlinear system with a dynamic subspace state space model,” Nonlinear Dynamics, vol. 75, no. 1-2, pp. 49–61, 2014. View at Publisher · View at Google Scholar · View at Scopus
  21. J. M. Ni, C. Shen, and F. Liu, “Estimation of the electromechanical characteristics of power systems based on a revised stochastic subspace method and the stabilization diagram,” Science China Technological Sciences, vol. 55, no. 6, pp. 1677–1687, 2012. View at Publisher · View at Google Scholar · View at Scopus
  22. T. Gustafsson, “Subspace identification using instrumental variable techniques,” Automatica, vol. 37, no. 12, pp. 2005–2010, 2001. View at Publisher · View at Google Scholar · View at Scopus
  23. B. Huang, S. X. Ding, and S. J. Qin, “Closed-loop subspace identification: an orthogonal projection approach,” Journal of Process Control, vol. 15, no. 1, pp. 53–66, 2005. View at Publisher · View at Google Scholar · View at Scopus
  24. J. Wang and S. J. Qin, “Closed-loop subspace identification using the parity space,” Automatica, vol. 42, no. 2, pp. 315–320, 2006. View at Publisher · View at Google Scholar · View at Scopus