Kalman Filter-Based RAIM for Reliable Aircraft Positioning with GPS and NavIC Constellations

This paper presents a novel Kalman filter (KF)-based receiver autonomous integrity monitoring (RAIM) algorithm for reliable aircraft positioning with global navigation satellite systems (GNSS). The presented method overcomes major limitations of the authors’ previous work, and uses two GNSS, namely, Navigation with Indian Constellation (NavIC) of India and the Global Positioning System (GPS). The algorithm is developed in the range domain and compared with two existing approaches—one each for the weighted least squares navigation filter and KF. Extensive simulations were carried out for an unmanned aircraft flight path over the Indian sub-continent for validation of the new approach. Although both existing methods outperform the new one, the work is significant for the following reasons. KF is an integral part of advanced navigation systems that can address frequent loss of GNSS signals (e.g., vector tracking and multi-sensor integration). Developing KF RAIM algorithms is essential to ensuring their reliability. KF solution separation (or position domain) RAIM offers good performance at the cost of high computational load. Presented range domain KF RAIM, on the other hand, offers satisfactory performance to a certain extent, eliminating a major issue of growing position error bounds over time. It requires moderate computational resources, and hence, shows promise for real-time implementations in avionics. Simulation results also indicate that addition of NavIC alongside GPS can substantially improve RAIM performance, particularly in poor geometries.


Introduction
In recent years, global navigation satellite systems (GNSS) have evolved into important infrastructure of modern society, with day-to-day activities increasingly relying on their positioning, navigation and timing services [1,2]. GNSS is now so deeply entrenched in everyday life that disruption or malfunction of its operations can lead to undesirable consequences. Among many other diverse applications, aircraft navigation in particular has immensely benefited from GNSS during different flight phases (e.g., cruising, terminal, and approach to landing) [3][4][5][6]. GNSS is also identified as a key sensor on-board unmanned aerial vehicles (UAV) [7][8][9]. In spite of ever expanding use, major concerns in aviation applications remain to be reliability and continuity of GNSS services. The reason is that a system anomaly/loss of signal in aviation can jeopardize safety, or cause legal or economic liabilities.
In order to guarantee reliability, a GNSS signal integrity monitor is integrated with avionics. Integrity is considered an essential performance metric in aviation. It is a measure of trustworthiness of a system [10]. An integrity monitor is designed to issue alarms to pilots within a prescribed time, when reliable GNSS performance cannot be ensured. To this end, it carries out continuous fault detection tests and raises an alarm in case a fault is detected. If no fault is identified, the monitor outputs be noted that constellation-wide faults due to a common cause are not considered in the current work. Measurements of a given epoch are also assumed to be uncorrelated with each other.
Extensive simulations over the primary coverage area of NavIC on the Indian sub-continent were performed for a 20-min UAV trajectory. With GPS and NavIC measurements, it is shown that the existing methods used as reference outperform the new one. However, the work is significant for the following reasons. KF is an integral part of advanced navigation systems, as noted earlier. Therefore, developing KF RAIM algorithms is essential to ensuring their reliability. The improvement of KF solution separation RAIM is obtained at the expense of high computational load. The presented range-based KF RAIM, on the other hand, offers satisfactory performance to a certain extent with moderate computational resources. It has potential for extensions to advanced methods such as vector receivers, where running parallel vector tracking loops for solution separation RAIM would be too complex and prohibitively expensive. Further, simulation results indicate that the addition of NavIC alongside GPS can substantially improve performance, particularly in poor geometries.

Overview of Existing KF RAIM Algorithms
The developed KF RAIM is compared with two existing integrity algorithms. The first one is range-based RAIM with WLS [20]. Although snapshot in nature, it is chosen for performance comparison because it works in the range domain, like the presented method. The second one, KF RAIM in the position domain, implements the solution separation approach with the EKF. Since the WLS RAIM algorithm is described in detail in [20], it is not repeated here. First, a brief overview of existing solution separation KF RAIM is provided. Following this, background of the existing range-based KF RAIM, which is extended further later in this paper, is presented for ease of understanding.

Solution Separation KF RAIM
First, error models and key integrity parameters are described. Although these are discussed in the context of solution separation method, these are used by all RAIM algorithms of this paper. The following error models are considered. Ionospheric error is assumed to be removed using dual frequency measurements. Broadcast satellite clock and ephemeris error, residual tropospheric delay (i.e., remaining error after removing the modeled component), noise and multipath in pseudorange measurements are modeled as a first order Gauss Markov (GM) process with 100 s time constant. Their standard deviations are taken from [44,45] for carrier smoothed code pseudorange measurements of GPS in L 1 and L 5 bands for airborne users. More realistic error models with different time constants for individual components will be considered in future work. Due to lack of availability of similar error models for NavIC, the same equations as those for GPS are assumed with frequency bands appropriately changed. Only white noise is modeled in the pseudorange rate measurements. Noise standard deviation corresponds to a C/N 0 of 45 dB-Hz.
The prior probability of narrow fault P F for both GPS and NavIC constellations is assumed to be 10 −4 /satellite/h. Over the years, P F is found to be less than this value for GPS [42]. For NavIC being a new constellation, many years of data are yet to be analyzed before arriving at a suitable probability. Constellation-wide faults are not considered in this work. The allocated probability of hazardously misleading information (HMI) (i.e., the probability that position error exceeds the PL but no alert is issued), P alloc HMI , is considered 10 −7 /h [46]. P alloc HMI for the vertical component, P alloc HMI, V is 9 × 10 −8 /h (i.e., 90% of P alloc HMI ). P alloc HMI for the horizontal component, P alloc HMI, H is 1 × 10 −8 /h. For the assumed P alloc HMI (or allocated integrity budget) and P F , monitoring up to two simultaneous independent satellite failures is sufficient, as discussed in [24]. This is because three or more simultaneous failures will have a probability less than the total integrity budget. The probability of three or more unmonitored faults P um is therefore subtracted from the integrity budget. The continuity budget allocated to disruptions due to false alert P FA is 3.33 × 10 −7 /sample [47]. P FA and P alloc HMI values chosen correspond to the en route phase of civil aviation. where h = 1 and 2 denote north and east components, respectively. NavIC constellation currently has seven operational satellites in orbit. Assuming a maximum of twelve GPS satellites, the sum of the preceding probabilities is slightly less than half of (P alloc HMI, H − 0.1P um ). For the vertical component of position, the allocated probabilities are P alloc, n f HMI, 3 = 3.6 × 10 −10 /h, P alloc, s f HMI, 3 = 3.6 × 10 −10 /h, P alloc, d f The sum of the above three probabilities is a little less than (P alloc HMI, V − 0.9P um ). The probabilities are not optimally chosen to minimize PLs due to higher computational load involved in optimization. As an alternative, a larger value (about 99% of the total probability) is allocated to dual-satellite fault mode instead of assigning one third of the sum to each of NF, SF and DF. This is to reduce PLs as the DF mode is found to result in higher PLs.
Position domain KF RAIM implemented in this paper adapts solution separation approach to EKF. It runs parallel EKFs-one for each fault mode plus an all-in-view filter for NF condition. Each EKF estimates errors in three coordinates of aircraft position and velocity in the Earth centered Earth Fixed (ECEF) frame, and receiver clock biases and drifts for GPS and NavIC from their respective a priori values. In addition, each of them has additional states (one for each pseudorange measurement) to model time correlated errors as first order GM processes. Key equations of KF solution separation RAIM are provided in Appendix A. Satellite position, velocity and modeled tropospheric delays are calculated only in the all-in-view filter and shared among all subset filters-one for each fault mode-to reduce the computational load. State error covariances are updated using a numerically stable forward recursive modified rank-one update algorithm [43]. Next, an overview of existing range-based KF RAIM is provided to lay foundation for the modified algorithm described later in this paper.

Existing Range-Based KF RAIM
In the fault detection method of existing range-based RAIM, three test statistics are formed, each with KF measurement innovations of a number of epochs [33]. The reason for three tests is described in detail in the reference. The first statistic, α 1 , is calculated with measurement innovations of M epochs including the current epoch t k . The second statistic, α 2 , has (N − M) epochs from k − N + 1 to k − M. The third test statistic, α 3 , is formed with innovations of all epochs before N epochs, whose contributions are included in α 3 before these terms are discarded (or removed from memory). At a time, only N epoch terms are retained in memory, but test statistics are formulated with innovations of all epochs. N and M are adaptively determined. Mathematically, the test statistics at t k , (α j, k ; j = 1, ..., 3), are given by where ∆ρ and ∆ρ denote EKF pseudorange and pseudorange rate innovation vectors, respectively, at epoch . m j and p j hold different values for different j, as mentioned before. I is the identity matrix of size n × n; and n is the number of visible satellites at an epoch. n can change with time. G and G rr are computed by WLS estimation. They are where W is the pseudorange measurement error covariance matrix at t th epoch. H is linearized pseudorange (or pseudorange rate) measurement model matrix. W rr is the pseudorange rate error covariance matrix. The error model accounts for unmodeled propagation delays by inflating noise covariance, and neglects error in broadcast satellite clock terms and ephemeris, and multipath. Measurement error between time epochs is thus assumed uncorrelated. Fault detection method considers WLS estimation whereas PL calculation is based on EKF position error. A WLS estimation-based fault detection test with EKF innovations allows simplified calculations of mean position error bounds, as discussed in the reference. With WLS estimation-based fault detection, each term of the test statistic under summation in Equation (3) and in the denominator of failure mode slope (FMS) would contribute fault of only the corresponding epoch, not any other epochs. This simplifies FMS calculation, and allows the formation of a separate FMS for mean position error bound pertaining to a test statistic. As a result, more than one test can be formulated. The noise statistics of αs can also be easily computed. It was proved in [31] that each term in Equation (3) generated with KF innovations has the same noise distribution as that of the WLS navigation filter. The threshold T th, j for an α j is therefore computed from the Chi square distribution with probability P FA /3 and appropriate degrees of freedom (DOF). A fault is declared if at least one of the test statistics crosses its threshold.
The bounds on horizontal and vertical mean position error under a fault, HPE U i and VPE U i , respectively, for a satellite i are derived as where HPE U i, j and VPE U i, j are bounds corresponding to α j (i.e., they consider the same epoch terms as those of α j ), with j = 1, 2, 3. µ error, U n (1 : 3) is a three-element vector of north, east and down mean position error bounds resulting from second order terms of nonlinear GNSS measurement equations.
It is important to note that the DOF of α 3 , and HPE U i, 3 and VPE U i, 3 grow with time. However, as HPE U i, 3 and VPE U i, 3 have very small contributions, VPE U i and HPE U i do not have the same growing trend over time. It is justified through simulation results that with three test statistics and mean position error bounds as calculated in the preceding equations, PL does not grow over time, unlike the PL with a single test statistic. Faults are also detected in a shorter time than that with a single test formed using all innovation terms from the beginning up to the current epoch. The reason is that the growing DOF of a single test statistic increases the fault detection threshold fast, resulting in late detection of faults. PLs are computed using the above-mentioned mean position error bounds under the assumption of a single satellite fault.
As noted before, SKF is used as navigation processor for range domain KF RAIM algorithm developed in this paper. Its formulation is briefly provided next.

Schmidt KF
In this section, the SKF formulation of KF is briefly outlined. Interest readers can refer to [43,48] for more details. The reason for using SKF will be explained in the next section while discussing the fault detection test of range-based KF RAIM. Although sub-optimal in nature, SKF accounts for the effects of time correlated errors without estimating additional states. It also has fast execution time.
It is desired to determine aircraft absolute position and velocity, and receiver clock biases and drifts for GPS and NavIC. Their true values and a priori estimate are denoted as κ andκ − , respectively. The difference between κ andκ − , ∆κ k is estimated. Thus, the SKF state vector at t k is where x, y and z are three position components of aircraft expressed in ECEF coordinate frame.ẋ,ẏ andż are velocity components in ECEF. b G clk and b N clk are the receiver clock biases for GPS and NavIC, respectively, andḃ G clk andḃ N clk are the receiver clock drifts for the two constellations. Broadcast ephemeris and clock error, residual tropospheric delay along the line of sight, multipath and noise in carrier smoothed code pseudorange measurements of a satellite are modeled as a first order GM process with time constant 100 s. Let σ i ura , σ i tr , σ i mp and σ i n be the standard deviations of broadcast ephemeris and clock error, residual tropospheric delay, multipath and noise, respectively, for ith satellite. Their elevation angle dependent models for airborne users are taken from [44,45] for GPS L 1 and L 5 band measurements and provided next.
where θ i is the elevation angle of ith GPS satellite in radians. σ ura is assumed to be 0.75. Thus, the resultant error is a GM process with 100 s time constant, and has a standard deviation of (σ i ura ) 2 + (σ i tr ) 2 + (σ i mp ) 2 + (σ i n ) 2 for ith GPS satellite. The same error models are used for NavIC, but the frequencies are changed to L 5 and S bands.
Pseudorange rate measurement error is modeled as white Gaussian noise with standard deviation corresponding to a C/N 0 of 45 dB-Hz. Thus, there are n GM processes, where n is the number of visible satellites at an epoch. In the SKF algorithm, states of GM processes are not estimated, but their effects are considered. If the vector p gm, k−1 has the states of n GM processes at t k−1 , then ∆κ k−1 and p gm, k−1 can be combined as X k−1 = [∆κ T k−1 p T gm, k−1 ] T . The process and linearized measurement models of X k and measurement update equations of SKF are described in Appendix B. Next, the range-based KF RAIM algorithm is presented.

Range-Based KF RAIM Algorithm
First, the fault detection method in the presence of time correlated errors is developed. Then, mean position error bounds and PLs are computed considering up to two simultaneous independent satellite failures.

Fault Detection Method
The test statistic described by Equation (3) earlier is modified to account for time correlated errors in pseudorange measurements, as explained next. Since no time correlated error is modeled in the pseudorange rate measurements, part of the test statistics with pseudorange rate innovations remains the same as before. It can be shown that where ∆ρ * T = W −1/2 ∆ρ T and G * = W −1/2 H (H T W −1 H ) −1 × H T W −1/2 where '×' is matrix multiplication. Using the idempotent property of (I − G * ) [18] With singular value decomposition, (I − G * ) can be written as L is a diagonal matrix. Its diagonal elements are ones and zeros. The number of zeros is equal to the dimension of the null space of (I − G * ). For a single constellation it is four whereas for dual constellations with a separate receiver clock bias state for each, it is five. Let the null space dimension be denoted as n null . It can be shown that (I − G * ) is also equal to where J(1 : n 1 , 1 : n 2 ) represents first n 1 rows and n 2 columns of a matrix J. Using Equations (12) and (13), and replacing each (I − G * ) in Equation (13) with Equation (15) and S (1 : n, 1 : n − n null ) with S , part of the test statistic formed with pseudorange innovations, ∑ where ∆ρ * stacks all ∆ρ * vectors with varying from the lower limit p j to the upper limit m j of Equation (3).S is a block matrix having non-zero sub-matrices (or blocks) along the diagonal, with the th block given by S . Θ is a square matrix with number of rows/columns equal to ∑ m j =p j n. It is introduced in the middle and defined later. When it is the identity matrix, Equation (16) reduces to Equation (3). For ease of representation, the time epoch k is dropped from the subscript. The mathematical model of ∆ρ * is given by H * is a block matrix with non-zero blocks along the diagonal. Its th block is W −1/2 H . H is defined after Equation (5). Only difference is that H has an additional column for receiver clock bias corresponding to the second constellation. ∆κ p stacks all [∆x ∆y ∆z ∆b G clk ∆b N clk ] T vectors for all noted after Equation (16). η 1 is measurement error. The diagonal terms of its covariance matrix, Ψ, are unity.S S T present on either side of Θ in Equation (16) removes H * ∆κ p from ∆ρ * , ensuring that the fault of only th epoch is injected through ∆ρ * . This enables formation of more than one test statistic, each with a certain number of epochs, as before. If a term from one epoch also injects faults of other epochs, formation of separate test statistics is difficult from the point of view of the corresponding FMS computation, as discussed earlier. This would be the case ifS S T were omitted from Equation (16).
whereS T ΘS (or Θ 1 ) has ∑ m j =p j (n − n null ) rows/columns. The covariance matrix Ψ 1 of error in ∆ρ , η 2 is given by where E is the expectation operator. Ψ is the covariance matrix of measurement error in ∆ρ * , η 1 . Measurement error comprises broadcast ephemeris and satellite clock error, unmodeled tropospheric delay, multipath and noise. It is assumed that the pseudorange measurement error of each visible satellite is a first order GM process with time constant 1/β and of zero mean under no faults.
Since time correlated errors are not estimated by the SKF, they remain unchanged in the pseudorange innovation terms. In this context, it should be noted that if time correlated errors were estimated by augmenting states in EKF, it would not be possible to formulate the RAIM algorithm described here. This is because if the augmented states for time correlated errors (one for each pseudorange innovation) were included in the WLS estimation of fault detection method, the resulting system would be under-determined. If the augmented states of EKF were not included in WLS estimation, faults of other epochs would be injected through the pseudorange innovation term of an epoch. This is because of possible wrong estimation of augmented states in EKF in the presence of a fault. This would complicate the calculation of FMS, and a separate FMS for each test statistic would not be possible. The reason for preferring a WLS-based fault detection approach was discussed in Section 3.2 and before Equation (18). Therefore, implementation of the SKF is essential to the formulation of the range-based KF RAIM algorithm presented in this paper. SKF is computationally efficient, and provides reasonably good performance as compared to that of the EKF that estimates additional states for time correlated errors [43]. Next, Ψ and Θ are formulated for the first test statistic with M epochs (see paragraph before Equation (3) for the definition of M).

Formulation of Ψ for First Test Statistic
First, Ψ is determined. Expanding η where the superscript within parentheses denotes time instants, and the argument indicates pseudorange measurement index (or satellite index) at a time instant. Total number of visible satellites, n can vary across epochs. By definition, where [0] is a matrix with zero elements. Its dimension is shown in the subscript. For any satellite i, at epochs and + q, with q > 0 where η is measurement error in ∆ρ (which stacks relevant ∆ρ terms), and σ denotes standard deviation of its individual element. Using GM model, one can write where ∆t = measurement update interval = 1 s. ζ ( ) i is white noise sample for satellite i and uncorrelated with η ( ) (i). Simplifying, Substituting this into Equation (21) The term can be written as It is a product of q ratios. Each ratio will have an accompanying exp (−β∆t) when exp (−β∆t) q is multiplied with the product. Using Equations (9)-(11), Figure 2 shows how the upper and lower bounds of exp (−β∆t) × σ ( ) (i) σ ( +1) (i) (=a 1 ) vary for GPS satellites with elevation angles θ ( ) at epoch . It is evident from the equations that maximum change in elevation angle between two successive epochs would maximize a 1 , which is also higher when θ ( +1) > θ ( ) . The maximum change in elevation angle between two epochs of 1 s interval is calculated in radians as follows.
where v s max is the maximum satellite speed, v u max is the maximum user speed (assumed Mach 1) and r su min is the minimum range between a satellite and user at an elevation angle. The maximum altitude of the user is 15 km. When θ ( +1) < θ ( ) , following the same procedure, lower bounds of a 1 are obtained. The term a 1 calculated for GPS is found to be greater and less than that of NavIC with respect to upper and lower bounds, respectively. This is because NavIC satellites are located in geosynchronous orbits. As shown in the figure, the largest value of a 1 , a max is 0.9911 at 5 • mask angle. The lowest value of a 1 , a min is 0.9890. A matrix Ψ max is modeled as The determinant of Ψ max can be easily calculated by the product of the determinants of all sub-matrices-each formed by grouping together the non-zero elements for a particular satellite. This is because measurements are assumed uncorrelated across satellites. Thus, with some rearrangements of rows and columns, it can be shown that Ψ max is a block diagonal matrix formed by these sub-matrices.
Its dimension is determined by the number of time instants when the satellite is visible within M epochs. If it is visible throughout, then the determinant of the sub-matrix is (1 − a 2 max ) (M−1) . It is greater than zero, and approaches zero as M increases, but theoretically is not equal to zero for a finite M. All leading principal minors are also positive. Thus, each sub-matrix is positive definite, and has all positive eigenvalues. Combining all sub-matrices, one can say that Ψ max is positive definite.
In case of Ψ, following Equations (26) and (27) , where a 11 , a 21 , and a 31 are values of a 1 at epochs t = t 1 , t 2 and t 3 , respectively. Its determinant is ∏ ). Using the same argument as earlier, Ψ is positive definite. SinceS has full column rank, Ψ 1 (=S T ΨS ) is also positive definite. This implies that there exists a matrix Ψ 2 such that For a 1 between [a min , a max ], it is evident from the mathematical expressions that the determinant where all diagonal terms of ∆Ψ i sub are zero, and off-diagonal terms are positive ∀ a 1 = a max . This is equivalent to perturbing Ψ i sub so as to make it move closer to singularity. Hence, the condition number of Ψ i sub, max (ratio of maximum to minimum eigenvalue for a symmetric matrix) would be greater. Figure 3 justifies this for different M by generating Ψ i sub 10 4 times for each M with uniformly random a 1 between [a min , a max ], and comparing its eigenvalues with those of Ψ i sub, max . From the figure, it can be concluded that the minimum (or maximum) eigenvalue of Ψ i sub, max is less (or greater) than the respective eigenvalue of Ψ i sub , when a 1 s are not equal to a max . Thus, the maximum and minimum eigenvalues of Ψ can be related to those of Ψ max as where λ J represents an eigenvalue of a matrix J and its superscript shows if it is maximum or minimum. Equation (30) will be useful while discussing the second test statistic.

Formulation of Θ for First Test Statistic
Next, Θ is determined so that the fault detection threshold of α j can be obtained. For ease of understanding, (α j ) 2 and the mathematical model of ∆ρ * are repeated here.
SinceS S T removes H * ∆κ p from ∆ρ * , the distribution of (α j ) 2 is determined by the measurement error η 1 . Hence, only η 1 of ∆ρ * will be considered in subsequent derivations. Using measurement error, Equation (31) results in Putting Reference [49] proves that the fault detection threshold of the preceding term can be obtained from the Chi square distribution if the maximum eigenvalue of Θ 2 is less than or equal to one. Θ is determined next to satisfy this. If λ Θ, 2 is an eigenvalue of Θ 2 = Ψ 2 T S T ΘS Ψ 2 , and v is the corresponding eigenvector, then Multiplying Ψ 2 on either side and denoting Ψ 2 v as v 1 Replacing Ψ 2 Ψ 2 T with Ψ 1 and S T ΘS with Θ 1 , it is clear that the eigenvalues of Ψ 1 Θ 1 (denoted by Ω later) are identical with those of Θ 2 . Therefore, in subsequent derivations, the eigenvalues of Ω will be considered.
The following form of Θ is proposed.
where c is a positive scale factor, which is determined such that the maximum eigenvalue of Ω is less than or equal to unity. From the preceding equation it is evident that for a positive value of c, Θ is a positive definite matrix. Since bothS and Ψ 2 have full column rank, Ψ 2 T S T ΘS Ψ 2 or Ω (=Ψ 1 Θ 1 ) have all positive eigenvalues. Next, it will be proved that there is a value of c, for which the maximum eigenvalue of Ω, λ max Ω is less than unity. If ||.|| 2 is the 2-norm of a matrix It should be noted that ||S || 2 is unity. Since both Ψ and Θ are symmetric, their maximum eigenvalues are equal to respective 2-norms. Hence, λ max Ψ = ||Ψ|| 2 and λ max Θ = ||Θ|| 2 . An expression of λ max Θ is derived next. If λ J represents an eigenvalue of a matrix J, then it can be easily shown that For matrix Θ = (cI + Ψ) −1 , λ Θ = 1/(c + λ Ψ ). Thus, the maximum eigenvalue of Θ, λ max where λ min Ψ is the minimum eigenvalue of Ψ. Therefore, If c is chosen equal to (λ max Ψ − λ min Ψ ), λ max Ω is less than unity. However, the above c results in a large FMS corresponding to the first test statistic. Since the FMS of the first test statistic significantly contributes to determining the PL, it should be reduced by selecting a small c. A computationally efficient algorithm to adaptively find Θ for α 1 at every epoch is provided in Appendix C. Ψ is recursively computed at every t k with Ψ i sub for all satellites. It should be noted that even after a satellite sets, its measurement innovations are used to form test statistics as long as they are within N epochs from the current epoch.

Formulation of Θ for Second Test Statistic
Ψ and Ψ max for the second test statistic, α 2 can be formulated the same way as before, but for (N − M) epochs. As the FMS corresponding to α 2 is generally small, Θ for it is calculated as where c is provided after Equation (40), but calculated with the minimum and maximum eigenvalues of Ψ max . In this case, Substituting for c and using Equation (30), it is evident that λ max Ω < 1. Θ can be pre-determined for various (N − M). Since Ψ max can be put in the form of a block diagonal matrix, it is sufficient to pre-calculate the inverse of (cI + Ψ i sub, max ) for various (N − M). Ψ i sub, max and its pre-calculated inverse are loaded into the program. With this, Θ is appropriately formed at every t k for all visible satellites over epochs t k−N+1 to t k−M . If a satellite q is visible for a shorter time during this interval, then part of Θ corresponding to Ψ q sub, max is recursively calculated starting from (cI + Ψ i sub, max ) −1 , following [50]. Satellite i is assumed to be visible for the longest duration over (N − M) epochs. Hence, λ min (or max) Ψ, max = minimum (or maximum) eigenvalue of Ψ i sub, max . Next, the final form of test statistics and their thresholds are discussed.

Test Statistics and Thresholds
In this work, the third test formulated with all discarded terms of epochs before t k−N+1 is not performed. The following subsection computes a bound on the contributions of the discarded terms to mean position error without requiring an FMS (or a corresponding test). Thus, Equation (3) is modified as (43) where j = 1, 2. SinceS S T is a block diagonal matrix with each block from an epoch, elements of the vectorS S T ∆ρ * corresponding to an epoch are calculated at that epoch and then stacked along with those of other epochs. This eliminates unnecessarily large matrix multiplications.
With the designed Ψ and Θ, fault detection thresholds, T 2 th, 1 for α 2 1 and T 2 th, 2 for α 2 2 , can be determined from the Chi square distribution with probability P FA /2 and DOF = ∑ where n null is the dimension of the null space of (I − G * ). The number two is multiplied for the second term of Equation (43). If either α 1 exceeds T th, 1 or α 2 exceeds T th, 2 or both happen, a fault is declared. Calculation of the mean position error bounds is detailed next.

Mean Position Error Bounds under Faults
Changes to mean position error bound computation mentioned in Section 3.2 are explained in this subsection. Since the contribution of nonlinear measurement model (see Equations (6) and (7) and discussion after them) was found to be very small, it is not considered in this paper. Next, the three terms under the summation of Equation (7) are discussed for the vertical position error (VPE). Terms for horizontal position error (HPE) are changed accordingly. The mean position error, µ k at time t k in the north-east-down (NED) frame is N terms corresponding to α 1 & α 2 (44) where µ k = C n e µ k (1 : 2 : 5), µ k is the mean error of all estimated states and (1 : 2 : 5) represents elements for the position states, B k−j = C n e B k−j (1 : 2 : 5, :), B k−j = ∏ k−j =k A K s, k−j−1 , A = (I − K s, C s, )F s, −1 . K s , F s and C s are Kalman gain, state transition and measurement model matrices of estimated SKF states, respectively. They are defined in Appendix B. K s, k = C n e K s, k (1 : 2 : 5, :), where J(1 : 2 : 5, :) has rows one, three and five of matrix J. C n e is the ECEF to NED frame coordinate transformation matrix. N terms corresponding to two test statistics and discarded terms are separately grouped with under braces. Fault is assumed to have started at an epoch m and denoted by vector f k−m+1 at t k . Up to two independent satellite faults are considered. As noted earlier, three or more faults are not monitored. Thus, for any two satellites s 1 and s 2 ,  [20] corresponding to the first test statistic α 1 with terms of M time epochs from k to k − M + 1 is given by where f w , the worst-case fault vector, maximizes the FMS. A pictorial depiction of FMS is provided in Figure 4. For a DF mode in satellites s 1 and s 2 , L is where '×' represents matrix multiplication operation. Although it is not evident, the number of visible satellites, n can vary across epochs. f w has the corresponding b 1 ,ḃ 1 , b 2 andḃ 2 terms from Equation (44). For a SF mode, elements for only satellite s 1 are considered.
It should be noted that for the HPE, first and second rows of B and K s are used to form L. f w is composed of the corresponding b 1 ,ḃ 1 , b 2 andḃ 2 terms. Λ of Equation (45) is given as where Θ 3 is formed with terms of W −1/2S S T ΘS S T W −1/2 of Equation (43) corresponding to faulty satellites of a fault mode. W is a block diagonal matrix whose th block is W , and ∆ρ * = W −1/2 ∆ρ. Θ 4 is a diagonal matrix with diagonal entries constituted from relevant terms of W rr, −1 (I − G rr ), = k (=p 1 ), k − 1, ..., k − M + 1 (=m 1 ), of Equation (43). (maximum FMS) α1, i is given by the square root of the maximum eigenvalue of LΛ −1 , λ max LΛ −1 . Since Θ 4 is a diagonal matrix, Λ −1 requires the inverse of an M × M matrix for SF cases and 2M × 2M matrix for DF cases. Thus, the upper bound of mean VPE corresponding to α 1 for fault mode i is where Γ = f T w Λ f w . It is determined next corresponding to the allocated probability of missed detection for fault mode i, P alloc MD, i . P alloc MD, i for SF and DF modes is P alloc MD, i = P alloc, s f HMI, 3 /(n max P F ) (for SF modes) P alloc MD, i = P alloc, d f HMI, 3 /(0.5n max (n max − 1)P 2 F ) (for DF modes) (49) where n max is the maximum number of visible satellites. Using Equation (35), and putting pseudorange rate measurement error, η 4 in place of pseudorange rate innovations, ∆ρ for the same reason as that discussed before Equation (33), Equation (43) for α 1 is written as Replacing W −1/2 rr, η 4, with η 5, , and W rr, Similar to Equation (15) (I − G rr, * ) = S rr (1 : n, 1 : n − n null )(S rr (1 : n, 1 : n − n null )) T Replacing S rr (1 : n, 1 : n − n null ) with S rr, , and substituting Equation (52) into Equation (51) is the matrix square root of Θ 2 , and stacking (S rr, ) T η 5, for all relevant values of in η With known error covariance of pseudorange rate measurements, η 7 has unit-variance independent elements. Hence, its covariance is identity matrix. Since η 3 also consists of unit-variance independent elements (see before Equation (34)), the covariance of η 6 is Θ 2 . It was shown in the previous subsection that the eigenvalues of Θ 2 lie between zero and one. The error covariance of [η 6 T η T 7 ] T (=η) denoted as Θ is a block diagonal matrix. The first block is Θ 2 , and the second block is identity matrix. Thus, the maximum eigenvalue of Θ, λ max Θ is unity. The minimum eigenvalue λ min Θ is greater than zero and equal to the minimum eigenvalue of Θ 2 . Although not considered in this paper, unknown error covariance of pseudorange rate measurements can be accommodated following [49]. If the mean ofη due to faults is represented asf , then the probability of missed detection P MD is given by where N stands for multi-variate Gaussian distribution, and V is the variable of integration. The square of magnitude off ,f Tf , is Γ, which is the same as f T w Λ f w under fault mode i (see after Equation (48)). Let the non-central Chi square distribution P ncx ((T * th, 1 ) 2 , DOF, Γ * ) with non-centrality parameter Γ * and for some threshold (T * th, 1 ) 2 have a probability P 1 , then P MD P 1 if [49] (T * th, If Γ * is obtained such that P ncx ((T * th, 1 ) 2 , DOF, Γ * ) = P alloc MD, i , DOF = number of rows/columns of Θ, and Γ is calculated from Equation (56) is generally found to be higher than M. Therefore, for the calculation of Θ 3 of Equation (47) In simulation studies with realistic error models for airborne users, this is not found to be a time intensive operation. This will be evident later with algorithm execution times. Second, Θ is pre-calculated using Ψ max instead of the actual Ψ to avoid forming a large Ψ over N − M epochs. λ min Θ of Equation (56) is also calculated with Ψ max in the place of Ψ. In order to ensure that the corresponding Γ is overbounded, the following is justified in Appendix D. For a given probability P 1 , Γ increases with decreasing λ min Θ , and the calculated λ min Θ is lower than the actual value obtained with Ψ.
The third bound VPE U 3 accounts for the contributions of all discarded terms before epoch t k−N+1 to the mean position error. Since it remains the same for all fault modes, the fault mode index i is dropped from its subscript. It is presented next.

3
As noted earlier, the third test with discarded terms has been eliminated in the current work. However, discarded terms may have an effect on current mean position error at t k , which should be bounded by VPE U 3 . In this context, a pertinent question arises as to which faults should be considered for the bound. Suppose at a time epoch t , f T w Λ f w corresponding to a fault detection test is equal to Γ, which is determined from P alloc MD, i under fault mode i. Hence, for any fault f under the same fault mode such that f T Λ f > Γ, the probability of missed detection is less than P alloc MD, i at t . Positinon error due to this fault is not bounded at t . It also need not be bounded at subsequent epochs. This is because if HMI occurs at t due to this fault with missed detection probability less than P alloc MD, i , then bounding its contribution to position error at subsequent epochs holds little meaning. It can be explained with the help of the relative frequency approach of the probability theory. Thus, VPE U 3 at t k bounds the contributions of faults with probability of missed detection more than the allocated probability corresponding to both α 1 and α 2 and all fault modes over all epochs from t 1 up to t k−N . It should be noted that VPE U 2, i , if obtained in a way similar to that of VPE U 3 , would be too large and inflate the PL to a great extent. Therefore, VPE U 2, i is not determined this way. VPE U 3 is [33] For the HPE, the factor 1 × 10 −6 is modified to √ 2 × 10 −6 to include north and east components. N max is the maximum of all values of N up to t k . N at an epoch is determined in the following way for an asymptotically stable filter. At an epoch t k , the last B matrix (see after Equation (44) for definition) with 2-norm larger than or equal to 1 × 10 −6 is identified. If the magnitude of N at t k−1 is more than the index of that B matrix, N is set to that index plus one. One is added to account for the new term at t k . All B matrices with indices higher than updated (N − 1) are not considered anymore. On the other hand, if the index is the same as N at t k−1 , N is only increased by one to include the new term at t k . All terms including B matrices from epochs before t k−N+1 are discarded or deleted from memory. Note that a small threshold 1 × 10 −6 is chosen to calculate N. A smaller value can also be chosen, but that will increase N. U max is an upper bound such that a discarded B matrix cannot have a 2-norm outside the range of U max × 10 −6 and zero after it is discarded. A computationally efficient algorithm to find U max at each epoch is provided in [33]. The denominator of Equation (57) results from an infinite geometric series with common ratio 0.1. Next, the maximum values of bias and bias rate,b max and˙b max , respectively from faults mentioned before Equation (57) are calculated differently in this paper as follows.b max = Γ max /s min ,˙b max = 2 Γ max /s min (58) Each term of the preceding equation is obtained next. Γ max is the maximum of all Γs calculated for SF and DF modes corresponding to α 1 and α 2 over all epochs up to t k . Using f T w Λ f w = Γ, and the fact that Θ 4 of Λ in Equation (47) is a diagonal matrix, s min is where the operator "min" denotes minimum over all terms provided below it or within parentheses next to it, and diag(.) has all diagonal elements of a square matrix. It should be noted that the multiplication of two in˙b max of Equation (58) accounts for DF modes (i.e., faults in any two pseudorange rate measurements). In order to determineb max , Θ 3 of Λ in Equation (47) is used. If f ρ represents faults in only pseudorange measurements for a fault mode, then one can write Multiplying and dividing f T ρ f ρ on the left-hand side Noting that the minimum value of It is apparent that forb max , the two need not be multiplied as f T ρ f ρ includes DF. Next, the PLs are determined.

Protection Levels (PLs)
The VPL is determined in this subsection without assuming statistical independence between test statistics and position error. The HPL can also be calculated the same way. The probabilities of HMI for the vertical position under NF, SF and DF are P Thus, VPL NF is obtained by setting P(( v > VPL NF )|NF) equal to P alloc, n f HMI, 3 . For an SF mode, P s f HMI, 3 is given by where n s f is the number of SF modes, VPL SF is the VPL under SF, and P F is the prior probability of a satellite fault defined in Section 3. Representing the probability of missed detection for all SF modes as P MD, s f and following the same approach as before Γs for α 1 and α 2 are determined such that P(α 1 T th, 1 |SF) ≤ P alloc MD, s f , and P(α 2 T th, 2 |SF) ≤ P alloc MD, s f , respectively (see after Equation (56)). P alloc MD, s f is equal to VPL SF is where the operator "max" denotes maximum over all terms provided below it or within parentheses next to it. P( v > γ s σ v ) = P alloc MD, s f . σ v is the error standard deviation of the estimated vertical position.
Similarly for DF modes, the probability of missed detection, P MD, d f is The corresponding Γs are also calculated for DF modes. The VPL for DF modes, VPL DF is

Simulation Studies
Primary coverage area of NavIC over the Indian sub-continent is chosen for simulation studies in this paper. A UAV is assumed to fly at an altitude of 184 m above the WGS84 reference ellipsoid at speed 10 m/s with a constant heading of −45 • for 20 min. UAV speed, altitude and flying time chosen are commensurate with those of small unmanned aircraft [51]. The starting point of the flight path is varied from 0 • to 40 • N latitudes and 60 • E to 105 • E longitudes in steps of 1 • . There are a total of 1886 simulations, each of duration 20 min. Every 50 simulations, the start time is increased by 7.5 h. At a time, maximum twelve GPS satellites are simulated due to channel constraints of the simulation platform. Dual frequency GPS (L 1 and L 5 ) and NavIC (L 5 and S) measurements are generated for the UAV flight path. Error models of carrier smoothed pseudorange measurements are discussed in Section 4 while describing SKF. White noise covariance of pseudorange rate measurements is obtained from the tracking error of a second order frequency locked loop (FLL) for 45 dB-Hz C/N 0 . The FLL is assumed to have a cross product discriminator, a coherent integration time of 10 ms and two-sided noise equivalent bandwidth of 4 Hz. Power spectral density of a temperature controlled crystal oscillator is considered for simulation of receiver clock error. Measurement update interval of both WLS and KF is one second.
VPLs and HPLs of the developed range-based KF RAIM algorithm are illustrated in Figure 5. Both GPS and NavIC constellations are considered. The geometric dilution of precision (GDOP) varies between 0.5 and 3.5. It is within 2.5 for 97.6% of the time as against 73.3% with GPS alone. For illustration purposes, PLs calculated every second are grouped into five bins, whose center points are indicated along the x-axis. Percentage of time along the y-axis is calculated as the ratio of the number of time instants at which the PL falls within a given range to the total number of time instants over all 1886 simulations times one hundred. KF VPL and HPL are within 40 m 97.16% and 99.35% of the time, respectively.
Evolution of PLs with time are also plotted in Figure 5 for specific GPS plus NavIC and GPS-only cases. For clarity, only two scenarios are selected. They have highest and lowest GDOPs corresponding to dual constellations. With the lowest GDOP, PLs are reduced by about 5 m to 15 m when both constellations are used as compared to PLs for GPS alone. For the highest GDOP, PLs are hundreds of thousands of meters with standalone GPS (having five-six visible satellites). This is shown separately in Figure 6 along with PLs of KF solution separation RAIM, with the assumption of a single satellite fault. Although not shown, a few other high GDOP cases are studied and found to have similar performance with GPS. Thus, it can be concluded that improvements obtained by using NavIC alongside GPS are more significant in poor geometries. M and N terms of KF range-based RAIM for a simulation are illustrated in Figure 7. They are of the same order in all simulations.   Having shown the PLs of range-based KF RAIM, its performance is compared against WLS range-based RAIM and KF solution separation RAIM using both constellations. Figure 8 depicts comparison results. Range-based KF RAIM has VPL within 3 to 11 m of that of range-based WLS RAIM 90.4% of the time and HPL within −1 to 11 m 98.82% of the time. Its VPL, on the other hand, is within −3 to 11 m of that of its solution separation counterpart 95.49% of the time, and HPL is within −3 to 11 m 97.28% of the time. Using a vertical alert limit of 50 m and horizontal alert limit of 40 m, availability of the three algorithms is shown in Table 1. Availability is calculated as the percentage of simulations, in which the PL is below the corresponding alert limit for the whole duration of 20 min. PL depends on both P FA and P HMI , for which no proper values have yet been fixed for UAV applications. Thus, availability can only be treated as another way of comparing algorithm performances. It is apparent from the table that the two existing algorithms offer higher availability than range-based KF RAIM. It is observed that its PL sometimes exhibits an instantaneous change when a satellite sets or rises. If such a change can be smoothed out by some method, its availability would increase to 98.46% in the vertical direction and 98.89% in the horizontal direction. This will be carried out in future work.
All 1886 simulations were run in MATLAB on a DELL core i7-8700 (CPU @ 3.2 GHz) computer with a Linux operating system, and with 16 GB RAM and 20 GB swap memory. WLS range-based RAIM took 11 h to complete all simulations. Range-based KF RAIM algorithm needed 3 days and 16 h whereas KF solution separation RAIM required 9 days. In order to see execution times of the algorithms on a low end computer, time needed by them for a simulation on an IBM Thinkpad laptop is also noted. The laptop has Linux operating system, core 2 duo processor (CPU @ 2 GHz), 1 GB RAM and 2 GB swap memory. The three algorithms needed 20 s, 600 s and 1521 s, respectively, for a 20-min (1200 s) run. Time includes that needed for running the navigation filter as well. A few other simulations also took similar times. MATLAB tic and toc functions are used to record execution times. Based on the above performance analyses, it can be concluded that WLS has the lowest execution time. However, since WLS is not preferred in advanced navigation methods, as noted earlier, KF-based RAIM algorithms are essential to ensuring reliability of such systems. Of the two KF algorithms, solution separation RAIM provides lower PLs, but has a much higher computational load. The computational burden can be reduced by fault grouping at the expense of inflated PLs [28].
Range-based KF RAIM, on the other hand, offers satisfactory performance to a certain extent with moderate computational resources, thereby having potential for real-time implementations. Its comparatively large PL is partially attributed to the fact that SKF is sub-optimal in nature for not estimating the GM states. SKF formulation is necessary for reasons discussed earlier. It allows the formation of more than one test statistic using innovation terms of a finite number of epochs and a separate FMS for each statistic. This prevents the PL from growing with time, which is considered a major challenge with KF RAIM in the range domain. There may be some scope of fine tuning the algorithm further for lower PLs, which will be attempted in the future.

Conclusions
A novel range domain KF RAIM algorithm is presented in this paper, building upon previous work of the authors. Major changes in the current algorithm include modeling time correlated errors of the pseudorange measurements and accounting for multi-satellite failures. The integrity monitoring performance of the developed algorithm is studied for a UAV trajectory with GPS and NavIC over the primary coverage area of the latter. It is shown that both WLS range-based and KF solution separation RAIM algorithms outperform the new method. However, WLS RAIM has limitations in that it is not ideal for advanced navigation systems. KF solution separation RAIM, on the other hand, offers lower PLs than the range-based KF algorithm, but with a much higher computational load. PLs of range-based KF RAIM are within 10 m of those of KF solution separation RAIM at least 95% of the time. Due to its low execution times, range-based KF RAIM shows promise for real-time implementations in avionics. It also has potential for extensions to advanced methods such as vector receivers, where running parallel vector loop filters for solution separation RAIM would be prohibitively expensive. Further, simulation results indicate that addition of NavIC alongside GPS can substantially improve RAIM performance, particularly in poor geometries.
Future extensions of the current work include extensive testing of RAIM for various geometries with GPS and NavIC over the latter's primary and secondary coverage areas for different dynamic scenarios. A suitable error model also needs to be developed for NavIC and validated with real satellite data. The error model will be used for experimental validation of the RAIM algorithms. In addition, it is important to take into account constellation-wide faults for integrity monitoring with multiple constellations. While this has been addressed in the literature for solution separation or ARAIM, a suitable methodology has to be devised for range-based RAIM. It is noted that the developed KF RAIM algorithm at times exhibits an instantaneous change in PL when satellite visibility changes. This should be smoothed out with appropriate modifications. Fine tuning of the algorithm will also be attempted to reduce PLs. Effects of different P FA and P HMI on the proposed method will be investigated. Moreover, the algorithm should be able to deal with unknown error covariances. Its performance in poor signal environments also warrants further studies. Extending the developed KF RAIM to sensor integration and vector tracking is another important area that needs to be explored in depth.
where h = 1, 2 and 3, denote north, east and down components, respectively. Q is the tail probability of zero mean unit normal distribution given by 1 √ 2π ∞ q e −q 2 /2 dq, and Q −1 is its inverse. "0" represents NF condition. i (or j) denotes a fault mode with the corresponding EKF using a subset of all measurements. The subset is determined by excluding one (or two) satellite(s) assumed faulty in the fault mode.
The threshold for a fault mode i (SF) or j (DF) is where N f ault is the total number of fault modes. PLs are calculated when test statistic | κ i/j n (h) − κ 0 n (h)| is less than respective threshold T i/j SS, h ∀ i, j. |.| represents the absolute value. κ i/j n (h) is the ith or jth subset filter estimate, with h = 1, 2 and 3 for north, east and down components of position, respectively. κ 0 n (h) denotes the corresponding all-in-view estimate. In this context, it should be noted that a change is made to Equation (A4). It is found that the (σ i/j h ) 2 − (σ 0 h ) 2 term of the equation results in the test statistics crossing their respective thresholds a few times even in no fault scenarios, which is not statistically consistent. Therefore, it is changed to (σ i/j h ) 2 + (σ 0 h ) 2 to overbound the distribution.

Appendix B. SKF Process and Measurement Models and Key Equations
The process and linearized measurement models of X k (=[∆κ T k−1 p T gm, k−1 ] T ) for SKF are where the subscript "s" is for estimated states. Subscript "p" represents GM states, which are not estimated in SKF and treated as "consider" parameters. F s is the state transition matrix of ∆κ and linear for position and velocity vectors expressed in ECEF [52]. F p = exp (−β∆t) I n×n , ∆t = measurement update interval = 1 s. β is the inverse of time constant of the GM process. It is 100 s (smoothing interval of carrier smoothed pseudorange measurements). w k−1 denotes process noise vector at t k−1 , which is distributed as Gaussian with zero mean and covariance Q. ξ k is the measurement noise vector assumed to be Gaussian distributed with zero mean and covariance W . Q and W are given by Q s is determined based on prior information about user dynamics and receiver clock power spectral density. Q p is a diagonal matrix. Its ith diagonal term is given by (1 − (exp (−β∆t) ) 2 ) × ((σ i ura ) 2 + (σ i tr ) 2 + (σ i mp ) 2 + (σ i n ) 2 ) so that the standard deviation of the ith GM state is equal to (σ i ura ) 2 + (σ i tr ) 2 + (σ i mp ) 2 + (σ i n ) 2 at steady state. In W , a small variance of 10 −6 is used for pseudorange measurements to avoid numerical issues with the SKF implementation. When a satellite sets or a new satellite comes into visibility, the variance is increased to 0.05 for 10 s to deal with transients. W rr is the pseudorange rate error covariance matrix whose diagonal terms are appropriately chosen for a C/N 0 of 45 dB-Hz. The measurement model matrix C s, k is defined in [52]. It has additional columns for receiver clock bias and drift of NavIC, in addition to those for GPS. C p, k = I n×n [0] n×n .
Its dimension is provided in the subscript. The a priori error covariance matrix Σ − and Kalman gain K can be partitioned for estimated states and "consider" parameters as it should be noted that zero rows of the Kalman gain matrix correspond to GM states p gm . The a posteriori error covariance is given by R k is given by (C k Σ − k C T + W ). A numerically stable forward recursive modified rank-one update algorithm to find Σ + k is implemented in this paper, following [43]. Its derivation can be found in [48]. K s, k is derived from Σ − k , and also recursively calculated. Time update steps are the same as those of an EKF. The measurement update equation of SKF is given by Since GM states are initialized to zero, they always remain zero in the SKF implementation.
Second, it is justified that the calculated λ min Θ for VPE U 2, i is lower than its actual value. From discussions before Equations (37) and (55), it is evident that the actual λ min Θ is equal to λ min Ω . Next, it is briefly shown that λ min ΨΘ > λ min Ψ, maxΘ , where Θ = (cI + Ψ max ) −1 for α 2 and VPE U 2, i . The minimum eigenvalue of (Ψ max Θ) is if Ψ max and Ψ are related as Ψ max = Ψ + ∆Ψ, then substituting for Θ and Ψ max , λ min ΨΘ can be shown to be where D = (cI + ∆Ψ)Ψ −1 , and the maximum eigenvalue of D is λ max D . For a 1 of Ψ between a min and a max , ∆Ψ has elements lower than 0.08. Next, Ψ sub is simulated five thousand times for each (N − M) (number of epochs for α 2 ) with uniformly random a 1 between a min and a max , and (N − M) is varied from 1 to 150. With the simulated Ψ sub Figure A1 justifies that Hence, By definition,S TS = I and ||S || 2 = 1. Next, it was found through simulations that the minimum eigenvalue ofS T ΨS S T ΘS is greater than that ofS T Ψ maxS S T ΘS . The MATLAB code to validate this is provided in the supplementary materials (see matlab_prog_min_lambda_theta.m). Thus, the multiplication ofS preserves the inequality of the preceding equation. The calculated λ min Θ for VPE U 2, i is the minimum eigenvalue ofS T Ψ maxS S T ΘS instead of the actual one (i.e., minimum eigenvalue ofS T ΨS S T ΘS = Ω). Therefore, the calculated λ min Θ is lower than its actual value, as desired.