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

Authors: Minghong She; Pengju Zhao
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

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.


---

## Body

## 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 [1–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–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 Section2 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 Figure1. The device model contains deterministic part P, and random parts are obtained by filtering white noise sequence ηy with the noise filter Hp. Therefore, the device model can be represented as
(1)yk=Pzuk+Hpzηyk.Figure 1
Schematic diagram of closed-loop system.In the formula,u∈ℛnu and y∈ℛny are input and output vectors, respectively, ηy∈ℛny 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:(2)ut−r2t=Czr1t−yt+Hczηut.In the formula, the variance matrix is the zero mean ofΔηu>0. The white noise sequence ηu∈ℛnu is filtered by Hc. Exogenous inputs of r1∈ℛny and r2∈ℛnu 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 r1 and r2. In Figure 1, the system scheme based on the following loop states recognition problem of this study: given the exogenous input of r1k and r2k, and the input and output sequences of uk and yk, 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 signalst∈ℛns and t∈Z are quasi steady processes satisfying the following two conditions:
(3)Est=mst,mst2≤C,∀t∈Z,Rssτ=limN→∞1N∑t=0N−τ−1Est+τstT.In the formula,Rssτ2≤C, ∀τ∈Z. E indicates expected operation. The function Rssτ: Z→ℛns×ns is called the autocorrelation function of st. Similarly, if ωt∈ℛnω is a quasi stationary signal, then the cross-correlation function Rswτ:Z→ℛns×nw of st and ωt is computed as
(4)Rswτ=limN→∞1N∑t=0N−τ−1Est+τωtTIf onlyN data is available, the estimates of the autocorrelation function and cross-correlation function can be calculated as
(5)R^ssτ=1N∑t=0N−τ−1st+τstT,R^swτ=1N∑t=0N−τ−1st+τωtT.At the same time, it is assumed that inN→∞ case, Rssτ and Rswτ 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 modelP shown in Figure 1 is rewritten as the following state space model form:
(6)xpk+1=Apxpk+Bpuk,yk=Cpxpk+Dpuk+vk.In the formula,xp∈ℛnp is the state vector of the device, and the system matrices are Ap∈ℛnp×np, Bp∈ℛnp×nu, Cp∈ℛny×np, and Dp∈ℛny×nu. If the exogenous signal r=r1 or r=r2 is related to the input u, the output y, and the noise v, the cross-correlation functions Ryrτ∈ℛny×ny, Rurτ∈ℛnu×ny, and Rvrτ∈ℛny×ny exist. The covariance Rxprτ∈ℛnp×ny of the states xp and rt is defined, and then the correlation function of the device can be represented as in the state space matrix.
(7)Rxprτ+1=ApRxprτ+BpRurτ,Ryrτ=CpRxprτ+DpRurτ+Rvrτ.It is assumed that the correlation functionsRyrτ and Rurτ are known within a certain interval τ∈τ0,τ1, so that R^yrτ is the estimate of the correlation function Ryrτ, and the definitions of R^urτ and R^vrτ are similar. Then, the correlation function shown in (7) can be rewritten as
(8)Rxprτ+1=ApRxprτ+BpR^urτ,R^yrτ=CpRxprτ+DpR^urτ+R^vrτ.Since the external input instrument signalsr are uncorrelated with noise of ηy and v, all elements of R^vrτ along with N→∞ will gradually tend to 0. Then, the correlation function (8) will converge to
(9)Rxprτ+1=ApRxprτ+BpR^urτ,R^yrτ=CpRxprτ+DpR^urτ.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.


## 2.1. Statement of the Problem and Notation

In the UAV servo closed-loop control system, there is an unknown device, as shown in Figure1. The device model contains deterministic part P, and random parts are obtained by filtering white noise sequence ηy with the noise filter Hp. Therefore, the device model can be represented as
(1)yk=Pzuk+Hpzηyk.Figure 1
Schematic diagram of closed-loop system.In the formula,u∈ℛnu and y∈ℛny are input and output vectors, respectively, ηy∈ℛny 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:(2)ut−r2t=Czr1t−yt+Hczηut.In the formula, the variance matrix is the zero mean ofΔηu>0. The white noise sequence ηu∈ℛnu is filtered by Hc. Exogenous inputs of r1∈ℛny and r2∈ℛnu 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 r1 and r2. In Figure 1, the system scheme based on the following loop states recognition problem of this study: given the exogenous input of r1k and r2k, and the input and output sequences of uk and yk, 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 signalst∈ℛns and t∈Z are quasi steady processes satisfying the following two conditions:
(3)Est=mst,mst2≤C,∀t∈Z,Rssτ=limN→∞1N∑t=0N−τ−1Est+τstT.In the formula,Rssτ2≤C, ∀τ∈Z. E indicates expected operation. The function Rssτ: Z→ℛns×ns is called the autocorrelation function of st. Similarly, if ωt∈ℛnω is a quasi stationary signal, then the cross-correlation function Rswτ:Z→ℛns×nw of st and ωt is computed as
(4)Rswτ=limN→∞1N∑t=0N−τ−1Est+τωtTIf onlyN data is available, the estimates of the autocorrelation function and cross-correlation function can be calculated as
(5)R^ssτ=1N∑t=0N−τ−1st+τstT,R^swτ=1N∑t=0N−τ−1st+τωtT.At the same time, it is assumed that inN→∞ case, Rssτ and Rswτ 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 modelP shown in Figure 1 is rewritten as the following state space model form:
(6)xpk+1=Apxpk+Bpuk,yk=Cpxpk+Dpuk+vk.In the formula,xp∈ℛnp is the state vector of the device, and the system matrices are Ap∈ℛnp×np, Bp∈ℛnp×nu, Cp∈ℛny×np, and Dp∈ℛny×nu. If the exogenous signal r=r1 or r=r2 is related to the input u, the output y, and the noise v, the cross-correlation functions Ryrτ∈ℛny×ny, Rurτ∈ℛnu×ny, and Rvrτ∈ℛny×ny exist. The covariance Rxprτ∈ℛnp×ny of the states xp and rt is defined, and then the correlation function of the device can be represented as in the state space matrix.
(7)Rxprτ+1=ApRxprτ+BpRurτ,Ryrτ=CpRxprτ+DpRurτ+Rvrτ.It is assumed that the correlation functionsRyrτ and Rurτ are known within a certain interval τ∈τ0,τ1, so that R^yrτ is the estimate of the correlation function Ryrτ, and the definitions of R^urτ and R^vrτ are similar. Then, the correlation function shown in (7) can be rewritten as
(8)Rxprτ+1=ApRxprτ+BpR^urτ,R^yrτ=CpRxprτ+DpR^urτ+R^vrτ.Since the external input instrument signalsr are uncorrelated with noise of ηy and v, all elements of R^vrτ along with N→∞ will gradually tend to 0. Then, the correlation function (8) will converge to
(9)Rxprτ+1=ApRxprτ+BpR^urτ,R^yrτ=CpRxprτ+DpR^urτ.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 ofR^τ0∣τi−1yr and R^τ0∣τi−1ur, including i rows and j columns, as defined below:
(10)R^τ0τi−1yr=R^yrτ0R^yrτ1⋯R^yrτj−1R^yrτ1R^yrτ2⋯R^yrτj⋮⋮⋱⋮R^yrτi−1R^yrτi⋯R^yrτj+i−2,R^τ0τi−1ur=R^urτ0R^urτ1⋯R^urτj−1R^urτ1R^urτ2⋯R^urτj⋮⋮⋱⋮R^urτi−1R^urτi⋯R^urτj+i−2.In the formula,R^τ0∣τi−1yr∈ℛnyi×nrj, R^τ0∣τi−1ur∈ℛnui×nrj, each element of R^τ0∣τi−1yr and R^τ0∣τi−1ur is the correlation function data, and i and j are user-defined indexes. According to (9), the above matrix satisfies the relation.
(11)R^τ0∣τi−1yr=ΓRτ0xr+T0∣i−1R^τ0∣τi−1ur.In the formula, vectorRr0xr consists of data from the state cross-correlation function.
(12)Rr0xr=Rxrτ0Rxrτ1⋯Rxrτj−1.The extended observation matrix and the lower triangular block Toeplitz matrix are defined as follows:(13)Γ=CpCpAp⋯CpApi−1T,(14)T0∣i−1=Dp0⋯0CpApDp⋯0⋮⋮⋱⋮CpApi−2BpCpApi−3Bp⋯Dp.By applying the one-step displacement process, the displacement equation corresponding to (11) can be defined as
(15)R^τ1∣τiyr=ΓApRτ0xr+T0∣iR^τ0∣τiur.In the formula, the matrixT0∣i can be added to the left of the T0∣i−1 with a column of zeros, as follows:
(16)T0∣i=0Dp0⋯00CpApDp⋯0⋮⋮⋮⋱⋮0CpApi−2BpCpApi−3Bp⋯Dp.Similarly, a zero line extension is added to the bottom of theR^τ0∣τi−1ur to obtain the representation of R^τ0∣τiur.


### 3.2. Deletion of Input Items Based on Null Space Projection

Although the range of the extended observable matrixΓ and the Toeplitz matrix T0∣i−1 is contained in the R^τ0∣τi−1yr of (11), the range of T0∣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 R^τ0∣τi−1ur can be computed, the projection matrix is first defined as follows:
(17)∏τ0∣τi−1ur=Inrj−R^τ0∣τi−1urTR^τ0∣τt−1urR^τ0∣τi−1urT†R^τ0∣τi−1ur.Theorem 1.
R^τ0∣τi−1yr, R^τ0∣τi−1ur, and T0∣i−1 are defined by the above formula, and ∏τ0∣τkur is defined according to (17) as follows:
(18)∏τ0∣τkur=Inrj−R^τ0∣τkurTR^τ0∣τkurR^τ0∣τkurT†R^τ0∣τkur.
Then,(19)T0∣i−1R^τ0∣τi−1ur∏τ0∣τkur=0.
In the formula,R^τ0∣τkur∈ℛnuk+1×nrj,τi−1≤τk≤τl.Proof.
Forτk=τi−1, (17) is brought into (19). For τi−1<τk, for example, τk=τi−1+s, item T0∣i−1 can be obtained by extending zero column to the right side of the matrix.
(20)T0∣i−1s=Dp0⋯00⋯0CpApDp⋯00⋯0⋮⋮⋱⋮⋮⋮⋮CpApi−2BpCpApi−3Bp⋯Dp0⋯0.
Similarly,R^τ0∣τi−1ur can obtain R^τ0∣τi−1+sur by adding zero line vectors. For all s≥0, T0∣i−1s is a fixed value, and then available
(21)T0∣i−1sR^τ0∣τi−1+sur=T0∣i−1R^τ0∣τi−1ur,T0∣i−1sR^τ0∣τi−1+sur∏τ0∣τkur=0.
The required proof results can be obtained by replacingT0∣i−1sR^τ0∣τi−1+sur with T0∣i−1R^τ0∣τi−1+sur. τk>τl is not allowed here to ensure that (17) satisfies the singularity.


### 3.3. System Dynamics Estimation

The direct result obtained in Theorem1 is that the same projection ∏τ0∣τiur removes the fixed items including T0∣i−1 and T1∣i from R^τ0∣τi−1yr and R^τ0∣τiyr by the right projection of the (11) and (15), and then the following equations can be obtained:
(22)R^τ0∣τi−1yr∏τ0∣τiur=ΓRτ0xr∏τ0∣τiur,(23)R^τ1∣τiyr∏τ0∣τiur=ΓApRτ0xr∏τ0∣τiur.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
(24)R^τ0∣τi−1yr∏τ0∣τiur=UnUs∑n00∑sVnTVsT.(25)Γ^=Un∑n1/2.According to (22) and (23), the minimization problem of the least squares solution is as follows:
(26)J=argminApApΓ†R^τ0∣τi−1yr∏τ0∣τiur−Γ†R^τ1∣τiyr∏τ0∣τiur.A^p=∑n−1/2UnTR^τ1∣τiyr∏τ0∣τiurVn∑n−1/2.The system matrixCp can be calculated by using the ny rows of the extended observation matrix.
(27)C^p=Γ^1:ny,:.Based on the estimation ofAp and Cp, the entire problem becomes linear in the unknown Bp and Dp. Then, the estimated form of the output yk is
(28)y^k=C^pA^px^0+ukT⊗InyvecD^p+∑t=0k−1ukT⊗C^pA^pk−t−1vecB^p.In the formula, the estimates ofBp and Dp can be obtained by solving the linear least squares method.
(29)minBp,Dp1N∑k=0N−1yk−φTkθ22.In the formula, where(30)φTk=ukT⊗Iny∑k=0s−1ukT⊗C^pA^pk−s−1,θ=vecD^pvecB^pT.


### 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, R^τ0∣τi−1yr, ∏τ0∣τiur, and R^τ1∣τiyr∏τ0∣τiur can be computed more efficiently by following LQ decomposition:
(31)R^τ0∣τiurR^τ0∣τi−1yr=L110L21L22Q1TQ2T,(32)R^τ0∣τiurR^τ1∣τiyr=L~110L~21L~22Q~1TQ~2T.Then the outputR^τ0∣τi−1yr can be calculated as
(33)R^τ0∣τi−1yr∏τ0∣τiur=ΓRτ0xr+T0∣i−1L11Q1T=L21Q1T+L22Q2T.Formula (33) is multiplied by Q2 to get L22=ΓRτ0xrQ2 and ImL22=ImΓ. The algorithm then performs L22Q2T and L~22Q~2T at R^τ0∣τi−1yr∏τ0∣τiur and R^τ1∣τi−1yr∏τ0∣τiur, respectively. In addition, the singular value decomposition is performed on L22Q2T 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 Figure2. First, the closed-loop prediction of the parameters of the Γ, L^, and G^ 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.Figure 2
Closed-loop subspace identification process.In Figure2, Γ^ 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.


## 3.1. Block Hankel Matrix and Related Data Estimation Equations

Construct the correlation function to estimate the partitioned Hankel matrices ofR^τ0∣τi−1yr and R^τ0∣τi−1ur, including i rows and j columns, as defined below:
(10)R^τ0τi−1yr=R^yrτ0R^yrτ1⋯R^yrτj−1R^yrτ1R^yrτ2⋯R^yrτj⋮⋮⋱⋮R^yrτi−1R^yrτi⋯R^yrτj+i−2,R^τ0τi−1ur=R^urτ0R^urτ1⋯R^urτj−1R^urτ1R^urτ2⋯R^urτj⋮⋮⋱⋮R^urτi−1R^urτi⋯R^urτj+i−2.In the formula,R^τ0∣τi−1yr∈ℛnyi×nrj, R^τ0∣τi−1ur∈ℛnui×nrj, each element of R^τ0∣τi−1yr and R^τ0∣τi−1ur is the correlation function data, and i and j are user-defined indexes. According to (9), the above matrix satisfies the relation.
(11)R^τ0∣τi−1yr=ΓRτ0xr+T0∣i−1R^τ0∣τi−1ur.In the formula, vectorRr0xr consists of data from the state cross-correlation function.
(12)Rr0xr=Rxrτ0Rxrτ1⋯Rxrτj−1.The extended observation matrix and the lower triangular block Toeplitz matrix are defined as follows:(13)Γ=CpCpAp⋯CpApi−1T,(14)T0∣i−1=Dp0⋯0CpApDp⋯0⋮⋮⋱⋮CpApi−2BpCpApi−3Bp⋯Dp.By applying the one-step displacement process, the displacement equation corresponding to (11) can be defined as
(15)R^τ1∣τiyr=ΓApRτ0xr+T0∣iR^τ0∣τiur.In the formula, the matrixT0∣i can be added to the left of the T0∣i−1 with a column of zeros, as follows:
(16)T0∣i=0Dp0⋯00CpApDp⋯0⋮⋮⋮⋱⋮0CpApi−2BpCpApi−3Bp⋯Dp.Similarly, a zero line extension is added to the bottom of theR^τ0∣τi−1ur to obtain the representation of R^τ0∣τiur.


## 3.2. Deletion of Input Items Based on Null Space Projection

Although the range of the extended observable matrixΓ and the Toeplitz matrix T0∣i−1 is contained in the R^τ0∣τi−1yr of (11), the range of T0∣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 R^τ0∣τi−1ur can be computed, the projection matrix is first defined as follows:
(17)∏τ0∣τi−1ur=Inrj−R^τ0∣τi−1urTR^τ0∣τt−1urR^τ0∣τi−1urT†R^τ0∣τi−1ur.Theorem 1.
R^τ0∣τi−1yr, R^τ0∣τi−1ur, and T0∣i−1 are defined by the above formula, and ∏τ0∣τkur is defined according to (17) as follows:
(18)∏τ0∣τkur=Inrj−R^τ0∣τkurTR^τ0∣τkurR^τ0∣τkurT†R^τ0∣τkur.
Then,(19)T0∣i−1R^τ0∣τi−1ur∏τ0∣τkur=0.
In the formula,R^τ0∣τkur∈ℛnuk+1×nrj,τi−1≤τk≤τl.Proof.
Forτk=τi−1, (17) is brought into (19). For τi−1<τk, for example, τk=τi−1+s, item T0∣i−1 can be obtained by extending zero column to the right side of the matrix.
(20)T0∣i−1s=Dp0⋯00⋯0CpApDp⋯00⋯0⋮⋮⋱⋮⋮⋮⋮CpApi−2BpCpApi−3Bp⋯Dp0⋯0.
Similarly,R^τ0∣τi−1ur can obtain R^τ0∣τi−1+sur by adding zero line vectors. For all s≥0, T0∣i−1s is a fixed value, and then available
(21)T0∣i−1sR^τ0∣τi−1+sur=T0∣i−1R^τ0∣τi−1ur,T0∣i−1sR^τ0∣τi−1+sur∏τ0∣τkur=0.
The required proof results can be obtained by replacingT0∣i−1sR^τ0∣τi−1+sur with T0∣i−1R^τ0∣τi−1+sur. τk>τl is not allowed here to ensure that (17) satisfies the singularity.


## 3.3. System Dynamics Estimation

The direct result obtained in Theorem1 is that the same projection ∏τ0∣τiur removes the fixed items including T0∣i−1 and T1∣i from R^τ0∣τi−1yr and R^τ0∣τiyr by the right projection of the (11) and (15), and then the following equations can be obtained:
(22)R^τ0∣τi−1yr∏τ0∣τiur=ΓRτ0xr∏τ0∣τiur,(23)R^τ1∣τiyr∏τ0∣τiur=ΓApRτ0xr∏τ0∣τiur.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
(24)R^τ0∣τi−1yr∏τ0∣τiur=UnUs∑n00∑sVnTVsT.(25)Γ^=Un∑n1/2.According to (22) and (23), the minimization problem of the least squares solution is as follows:
(26)J=argminApApΓ†R^τ0∣τi−1yr∏τ0∣τiur−Γ†R^τ1∣τiyr∏τ0∣τiur.A^p=∑n−1/2UnTR^τ1∣τiyr∏τ0∣τiurVn∑n−1/2.The system matrixCp can be calculated by using the ny rows of the extended observation matrix.
(27)C^p=Γ^1:ny,:.Based on the estimation ofAp and Cp, the entire problem becomes linear in the unknown Bp and Dp. Then, the estimated form of the output yk is
(28)y^k=C^pA^px^0+ukT⊗InyvecD^p+∑t=0k−1ukT⊗C^pA^pk−t−1vecB^p.In the formula, the estimates ofBp and Dp can be obtained by solving the linear least squares method.
(29)minBp,Dp1N∑k=0N−1yk−φTkθ22.In the formula, where(30)φTk=ukT⊗Iny∑k=0s−1ukT⊗C^pA^pk−s−1,θ=vecD^pvecB^pT.


## 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, R^τ0∣τi−1yr, ∏τ0∣τiur, and R^τ1∣τiyr∏τ0∣τiur can be computed more efficiently by following LQ decomposition:
(31)R^τ0∣τiurR^τ0∣τi−1yr=L110L21L22Q1TQ2T,(32)R^τ0∣τiurR^τ1∣τiyr=L~110L~21L~22Q~1TQ~2T.Then the outputR^τ0∣τi−1yr can be calculated as
(33)R^τ0∣τi−1yr∏τ0∣τiur=ΓRτ0xr+T0∣i−1L11Q1T=L21Q1T+L22Q2T.Formula (33) is multiplied by Q2 to get L22=ΓRτ0xrQ2 and ImL22=ImΓ. The algorithm then performs L22Q2T and L~22Q~2T at R^τ0∣τi−1yr∏τ0∣τiur and R^τ1∣τi−1yr∏τ0∣τiur, respectively. In addition, the singular value decomposition is performed on L22Q2T 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 Figure2. First, the closed-loop prediction of the parameters of the Γ, L^, and G^ 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.Figure 2
Closed-loop subspace identification process.In Figure2, Γ^ 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.


## 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
(34)dXdt=AX+DU,Y=CX+BU.In the formula,X=I¯,ϕ¯ is the inductance matrix, U=V¯ is the input matrix, and Y=I¯ is the output matrix. And there are
(35)A=−1στs+1−σστrLmσLsLrτr+jLmLsLrpωrLmτr−1τr+jpωr,B=1στs0T,C=1,0,D=0.In the formula,σ=1−Lm2/LsLr, τs=Ls/Rs, and τr=Lr/Rr. Lm is excitation inductance, Lr is rotor inductance, Ls is stator inductance, Rs is stator resistance, and Rr 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:Rs=10Ω, Rr=6.58Ω, Lr=0.31H, Ls=0.31H, and ωr=157.08rad/s.The sampling period used is equal to 10−4s. According to the experimental parameters, the following induction motor model parameter matrices can be obtained:
(36)A=−200.32247.02+j3656.165.73−21.23+j314.16,The number of conditions can be calculated as(37)condA=322.47.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 Figures3–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.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 Figures3–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 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 Figures7–9.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.Figures7–9 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 Figure10. 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.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 timet=8s, the motor speed changes from 750r/min to 1200r/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
Convergence of identification results.
Time (s)
Identification error
Rotor resistance/Ω
LMAg/H
9
5.363
8.351
13
2.585
5.262
17
0.822
2.185
21
0.287
0.873Table1 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=8s changes step by step. Here, we start recording the identification error from t=9s. According to the experimental results, the identification error of the algorithm decreases as time goes by, which shows that the algorithm has better convergence.


## 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
(34)dXdt=AX+DU,Y=CX+BU.In the formula,X=I¯,ϕ¯ is the inductance matrix, U=V¯ is the input matrix, and Y=I¯ is the output matrix. And there are
(35)A=−1στs+1−σστrLmσLsLrτr+jLmLsLrpωrLmτr−1τr+jpωr,B=1στs0T,C=1,0,D=0.In the formula,σ=1−Lm2/LsLr, τs=Ls/Rs, and τr=Lr/Rr. Lm is excitation inductance, Lr is rotor inductance, Ls is stator inductance, Rs is stator resistance, and Rr 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:Rs=10Ω, Rr=6.58Ω, Lr=0.31H, Ls=0.31H, and ωr=157.08rad/s.The sampling period used is equal to 10−4s. According to the experimental parameters, the following induction motor model parameter matrices can be obtained:
(36)A=−200.32247.02+j3656.165.73−21.23+j314.16,The number of conditions can be calculated as(37)condA=322.47.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 Figures3–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.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 Figures3–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 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 Figures7–9.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.Figures7–9 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 Figure10. 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.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 timet=8s, the motor speed changes from 750r/min to 1200r/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
Convergence of identification results.
Time (s)
Identification error
Rotor resistance/Ω
LMAg/H
9
5.363
8.351
13
2.585
5.262
17
0.822
2.185
21
0.287
0.873Table1 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=8s changes step by step. Here, we start recording the identification error from t=9s. 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.



---

*Source: 2615203-2018-04-05.xml*