Optimal input design for multi UAVs formation anomaly detection

ISA Transactions 91 (2019) 157–165

Contents lists available at ScienceDirect

ISA Transactions journal homepage: www.elsevier.com/locate/isatrans

Research article

Optimal input design for multi UAVs formation anomaly detection Hong Wang-jian Universidad de Monterrey, Ave. Morones Prieto 4500 Pte., Jesú,s M. Garza, San Pedro Garza García, N.L., Mexico School of Engineering and Sciences, Tecnologico de Monterrey, Monterrey, Mexico

highlights • Dynamic programming strategy is proposed to consider a hypothesis testing problem in anomaly detection, further the main idea of using dynamic programming strategy is to compute the risk of decision corresponding to each choice policy.

• Optimal input signal design for statistical description of noise is given more specifically. • The variances of unknown parameters are derived corresponding to least squares method and sparse optimization algorithm respectively.

article

info

Article history: Received 6 September 2018 Received in revised form 24 December 2018 Accepted 24 January 2019 Available online 4 February 2019 Keywords: Multi UAVs formation Anomaly detection Optimal input design Dynamic programming

a b s t r a c t As the input signal must be informative enough, so that the resulting dataset is enough informative to excite the identification experiment for multi UAVs formation anomaly detection. Based on our previous work on multi UAVs formation anomaly detection, the optimal input signals are designed for two identification strategies, i.e. least squares estimation and improved sparse estimation. Using the variance of the asymptotic distribution corresponding to the unknown parameters, the common trace operation is chosen to construct one numerical optimization problem, whose solution is corresponded to the optimal power spectral. After giving the detailed minimization process, we see that the power spectral corresponding to the optimal input signal is a constant. In addition, for the sake of completeness, one dynamic programming technique in multi UAVs formation anomaly detection is added to complete our early research. Finally, one numerical example illustrates the effectiveness of our proposed theories. © 2019 ISA. Published by Elsevier Ltd. All rights reserved.

1. Introduction Unmanned aerial vehicles (UAVs) are used as an efficient detection, attack and defense devices, and gained great superiority in our military field in these recent years. To improve their detection and attack efficiency, and extend the applied range, one novel concept of multi UAVs formation is proposed in America. Multi UAVs formation means that many UAVs fly in a special flight mode throughout the entire task completion. The aim of the formation flight is to maintain the relative position and attitude corresponding to each UAV, thus expanding the vision field, achieving the full ground or controlling a target, and improving the integration efficiency. When formation task or battlefield environment change with time increases, the entire multi UAVs formation will also adjust adaptively. Research regarding on multi UAVs formation has primarily concentrated on formation control, formation aerodynamic coupling and formation communication positioning systems [1]. These three categories can be further divided into trajectory planning, trajectory maintenance, trajectory control and autonomous reconfiguration. However little research has been reported on the E-mail address: [email protected]. https://doi.org/10.1016/j.isatra.2019.01.027 0019-0578/© 2019 ISA. Published by Elsevier Ltd. All rights reserved.

detection of multi UAVs formation anomalies. To the best of the author’s knowledge, the problem of anomaly detection for multi UAVs formation flight is firstly put forth by the author [2]. This problem lies in the judgment of the behavior and characteristics of a certain number of subsystems, which differ from most other subsystems in large scale systems. The idea of anomaly detection is derived from statistical signal processing, used to solve the decision or hypothesis testing problem. After applying some probability inequalities to compute probability levels for each hypothesis, the hypothesis can be determined by comparing the respective probability level. Some statistical hypothesis testing methods are introduced in [3], including the cluster method, classification method, nearest neighbor method, information entropy method and the spectral analysis method, etc. The cluster method is used for image compression sensing [4], further J H Braslavsky [5] considers the credibility of applying minimum variance control in Gaussian communication channels and the information entropy method is applied to determine the probability of failure with the application of that minimum variance control strategy. In O C Imer and S Yukgel [6], the spectral analysis method is proposed to judge the hypothesis testing probability regarding to system stability when data drop information is added to the communication channel. In V Gupta and D Spanos [7], the optimal control of dynamical

158

H. Wang-jian / ISA Transactions 91 (2019) 157–165

systems over unreliable communication links is considered, it is pointed out that when the number of hypothesis is greater than two, it becomes difficult to solve multi hypothesis testing problem, and the number of hypothesis will increase with the number of systems increase [8]. Therefore the process of anomaly detection is proposed here to use multi hypothesis testing to determine the statistical probability, which is employed to judge whether one certain of system is abnormal. The anomaly detection more widely exists in the fault identification and credibility. Here firstly let us give an explicit description of this problem of multi UAVs formation anomaly detection. Assume the formation flight consists of N UAVs, and these N UAVs show similar dynamics. Much data is collected through ground control station to supervise these N UAVs and enhance the flight safety of formation. Classical methods used to detect anomaly need the knowledge about the number of the abnormal UAVs k, then from statistical signal processing theory [9], anomaly detection is achieved by one multi hypothesis testing in which k UAVs are chosen from a set of N UAVs, leading to a total of CNk possible hypotheses. Through comparing observed value and its predictive value, statistical testing operation is used to judge the hypothesis. But here in order to avoid classical statistical knowledge, we propose to apply the system identification method to achieve the goal without any knowledge of these N UAVs. This proposed system identification method uses only the input–output data to achieve the anomaly detection of multi UAVs formation. From flight dynamic analysis and that each nonlinear function can be approximated by one appropriate linear function, one unknown linear relation is determined to describe the measurable quantities of the N UAVs. After chosen a certain period of time as the total number of observations, the linear input–output relations corresponding to N UAVs are obtained, as nonlinear forms can be rewritten as their corresponding linear regression models. These unknown linear relations show the relationship among the angle of attack, velocity, dynamic pressure, elevator deflection and the stabilizer deflection angle [10]. The anomaly detection problem means that the difference of each UAV corresponds to the difference of one unknown parameter in the unknown linear relation. By the use of this difference, the anomaly detection problem is changed to identify these N unknown parameters using only input–output data, as each unknown parameter corresponds to each UAV. To solve this identification problem, the common maximum likelihood estimation is used to identify these N unknown parameters. Generally when the number of abnormal UAVs is known as k in prior, then these k parameters with large errors among N parameters are regarded as abnormal states, and other N − k states are normal states. The above anomaly detection process is similar to data clustering. Due to that N UAVs have similar dynamics, then all N unknown parameters must be equal or nearly equal. After conducting cluster analysis with respect to these identified N parameters, the parameters with neglected errors in the admissible range are clustered as first class, while another second class includes the parameters with large errors. This second class can be regarded as the abnormal states of UAVs. As the concept of anomaly detection in multi UAVs formation is proposed originally by the author, the interested reader can refer to two main published papers [2,11] and one recent paper [12]. Paper [2] constructs one sparse optimization problem to identify a linear unknown parameter vector, further the optimum necessary condition and the fast gradient algorithm are respectively proposed to estimate the optimum values. Furthermore in [11], one unknown nonlinear relation is established to describe the measurable quantities, then after expanding it by the orthonormal basis functions, one bias compensated least squares method is applied to solve the unknown parameter vector. And our recent paper [12] puts forth one projection algorithm with dead zone to identify these unknown

parameters in the presence of unknown but bounded noise. Also in that paper, dynamic programming technique is introduced to balance the desire for lower present cost with the undesirability of high future cost in determining the anomaly detector. The advantage of introducing dynamic programming technique is that the cost of collecting new observations and the higher probability of accepting the wrong hypothesis can be compensated. As the problem of anomaly detection for multi UAVs is not considered yet in research field, except above mentioned three papers, so it is very necessary to continue this topic deeply based on our above three papers. Here in this paper, we continue the research on the anomaly detection in multi UAVs formation deeply on the basis of our three papers, where experimental data are usually obtained by means of measurement procedures that are known to be affected by external noise. It is well known that designing one identification experiment includes many choices. Firstly it is essential to obtain some information about the considered system, such as input and output signal. But in some cases, it is not clear which signal can be used as input signal. This tells that one has to decide which signal should be used to excite the considered system during the whole identification or estimation experiment. Generally the chosen input would be informative enough, so that the collected dataset is enough informative to allow a good identification procedure. From the classical system identification theory point of view, the problem of designing optimal input signal is transformed into one numerical optimization problem, whose performance function is constructed by one engineer. And the decision variable is chosen as the cross correlation function or power spectral corresponding to the optimal input signal. Then through solving this numerical optimization problem, the optimal power spectral is derived. Most of the results available from references are based on a statistical description of the uncertainty or noise affecting the measurement data. As the least squares method and maximum likelihood method are suited for identification accuracy under one condition that the considered noise may be a zero mean random signal, and the noise is statistically independent on the input. But zero mean random white noise is an ideal case, it will not exist in engineering. In additional, deriving statistical properties of noise is often very difficult in practice as it is usually not possible to measure noise directly. A worthwhile alternative to the statistical description is the so-called bounded error characterization where measurement noises are assumed to be unknown but bounded. Such a description may be more appropriate in many cases where a priori statistical information is not available. Based on the unknown but bounded description of the measurement noise, one projection algorithm with dead zone is proposed to deal with this general case in the presence of bounded noise. The first case of statistical description of noise is considered in [2] and [11], where least squares method and sparse optimization algorithm are applied to identify the unknown parameters. Due to the performance function in designing the optimal input signal is dependent on the variance of unknown parameters, so a deep asymptotic analysis of unknown parameters is derived explicitly. According to the asymptotic analysis or the variance of asymptotic distribution, the common trace operation is chosen as our performance function. Combining this trace operation and one unknown decision variable, we construct one numerical optimization problem to design the optimal input signal, whose optimal power spectral can be obtained. From these results, we see that whatever least squares method, sparse optimization algorithm or any other identification algorithms in the statistical description of noise, the optimal input signal is determined with a constant power spectrum. Moreover for the sake of completeness, one dynamic programming technique in anomaly detection is added to complete our paper [12], where we only mentioned dynamic programming technique preliminary.

H. Wang-jian / ISA Transactions 91 (2019) 157–165

159

velocity, dynamic pressure, thrust acting on UAV and other flight operation qualities. Then we write it as follows. y (t ) = F (u (t )) + e (t )

Fig. 1. Structure of multi UAVs formation flight system.

(1)

where y (t ) is a vector of six accelerations(three linear accelerations at the aircraft center-of-gravity and three rotational accelerations), u (t ) is a vector of inputs, e (t ) is the white noise from flight condition variables and turbulence, t is the time instant. The nonlinear relation F (u (t )) describes the dependence of the accelerations on vehicle mass, rotational inertia, forces, and moments acting on the vehicle. The components of the input include a vector of flight surface angles s, a vector of velocities and rotational rates v , dynamic pressure q and thrust τ . The components of the input u (t ) and accelerations y (t ) are available in the aircraft avionics systems and are collected as part of flight operations quality assurance datasets. For a single flight, we might have a number of sensors and their recorded values during flight. For example, an airspeed sensor for a single flight would contain the speed of the aircraft during the entire flight. The dataset collected during the flight can be represented as

{u (t ) , y (t )}M t =1

(2)

where M is the number of samples, u (t ) and y (t ) are the input– output at time sample t respectively. As we consider the multi UAVs, we extend the sample database of single UAV to the sample database of multi UAVs. The database consists of a collection of these flights. That is. Fig. 2. The simplified structure.

So here we give a detailed idea of dynamic programming technique that is proposed to consider a hypothesis testing problem in anomaly detection. In this dynamic programming technique, a decision maker can make observation, at a risk of cost each, relating to two hypotheses. Given a new measurement data, one can either accept one of hypotheses or delay the decision, pay the risk of cost, and obtain a measurement data. The reason of introducing dynamic programming into anomaly detection is that the risk of cost is more actual than theory, and the optimality of the decision policy can be guaranteed. This paper is organized as follows. In Section 2, the problem of anomaly detection for multi UAVs formulation is formulated. In Section 3, dynamic programming strategy is proposed to consider a hypothesis testing problem in anomaly detection, further the main idea of using dynamic programming strategy is to compute the risk of decision corresponding to each choice policy, and determine those abnormal UAVs with lower risk. Optimal input signal design for statistical description of noise is given more specifically in Section 4, where the variances of unknown parameters are derived corresponding to least squares method and sparse optimization algorithm respectively. In Section 5, one numerical example illustrates the effectiveness of our proposed theories. Section 6 ends the paper with final conclusion and points our other ongoing research.

{

{ui (t ) , yi (t )}M t =1

}N

(3)

i=1

where i indexes the flight in the database, N is the total number of UAVs, ui (t ) and yi (t ) obey the nonlinear relation (1). That is yi (t ) = Fi (ui (t )) + ei (t ) , i = 1 . . . N , t = 1 . . . M

(4)

As the nonlinear relation is unknown, it is difficult to identify it. The commonly used identification method is to expand the nonlinear relation as one basic functions expansion form. Fi (ui (t )) =

n ∑

θi0,l ϕi,l (ui (t ))

(5)

l=1

} N ,n

where ϕi,l : R → R i=1,l=1 are a-priori chosen functions belonging to the canonical polynomial basis in function space and the most common form of basic functions are the radial basis function kernel in neural networks. From mathematical approximation theory, when n takes large, the basis functions expansion form can approximate the original nonlinear relation very well and the approximation { } accuracy can be attained at any specified requirements. θi0,l are the unknown parameters, superscript 0 denotes the true parameters. The basis functions in (5) depend only on input variable ui (t ), but not on output variable yi (t ). Substituting (5) into (4), we obtain:

{

yi (t ) =

n ∑

θi0,l ϕi,l (ui (t )) + ei (t ) , i = 1 . . . N , t = 1 . . . M

(6)

l=1

2. Problem description

To rewrite (6) as the matrix form, construct the following vectors

Assume that formation flight system in Fig. 1 and its simplified structure in Fig. 2 include N UAVs, and all UAVs exhibit similar dynamics. Using the stress analysis process of flight control systems [13], one nonlinear dynamical model of N UAVs is established. This nonlinear relation indicates that when force is applied to one UAV, the acceleration of each UAV is a function with respect to some certain variables, for example aircraft control surface position,

[ ϕi (ui (t )) = ϕi,1 (ui (t )) ϕi,2 (ui (t )) [ ]T θi0 = θi0,1 θi0,2 · · · θi0,n

···

]T ϕi,n (ui (t ))

Then relation (6) can be formulated as: yi (t ) = ϕiT (ui (t )) θi0 + ei (t )

(7)

Based on system identification theory, we see that there exists one parameter vector θi0 such that ui (t ) and yi (t ) have the linear

160

H. Wang-jian / ISA Transactions 91 (2019) 157–165

regression form. The parameterized model with respect to Eq. (7) is that. yi (t ) = ϕiT (ui (t )) θi + ei (t ) = ϕi (t )T θi + ei (t ) , i = 1, 2 . . . N (8) where t is time instant, i denotes one UAV, yi (t ) ∈ R is the observed output at time instant t, ϕi (t ) is one regression vector at time instant t. θi is an unknown parameter, ei (t ) is the external noise. Consider each UAV, {yi (t ) , ϕi (t )}M t =1 denotes one dataset collected by ground control station at time instant t, and M is the total number of time instant t. For one single flight, a number of sensors record values during flight, for example an airspeed sensor for a single flight contains the speed of the aircraft over the entire flight. These recorded values are all included in dataset {yi (t ) , ϕi (t )}M t =1 . Among these N UAVs, for i = 1, 2 . . . N, as all UAVs have similar dynamics and each UAV corresponds to each parameter θi , then all of these N parameters must be equal to each other, i.e. it holds that.

θ1 = θ2 = · · · = θN = θ0

(9)

where θ0 is a true or nominated parameter value, if Eq. (2) holds, then all N UAVs are in normal states. On the contrary, if one parameter θj deviates from the nominated parameter θ0 , then UAV j is deemed as abnormal state. Due to the number of abnormal UAVs k is not known in prior, the total number of possible choice of choosing k UAVs from a set of N UAVs leads to CNk . N!

CNk =

(N − k)!k!

From above description about anomaly detection problem, the main contribution of applying dynamic programming technique is to compute the risk of decision corresponding to each choice policy by using probability inequalities, and determine which k UAVs are in abnormal states with lower risk. From Eqs. (8) and (9) given above, we see that after all unknown parameter vectors θi , i = 1, 2 . . . N are identified, that special nominated parameter value θ0 can be given as the average value of these N parameters.

θ0 =

θ1 + θ2 + . . . θN N

Also those UAVs with parameters deviate from the nominated parameter θ0 are abnormal. The above anomaly detection problem in multi UAVs formation can be solved by the maximum likelihood method, where noise ei (t ) in Eq. (8) is assumed to be one Gaussian white noise with zero mean and unknown variance. When to identify the unknown parameter vectors θˆi , the classical least squares estimation θˆi is defined as.

( T )−1 T φi φi φi Yi 1 θˆi = arg min |Yi − φi θi |2 = M M M θ [ Mi ]−1 M ∑ ∑ = ϕi (t )T ϕi (t ) ϕi (t )T yi (t ) t =1

(10)

t =1

where two vectors in Eq. (10) are given respectively as

{

Yi = yi (1)

[ [

φi = ϕi (1)

]T · · · yi (M ) ]T ϕi (2) · · · ϕi (M )

yi (2)

(11)

To improve the convergence rate and guarantee that the least squares estimation θˆi is unbiased, we always add sparseness property into Eq. (8), and then an improved sparse estimation θˆi,λ is obtained

θˆi,λ =

[ M ∑ t =1

ϕi (t )T ϕi (t ) + λI

]−1 [ M ∑ t =1

] ϕi (t )T yi (t ) + λθ0

(12)

where λ is a regularization factor, it is used to reduce the number of outliers. In practice, this regularization factor λ can be chosen as λ ∈ (0, 1). The detailed analysis about above least squares estimation θˆi and that improved sparse estimation θˆi,λ can be seen our two papers [2,11]. These two estimations θˆi and θˆi,λ will be beneficial for our later designing optimal input signal, because we need to compute their variances. 3. Dynamic programming techniques in anomaly detection As in our paper [12], we only have an idea of applying dynamic programming technique into our designed anomaly detector, but do not give a detailed procedure about how to apply. So here in this section, we study this subject to remedy the deficiency in paper [12]. Dynamic programming deals with situations where decisions are made in stages. The outcome of each decision may not be fully predictable but can be anticipated to some extent before the next decision is made. A key aspect of such situations is that decisions cannot be viewed in isolation since one must balance the desire for lower present cost with the undesirability of high future cost. The dynamic programming technique captures this tradeoff. At each stage, it ranks decisions based on the sum of the present cost and the expected future cost, assuming optimal decision making for subsequent stage. Here in this section, dynamic programming technique is introduced in the anomaly detection process in multi UAVs formation. In our paper [2], a new anomaly detector is summarized as a hypothesis testing problem. Given the residual εi (t ) at time instant t, to achieve the goal of anomaly detection, we compute the residual of each UAV by using identified model.

εi (t ) = yi (t ) − ϕiT (ui (t )) θi

(13)

The idea of anomaly detector based on residual is from statistical process control theory, it is often used to detect the signal’s persistent mutants. Based on Eq. (16) and parameter estimation θˆi , the variance of residual εi is defined as. W =

1 K −1

×

∑

∑

yi (t ) yTi (t ) + θi

i,t

∑

ϕi (ui (t )) ϕiT (ui (t )) θiT − 2θi

i,t

ϕ (ui (t )) T i

yTi

(t )

(14)

i ,t

where constant K = MN denotes the total number of the sample data. Given residual εi (t ) at time instant t, the detector D0 is defined. D0 (εi (t )) =

{

1 0

2

if W −1/2 εi (t ) ≥ T0 other w ise



(15)

where 1 corresponds to the abnormal state, and 0 is the normal state, T0 is the chosen threshold. But in above anomaly detector D0 , new observations are neglected and a decision may pay a cost or risk to accept one of these two hypotheses. So there exists a compromise between the cost of new observations and the higher probability of accepting the wrong hypothesis. Let yi (1) , yi (2) . . . yi (M ) be the sequence of observations. Based on (1) and (2), the mean value corresponding to each parameter θi is nominated value θ0 or its identified value θˆi . Then the probability distribution of observation yi (t ) is either f0 or f1 , and that we are trying to decide on one of these.

⎧ [ ( )2 ] ⎪ yi (t ) − ϕi (t )T θ0 1 ⎪ ⎪ ⎪ f0 = √ exp − ⎪ ⎪ 2σ 2π ⎨ ⎡ ( )2 ⎤ T ˆ y t − ϕ t θ ( ) ( ) ⎪ i i i ⎪ 1 ⎢ ⎥ ⎪ ⎪ f1 = √ exp ⎣− ⎦ ⎪ ⎪ 2σ ⎩ 2π

(16)

H. Wang-jian / ISA Transactions 91 (2019) 157–165

where σ is the variance of noise ei (t ), here for the sake of simplicity, we give a statistical description on noise ei (t ) not its bounded characterization. In case of bounded noise, dynamic programming technique in anomaly detection is the topic of our next paper, where game is introduced to solve the case of bounded noise. For any observation yi (t ), f0 and f1 denote the probabilities of yi (t ), when f0 and f1 are true distributions, respectively. At time instant t, using observations yi (1) , yi (2) . . . yi (t ), we may either stop collecting new observation and accept either f0 or f1 , or we may collect new observations at cost C > 0. If we stop collecting and make a choice, then we obtain zero cost when our choice is correct, and cost L0 and L1 if we choose incorrectly f0 and f1 . Assume a priori probability that the true distribution is f0 be p. Then we encounter one hypothesis testing problem with two states.

{

x0 x1

true distribution is f0 true distribution is f1

(17)

From convex optimization theory [14], function At (pt ) is concave, and satisfy for all t, and p ∈ [0, 1].

{

At (0) = At (1) = 0 At −1 (p) ≤ At (p)

pt = P (xt = x0 /yi (1) , yi (2) , . . . , yi (t ))

Using the concave property of function At (pt ), for all p ∈ [0, 1], we have. JM −1 (p) ≤ min [(1 − p) L0 , pL1 ] = JM (p) Applying the stationary of the system and the monotonicity property [15], we obtain. Jt (p) ≤ Jt +1 (p)

⎧ ⎪ pf0 (yi (1)) ⎪ ⎪ ⎪ p0 = ⎪ ⎪ pf0 (yi (1)) + (1 − p) f1 (yi (1)) ⎪ ⎨ pt f0 (yi (t + 1)) pt +1 = , ⎪ ⎪ ⎪ pt f0 (yi (t + 1)) + (1 − pt ) f1 (yi (t + 1)) ⎪ ⎪ ⎪ ⎪ ⎩ t = 1, 2 . . . M − 1

At +1 (p) ≤ At (p)

(19)

(20)

[

(1 − pt ) L0 , C + Eyi (t +1) )}] ( { pt f0 (yi (t + 1)) Jt +1 pt f0 (yi (t + 1)) + (1 − pt ) f1 (yi (t + 1))

L0

L1

Then as M → ∞, the sequences {αt } , {βt } converge to their limit values α, β respectively, and the optimal policy is approximated by the stationary policy.

⎧ ⎨accept f0

if pt ≥ α if pt < β new obser v ations

if

β < pt < α

⎧ pf (y (1)) f0 (yi (2)) . . . f0 (yi (i)) ⎪ ⎨pt = 0 i ∗ ∗ = pf0 (yi (1)) f0 (yi (2)) . . . f0 (yi (i)) + (1 − p) f1 (yi (1)) ⎪ ⎩ f1 (yi (2)) . . . f1 (yi (i)) accept f0 accept f1 continue

Jt (pt ) = min [(1 − pt ) L0 , pt L1 , C + At (pt )] pt f0 (yi (t + 1))

)}

pt f0 (yi (t + 1)) + (1 − pt ) f1 (yi (t + 1))

(22) An optimal policy for the last time instant is obtained from the minimization.

L0 + L1

βt < pt < αt

if

(24)

From Eq. (42), the stationary policy can be reformulated as.

For t = 1, 2, . . . , M, we have

L0

if pt ≥ αt if pt < βt new obser v ations

C ⎪ ⎩· · · ≥ βt +1 ≥ βt ≥ βt −1 ≥ · · · ≥

{

⎧ ⎪ ⎪ ⎨

γ =

L0 L1 L0 + L1

⎧ C ⎪ ⎨· · · ≤ αt +1 ≤ αt ≤ αt −1 ≤ · · · ≤ 1 −

p (yi (t + 1)) = pt f0 (yi (t + 1)) + (1 − pt ) f1 (yi (t + 1))

where γ is determined as the relation.

<

Now the conditional probability pt is obtained as.

where expectation operation over yi (t + 1) is taken with respect to the probability distribution.

pM ≥ γ pM < γ

L0 + L1

βt L1 = C + At (βt ) (1 − αt ) L0 = C + At (αt )

accept f1 ⎩ continue

(21)

if if

)

We have monotonicity inequalities.

where (1 − pM ) L0 is the expected cost for accepting f0 , and pM L1 is the expected cost for accepting f1 . Using Eq. (37), the optimal expected cost for time instant t is given as.

accept f0 accept f1

L0

Then an optimal policy for each time instant t is given as. accept f0 accept f1 continue

{

{

Then the optimal expected cost for the last observation yi (M ) JM (pM ) = min [(1 − pM ) L0 , pM L1 ]

( C + AM −1

where scalars αt , βt are determined as.

is.

{

∀t and p ∈ [0, 1]

Using Eq. (40), for all t and p ∈ [0, 1], we have

(18)

From probability theory [13], that conditional probability pt can be generated recursively as.

{ ( ⎪ ⎪ ⎩At (pt ) = Eyi (t +1) Jt +1

(23)

From Eq. (41), we obtain that if

One conditional probability from dynamic programming algorithm is given as.

Jt (pt ) = min

161

if Rt ≥ A if Rt < B new obser v ations

if

B < Rt < A

(25)

where

⎧ ⎪ (1 − p) β (1 − p) α ⎪ ⎪ ) ,B = ( ⎨A = p (1 − α) p 1−β ⎪ ⎪ f (y (1)) f0 (yi (2)) . . . f0 (yi (i)) ⎪ ⎩Rt = 0 i f1 (yi (1)) f1 (yi (2)) . . . f1 (yi (i)) And Rt can be generated by means of its recursive form. Rt +1 =

f0 (yi (t + 1)) f1 (yi (t + 1))

Rt

The optimal policy (25) is known as the sequential probability ratio test and used to replace that anomaly detector in [2]. The

162

H. Wang-jian / ISA Transactions 91 (2019) 157–165

advantage of our stationary policy is that here the cost of collecting new observations and the higher probability of accepting the wrong hypothesis are all taken into account simultaneously. The computational complexity of the proposed dynamic programming algorithm is known as ‘‘real arithmetic model of computations’’ as opposed to ‘‘integer arithmetic model’’. We assume that the computations are carried out by an idealized version of the usual computer which is capable to store countably many reals and can perform with them the standard exact real arithmetic operations.

From Eq. (30), the least squares estimation θˆi is unbiased. When to design optimal input signal, we continue to compute the variance of the asymptotic distribution.

[ ]2 = E θˆi − θ0 [ M ]−1 [ M ] ∑ ∑ T T = ϕi (t ) ϕi (t ) ϕi (t ) ei (t )

( )

cov θˆi

t =1

[ ×

M ∑

]−1 ]T [t =M1 ∑ T T ϕi (t ) ϕi (t ) ϕi (t ) ei (t ) t =1

[ M t =1 ]−1 ∑ T = ϕi (t ) ϕi (t )

4. Optimal input design for statistical noise According(to linear ) regression model (8) and its two parameter

t =1

estimations θˆi , θˆi,λ , now in this section we turn to consider the

where we use the following property of white noise.

main contribution of this paper. Due to the important procedure in designing optimal input signal for multi UAVs formation anomaly detection is to compute the variances of the asymptotic distribution, which ( are used ) to measure the accuracy about parameter esti-

E ei (t ) ei (s)T =

mations θˆi , θˆi,λ . As there are two kinds of parameter estimations

( ) θˆi , θˆi,λ , so we consider them separately as the following two

(31)

[

]

{

1 0

t=s t ̸= s

(32)

Based on variance matrix (31), its trace operation can be denoted as that.

[ M ]−1 ( )] ∑ T trace cov θˆi = trace ϕi (t ) ϕi (t ) [

(33)

t =1

cases.

If regressor vector ϕi (t ) is chosen as follows. ui (t ) ϕi (t ) = ui (t − 1)

]

[

4.1. The first case for least squares estimation

Then the above trace operation will be simplified as. Rewrite the linear regression model (8) to describe the linear relation for multi UAVs formation anomaly detection. yi (t ) = ϕi (t )T θi + ei (t ) , i = 1, 2 . . . N ,

t = 1, 2 . . . m

(26)

where the regressor vector ϕi (t ) is dependent on input signal ui (t ) or yi (t ) simultaneously, i.e.

{ ϕ (u (t )) ϕi (t ) = i i ϕi (yi (t ) , ui (t ))

(27)

[ M [ ]−1 ] ( )] ∑ [ ]T ui (t ) ˆ ui (t ) ui (t − 1) trace cov θi = trace ui (t − 1) t =1 ⎤−1 ⎡ M M ∑ ∑ 2 ui (t − 1) ui (t )⎥ ui (t ) ⎢ ⎥ ⎢ t =1 t =1 ⎥ = trace ⎢ M M ⎥ ⎢∑ ∑ ⎦ ⎣ 2 ui (t − 1) ui (t − 1) ui (t ) [

t =1

Let us continue to compute the least squares parameter estimation θˆi , and then it holds.

θˆi =

[ M ∑

]−1 [ ϕi (t )T ϕi (t )

M ∑

t =1

= θ0 +

[ M ∑

t =1 −1

] ϕi (t ) ϕi (t ) T

] ϕi (t )T yi (t )

[ M ∑

t =1

]

(28)

[ M ∑

ϕi (t ) ei (t )

]−1 [ M ] ∑ T ϕi (t ) ϕi (t ) ϕi (t ) ei (t ) T

(29)

t =1

As {ei (t )}M t =1 is an independent and identically distributed white noise, after taking limit operation on Eq. (29), we have.

( ) lim θˆi − θ0 = lim

M →∞

[

M →∞

×

M ∑

]−1

t =1

[ M ]−1 ( )] ∑ T trace cov θˆi = trace ϕi (t ) ϕi (t ) = [

ϕi (t )T ϕi (t )

2 M σu2

(35)

where σu2 is the variance function of the input signal, and its power spectral is φu (w). Applying the equivalence between variance function and its power spectral, the problem of optimal input signal design is formulated as the following Theorem 1. Theorem 1. The problem of optimal input signal design in multi UAVs formation anomaly detection can be summarized as the following minimization problem with finite input signal energy constrained condition.

⎧ ∫ ⎪ ⎪ ⎨ min φu (w)

2π 0

2 M φ∫ u (w) 2π

⎪ ⎪ ⎩subject to

dw (36)

φu (w) dw ≤ L2

0

The power spectral of the optimal input signal is given that.

t =1

[ M ∑

Using random process theory, we write the trace operation as that.

t =1

t =1

t =1

(34)

T

where we substitute equation (7) in above equation (28). To measure the bias error between the least squares estimation θˆi and its true value θ0 , we compute the difference.

θˆi − θ0 =

t =1

] ϕi (t ) ei (t ) = 0 T

(30)

φuopt (w) = ∫ 2π 0

L2 2 d M

w

=

L2 2 M

× 2π

=

ML2 4π

(37)

H. Wang-jian / ISA Transactions 91 (2019) 157–165

[

From Eq. (31), we see that, the power spectral of the optimal input signal is a constant under the condition that the least squares method is used to identify the unknown parameter vector.

×

163

Mui (t )2

∗

θˆi,λ =

[ M ∑

ϕi (t )T ϕi (t ) + λI

]−1 [ M ∑

t =1

∗

θˆi,λ =

[

ϕi (t )T yi (t ) + λθ0

(45)

η=

= θ0 +

[ M ∑

ϕi (t )T ϕi (t ) + λI

t =1

(46)

Theorem 2. The problem of optimal input signal design in multi UAVs formation anomaly detection can be summarized as the following minimization problem with finite input signal energy constrained condition. 2M η

⎧ ⎪ ⎪ ⎨min η

ϕi (t )T ei (t )

]−t =1 1 M ∑

φu (w) dw

Then the problem of optimal input signal design in multi UAVs formation anomaly detection can be reformulated as Theorem 2.

[ M t =1 ]−1t =[1 M ] ∑ ∑ T T = ϕi (t ) ϕi (t ) + λI ϕi (t ) ϕi (t ) + λI θ0 ]−1 t =M1 ∑

2π

∫ 0

ϕi (t ) ϕi (t ) + λI ] [ M M ∑ ∑ T T × ϕi (t ) ϕi (t ) θ0 + ϕi (t ) ei (t ) + λθ0

t =1

∫ 2π 2M 0 φu (w) dw = ]2 [ ∫ )2 2π Mui (t )2 + λ M φu (w) dw + λ 2Mui (t )2

0

T

ϕi (t )T ϕi (t ) + λI

= (

To simplify notation, we set the decision variable as that.

t =1

+

)]

(38)

]−1

[ M ∑

(

trace cov θˆi,λ

]

t =1

t =1

(44)

Taking the trace operation, it holds that.

Substituting Eq. (7) into (38), it holds that.

[ M ∑

]−1 ∗ 2 Mui (t ) + λ

Mui (t )2 + λ

4.2. The second case for improved sparse estimation As one regularization factor λ is introduced in that improved sparse estimation, then it will increase the complexity in computing the variance of the asymptotic distribution. Similarly we rewrite the improved sparse estimation θˆi,λ .

∗ Mui (t )2

[ ×

]

[M η + λ]2 ∫ ⎪ ⎪ ⎩subject to η =

ϕi (t )T ei (t )

t =1

2π

(47)

φu (w) dw ≤ L3

0

where L3 is the upper bound of the finite input signal energy, (39) The power spectral of the optimal input signal is given that.

Then computing the difference, we have.

θˆi,λ − θ0 =

[ M ∑

]−1 ϕi (t ) ϕi (t ) + λI T

t =1

M ∑

ϕi (t ) ei (t ) T

(40)

t =1

Furthermore computing the variance of the asymptotic distribution, it holds that.

(

cov θˆi,λ

]2 = E θˆi,λ − θ0 ] [ M ]−1 [ M ∑ ∑ T T = ϕi (t ) ϕi (t ) + λI ϕi (t ) ϕi (t )

)

[

t =1

[ ×

]−1t =1

M

∑

ui ( t ) ui (t − 1)

Proof. Due to 2M η

η

M ∑

ϕi (t )T ϕi (t ) + λI

t =1

ui (t − 1) ui (t )

ui (t − 1) ui (t )

u2i (t − 1)

M ∑ ui (t )2 + λ ⎢ ⎢ t =1 = ⎢ M ⎢∑ ⎣ ui (t − 1) ui (t ) t =1

cov θˆi,λ

)

[ =

Mui (t )2 + λ

∗

2M

η + 2M λ +

λ2 η

]

M ∑

⎤ ui (t − 1) ui (t )⎥

t =1 M ∑

ui (t − 1)2 + λ

⎥ ⎥ ⎥ ⎦

(43)

t =1

]−1 ∗ Mui (t )2 + λ

2M η [M η + λ]2

= arg min [ η

M2

2M

η + 2M λ +

λ2 η

] (50)

Applying geometric inequality, it holds that. (42)

Substituting Eq. (43) into the variance of the asymptotic distribution, as that.

(

M2

[ [ ] ] λ2 λ2 = arg max M 2 η + 2M λ + = arg max M 2 η + η η η η

]

u2i (t )

⎡

]

]=[

(49)

arg min

Then we have the following relations, for example.

[

M 2 η2 + 2M λη + λ2

Then

]

[

2M η

= [

(41)

ϕi (t )T ϕi (t ) + λI

ϕi (t )T ϕi (t ) =

(48)

⎪ ⎩φ opt (w) = L3 u 2π

[M η + λ]2

Also when regressor vector ϕi (t ) is the form of that.

ϕi (t ) =

φuopt (w) dw = L3

0

t =1

[

2π

⎧ ∫ ⎪ ⎨η =

√ 2 λ λ2 M 2η + ≥ 2 M 2η × = 2M η η η

(51)

The equality holds under the condition that. M 2η =

λ2 λ ⇒η= η M

(52)

But inequality (51) signified the minimization operation, not the maximization operation. In order to solve maximization operation, we can use the mathematical analysis theory.

164

H. Wang-jian / ISA Transactions 91 (2019) 157–165

Fig. 3. The information flow in formation flight.

Firstly by differentiation with respect to η and by setting the derivative equal to zero, it holds that.

λ2 M η+ η

[

2

]′

= M2 − ⇓

η=

λ2 =0 η2

Fig. 4. Parameter estimation errors among parameter estimations and their nominated value.

(53)

λ M

λ

[ Then η2 ]= M is a point of inflexion on performance function M 2 η + λη , it means that if Mλ ≤ η ≤ L3 , then it holds that. M2 −

λ2 ≥0 η2 λ

If 0 ≤ η ≤ M2 −

(54) M

, then it holds that.

λ2 ≤0 η2

(55)

[ These two ] Eqs. (54) and (55) tell us that performance function ] [ 2 M 2 η + λη is an increasing function in interval η ∈ Mλ L3 , [ ] otherwise it is a decreasing function in interval η ∈ 0 Mλ . So the maximization operation is taken as the set point in interval [ ] η ∈ λ ML3 , i.e.

ηopt = L3

(56)

Substituting Eq. (56) into the definition, it holds that.

⎧ ∫ ⎪ opt ⎪ η = ⎪ ⎨

2π

0

⎪ ⎪ ⎪ ⎩

Fig. 5. Parameter estimations.

φ

opt u

φuopt (w) dw = L3 (57)

⇓ (w) =

L3

2π Combining two optimal power spectrums in Eqs. (48) and (57), we conclude that whatever least squares method or sparse optimization algorithm is applied, the optimal input signal, used to excite the whole identification experiment, may be a constant. 5. Numerical example Consider the multi UAVs formation flight process and the flight performance data corresponding to each UAV can be collected in ground control station to construct database. One linear unknown relation can be set up by using flight data. The whole flow of information with respect to the ground control station is seen in Fig. 3. In Fig. 3, the formation flight data are stored in database, then one linear regression model is constructed through the database. This model uses different unknown parameters to represent different UAVs. Based on this linear regression model with respect

to unknown parameters, we apply the projection algorithm with dead zone to solve all unknown parameters. Assume that there are nine UAVs, where three UAVs are known as abnormal in prior. The number of collected measurements is 500, i.e. in equation (4), it holds that. N = 9, k = 3, M = 500 According to the ith UAV at time instant t, yi (t ) denotes the angle of attack, ϕi (t ) is one regression vector consisted by velocity, dynamic pressure, elevator deflection and the stabilizer deflection angle, and measurement noise ei (t ) is a bounded noise, where the upper bound is 0.5, i.e.

|ei (t )| ≤ 0.5 Further let the mean and variance corresponding to nominated parameter value θ0 in normal formation flight be that. 0.8 0.04 ⎢ −3 ⎥ ∑ ⎢0.12 θ0 = ⎣ , =⎣ −0.6⎦ 0.02 0 0.5 0.02

⎡

⎤

⎡

0.12 0.08 0.09 0.2

0.02 0.08 0.03 0.1

0.02 0.01⎥ 0 .1 ⎦ 0.05

⎤

H. Wang-jian / ISA Transactions 91 (2019) 157–165

165

Fig. 6. The optimal input signal.

Then the nominated parameter θ0 is a stochastic variable.

( θ0 ∼ N θ 0 ,

) ∑ 0

That regression vector ϕi (t ) is also a stochastic variable.

)

( ϕi (t ) ∼ N ϕ,

∑ Acknowledgment

ϕ

where 0.9 0.25 ⎢−1.2⎥ ∑ ⎢−0.02 ϕ=⎣ , =⎣ −2.8⎦ 0.12 ϕ 0.7 −0.04

⎡

⎤

in deriving the variance, we need an assumption about the independent and identically distributed white noise. But white noise is an ideal case, as it will not exist in any engineering. So to relax the strict description on noise, we should design the general optimal input signal design in the presence of unknown but bounded noise. This new subject aims at minimizing the worst case identification error with respect to the chosen input signal.

⎡

−0.02 0.45 0.03 −0.52

0.12 0.03 0 .9 −1.2

⎤ −0.04 −0.05⎥ −1.2 ⎦

This work is partially supported by the Grants from the National Science Foundation of China (No. 61364014) and (No. jxxjb18020). Conflict of interest

3

Substituting above data information into Eq. (28) and applying our least squares algorithm to solve these nine unknown parameters, the parameter estimation errors among these nine parameter estimators {θi }Ni=1 and their nominated value θ0 is plotted in Fig. 4. From Fig. 4, we see that the parameter error between the first parameter θ1 and nominated value θ0 converges to zero, this result also holds for other parameters such as θ3 , θ4 , θ5 , θ7 , θ8 . Then those six UAVs corresponding to these six parameters are normal. Meanwhile those three parameters (θ2 , θ6 , θ9 ) deviate from nominated value θ0 , it means that these three UAVs are regarded as abnormal state. In Fig. 5, the parameter estimations corresponding to the first UAV is given. From the iterative curves in Fig. 5, we see that these identified parameters will all converge to their own true values, i.e. the estimation errors between parameter estimator and their nominated value will be zero. Then we can conclude that the first UAV is normal. The similar analysis can be expanded for other UAVs in the whole formation flight. Here the input signal is chosen as a binary signal, its amplitude is ±2, and its length is 10S. We use one optimization method to solve that minimization problem with constrained condition, then the obtained input signal is shown in Fig. 6 after twenty iterations. From Fig. 6, we see that the optimal input signal flutters around the zero amplitude, then it coincides with the binary signal. 6. Conclusion This paper continues to design the optimal input signal for multi UAVs formation anomaly detection, in order to make the designed input signal excite the whole identification experiment. According to two identification strategies for anomaly detection, we compute the variance of the asymptotic distribution, and choose the trace operation to construct one numerical optimization problem, whose solution is corresponded to the optimal power spectral. But

The author declares that there is no conflict of interests regarding the publication of this paper. References [1] Karimoddini A, Lin H. Hybrid three dimensional formation control for unmanned helicopters. Automatica 2013;49(2):424–33. [2] Jian-hong Wang. Sparse optimization algorithm in multi UAVs formation anomaly detection. Electron Opt Control 2015;22(8):1–7. [3] Ljung L. System identification: Theory for the user. Prentice Hall, Upper Saddle River: Prentice Hall; 1999. [4] Donoho DL. Compressed sensing. IEEE Trans Inform Theory 2006;52(4):1289– 306. [5] Braslavsky JH. Minimum variance control over a gaussian communication channel. IEEE Trans Automat Control 2011;56(8):744–56. [6] Imer OC, Yuksel S. Optimal control of dynamical systems over unreliable communication links. Automatica 2006;42(9):1429–40. [7] Gupta V, Spanos D. Optimal LQG control across a packet dropping link. Systems Control Lett 2007;56(6):439–46. [8] You Keyou, Fu Minyue. Mean square stability for kalman filtering with Markovian packet losses. Automatica 2011;47(12):2647–57. [9] Silva EI, Goodwin GC. Control system design subject to SNR constraints. Automatica 2010;46(2):299–310. [10] Polyanskiy Y, Poor HV. Channel coding rate in the finite blocklength regime. IEEE Trans Inform Theory 2010;56(5):2307–59. [11] Jian-hong Wang. Bias compensated estimation in multi UAVs formation anomaly detection. J Control Syst Eng 2016;4(1):40–50. [12] Jian-hong Wang. Dynamic programming technique in multi UAVs formation anomaly detection. Int J Innovative Comput Inf Control 2018;14(6):1977–81. [13] Pintelon R, Schoukens J. System identification: a frequency domain approach. New York: IEEE Press; 2001. [14] Boyd S, Vandenberghe L. Convex optimization. Cambridge University Press; 2004. [15] Berstekas DP, Tsitsiklis JN. Parallel and distributed computation: Numerical methods. Athema Scientific Press; 1997.