Kalman Filtering with Equality and Inequality State Constraints
aa r X i v : . [ m a t h . O C ] S e p Report no. 07/18
Kalman Filtering with Equality and Inequality StateConstraints
Nachi Gupta Raphael Hauser
Oxford University Computing Laboratory, Numerical Analysis Group,Wolfson Building, Parks Road, Oxford OX1 3QD, U.K.
Both constrained and unconstrained optimization problems regularly appearin recursive tracking problems engineers currently address – however, constraintsare rarely exploited for these applications. We define the Kalman Filter anddiscuss two different approaches to incorporating constraints. Each of these ap-proaches are first applied to equality constraints and then extended to inequalityconstraints. We discuss methods for dealing with nonlinear constraints and forconstraining the state prediction. Finally, some experiments are provided to indi-cate the usefulness of such methods.
Key words and phrases:
Constrained optimization, Kalman filtering, Nonlinear filters,Optimization methods, Quadratic programming, State estimation
The first author would like to thank the Clarendon Bursary for financial support.Oxford University Computing LaboratoryNumerical Analysis GroupWolfson BuildingParks RoadOxford, England OX1 3QD
E-mail: [email protected]
February, 2008
Introduction
Kalman Filtering [8] is a method to make real-time predictions for systems with some knowndynamics. Traditionally, problems requiring Kalman Filtering have been complex and non-linear. Many advances have been made in the direction of dealing with nonlinearities (e.g.,Extended Kalman Filter [1], Unscented Kalman Filter [7]). These problems also tend to haveinherent state space equality constraints (e.g., a fixed speed for a robotic arm) and state space inequality constraints (e.g., maximum attainable speed of a motor). In the past, less interesthas been generated towards constrained Kalman Filtering, partly because constraints can bedifficult to model. As a result, constraints are often neglected in standard Kalman Filteringapplications.The extension to Kalman Filtering with known equality constraints on the state space isdiscussed in [5, 11, 13, 14, 16]. In this paper, we discuss two distinct methods to incorporateconstraints into a Kalman Filter. Initially, we discuss these in the framework of equality con-straints. The first method, projecting the updated state estimate onto the constrained region,appears with some discussion in [5, 11]. We propose another method, which is to restrict theoptimal Kalman Gain so the updated state estimate will not violate the constraint. With somealgebraic manipulation, the second method is shown to be a special case of the first method.We extend both of these concepts to Kalman Filtering with inequality constraints in thestate space. This generalization for the first approach was discussed in [12]. Constraining theoptimal Kalman Gain was briefly discussed in [10]. Further, we will also make the extensionto incorporating state space constraints in Kalman Filter predictions.Analogous to the way a Kalman Filter can be extended to solve problems containing non-linearities in the dynamics using an Extended Kalman Filter by linearizing locally (or by usingan Unscented Kalman Filter), linear inequality constrained filtering can similarly be extendedto problems with nonlinear constraints by linearizing locally (or by way of another schemelike an Unscented Kalman Filter). The accuracy achieved by methods dealing with nonlinearconstraints will naturally depend on the structure and curvature of the nonlinear function it-self. In the two experiments we provide, we look at incorporating inequality constraints to atracking problem with nonlinear dynamics.
A discrete-time Kalman Filter [8] attempts to find the best running estimate for a recursivesystem governed by the following model : x k = F k , k − x k − + u k , k − , u k , k − ∼ N (cid:0) , Q k , k − (cid:1) (2.1) z k = H k x k + v k , v k ∼ N ( , R k ) (2.2) The similar extension for the method of [16] was made in [6]. The subscript k on a variable stands for the k -th time step, the mathematical notation N ( m , S ) denotes anormally distributed random vector with mean m and covariance S , and all vectors in this paper are columnvectors (unless we are explicitly taking the transpose of the vector). x k is an n -vector that represents the true state of the underlying system and F k , k − is an n × n matrix that describes the transition dynamics of the system from x k − to x k . Themeasurement made by the observer is an m -vector z k , and H k is an m × n matrix that trans-forms a vector from the state space into the appropriate vector in the measurement space. Thenoise terms u k , k − (an n -vector) and v k (an m -vector) encompass known and unknown errorsin F k , k − and H k and are normally distributed with mean 0 and covariances given by n × n matrix Q k , k − and m × m matrix R k , respectively. At each iteration, the Kalman Filter makesa state prediction for x k , denoted ˆ x k | k − . We use the notation k | k − k − k .The state prediction error ˜ x k | k − is defined as the difference between the true state and thestate prediction, as below. ˜ x k | k − = x k − ˆ x k | k − (2.3)The covariance structure for the expected error on the state prediction is defined as theexpectation of the outer product of the state prediction error. We call this covariance structurethe error covariance prediction and denote it P k | k − . P k | k − = E h(cid:0) ˜ x k | k − (cid:1) (cid:0) ˜ x k | k − (cid:1) ′ i (2.4)The filter will also provide an updated state estimate for x k , given all the measurementsprovided up to and including time step k . We denote these estimates by ˆ x k | k . We similarlydefine the state estimate error ˜ x k | k as below.˜ x k | k = x k − ˆ x k | k (2.5)The expectation of the outer product of the state estimate error represents the covariancestructure of the expected errors on the state estimate, which we call the updated error covari-ance and denote P k | k . P k | k = E h(cid:0) ˜ x k | k (cid:1) (cid:0) ˜ x k | k (cid:1) ′ i (2.6)At time-step k , we can make a prediction for the underlying state of the system by allowingthe state to transition forward using our model for the dynamics and noting that E (cid:2) u k , k − (cid:3) = x k | k − = F k , k − ˆ x k − | k − (2.7)If we expand the expectation in Equation (2.4), we have the following equation for theerror covariance prediction. P k | k − = F k , k − P k − | k − F ′ k , k − + Q k , k − (2.8)We can transform our state prediction into the measurement space, which is a predictionfor the measurement we now expect to observe. We use the prime notation on a vector or a matrix to denote its transpose throughout this paper. z k | k − = H k ˆ x k | k − (2.9)The difference between the observed measurement and our predicted measurement is themeasurement residual, which we are hoping to minimize in this algorithm. n k = z k − ˆ z k | k − (2.10)We can also calculate the associated covariance for the measurement residual, which is theexpectation of the outer product of the measurement residual with itself, E (cid:2) n k n ′ k (cid:3) . We call thisthe measurement residual covariance. S k = H k P k | k − H ′ k + R k (2.11)We can now define our updated state estimate as our prediction plus some perturbation,which is given by a weighting factor times the measurement residual. The weighting factor,called the Kalman Gain, will be discussed below.ˆ x k | k = ˆ x k | k − + K k n k (2.12)Naturally, we can also calculate the updated error covariance by expanding the outer prod-uct in Equation (2.6). P k | k = ( I − K k H k ) P k | k − ( I − K k H k ) ′ + K k R k K ′ k (2.13)Now we would like to find the Kalman Gain K k , which minimizes the mean square stateestimate error, E h(cid:12)(cid:12) ˜ x k | k (cid:12)(cid:12) i . This is the same as minimizing the trace of the updated errorcovariance matrix above. After some calculus, we find the optimal gain that achieves this,written below. K k = P k | k − H ′ k S − k (2.14)The covariance matrices in the Kalman Filter provide us with a measure for uncertaintyin our predictions and updated state estimate. This is a very important feature for the variousapplications of filtering since we then know how much to trust our predictions and estimates.Also, since the method is recursive, we need to provide an initial covariance that is largeenough to contain the initial state to ensure comprehensible performance. For a more detaileddiscussion of Kalman Filtering, we refer the reader to the following book [1]. The I in Equation (2.13) represents the n × n identity matrix. Throughout this paper, we use I to denote thesame matrix, except in Appendix A, where I is the appropriately sized identity matrix. Note that v ′ v = trace [ vv ′ ] for some vector v . We could also minimize the mean square state estimate error in the N norm, where N is a positive definiteand symmetric weighting matrix. In the N norm, the optimal gain would be K Nk = N K k . Equality Constrained Kalman Filtering
A number of approaches have been proposed for solving the equality constrained KalmanFiltering problem [5, 11, 13, 14, 16]. In this paper, we show two different methods. The firstmethod will restrict the state at each iteration to lie in the equality constrained space. Thesecond method will start with a constrained prediction, and restrict the Kalman Gain so thatthe estimate will lie in the constrained space. Our equality constraints in this paper will bedefined as below, where A is a q × n matrix, b a q -vector, and x k , the state, is a n -vector. Ax k = b (3.1)So we would like our updated state estimate to satisfy the constraint at each iteration, asbelow. A ˆ x k | k = b (3.2)Similarly, we may also like the state prediction to be constrained, which would allow abetter forecast for the system. A ˆ x k | k − = b (3.3)In the following subsections, we will discuss methods for constraining the updated stateestimate. In Section 4, we will extend these concepts and formulations to the inequality con-strained case, and in Section 6, we will address the problem of constraining the prediction, aswell. We can solve the following minimization problem for a given time-step k , where ˆ x Pk | k is theconstrained estimate, W k is any positive definite symmetric weighting matrix, and ˆ x k | k is theunconstrained Kalman Filter updated estimate.ˆ x Pk | k = arg min x ∈ R n n(cid:0) x − ˆ x k | k (cid:1) ′ W k (cid:0) x − ˆ x k | k (cid:1) : Ax = b o (3.4)The best constrained estimate is then given byˆ x Pk | k = ˆ x k | k − W − k A ′ (cid:0) AW − k A ′ (cid:1) − (cid:0) A ˆ x k | k − b (cid:1) (3.5)To find the updated error covariance matrix of the equality constrained filter, we first definethe matrix ¡ below. ¡ = W − k A ′ (cid:0) AW − k A ′ (cid:1) − (3.6)Equation (3.5) can then be re-written as following. A and b can be different for different k . We don’t subscript each A and b to avoid confusion. Note that ¡ A is a projection matrix, as is ( I − ¡ A ) , by definition. If A is poorly conditioned, we can use a QRfactorization to avoid squaring the condition number. x Pk | k = ˆ x k | k − ¡ (cid:0) A ˆ x k | k − b (cid:1) (3.7)We can find a reduced form for x k − ˆ x Pk | k as below. x k − ˆ x Pk | k = x k − ˆ x k | k + ¡ (cid:0) A ˆ x k | k − b − ( Ax k − b ) (cid:1) (3.8a) = x k − ˆ x k | k + ¡ (cid:0) A ˆ x k | k − Ax k (cid:1) (3.8b) = − ( I − ¡ A ) (cid:0) ˆ x k | k − x k (cid:1) (3.8c)Using the definition of the error covariance matrix, we arrive at the following expression. P Pk | k = E (cid:20)(cid:16) x k − ˆ x Pk | k (cid:17) (cid:16) x k − ˆ x Pk | k (cid:17) ′ (cid:21) (3.9a) = E h ( I − ¡ A ) (cid:0) ˆ x k | k − x k (cid:1) (cid:0) ˆ x k | k − x k (cid:1) ′ ( I − ¡ A ) ′ i (3.9b) = ( I − ¡ A ) P k | k ( I − ¡ A ) ′ (3.9c) = P k | k − ¡ AP k | k − P k | k A ′ ¡ ′ + ¡ AP k | k A ′ ¡ ′ (3.9d) = P k | k − ¡ AP k | k (3.9e) = ( I − ¡ A ) P k | k (3.9f)It can be shown that choosing W k = P − k | k results in the smallest updated error covariance.This also provides a measure of the information in the state at k . Alternatively, we can expand the updated state estimate term in Equation (3.2) using Equation(2.12). A (cid:0) ˆ x k | k − + K k n k (cid:1) = b (3.10)Then, we can choose a Kalman Gain K Rk , that forces the updated state estimate to be inthe constrained space. In the unconstrained case, we chose the optimal Kalman Gain K k , bysolving the minimization problem below which yields Equation (2.14). K k = arg min K ∈ R n × m trace (cid:2) ( I − KH k ) P k | k − ( I − KH k ) ′ + KR k K ′ (cid:3) (3.11) If M and N are covariance matrices, we say N is smaller than M if M − N is positive semidefinite. Anotherformulation for incorporating equality constraints into a Kalman Filter is by observing the constraints as pseudo-measurements [14, 16]. When W k is chosen to be P − k | k , both of these methods are mathematically equivalent [5].Also, a more numerically stable form of Equation (3.9) with discussion is provided in [5]. K Rk that satisfies the constrained optimization problem writtenbelow for a given time-step k . K Rk = arg min K ∈ R n × m trace (cid:2) ( I − KH k ) P k | k − ( I − KH k ) ′ + KR k K ′ (cid:3) s.t. A (cid:0) ˆ x k | k − + K n k (cid:1) = b (3.12)We will solve this problem using the method of Lagrange Multipliers. First, we take thesteps below, using the vec notation (column stacking matrices so they appear as long vectors,see Appendix A) to convert all appearances of K in Equation (4.8) into long vectors. Let usbegin by expanding the following term. trace (cid:2) ( I − KH k ) P k | k − ( I − KH k ) ′ + KR k K ′ (cid:3) = trace (cid:2) P k | k − − KH k P k | k − − P k | k − H ′ k K ′ + KH k P k | k − H ′ k K ′ + KR k K ′ (cid:3) (2.11) = trace (cid:2) P k | k − − KH k P k | k − − P k | k − H ′ k K ′ + KS k K ′ (cid:3) = trace (cid:2) P k | k − (cid:3) − trace (cid:2) KH k P k | k − (cid:3) − trace (cid:2) P k | k − H ′ k K ′ (cid:3) + trace (cid:2) KS k K ′ (cid:3) (3.13a)We now expand the last three terms in Equation (3.13a) one at a time. trace (cid:2) KH k P k | k − (cid:3) (A.9) = vec h(cid:0) H k P k | k − (cid:1) ′ i ′ vec [ K ]= vec (cid:2) P k | k − H ′ k (cid:3) ′ vec [ K ] (3.14)trace (cid:2) P k | k − H ′ k K ′ (cid:3) (A.9) = vec [ K ] ′ vec (cid:2) P k | k − H ′ k (cid:3) (3.15)trace (cid:2) KS k K ′ (cid:3) (A.9) = vec [ K ] ′ vec [ KS k ] (A.7) = vec [ K ] ′ ( S ⊗ I ) vec [ K ] (3.16)Remembering that trace (cid:2) P k | k − (cid:3) is constant, our objective function can be written as be-low. vec [ K ] ′ ( I ⊗ S k ) vec (cid:2) K ′ (cid:3) − vec (cid:2) P k | k − H ′ k (cid:3) ′ vec [ K ] − vec [ K ] ′ vec (cid:2) P k | k − H ′ k (cid:3) (3.17)Using Equation (A.8) on the equality constraints, our minimization problem is the follow-ing. Throughout this paper, a number in parentheses above an equals sign means we made use of this equationnumber. We use the symmetry of P k | k − in Equation (3.14) and the symmetry of S k in Equation (3.16). Rk = arg min K ∈ R n × m vec [ K ] ′ ( S k ⊗ I ) vec [ K ] − vec (cid:2) P k | k − H ′ k (cid:3) ′ vec [ K ] − vec [ K ] ′ vec (cid:2) P k | k − H ′ k (cid:3) s.t. (cid:0) n ′ k ⊗ A (cid:1) vec [ K ] = b − A ˆ x k | k − (3.18)Further, we simplify this problem so the minimization problem has only one quadraticterm. We complete the square as follows. We want to find the unknown variable m which willcancel the linear term. Let the quadratic term appear as follows. Note that the non-“vec [ K ] "term is dropped as is is irrelevant for the minimization problem. ( vec [ K ] + m ) ′ ( S k ⊗ I ) ( vec [ K ] + m ) (3.19)The linear term in the expansion above is the following.vec [ K ] ′ ( S k ⊗ I ) m + m ′ ( S k ⊗ I ) vec [ K ] (3.20)So we require that the two equations below hold. ( S k ⊗ I ) m = − vec (cid:2) P k | k − H ′ k (cid:3) m ′ ( S k ⊗ I ) = − vec (cid:2) P k | k − H ′ k (cid:3) ′ (3.21)This leads to the following value for m . m (A.3) = − (cid:0) S − k ⊗ I (cid:1) vec (cid:2) P k | k − H ′ k (cid:3) (A.8) = − vec (cid:2) P k | k − H ′ k S − k (cid:3) (2.14) = − vec [ K k ] (3.22)Using Equation (A.6), our quadratic term in the minimization problem becomes the fol-lowing. ( vec [ K − K k ]) ′ ( S k ⊗ I ) ( vec [ K − K k ]) (3.23)Let l = vec [ K − K k ] . Then our minimization problem becomes the following. K Rk = arg min l ∈ R mn l ′ ( S k ⊗ I ) l s.t. (cid:0) n ′ k ⊗ A (cid:1) ( l + vec [ K k ]) = b − A ˆ x k | k − (3.24)We can then re-write the constraint taking the vec [ K k ] term to the other side as below. (cid:0) n ′ k ⊗ A (cid:1) l = b − A ˆ x k | k − − (cid:0) n ′ k ⊗ A (cid:1) vec [ K k ] (A.8) = b − A ˆ x k | k − − vec [ AK k n k ]= b − A ˆ x k | k − − AK k n k (2.12) = b − A ˆ x k | k (3.25)9his results in the following simplified form. K Rk = arg min l ∈ R mn l ′ ( S k ⊗ I ) l s.t. (cid:0) n ′ k ⊗ A (cid:1) l = b − A ˆ x k | k (3.26)We form the Lagrangian L , where we introduce q Lagrange Multipliers in vector l = (cid:0) l , l , . . . , l q (cid:1) ′ L = l ′ ( S k ⊗ I ) l − l ′ (cid:2)(cid:0) n ′ k ⊗ A (cid:1) l − b + A ˆ x k | k (cid:3) (3.27)We take the partial derivative with respect to l . ¶ L ¶ l = l ′ ( S k ⊗ I ) − l ′ (cid:0) n ′ k ⊗ A (cid:1) (3.28)Similarly we can take the partial derivative with respect to the vector l . ¶ L ¶l = (cid:0) n ′ k ⊗ A (cid:1) l − b + A ˆ x k | k (3.29)When both of these derivatives are set equal to the appropriate size zero vector, we havethe solution to the system. Taking the transpose of Equation (3.28), we can write this systemas Mn = p with the following block definitions for M , n , and p . M = (cid:20) S k ⊗ I n k ⊗ A ′ n ′ k ⊗ A [ q × q ] (cid:21) (3.30) n = (cid:20) l l (cid:21) (3.31) p = (cid:20) [ mn × ] b − A ˆ x k | k (cid:21) (3.32)We solve this system for vector n in Appendix C. The solution for l is pasted below. (cid:16)h S − k n k (cid:0) n ′ k S − k n k (cid:1) − i ⊗ h A ′ (cid:0) AA ′ (cid:1) − i(cid:17) (cid:0) b − A ˆ x k | k (cid:1) (3.33)Bearing in mind that b − A ˆ x k | k = vec (cid:2) b − A ˆ x k | k (cid:3) , we can use Equation (A.8) to re-write l as below. vec h A ′ (cid:0) AA ′ (cid:1) − (cid:0) b − A ˆ x k | k (cid:1) (cid:0) n ′ k S − k n k (cid:1) − n ′ k S − k i (3.34)The resulting matrix inside the vec operation is then an n by m matrix. Remembering thedefinition for l , we notice that K − K k results in an n by m matrix also. Since both of thecomponents inside the vec operation result in matrices of the same size, we can safely remove We used the symmetry of ( S k ⊗ I ) here. Here we used the symmetry of S − k and (cid:0) n ′ k S − k n k (cid:1) − (the latter of which is actually just a scalar). K Rk . K k − A ′ (cid:0) AA ′ (cid:1) − (cid:0) A ˆ x k | k − b (cid:1) (cid:0) n ′ k S − k n k (cid:1) − n ′ k S − k (3.35)If we now substitute this Kalman Gain into Equation (2.12) to find the constrained updatedstate estimate, we end up with the following.ˆ x Rk | k = ˆ x k | k − A ′ (cid:0) AA ′ (cid:1) − (cid:0) A ˆ x k | k − b (cid:1) (3.36)This is of course equivalent to the result of Equation (3.5) with the weighting matrix W k chosen as the identity matrix. The error covariance for this estimate is given by Equation(3.9). In the more general case of this problem, we may encounter equality and inequality constraints,as given below. Ax k = bCx k ≤ d (4.1)So we would like our updated state estimate to satisfy the constraint at each iteration, asbelow. A ˆ x k | k = bC ˆ x k | k ≤ d (4.2)Similarly, we may also like the state prediction to be constrained, which would allow abetter forecast for the system. A ˆ x k | k − = bC ˆ x k | k − ≤ d (4.3)We will present two analogous methods to those presented for the equality constrainedcase. In the first method, we will run the unconstrained filter, and at each iteration constrainthe updated state estimate to lie in the constrained space. In the second method, we willfind a Kalman Gain ˇ K Rk such that the the updated state estimate will be forced to lie in theconstrained space. In both methods, we will no longer be able to find an analytic solution asbefore. Instead, we use numerical methods. We can use the unconstrained or constrained Kalman Gain to find this error covariance matrix. Since theconstrained Kalman Gain is suboptimal for the unconstrained problem, before projecting onto the constrainedspace, the constrained covariance will be different from the unconstrained covariance. However, the differencelies exactly in the space orthogonal to which the covariance is projected onto by Equation (3.9). The proof isomitted for brevity. C and d can be different for different k . We don’t subscript each C and d to simplify notation. .1 By Projecting the Unconstrained Estimate Given the best unconstrained estimate, we could solve the following minimization problemfor a given time-step k , where ˇ x Pk | k is the inequality constrained estimate and W k is any positivedefinite symmetric weighting matrix.ˇ x Pk | k = arg min x (cid:0) x − ˆ x k | k (cid:1) ′ W k (cid:0) x − ˆ x k | k (cid:1) s.t. Ax = bCx ≤ d (4.4)For solving this inequality constrained optimization problem, we can use a variety of stan-dard methods, or even an out-of-the-box solver, like fmincon in Matlab. Here we use an ac-tive set method [4]. This is a common method for dealing with inequality constraints, wherewe treat a subset of the constraints (called the active set) as additional equality constraints.We ignore any inactive constraints when solving our optimization problem. After solving theproblem, we check if our solution lies in the space given by the inequality constraints. If itdoesn’t, we start from the solution in our previous iteration and move in the direction of thenew solution until we hit a set of constraints. For each iteration, the active set is made up ofthose inequality constraints with non-zero Lagrange Multipliers.We first find the best estimate (using Equation (3.5) for the equality constrained problemwith the equality constraints given in Equation (4.1) plus the active set of inequality constraints.Let us call the solution to this ˇ x P ∗ k | k , j since we have not yet checked if the solution lies in theinequality constrained space. In order to check this, we find the vector that we moved alongto reach ˇ x P ∗ k | k , j . This is given by the following. s = ˇ x P ∗ k | k , j − ˇ x Pk | k , j − (4.5)We now iterate through each of our inequality constraints, to check if they are satisfied. Ifthey are all satisfied, we choose t max =
1. If they are not, we choose the largest value of t max such that ˆ x k | k , j − + t max s lies in the inequality constrained space. We choose our estimate tobe ˇ x Pk | k , j = ˇ x Pk | k , j − + t max s (4.6)If we find the solution has converged within a pre-specified error, or we have reached apre-specified maximum number of iterations, we choose this as the updated state estimate toour inequality constrained problem, denoted ˇ x Pk | k . If we would like to take a further iterationon j , we check the Lagrange Multipliers at this new solution to determine the new active set. We then repeat by finding the best estimate for the equality constrained problem includingthe new active set as additional equality constraints. Since this is a Quadratic Programmingproblem, each step of j guarantees the same estimate or a better estimate. For the inequality constrained filter, we allow multiple iterations within each step. The j subscript indexesthese further iterations. The previous active set is not relevant. (cid:16) ˇ x Pk | k , j − ˇ x Pk | k , j − (cid:17) (cid:16) ˇ x Pk | k , j − ˇ x Pk | k , j − (cid:17) ′ (4.7)This is a measure of our convergence error and should typically be small relative to theunconstrained error covariance. We can then use Equation (3.9) to project the covariancematrix onto the constrained subspace, but we only use the defined equality constraints. We donot incorporate any constraints in the active set when computing Equation (3.9) since thesestill represent inequality constraints on the state. Ideally we would project the error covariancematrix into the inequality constrained subspace, but this projection is not trivial. We could solve this problem by restricting the optimal Kalman gain also, as we did for equalityconstraints previously. We seek the optimal K k that satisfies the constrained optimizationproblem written below for a given time-step k .ˇ K Rk = arg min K ∈ R n × m trace (cid:2) ( I − KH k ) P k | k − ( I − KH k ) ′ + KR k K ′ (cid:3) s.t. A (cid:0) ˆ x k | k − + K k n k (cid:1) = bC (cid:0) ˆ x k | k − + K k n k (cid:1) ≤ d (4.8)Again, we can solve this problem using any inequality constrained optimization method(e.g., fmincon in Matlab or the active set method used previously). Here we solved theoptimization problem using SDPT3, a Matlab package for solving semidefinite programmingproblems [15]. When calculating the covariance matrix for the inequality constrained estimate,we use the restricted Kalman Gain. Again, we can add on the safety term for the convergenceerror, by taking the outer product of the difference between the updated state estimates cal-culated by the restricted Kalman Gain for the last two iterations of SDPT3. This covariancematrix is then projected onto the subspace as in Equation (3.9) using the equality constraintsonly. Thus far, in the Kalman Filter we have dealt with linear models and constraints. A numberof methods have been proposed to handle nonlinear models (e.g., Extended Kalman Filter [1],Unscented Kalman Filter [7]). In this paper, we will focus on the most widely used of these, theExtended Kalman Filter. Let’s re-write the discrete unconstrained Kalman Filtering problemfrom Equations (2.1) and (2.2) below, incorporating nonlinear models. x k = f k , k − ( x k − ) + u k , k − , u k , k − ∼ N (cid:0) , Q k , k − (cid:1) (5.1) z k = h k ( x k ) + v k , v k ∼ N ( , R k ) (5.2)13n the above equations, we see that the transition matrix F k , k − has been replaced by thenonlinear vector-valued function f k , k − ( · ) , and similarly, the matrix H k , which transforms avector from the state space into the measurement space, has been replaced by the nonlinearvector-valued function h k ( · ) . The method proposed by the Extended Kalman Filter is to lin-earize the nonlinearities about the current state prediction (or estimate). That is, we choose F k , k − as the Jacobian of f k , k − evaluated at ˆ x k − | k − , and H k as the Jacobian of h k evaluated atˆ x k | k − and proceed as in the linear Kalman Filter of Section 2. Numerical accuracy of thesemethods tends to depend heavily on the nonlinear functions. If we have linear constraintsbut a nonlinear f k , k − ( · ) and h k ( · ) , we can adapt the Extended Kalman Filter to fit into theframework of the methods described thus far. Since equality and inequality constraints we model are often times nonlinear, it is important tomake the extension to nonlinear equality and inequality constrained Kalman Filtering for themethods discussed thus far. Without loss of generality, our discussion here will pertain only tononlinear inequality constraints. We can follow the same steps for equality constraints. Wereplace the linear inequality constraint on the state space by the following nonlinear inequal-ity constraint c ( x k ) = d , where c ( · ) is a vector-valued function. We can then linearize ourconstraint, c ( x k ) = d , about the current state prediction ˆ x k | k − , which gives us the following. c (cid:0) ˆ x k | k − (cid:1) + C (cid:0) x k − ˆ x k | k − (cid:1) / d (5.3)Here C is defined as the Jacobian of c evaluated at ˆ x k | k − . This indicates then, that thenonlinear constraint we would like to model can be approximated by the following linearconstraint Cx k / d + C ˆ x k | k − − c (cid:0) ˆ x k | k − (cid:1) (5.4)This constraint can be written as ˜ Cx k ≤ ˜ d , which is an approximation to the nonlinearinequality constraint. It is now in a form that can be used by the methods described thus far.The nonlinearities in both the constraints and the models, f k , k − ( · ) and h k ( · ) , could havebeen linearized using a number of different methods (e.g., a derivative-free method, a higherorder Taylor approximation). Also an iterative method could be used as in the Iterated Ex-tended Kalman Filter [1]. We can also do a midpoint approximation to find F k , k − by evaluating the Jacobian at (cid:0) ˆ x k − | k − + ˆ x k | k − (cid:1) / We replace the ‘ ≤ ’ sign with an ‘ = ’ sign and the ‘ / ’ with an ‘ ≈ ’ sign. This method is how the Extended Kalman Filter linearizes nonlinear functions for f k , k − ( · ) and h k ( · ) . Hereˆ x k | k − can be the state prediction of any of the constrained filters presented thus far and does not necessarily relateto the unconstrained state prediction. Constraining the State Prediction
We haven’t yet discussed whether the state prediction (Equation (2.7)) also should be con-strained. Forcing the constraints should provide a better prediction (which is used for fore-casting in the Kalman Filter). Ideally, the transition matrix F k , k − will take an updated stateestimate satisfying the constraints at time k − k . Of course this may not be the case. In fact, the constraints may depend on theupdated state estimate, which would be the case for nonlinear constraints. On the downside,constraining the state prediction increases computational cost per iteration.We propose three methods for dealing with the problem of constraining the state predic-tion. The first method is to project the matrix F k , k − onto the constrained space. This is onlypossible for the equality constraints, as there is no trivial way to project F k , k − to an inequal-ity constrained space. We can use the same projector as in Equation (3.9f) so we have thefollowing. F Pk , k − = ( I − ¡ A ) F k , k − (6.1)Under the assumption that we have constrained our updated state estimate, this new transi-tion matrix will make a prediction that will keep the estimate in the equality constrained space.Alternatively, if we weaken this assumption, i.e., we are not constraining the updated stateestimate, we could solve the minimization problem below (analogous to Equation (3.4)). Wecan also incorporate inequality constraints now.ˇ x Pk | k − = arg min x (cid:0) x − ˆ x k | k − (cid:1) ′ W k (cid:0) x − ˆ x k | k − (cid:1) s.t. Ax = bCx ≤ d (6.2)We can constrain the covariance matrix here also, in a similar fashion to the method de-scribed in Section 4.1. The third method is to add to the constrained problem the additionalconstraints below, which ensure that the chosen estimate will produce a prediction at the nextiteration that is also constrained. A k + F k + , k x k = b k + C k + F k + , k x k ≤ d k + (6.3)If A k + , b k + , C k + or d k + depend on the estimate (e.g., if we are linearizing nonlinearfunctions a ( · ) or b ( · ) , we can use an iterative method, which would resolve A k + and b k + using the current best updated state estimate (or prediction), re-calculate the best estimateusing A k + and b k + , and so forth until we are satisfied with the convergence. This methodwould be preferred since it looks ahead one time-step to choose a better estimate for the currentiteration. However, it can be far more expensive computationally. In these three methods, the symmetric weighting matrix W k can be different. The resulting ¡ can conse-quently also be different. Further, we can add constraints for some arbitrary n time-steps ahead. Experiments
We provide two related experiments here. We have a car driving along a straight road withthickness 2 meters. The driver of the car traces a noisy sine curve (with the noise lying only inthe frequency domain). The car is tagged with a device that transmits the location within someknown error. We would like to track the position of the car. In the first experiment, we filterover the noisy data with the knowledge that the underlying function is a noisy sine curve. Theinequality constrained methods will constrain the estimates to only take values in the interval [ − , ] . In the second experiment, we do not use the knowledge that the underlying curve is asine curve. Instead we attempt to recover the true data using an autoregressive model of order6 [3]. We do, however, assume our unknown function only takes values in the interval [ − , ] ,and we can again enforce these constraints when using the inequality constrained filter.The driver’s path is generated using the nonlinear stochastic process given by Equation(5.1). We start with the following initial point. x = (cid:20) (cid:21) (7.1)Our vector-valued transition function will depend on a discretization parameter T and canbe expressed as below. Here, we choose T to be p /
10, and we run the experiment from aninitial time of 0 to a final time of 10 p . f k , k − = (cid:20) ( x k − ) + T sin (( x k − ) + T ) (cid:21) (7.2)And for the process noise we choose the following. Q k , k − = (cid:20) .
00 0 m (cid:21) (7.3)The driver’s path is drawn out by the second element of the vector x k – the first elementacts as an underlying state to generate the second element, which also allows a natural methodto add noise in the frequency domain of the sine curve while keeping the process recursivelygenerated. To create the measurements, we use the model from Equation (2.2), where H k is the squareidentity matrix of dimension 2. We choose R k as below to noise the data. This considerablymasks the true underlying data as can be seen in Fig. 1. R k = (cid:20)
10 m
00 10 m (cid:21) (7.4) The figure only shows the noisy sine curve, which is the second element of the measurement vector. Thefirst element, which is a noisy straight line, isn’t plotted. ( m e t e r s ) Noisy MeasurementsTime (seconds)
Figure 1:
We take our sine curve, which is already noisy in the frequency domain (due to process noise),and add measurement noise. The underlying sine curve is significantly masked.
For the initial point of our filters, we choose the following point, which is different fromthe true initial point given in Equation (7.1).ˆ x | = (cid:20) (cid:21) (7.5)Our initial covariance is given as below. . P | = (cid:20) . . (cid:21) (7.6) Nonzero off-diagonal elements in the initial covariance matrix often help the filter converge more quickly
17n the filtering, we use the information that the underlying function is a sine curve, and ourtransition function f k , k − changes to reflect a recursion in the second element of x k – now wewill add on discretized pieces of a sine curve to our previous estimate. The function is givenexplicitly below. f k , k − = (cid:20) ( x k − ) + T ( x k − ) + sin (( x k − ) + T ) − sin (( x k − ) ) (cid:21) (7.7)For the Extended Kalman Filter formulation, we will also require the Jacobian of thismatrix denoted F k , k − , which is given below. F k , k − = (cid:20) (( x k − ) + T ) − cos (( x k − ) ) (cid:21) (7.8)The process noise Q k , k − , given below, is chosen similar to the noise used in generatingthe simulation, but is slightly larger to encompass both the noise in our above model and toprevent divergence due to numerical roundoff errors. The measurement noise R k is chosen thesame as in Equation (7.4). Q k , k − = (cid:20) .
00 0 . (cid:21) (7.9)The inequality constraints we enforce can be expressed using the notation throughout thechapter, with C and d as given below. C = (cid:20) − (cid:21) (7.10) d = (cid:20) (cid:21) (7.11)These constraints force the second element of the estimate x k | k (the sine portion) to lie inthe interval [ − , ] . We do not have any equality constraints in this experiment. We run theunconstrained Kalman Filter and both of the constrained methods discussed previously. A plotof the true position and estimates is given in Fig. 2. Notice that both constrained methods forcethe estimate to lie within the constrained space, while the unconstrained method can violatethe constraints. In the previous experiment, we used the knowledge that the underlying function was a noisysine curve. If this is not known, we face a significantly harder estimation problem. Let usassume nothing about the underlying function except that it must take values in the interval [ − , ] . A good model for estimating such an unknown function could be an autoregressivemodel. We can compare the unconstrained filter to the two constrained methods again usingthese assumption and an autoregressive model of order 6, or AR(6) as it is more commonlyreferred to. 18 ( m e t e r s ) True Position and EstimatesTime (seconds)
True PositionUnconstrainedActive Set MethodConstrained Kalman Gain
Figure 2:
We show our true underlying state, which is a sine curve noised in the frequency domain,along with the estimates from the unconstrained Kalman Filter, and both of our inequality constrainedmodifications. We also plotted dotted horizontal lines at the values -1 and 1. Both inequality con-strained methods do not allow the estimate to leave the constrained space.
In the previous example, we used a large measurement noise R k to emphasize the gainachieved by using the constraint information. Such a large R k is probably not very realistic,and when using an autoregressive model, it will be hard to track such a noisy signal. Togenerate the measurements, we again use Equation (2.2), this time with H k and R k as givenbelow. H k = (cid:2) (cid:3) (7.12) R k = (cid:2) . (cid:3) (7.13)19ur state will now be defined using the following 13-vector, in which the first element isthe current estimate, the next five elements are lags, the six elements afterwards are coefficientson the current estimate and the lags, and the last element is a constant term.ˆ x k | k = (cid:2) y k y k − · · · y k − a a · · · a (cid:3) ′ (7.14)Our matrix H k in the filter is a row vector with the first element 1, and all the rest as 0, so y k | k − is actually our prediction ˆ z k | k − in the filter, describing where we believe the expectedvalue of the next point in the time-series to lie. For the initial state, we choose a vector of allzeros, except the first and seventh element, which we choose as 1. This choice for the initialconditions leads to the first prediction on the time series being 1, which is incorrect as the trueunderlying state has expectation 0. For the initial covariance, we choose I [ × ] and add 0 . The transition function f k , k − for the AR(6) model is givenbelow. min ( , max ( − , a y k − + · · · + a y k − + a )) min ( , max ( − , y k − )) min ( , max ( − , y k − )) min ( , max ( − , y k − )) min ( , max ( − , y k − )) min ( , max ( − , y k − )) a a ... a a (7.15)Putting this into recursive notation, we have the following. min ( , max ( − , ( x k − ) ( x k − ) + · · · + ( x k − ) )) min ( , max ( − , ( x k − ) )) min ( , max ( − , ( x k − ) )) min ( , max ( − , ( x k − ) )) min ( , max ( − , ( x k − ) )) min ( , max ( − , ( x k − ) ))( x k − ) ( x k − ) ... ( x k − ) ( x k − ) (7.16)The Jacobian of f k , k − is given below. We ignore the min ( · ) and max ( · ) operators sincethe derivative is not continuous across them, and we can reach the bounds by numerical error. The bracket subscript notation is used through the remainder of this paper to indicate the size of zero matricesand identity matrices. ( x k − ) · · · ( x k − ) I [ × ] [ × ] ( x k − ) · · · ( x k − ) [ × ] [ × ] I [ × ] (7.17)For the process noise, we choose Q k , k − to be a diagonal matrix with the first entry as 0.1and all remaining entries as 10 − since we know the prediction phase of the autoregressivemodel very well. The inequality constraints we enforce can be expressed using the notationthroughout the chapter, with C as given below and d as a 12-vector of ones. C = " I [ × ] − I [ × ] [ × ] (7.18)These constraints force the current estimate and all of the lags to take values in the range [ − , ] . As an added feature of this filter, we are also estimating the lags at each iteration usingmore information although we don’t use it – this is a fixed interval smoothing. In Fig. 3, weplot the noisy measurements, true underlying state, and the filter estimates. Notice again thatthe constrained methods keep the estimates in the constrained space. Visually, we can see theimprovement particularly near the edges of the constrained space. We’ve provided two different formulations for including constraints into a Kalman Filter. Inthe equality constrained framework, these formulations have analytic formulas, one of whichis a special case of the other. In the inequality constrained case, we’ve shown two numericalmethods for constraining the estimate. We also discussed how to constrain the state predictionand how to handle nonlinearities. Our two examples show that these methods ensure theestimate lies in the constrained space, which provides a better estimate structure.
Appendix A Kron and Vec
In this appendix, we provide some definitions used earlier in the chapter. Given matrix A ∈ R m × n and B ∈ R p × q , we can define the right Kronecker product as below. ( A ⊗ B ) = a , B · · · a , n B ... . . . ... a m , B · · · a m , n B (A.1) The indices m , n , p , and q and all matrix definitions are independent of any used earlier. Also, the subscriptnotation a , n denotes the element in the first row and n -th column of A , and so forth. ( m e t e r s ) Measurement, True Position, and EstimatesTime (seconds)
Noisy MeasurementTrue PositionUnconstrainedActive Set MethodConstrained Kalman Gain
Figure 3:
We show our true underlying state, which is a sine curve noised in the frequency domain, thenoised measurements, and the estimates from the unconstrained and both inequality constrained filters.We also plotted dotted horizontal lines at the values -1 and 1. Both inequality constrained methods donot allow the estimate to leave the constrained space.
Given appropriately sized matrices A , B , C , and D such that all operations below are well-defined, we have the following equalities. ( A ⊗ B ) ′ = (cid:0) A ′ ⊗ B ′ (cid:1) (A.2) ( A ⊗ B ) − = (cid:0) A − ⊗ B − (cid:1) (A.3) ( A ⊗ B ) ( C ⊗ D ) = ( AC ⊗ BD ) (A.4)22e can also define the vectorization of an [ m × n ] matrix A , which is a linear transformationon a matrix that stacks the columns iteratively to form a long vector of size [ mn × ] , as below.vec [ A ] = a , ... a m , a , ... a m , ... a , n ... a m , n (A.5)Using the vec operator, we can state the trivial definition below.vec [ A + B ] = vec [ A ] + vec [ B ] (A.6)Combining the vec operator with the Kronecker product, we have the following.vec [ AB ] = (cid:0) B ′ ⊗ I (cid:1) vec [ A ] (A.7)vec [ ABC ] = (cid:0) C ′ ⊗ A (cid:1) vec [ B ] (A.8)We can express the trace of a product of matrices as below.trace [ AB ] = vec (cid:2) B ′ (cid:3) ′ vec [ A ] (A.9)trace [ ABC ] = vec [ B ] ′ ( I ⊗ C ) vec [ A ] (A.10a) = vec [ A ] ′ ( I ⊗ B ) vec [ C ] (A.10b) = vec [ A ] ′ ( C ⊗ I ) vec [ B ] (A.10c)For more information, please see [9]. Appendix B Analytic Block Representation for the inverseof a Saddle Point Matrix M S is a saddle point matrix if it has the block form below. The subscript S notation is used to differentiate these matrices from any matrices defined earlier. S = (cid:20) A S B ′ S B S − C S (cid:21) (B.1)In the case that A S is nonsingular and the Schur complement J S = − (cid:0) C S + B S A − S B ′ S (cid:1) isalso nonsingular in the above equation, it is known that the inverse of this saddle point matrixcan be expressed analytically by the following equation (see e.g., [2]). M − S = (cid:20) A − S + A − S B ′ S J − S B S A − S − A − S B ′ S J − S − J − S B S A − S J − S (cid:21) (B.2) Appendix C Solution to the system Mn = p Here we solve the system Mn = p from Equations (3.30), (3.31), and (3.32), re-stated below,for vector n . (cid:20) S k ⊗ I n k ⊗ A ′ n ′ k ⊗ A [ q × q ] (cid:21) (cid:20) l l (cid:21) = (cid:20) [ mn × ] b − A ˆ x k | k (cid:21) (C.1) M is a saddle point matrix with the following equations to fit the block structure of Equa-tion (B.1). A S = S k ⊗ I (C.2) B S = n ′ k ⊗ A (C.3) C S = [ q × q ] (C.4)We can calculate the term A − S B ′ S . A − S B ′ S = [ ( S k ⊗ I )] − (cid:0) n ′ k ⊗ A (cid:1) ′ (C.5a) (A.2)(A.3) = (cid:0) S − k ⊗ I (cid:1) (cid:0) n k ⊗ A ′ (cid:1) (C.5b) (A.4) = (cid:0) S − k n k (cid:1) ⊗ A ′ (C.5c)And as a result we have the following for J S . J S = − (cid:0) n ′ k ⊗ A (cid:1) (cid:2)(cid:0) S − k n k (cid:1) ⊗ A ′ (cid:3) (C.6a) (A.4) = − (cid:0) n ′ k S − k n k (cid:1) ⊗ (cid:0) AA ′ (cid:1) (C.6b) We use Equation (A.2) with B ′ S to arrive at the same term for B s in Equation (C.1). − S is then, as below. J − S = − (cid:2)(cid:0) n ′ k S − k n k (cid:1) ⊗ (cid:0) AA ′ (cid:1)(cid:3) − (C.7a) (A.3) = − (cid:0) n ′ k S − k n k (cid:1) − ⊗ (cid:0) AA ′ (cid:1) − (C.7b)For the upper right block of M − , we then have the following expression. A − S B ′ S J − S = (cid:2)(cid:0) S − k n k (cid:1) ⊗ A ′ (cid:3) h(cid:0) n ′ k S − k n k (cid:1) − ⊗ (cid:0) AA ′ (cid:1) − i (C.8a) (A.4) = h S − k n k (cid:0) n ′ k S − k n k (cid:1) − i ⊗ h A ′ (cid:0) AA ′ (cid:1) − i (C.8b)Since the first block element of p is a vector of zeros, we can solve for n to arrive at thefollowing solution for l . (cid:16)h S − k n k (cid:0) n ′ k S − k n k (cid:1) − i ⊗ h A ′ (cid:0) AA ′ (cid:1) − i(cid:17) (cid:0) b − A ˆ x k | k (cid:1) (C.9)The vector of Lagrange Multipliers l is given below. − h(cid:0) n ′ k S − k n k (cid:1) − ⊗ (cid:0) AA ′ (cid:1) − i (cid:0) b − A ˆ x k | k (cid:1) (C.10) References [1] Y. B AR -S HALOM , X. R. L I , AND
T. K
IRUBARAJAN , Estimation with Applications toTracking and Navigation , John Wiley and Sons, Inc., 2001.[2] M. B
ENZI , G. H. G
OLUB , AND
J. L
IESEN , Numerical solution of saddle point problems ,Acta Numerica, 14 (2005), pp. 1–137.[3] G. E. P. B
OX AND
G. M. J
ENKINS , Time Series Analysis. Forecasting and Control(Revised Edition) , Oakland: Holden-Day, (1976).[4] R. F
LETCHER , Practical methods of optimization. Vol. 2: Constrained Optimiza-tion , John Wiley & Sons Ltd., Chichester, 1981. Constrained optimization, A Wiley-Interscience Publication.[5] N. G
UPTA , Kalman filtering in the presence of state space equality constraints , inIEEE Proceedings of the 26th Chinese Control Conference, July 2007, arXiv:physics.ao-ph/0705.4563, Oxford na-tr:07/14.[6] N. G
UPTA , R. H
AUSER , AND
N. F. J
OHNSON , Forecasting financial time-series usingan artificial market model , in Proceedings of the 10th Annual Workshop on EconomicHeterogeneous Interacting Agents, June 2005, Oxford na-tr:05/09.257] S. J. J
ULIER AND
J. K. U
HLMANN , A new extension of the kalman filter to non-linear systems , in Proceedings of AeroSense: The 11th International Symposium onAerospace/Defence Sensing, Simulation and Controls, vol. 3, 1997, pp. 182–193.[8] R. E. K
ALMAN , A new approach to linear filtering and prediction problems , Transac-tions of the ASME–Journal of Basic Engineering, 82 (1960), pp. 35–45.[9] P. L
ANCASTER AND
M. T
ISMENETSKY , The Theory of Matrices: With Applications ,Academic Press Canada, 1985.[10] A. G. Q
URESHI , Constrained kalman filtering for image restoration , in Proceedings ofthe International Conference on Acoustics, Speech, and Signal Processing, vol. 3, 1989,pp. 1405 – 1408.[11] D. S
IMON AND
T. L. C
HIA , Kalman filtering with state equality constraints , IEEETransactions on Aerospace and Electronic Systems, 38 (2002), pp. 128–136.[12] D. S
IMON AND
D. L. S
IMON , Aircraft turbofan engine health estimation using con-strained kalman filtering , Journal of Engineering for Gas Turbines and Power, 127(2005), p. 323.[13] T. L. S
ONG , J. Y. A HN , AND
C. P
ARK , Suboptimal filter design with pseudomeasure-ments for target tracking , IEEE Transactions on Aerospace and Electronic Systems, 24(1988), pp. 28–39.[14] M. T
AHK AND
J. L. S
PEYER , Target tracking problems subject to kinematic constraints ,Proceedings of the 27th IEEE Conference on Decision and Control, (1988), pp. 1058–1059.[15] K. C. T OH , M. J. T ODD , AND
R. T
UTUNCU , SDPT3 — a Matlab software packagefor semidefinite programming , Optimization Methods and Software, 11 (1999), pp. 545–581.[16] L. S. W
ANG , Y. T. C
HIANG , AND
F. R. C