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

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.


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 [1][2][3][4].
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 [6][7][8][9][10][11], 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.

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 P, and random parts are obtained by filtering white noise sequence η y with the noise filter H p .Therefore, the device model can be represented as In the formula, u ∈ ℛ n u and y ∈ ℛ n y are input and output vectors, respectively, η y ∈ ℛ n y is the zero mean and covariance matrix, and Δ η y > 0 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 Δ η u > 0. The white noise sequence η u ∈ ℛ n u is filtered by H c .Exogenous inputs of r 1 ∈ ℛ n y and r 2 ∈ ℛ n u meet persistent excitation conditions and are independent of white noise η u and η y .The internal signals ω and v can be represented by η u and η y , which are statistically independent quantities of arbitrary colors and known external sequences r 1 and r 2 .In Figure 1, the system scheme based on the following loop states recognition problem of this study: given the exogenous input of r 1 k and r 2 k , and the input and output sequences of u k and y k , 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 s t ∈ ℛ n s and t ∈ Z are quasi steady processes satisfying the following two conditions: In the formula, R ss τ 2 ≤ C, ∀τ ∈ Z. E indicates expected operation.The function R ss τ : Z → ℛ n s ×n s is called the autocorrelation function of s t .Similarly, if ω t ∈ ℛ n ω is a quasi stationary signal, then the crosscorrelation function R sw τ : Z → ℛ n s ×n w of s t and ω t is computed as If only N 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 N → ∞ case, R ss τ and R sw τ are convergents, respectively.

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 P shown in Figure 1 is rewritten as the following state space model form: In the formula, x p ∈ ℛ n p is the state vector of the device, and the system matrices are A p ∈ ℛ n p ×n p , B p ∈ ℛ n p ×n u , C p ∈ ℛ n y ×n p , and D p ∈ ℛ n y ×n u .If the exogenous signal r = r 1 or r = r 2 is related to the input u, the output y, and the noise v, the cross-correlation functions R yr τ ∈ ℛ n y ×n y , R ur τ ∈ ℛ n u ×n y , and R vr τ ∈ ℛ n y ×n y exist.The covariance R x p r τ ∈ ℛ n p ×n y of the states x p and r t 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 R yr τ and R ur τ are known within a certain interval τ ∈ τ 0 , τ 1 , so that Ryr τ is the estimate of the correlation function R yr τ , and the definitions of Rur τ and Rvr τ are similar.Then, the correlation function shown in (7) can be rewritten as Since the external input instrument signals r are uncorrelated with noise of η y and v, all elements of Rvr τ along with N → ∞ 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.

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.

Block Hankel Matrix and Related Data Estimation
Equations.Construct the correlation function to estimate the partitioned Hankel matrices of Ryr τ 0 |τ i−1 and Rur τ 0 |τ i−1 , including i rows and j columns, as defined below: In the formula, Ryr τ 0 |τ i−1 ∈ ℛ n y i×n r j , Rur τ 0 |τ i−1 ∈ ℛ n u i×n r j , each element of Ryr τ 0 |τ i−1 and Rur τ 0 |τ i−1 is the correlation function data, and i and j are user-defined indexes.According to (9), the above matrix satisfies the relation.

Ryr
In the formula, vector R xr r0 consists of data from the state cross-correlation function.
Figure 1: Schematic diagram of closed-loop system.

International Journal of Aerospace Engineering
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 T 0|i can be added to the left of the T 0|i−1 with a column of zeros, as follows: Similarly, a zero line extension is added to the bottom of the Rur τ 0 |τ i−1 to obtain the representation of Rur τ 0 |τ i .

Deletion of Input Items Based on Null Space Projection.
Although the range of the extended observable matrix Γ and the Toeplitz matrix T 0|i−1 is contained in the Ryr τ 0 |τ i−1 of (11), the range of T 0|i−1 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 Rur τ 0 |τ i−1 can be computed, the projection matrix is first defined as follows: , and T 0|i−1 are defined by the above formula, and ∏ ur τ 0 |τ k is defined according to (17) as follows: In the formula, Proof.For τ k = τ i−1 , ( 17) is brought into (19).For τ i−1 < τ k , for example, τ k = τ i−1+s , item T 0|i−1 can be obtained by extending zero column to the right side of the matrix.
Similarly, Rur τ 0 |τ i−1 can obtain Rur τ 0 |τ i−1+s by adding zero line vectors.For all s ≥ 0, T s 0|i−1 is a fixed value, and then available The required proof results can be obtained by replacing T s 0|i−1

System Dynamics Estimation.
The direct result obtained in Theorem 1 is that the same projection ∏ ur τ 0 |τ i removes the fixed items including T 0|i−1 and T 1|i from Ryr τ 0 |τ i−1 and Ryr τ 0 |τ i by the right projection of the ( 11) and ( 15), and then the following equations can be obtained: Ryr Ryr 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 Ryr According to (22) and ( 23), the minimization problem of the least squares solution is as follows: 4 International Journal of Aerospace Engineering

J = arg min
The system matrix C p can be calculated by using the n y rows of the extended observation matrix.
Based on the estimation of A p and C p , the entire problem becomes linear in the unknown B p and D p .Then, the estimated form of the output y k is In the formula, the estimates of B p and D p can be obtained by solving the linear least squares method.

29
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, Ryr τ 0 |τ i−1 , ∏ ur τ 0 |τ i , and Ryr τ 1 |τ i ∏ ur τ 0 |τ i can be computed more efficiently by following LQ decomposition: Rur Then the output Ryr τ 0 |τ i−1 can be calculated as Ryr Formula (33) is multiplied by In addition, the singular value decomposition is performed on L 22 Q T 2 to achieve the estimation of the extended observation matrix Γ.

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 Γ, L, and Ĝ matrices is obtained based on the subspace projection method.Then, based on the matrix operation and singular value decomposition process, the values of X, A, C, and R are obtained.Finally, the system parameter matrices B and D are solved by the least squares method.
In Figure 2, Γ can be obtained according to (25), and B can be calculated according to the (31-32).X is the state vector matrix.A, B, C, and D are the state space matrices (for calculation process, see Section 2.3-system dynamic estimation part).T and R can be calculated according to (14) and (33), respectively.

Experimental Analysis
In the formula, σ = 1 − L 2 m /L s L r , τ s = L s /R s , and τ r = L r /R r .L m is excitation inductance, L r is rotor inductance, L s is stator inductance, R s is stator resistance, and R r 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: R s = 10 Ω, R r = 6 58 Ω, L r = 0 31 H, L s = 0 31 H, and ω r = 157 08 rad/s.The sampling period used is equal to 10 The number of conditions can be calculated as cond A = 322 47 37 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 3-5, 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.
As shown in Figures 3-5, 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 International Journal of Aerospace Engineering 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_N = 250 MW motor rated power, rated power factor cos_ N = 0 9, rated voltage U_N = 15 75 kV, rated speed of n_N = 250 r/min D_a1 = 8500 mm, the outside diameter of the stator core, stator core diameter D_i1 = 7500 mm, slot number Z = 360, the stator slot width b_s = 24 7 mm, damping bar diameter d_B = 25 mm, damping number n_B = 7, and air cooling mode.
Because the model belongs to the serial working mode, the proposed algorithm can realize simultaneous identification of magnetizing 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 t = 8 s, the motor speed changes from 750 r/min to 1200 r/min 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 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 t = 8 s changes step by step.Here, we start recording the identification error from t = 9 s.According to the experimental results, the identification error of the algorithm decreases as time goes by, which shows that the algorithm has better convergence.

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.

Table 1 :
Convergence of identification results.