State estimation method using median of multiple candidates for observation signals including outliers

The state estimation problem for systems in which observation outputs include outliers is addressed herein. When the observation output has outliers, the accuracy of the state estimation is dramatically worse. To overcome this problem, a novel observer structure using multiple candidates of the estimated state is proposed. First, multiple candidates of the estimated state are created; each candidate uses the sensing output value of a different detection timing. If outliers occur infrequently, eliminating candidates affected by outliers can prevent deterioration in estimation accuracy. Our proposed observer select one from the obtained candidates of the estimated state using a median or a weighted median operation. Through the median operation, the estimated state that does not use the outlier value is selected from these candidates. In addition, a method is provided to design the observer gains of these estimated state candidates based on a reachable set of the estimated state error, using Lyapunov-based inequalities. The effectiveness of the proposed observer is illustrated using numerical examples.


Introduction
In control engineering, state feedback control is one of the most important methods for control systems. In state feedback control, the sum of all the state values weighted at each time is calculated, and the calculated values are applied to the plant as an input signal. In general, it is difficult to measure all state values with sensors; therefore, the state observer is also used in the implementation of control systems [1][2][3][4][5][6]. The state observer is a system that provides an estimate of the internal state of a given real system from the input and measured output of the real system. It is used in various situations, such as when the sensor cannot be used or when not all the internal states can be observed because of cost problems. By including the model of the plant in the state observer, one can seek the estimated state x f of the internal state x by simulating the action using the information of input u and output y. Although state observers have been studied for a long time, various research fields involving state observers have progressed in recent years, such as vehicle-motion control systems [8,9], networked control [10] systems and cyber-physical systems. The state estimation by the state observer and the state estimation technique using a Kalman filter can be applied to the case of linear and nonlinear systems [4]. The appropriate estimation of the state quantity of a plant plays an important role in control system design. Various state observer design methods have been developed to achieve satisfactory estimation performance [15][16][17].
In general, when a state observer is used, it is assumed that noise and disturbance exist. It is known that the estimated accuracy is significantly degraded when problems related to outliers [11][12][13][14]19] and packet losses [21][22][23] occur. In the case of using a noncontact sensor, such as a visual sensor, outliers can sometimes occur. Because of the influence of dramatic changes in luminance, communication obstruction by moving obstacles, etc., there are some cases where the sensor temporarily detects a value that is clearly different from the intended signal to be observed. In such a case, even if the dramatic change in the value obtained by the sensor is temporary, because the influence on the control system increases in real-time control, it is necessary to remove the influence of the outliers appropriately. Furthermore, in networked control [21][22][23], packet loss is an important issue, as well as signal quantization and time delay. Packet loss can also be regarded as an outlier. Missing data occur when the data traffic congestion suddenly breaks because of the limited capacity of the information transmission. In recent years, research related to information security [24,25] has been advancing, and, considering the possibility of outliers arising because of data tampering, a theoretical construction regarding the removal of the influence of outliers is an important element.
State estimation problems caused by outliers have been studied previously [11,12,22,23]. A Kalman filter and H ∞ filter based on the statistical nature and frequency characteristics of outliers were designed [11,22]. Moreover, research regarding control systems, including packet loss, has been conducted. The performance of outlier values applied to a Kalman filter was analysed [23]. In previous work [13], the covariance matrix of the observation noise of the Kalman filter was updated by utilizing the observation error information from outliers. The design can also be implemented with a policy to stop updating when the observation error is large by increasing the covariance or positively updating if the error is small. However, there is a limit to the correspondence of utilizing statistical properties and covariance. In particular, it is difficult to eliminate the effect of outliers that have the possibility of changing to unexpectedly large values because of the dynamics arising in the linear filter and the accompanying signal processing.
In other research, however [13], robust state observers that limit the maximum value of the error of observation output were considered. These methods correspond to setting an absolute upper limit of the prediction error using dead zone components or a threshold. Consequently, it is possible to handle the influence items in a limited manner, although it is impossible to eliminate the influence of outliers. A moving horizon approach to remove the effect of outliers was proposed for systems with a single output [14]. In addition, a nonlinear/logical operation method that ignores outliers entirely by using a majority decision has been proposed [19]. In other work [19], the state observer eliminated the influence of outliers on state estimates by performing majority decision operations on multiple observed signals. Furthermore, a state observer that eliminates the influence of outliers by calculating thresholds by using multiple observable sensor pairs and removing signals exceeding the threshold was proposed [20]. Even if an enormous value is applied as an outlier, these methods can effectively perform state estimation because they include a structure for ignoring that value. However, the class of the target systems is limited because it is assumed that a plurality of output signals can be observed. Moreover, it cannot be applied in cases where all values are outliers or packet losses at the same time.
In this paper, a novel state observer is proposed that eliminates the influence arising from outliers at a specific time. As an alternative to methods where states are derived from multiple outputs [20], a method of state estimation is constructed that selects multiple candidates of estimated values in the direction of the time axis and then selects from the candidates. In general, the order of the state observers is equal to or smaller than that of the plants. In this study, the aim is to eliminate the influence of higher-order outliers. By deriving candidates of multiple state estimates and performing nonlinear processing to adopt candidates that are not influenced by the outliers in the values of candidates, a novel state observer that has the robustness needed to overcome outliers is constructed. As a specific nonlinear processing method, the median is used. The values from the candidates obtained by the median, which are not affected by the outliers, are selected. The proposed state observer is called the "median of candidate vectors" (MCV) observer.
This paper is structured as follows. In Section 2, the principle of the state observer and the median calculation method as a mathematical preparation are explained. In Section 3, the MCV observer is proposed as a state observer that is robust to outliers. First, the derivation of the candidate value of the state quantity estimate at time k + 1 using the state estimate values from step 0 to N − 1 is described. The estimate at time k + 1 is determined by applying the median operation to the candidate value. The state observer algorithm is also described. In Section 4, the matrix inequality condition for designing the observer gains based on a robust invariant ellipsoid is shown. In Section 5, numerical simulation that confirmed that the state can be estimated with high accuracy, even if outliers occur in the MCV observer, is described. Finally, this research is summarized in Section 6, and the scope for future work is discussed.
The basic idea of the algorithm of the MCV observer was presented at the SICE Annual Conference [31] and Trans. of SICE [32]. The observer gains were determined based on the linear matrix inequalities (LMIs) for evaluating an invariant set for a peak-bounded disturbance in [32] (Japanese article). This report is an advanced version of the presented work to generalize the multiple-input multiple-output setting. Moreover, weighted median is addressed and corresponding numerical results are discussed in this paper.

Plant
This study addressed a discrete-time linear timeinvariant system. The current time is taken as k, and the dynamic characteristics of a plant are given by the state equation as follows: Matrices A, B u , B d , C, and D are suitably defined matrices that correspond to the square matrix and the inputs and outputs expressing the dynamic characteristics of plant, and the order of the plant is m. A state of the plant is defined as x p (k) ∈ R m , the control input as u(k) ∈ R m i1 , and the observed output as y(k) ∈ R m o1 . Signals d(k) ∈ R m i2 and w(k) ∈ R m o2 are the process noise and the sensing noise, respectively. The observability of the plant is assumed.
Here, it is assumed that the initial state x p (0) is given as x p (0) = 0 and that the control input, process noise and sensing noise before time k = 0 are given as u(k) = 0, d(k) = 0 and w(k) = 0, respectively. It is assumed that the maximum absolute values of each element in d(k) and w(k) are less than or equal to 1. In another way, d(k) ∞ ≤ 1 and w(k) ∞ ≤ 1 hold.
The term y out (k) characterizes the outliers. The term y out (k) = 0 holds in ordinary times. When an outlier occurs at k = k 1 , y out (k 1 ) = 0 holds, and the magnitude of |y out (k 1 )| has no size limit. It is considered that the outlier values are drastically larger than the noise. In particular, if y out (k) = 0, ∀k holds, it means an observed output with no outliers. One can express the packet loss such that To distinguish the outliers from observations and noise, it is assumed here that the absolute value of y out (k) is sufficiently large compared with other signals.

Observation frequency of outliers in observed output
The detailed formulation of y out (k) is presented in this section. In this study, outliers are supposed to occur instantly. In other words, it is considered that outliers occur for a short period. Two integers F 1 and F 2 are used to characterize the outliers in the observed output. For any time interval [k, . . . , k + F 1 − 1] with time k, the maximum number of outliers appearing in the time interval is F 2 . In other words, F 1 determines the time interval and F 2 characterizes the outlier occurring frequency in the time interval. For example, in the case where an outlier occurs every 10 periods in the observed output, i.e. a periodically observed outlier. Then, F 1 = 10 and F 2 = 1 are set to characterize the outliers. Next, F 1 and F 2 are set based on the frequency of outliers in the sensor output at the time of implementation. Figure 1 shows an example of the observed output with outliers. Large periodic outliers for two output signals. 30-periodic outliers (red mark) and 100periodic outliers (blue mark) occur for output signals. In this example, one can set F 1 = 30 and F 2 = 2, and the number of maximum outlier occurrence times (the outliers occur for output 1 or output 2) in the time range [k, . . . , k + 29] is one for any range with k.
In this study, it is assumed that the following condition is satisfied for F 1 and F 2 : In various settings for outliers in the sensor outputs, the outliers occur sparsely. In such a case, the former assumption can be easily satisfied.

Median operation
Because the calculation cost of the median operation is small in many computer devices, the median operation can be used in a wide range of research fields. For example, it is used as a filter in image processing and is known to be effective in the removal of salt-and-pepper noise from original images. The median is a value located in the centre of the whole data if the data are arranged by order of magnitude. When the number of values is even, it can be obtained by taking the average value of two numbers located at the centre. For data ξ(i), i = 1, . . . , n when the values are sorted in ascending order is ξ(1),ξ(2), . . . ,x(n). The median value ξ m is obtained as follows: The median operator is denoted as MED [·], and the median value is given by the following equation: The weighted median filter was first introduced as a generalization of the median filter (5). Nonnegative integer weights are assigned to each value. For the data x(i), i = 1, . . . , n, the output ξ wm of the weighted median filter of span n associated with the weights is given by The symbol denotes duplication, i.e.
When the weight w i = 1 for all i, (7) becomes the standard median filter.

Outline
A state observer is proposed for improving the robustness against outliers ( Figure 2). In this section, candidate state vectors are described, and the median is used as the method for obtaining estimation values from them. Moreover, the implementation of the algorithm is described.

Candidate of estimated state using feed forward estimation
The usual state observer uses y(k) and u(k) to estimate one step after the state of time k. Therefore, for the case in which the observation signal y(k) at time k is replaced by an outlier, a proper estimation cannot be made. Thus if one uses the estimated values from the N step precondition as candidates and selects the appropriate state from those candidates, the estimated state x, which is not affected by outliers, can be generated.
First, candidates of state estimates are described. In the MCV observer, the estimated value at time k − i is used to estimate the state of time k + 1. Herein, the estimated value using a previous estimated state x(k − i) is denoted asx i (k + 1).
wherex i (k + 1) represents the estimated value of the state at time k + 1, and, for each vector, only the value of y(k − i) is used as the observed output. L i is an observer gain for obtainingx i (k + 1). The observed output at the other time (k − j, j = i) is not used. For the state estimated by the MCV observer not to be affected by outliers, it is necessary to place constraints on the number of occurrences of outliers within a certain interval. In the above description, N is selected as the number of candidates by the expression (9). The number of candidates N is set to satisfy the following equation using F 1 and F 2 in Section 2.2.
The case where F 1 = 10 and F 2 = 1 is considered as an example. Then, one can select N = 3. First, if i = 0 is calculated using (9),x 0 is given by the following equation: This is the same structure as the usual estimated structure of the Luenberger observer. In addition, for i = 1, 2,x is given. Here, each estimated value in i = 1, 2 transits the state in a feed forward manner with respect to the normal state estimation value.
As an example, the case of using three state estimate candidates given by (11)-(13) is considered. One can correctly detect, with y(k − 1) and y(k − 2) among these three candidates, y(k) = 0 by outliers. Here, it is assumed that the estimated value ofx 0 (k + 1) contains a large error. However,x 1 (k + 1) andx 2 (k + 1) are estimated precisely. The error contained in the given estimate is expected to be small. If a majority vote is taken among the three candidates, it is possible to estimate the correct state. The estimation is implemented by applying this consideration.
In this study, the median is used for estimation based on candidate values. As with the majority vote, it is expected that, among the candidates of state estimates aligned in the order of norms, outliers' influence will not go to the middle of the aligned ones. If an outlier occurs at time k, the influence of outliers appears inx 0 in the estimate of time k + 1. However, estimation can be performed correctly withx 1 ,x 2 . The influence of the outlier appears inx 1 at the next estimation time k + 2, and the influence of outliers inx 0 ,x 2 is not accepted. Thus, because only one ofx 1 ,x 2 andx 3 is affected by outliers, by taking these median values, the influence of outliers can be eliminated.
Setting up an expression for i = 0, 1, 2 . . . , N − 1 makes it possible to estimate the state via the median using the N candidate vectors of the estimated state.

end loop
From the above results, the proposed state estimation method is summarized as follows.

Algorithm of MCV observer
The configuration procedure of the proposed MCV observer based on the above expansion is shown below. The number of state estimate candidates N is determined using F 1 and F 2 .
To evaluate whetherx i (k) is affected by an outlier, the following output y e (i) is defined.
Here, y e (i) can be calculated because x(k − i − 1) and u(k − i − 1) are known values, and y(k − i) is the observed output. As an alternative expression, the signal y e (i) can be expressed as follows: where e is the error between x p (k) and x(k), i.e. e(k) = x p (k) − x(k) holds. If one focuses on the right-hand side of (15), the influence of y out (k − i) is greater than that of y out (k − i − 1), w(k − i) and d(k − i − 1) when an outlier occurs. Therefore, it is expected that one can judge whether the observed output y(k − i) includes an outlier by focusing on the magnitude of y e (i) of (14). Because N is chosen to be 2F 2 + 1 ≤ N, the number of candidates in y e (i) that are affected by outliers, i.e. y out = 0, is less than F 2 . This shows that for i * , which is the median of y e (i)(i = 0, . . . , N − 1), y e (i * ) does not contain an outlier, i.e. y out = 0. Therefore, the MCV observer in this study uses the signal y e (i) for selecting 1 from x i (k + 1) of candidate state estimates.
Using (9) and (14) with weighting vector T ∈ R 1×m o1 , the algorithm of the MCV observer is given as follows.
The proposed MCV observer algorithm is simple. Note that weight vector T should be selected not to disappear the outliers. The number of candidates affected by the outlier is smaller than F 2 . Therefore, the candidate selected by the median operation is not affected by Algorithm 2 Estimation structure of MCV observer with weighted median filter loop [1.] Calculatex i (k + 1) for all i = 0, · · · , N − 1 by (9) and calculate y e (i) in (14). [2.] Calculate the median to find the following value i * .
[3.] x(k + 1) =x i * (k + 1) and move to the next step. end loop outliers, because of (10). The MCV observer guarantees that the estimated state is not affected by outliers.
It is also possible to use a weighted median operation instead of the median operation in Algorithm 1. To use weighted median it is required to choice appropriate W for selecting suitable estimated state from candidates. Condition of W for the limited case with F 2 = 1 is as follows: As an example of the weight W for N = 4 is as follows: We can find that (17) satisfy the condition in (16). Median value is not affected by the outlier value for any time. Then, the modified algorithm of Algorithm 1 using median filter is given as follows.
Discussion about weight design of W is shown later part of this paper.

Observer gain design using robust invariant set
In this section, observer gains L i in (9) is designed based on the maximum value of an evaluation output. The following error system is obtained In addition to this equation, the following evaluation output is dealt with to consider the noise rejection performance by the MCV observer ( Figure 2).
E is an evaluation vector that characterizes the evaluation signal. When E = I in (20) is selected, the state error is evaluated by z(k). A design method for the observer gains L i was developed, focusing on the maximum absolute value of the signal z e (k).
An error system (19) designed based on the results of a robust invariant set for peak-bounded disturbance reported elsewhere [20,[26][27][28][29][30] is analysed. The robust invariant set in this method is based on an ellipsoid using a matrix inequality condition. A theorem from previous work [26] is used to find a common invariant set for all state candidatesx i . By setting up a common invariant set, the stateê i of the error system stays in the invariant set for any choice of i. Thus if one can find a common invariant set for all candidatesx i , the obtained set is also invariant to the state error of the MCV observer.
First, the gain design is considered using the state set, satisfying the following equation: The state error e(k) includes an ellipsoid characterized by a positive definite matrix P. First, it is assumed that the following equation holds for a certain k.
In this case, it is necessary to consider the conditions for the following equation to hold.
If (23) holds for all i = 0, . . . , N − 1, because one state candidate is selected fromê i (k + 1), then the following inequality holds, and (21) is a robust invariant set of the MCV observer.
The following lemma is given based on the theorem in [26] to obtain the analytical condition of the robust invariant set using common matrix P.

Lemma 4.1:
Assuming that L i is given for all i and that the noise d(k) ∞ ≤ 1 and w(k) ∞ ≤ 1 given. If the following conditions: are satisfied for all i, the robust invariant set is given by whereĀ i andB i are given as follows: = m i2 (i + 1) + m o2 (28) The inequality condition of Lemma 4.1 is attributed to LMI with respect to P > 0 if the scalar parameters α i are fixed. The following inequality conditions are further considered for z e (k) defined in (20).
The symbol γ is a scalar positive variable. The necessary and sufficient conditions for this inequality condition to hold for e i (k) are given as follows: If z e (k) is obtained for z e (k) defined in (20), then the following equation holds by (29) and (26).
The evaluation signal z e (k) is less than or equal to γ for all k by the inequality (31). Thus the following theorem is derived from the fact that the evaluation signal for a given MCV observer gain, which satisfies the conditions of (30) and Lemma 4.1, is characterized by γ .

Theorem 4.2:
Assuming that d(k) ∞ ≤ 1 and w(k) ∞ ≤ 1 are satisfied. For the given L i , the following inequality conditions are considered: Given P > 0, which satisfies the inequality condition for all i, the following equation holds for the evaluation signal z e of the system using the MCV observer of Algorithm 1.

e(k) T Pe(k) ≤ 1 (35)
Furthermore, the first inequality condition (32) holds. The condition (29) is satisfied, and z e (k) is a scalar; thus one obtains the following equation: From the above, |z e (k)| ≤ γ holds, and Theorem 4.2 is satisfied.
Here, one can fix α i , minimize γ as an LMI problem and can find the α i that gives the optimal value of γ . In this framework, it is important to determine whether the given parameters α i satisfy the condition α i ∈ [0, 1 − ρ(Ā i ) 2 ]. However, finding these parameters is not difficult because the search ranges are limited by [0, 1]. A method for determining the parameters is also discussed in [26].
The next step is to develop a design problem for the observer gains based on the analysis method characterized in Theorem 4.2. The Schur complement is applied to the inequality constraint of (32) in Theorem 4.2. Furthermore, Y i := −PL i are regarded as the design variables instead of L i . By applying the transformation, one can solve the following gain design problem.

Problem 4.1:
Find the minimum γ with P, Y i for which γ , P and Y i satisfy the following matrix inequality conditions for all i.
In this case, the matrix inequality of Problem 4.1, as well as the theorem giving the analysis conditions, can be considered as an LMI problem if the scalar parameters α i are fixed. Therefore, as in the analysis, one can fix α i and minimize γ as an LMI problem. Parameters α i are searched by hyperplane search that gives the optimal value γ . By solving Problem 4.1, optimal parameters P and Y i which satisfying |z e (k)| ≤ γ are obtained. Then, the optimal observer gains L i are given by L i = −P −1 Y i .

Evaluation of state candidates for weight design
A method about weight design of W is discussed in this section. The observer gain L i for all candidate (9) is assumed to be given in this section. When we estimate state using Algorithm 2, it is required to give appropriate setting about W. One approach for give W is introduced based on the robust invariant set. In Section 4.1, the observer gain L i is designed using common Lyapunov matrix P. On the other hand, the performance of each estimated state candidatex i is analysed using an individually assigned Lyapunov matrix P i . It can be assumed that small upper bound rather than γ is obtained when a certainx i is selected by median filter for all k. For evaluating potential performance of each candidatex i , we use γ i and analyse minimum γ i which satisfy the following equation: The inequality |Ee i (k)| ≤ γ i is satisfied by (37). The performance of each candidatex i can be estimated by solving the following problem.

Problem 4.2:
Find the minimum γ i with P i and α i for which γ i , P i and α i satisfy the following matrix inequality conditions: When L i is designed by solving Problem 4.2, the following equation is satisfied.
It can be regarded the estimated statex i is better when γ i is small. Therefore, it is useful to assign large value w i when γ i is smaller. It is expected that γ i is smaller when i is small value because the candidatex i with small number i is few affected by the process noise rather than the case with the candidatex i with large number i. One approach is to consider the ratio value of 1/γ i for assigning w i . To define w i , the condition (16) should also be satisfied.

Verification of the effect of MCV observer
The numerical example of the MCV observer is shown in this section. The plant dynamics are given as A test input signal for the plant in the simulation is as follows: The state estimation is performed. Problem 4.1 is solved by solving a semidefinite problem. In this section, the observed output shown in Figure 1 is used in the examples. F 1 = 5 and F 2 = 1 are given based on Figure 1 and Note that optimal variables α i are obtained by hyperplane search. Figure 3 shows the state of the plant and the estimated state by the MCV observer with gains (43). To compare the result, an estimated state by the well-known Luenberger observer (11) is shown in Figure 4. Dashed lines are the state of the plant, and red lines are the estimated state. It can be seen that the outliers in Figure 1 do not affect the estimated state in the MCV observer compared to the Luenberger observer. The optimal performance by (43) is characterized by γ = 1.4064, and the following result is obtained using the designed gains: |z e (k)| ≤ 1.4064 (44) Figure 5 shows the state error between the plant state and the estimated state. The effect of outliers in Figure 1 is not presented in Figures 3 and 5. Figure 6 shows the performance index signal z e (k). One can find that (44) is satisfied for all k from Figure 6. Then, the state estimation performance is compared by changing the candidate number N. In Table 1, γ is obtained by solving Problem 4.1 for each N. We can   find that when N increasing, γ becomes larger, and it is better to select small N under the condition (10).

Weighted median
In this section, effect of the weighted median is discussed using some examples. The same plant in Section 5.1 is used in this example, and N = 4 is selected which satisfy the condition (10). The solution   γ by solving Problem 4.1 is presented in Table 1 with N = 4. The performance of each candidatex i is estimated by solving Problem 4.2 and the results are shown in Table 2. We can find that performance with small i is better. Then, the selection frequency of each state candidate is presented using Algorithms 1 and 2. Algorithm 2 is MCV observer with the weighted median. (17) is selected as the weight W. The time series of the selection of the estimated state candidate number i is shown in Figures 7 and 8. Figure 7 is for Algorithm 1, and Figure 8 is for Algorithm 2 with (17). In addition, the number of times selected for each candidate is shown in Table 3. It can be seen that by using weight W, the selection frequency ofx i can be changed. For example, the selection frequency ofx 0 is the larger in case weight W is used. Candidates with larger values of w i are selected more often from Table 3. Figures 9 and 10 show its obtained evaluation signal z e (k) for each cases. We can find that it has potential to achieve better performance by MCV observer with weighted median because the amplitude of z e (k) with weighted median is slightly small.

Conclusion
This paper proposed an MCV observer that yields state estimation values independently as candidates and uses the median of the estimated states. Because multiple candidates are provided by the observed output y at independent sampling times, it is possible to obtain estimated states with robustness against outliers using a simple algorithm. This fact was verified using numerical examples, and the proposed method was shown to be effective. Moreover, the weighted median was addressed and corresponding numerical results were discussed in this paper. A design method for observer gains based on the robust invariant set in Problem 4.1 was also proposed. The design problem of the observer gains was formulated based on an LMI optimization problem.
In this case, all candidates of the state estimates are handled evenly. However, when there is no influence of outliers, the estimation accuracy is higher for those using state estimates closer to the current time.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work was supported by Japan Society for the Promotion of Science.

Notes on contributors
Hiroshi Okajima received his M.E. and Ph.D. degrees from Osaka University, Japan, in 2004 and 2007, respectively. He is presently an associate professor at Kumamoto University, Japan. His research interests include tracking control, analysis of non-minimum-phase systems, and data quantization for networked systems. He is a member of ISCIE and IEEE.