Dynamical systems: Luenberger observers, how to calculate the gain matrix and apply them?

Question

Summary

For the dynamical system given below, which defines a Luenberger observer, how does one explicitly obtain a value for the observer gain, $L$?

$\hat{x}(k+1) = A \hat{x}(k) + B u(k) + L[y(k) - \hat{y}(k)]$

$\hat{y}(k) = C \hat{x}(k)$

Are my assumptions and understanding in parts 1. - 5. below correct?

Context

Part of a colleague's work is a review on fault detection techniques. They wish to give the reader some insight into when it is appropriate to consider certain fault detection methods, and to give a high-level explanation as to why choices in certain papers worked well or were, in other cases, ill-informed. They asked me if I could provide some mathematical insight in some areas.

One of the most basic model-based methods for doing this work uses state observers.

From Wikipedia:

The state of a linear, time-invariant physical discrete-time system is assumed to satisfy

$x(k+1) = A x(k) + B u(k)$

$y(k) = C x(k) + D u(k)$

In order to make my question easier to answer I want to establish what the pertinent variables are and what they represent. We also make some modifications to the notation and a simplifying assumption:

• $x$ represents a state variable that, in the application domain, cannot be directly observed or measured.
• For known-input observers $u$ represents a state variable for which all of the values of $u(k)$ are known.$^1$
• $y$ represents a state variable that, in the application domain, is an 'output' measured at each $k$. This, practically speaking means that we know all the values of $y(k)$.
• We use a carat to denote estimated variables.
• For the application in question it is often assumed that $D=0$.

We make a distinction between the true value of the variables and the values that we estimate during our calculation:

$\hat{x}(k+1) = A \hat{x}(k) + B u(k)$

$\hat{y}(k) = C \hat{x}(k)$

$u$ is known without error at each time step, and in a practical sense we measure it to determine this value. We can only know $x$ without error if $x(0)$ is known without error, and the state space matrices $A$, $B$, and $C$, are known without error. In practice this never happens. Hence the above form of the equations which acknowledge that as we calculate $x$ and $y$, we will only have approximations of the true values; $\hat{x}$ and $\hat{y}$.

Our goal is to compare the estimated value $\hat{y}$ to the known value $y$. If the two values differ then we know that a measurement error has occurred (that is, the value of $y$ which we assumed was a correct representation of the system's behaviour is actually incorrect due to a measurement error), or that the system behaviour has changed unexpectedly (the values in $A$, $B$, and $C$ have changed, new terms have been added, the system has become non-linear or completely different).

The method for deciding that the values are sufficiently different for sufficiently long is not important, I don't think.

$^1$There is a class of observers, unknown-input observers, for which this is not the case, but that is not what is being asked about here.

My Understanding, Assumptions, and Questions

So, the first method I want to get a handle on is the Luenberger observer. The Luenberger observer seems to have been developed for control purposes as a way to make control systems tolerent to noise and changes in uncontrolled inputs, but the literature on fault detection definitely describes them as useful for detecting certain faults. In the Luenberger observer we have:

$\hat{x}(k+1) = A \hat{x}(k) + B u(k) + L[y(k) - \hat{y}(k)]$

$\hat{y}(k) = C \hat{x}(k)$

Again, from Wikipedia:

The observer is called asymptotically stable if the observer error $e(k) = \hat{x}(k) - x(k)$ converges to zero when $k \rightarrow > \infty$. For a Luenberger observer, the observer error satisfies $e(k+1) = (A - LC) e(k)$. The Luenberger observer for this discrete-time system is therefore asymptotically stable when the matrix $A-LC$ has all the eigenvalues inside the unit circle.

So the understanding/intuition I have here is as follows, I'm relatively confident in this but would appreciate it if anyone has corrections or extensions:

1. We have a dynamical system, $x(k+1) = A x(k) + B u(k)$. The 'visible' (measurable) outputs of this system are given by $y(k) = C x(k)$.
2. Our model of the system is not perfect, so we compute an estimated $\hat{y}(k)$, and we compare it to the measured value of $y(k)$.
3. If $\hat{y}(k)$ and $y(k)$ are sufficiently different for some definition of 'sufficiently different' we can conclude that either our measurement of $y(k)$ was erroneous, or that the model which was used to calculate $\hat{y}(k)$ became intolerably inaccurate at some point. We use this as a basis to detect different kinds of faults.
4. By adding the term $L[y(k) - \hat{y}(k)]$ to the first equation we have introduced a kind of feedback loop into the model. If either of the error types occur and $y(k) - \hat{y}(k) \neq \vec{0}$ then this will affect our calculation of $x(k+1)$, which in turn can (but might not, depending on L) affect our calculation of $\hat{y}(k+1)$.

There is a lot in the literature about properties of $L$, and inequalities or equalities that $L$ must satisfy. Wikipedia gives the example $A-LC$ has eigenvalues inside the unit circle. I have found many other relations for different situations and different forms of dynamical systems, but I have thus far failed to find any explanation of how to actually choose $L$.

5. I assume that $L$ must be chosen such that 'expected' problems (such as noise or changes in $u$) in the model are removed or compensated for, so that a family of $y(k) - \hat{y}(k)$ will produce a $y(k+1) - \hat{y}(k+1)$ within this family such that $|y(k+1) - \hat{y}(k+1)| < |y(k) - \hat{y}(k)|$.
   That is, for some family of residuals this system should remain stable. This is desireable from a control perspective and it translates relatively well to fault detection since we wish to be tolerant of noise and changes in system inputs.
   Any remaining $y(k) - \hat{y}(k)$ will cause successive values of $|y(k') - \hat{y}(k')|$ for $k'>k$ to be larger and larger. This would make 'expected' deviation from the model unimportant and magnify error since our estimated values $\hat{x}$ and $\hat{y}$ would likely diverge from normal. From a control perspective this would be unfortunate but completely expected behaviour given that the practical equivalent of this is a sensor failure or process fault. From the perspective of fault detection we can use this divergence to infer the existence of the failure.

I base this assumption on nothing other than how I would use/implement this idea, and have no idea if it is correct.

Are assumptions 1. - 5. correct?

With respect to 5. in particular, if it is incorrect, what is the correct intuition?

How does one actually choose an $L$ that has the desired behaviour outlined in 5. or a version of 5. that presents the correct intuition, if 5. is incorrect?

Answer 1

I'll address your five listed assumptions/questions, and hopefully it will give you a general sense of clarification.

1. Yes, when modeling a dynamical system using a linear time-invariant state-space model, the equation $y(k) = Cx(k)$ represents the "output" equation. In other words, the overall dynamics of the system are encapsulated in the state variables $x$, and the measurements of the system dynamics are captured in $y$. As an example, consider a simple pendulum oscillating back and forth about a pivot point. If you write the equations of motion for this system using Newton's second law, you'll find that the system dynamics are described by a second-order differential equation. Therefore, if you linearize the system (see, for example, Jacobian linearization), then the state variable $x$ used to encode the dynamics of the pendulum will be a $2\times 1$ vector. Typically, for single degree-of-freedom mechanical systems of this sort, one state variable, e.g., $x_1$, represents the "position" of the system (the angular position of the pendulum in this example). The other state variable in these types of mechanical systems, $x_2$, typically represents the "velocity" of the system (the rotational velocity of the pendulum in this example). Hence, if our pendulum system only has a "sensor" measuring the angular position of the pendulum, e.g., a rotary encoder, then the measured output would be the scalar values of position. In equations, this would be represented by $y(k) = Cx(k) = \begin{bmatrix}1 & 0 \end{bmatrix}\begin{bmatrix} x_1(k) \\ x_2(k)\end{bmatrix} = x_1(k)$ in the case that the state variable $x_1$ represents the angular position. If we have sensors measuring both position and velocity, then $C=I_2$, the $2\times 2$ identity matrix, since our output has two components: $y(k) = x(k) = (x_1(k),x_2(k))$. In the odd case that our system only has one sensor, and the sensor measures the average of the position and the velocity, then the output equation would read $y(k) = \begin{bmatrix}1/2 & 1/2 \end{bmatrix}\begin{bmatrix}x_1(k) \\ x_2(k) \end{bmatrix} = \frac{1}{2}(x_1(k)+x_2(k))$. The bottom line is this: the state variable $x$ contains all the "dynamical information" of the system, whereas the output variable $y$ is the information that you measure. If you have a lot of sensors, then the measurement matrix $C$ is relatively dense compared to the size of the state variable $x$, meaning that "most" of the dynamical information in $x$ is captured by your measurements in $y$. To be more precise, this condition is known as observability. In particular, the system is observable if and only if the observability matrix is rank-$n$, i.e., $$\text{rank}\begin{bmatrix}C \\ CA \\ \vdots \\ CA^{n-1}\end{bmatrix} = n.$$ In the case that the system is observable, you can actually exactly recover the state trajectory $\{x(k)\}_{k=0}^n$ from $n$ measurements $\{y(k)\}_{k=0}^n$. The condition that the pair $(A,C)$ is observable is very important when designing state observers (see below points).

2. Yes, one way to think of the Luenberger observer is as an estimate of the true system state in the presence of measurement and state disturbance/noise. Another way of thinking of the Luenberger observer is as a real-time implementable estimate of the system state based on measured data. This approach is in comparison to the offline approach, where, for an observable system, you would take at least $n$ measurements and solve an appropriate set of linear equations using a pseudoinverse in order to recover the exact system state during those previous times. Since the Luenberger observer converges toward exact state estimation for observable systems, it provides a useful state estimate even during the transient period before $n$ measurements are taken. Clearly, the offline version of state recovery is not as useful in real-time control applications as the online Luenberger state observer.

3. I'm not an expert in the area of fault detection, by any means. Therefore, I am not sure what type of measurements/deviations can be used to constitute a positive fault detection. In the case that $y$ is some kind of measurement of "fault" activity, and $\hat{y}$ is an estimate of this measurement based on the state estimate $\hat{x}$, then it seems reasonable to me that an unexpected change in fault activity would cause $y$ to deviate from its nominal operating point, and that the estimate $\hat{y}$ will lag in detecting such change. Again, I have no domain knowledge in this field, so I am not sure exactly what measurements are being taken, and don't really know how the resulting state estimation algorithm will respond in the case of fault activity.

4. Correct. When the estimated output, $\hat{y}$, deviates from the true measurement, $y$, the goal of the state observer is to try to adjust our state estimate $\hat{x}$ so that the estimated output better matches the measured output. So yes, the term $L(y(k)-\hat{y}(k))$ acts as a feedback term in the dynamics for the state estimate $\hat{x}$. In particular, note that when $y(k)=\hat{y}(k)$, our estimates are performing well with respect to our measurements, and therefore the observer dynamics become $\hat{x}(k+1) = A\hat{x}(k)+Bu(k)$. This shows that, at least for the time being, we "trust" our state estimate, and we expect the next best state estimate to simply be the state computed using the system's model and current estimated state.

5. Your intuition is correct. The primary goal in designing/choosing the observer gain $L$ is to ensure the stability of our estimates. A bit of rearranging of the state equations shows that $$e(k+1) = (A-LC)e(k),$$ where $e(k) = x(k)-\hat{x}(k)$ is the state estimation error. Therefore, if the closed-loop matrix $A-LC$ has eigenvalues with magnitude less that one (eigenvalues in the unit circle), then $\lim_{k\to\infty}e(k)=0$. This also implies that $y(k)-\hat{y}(k)\to 0$. Therefore, when designing $L$, we should choose it to ensure $A-LC$ has all eigenvalues in the unit circle. (Note that for continuous time systems, this requirement changes to having all eigenvalues in the left-half plane.) So now to your other question... how do we actually choose $L$? Well, remember the observability condition on the pair $(A,C)$ we discussed earlier? In the case that the system is observable, then the eigenvalues of $A-LC$ can be placed anywhere you'd like (with the restriction that complex eigenvalues come in complex conjugate pairs)! This method is called pole placement. In MATLAB, you can compute $L$ for your desired eigenvalues using the place command. You now might ask: which eigenvalues should we choose? Well, this is where control theory turns into control engineering... good pole placement requires a blend of experience, intuition, and black magic. In general, you'd like the poles to correspond to fast exponential decay with little to no amount of oscillation. For continuous time systems, these "good" poles are found far left from the origin of the complex plane, and relatively close to the real axis, since exponential modes decay faster and with less oscillation in these regions. For discrete time systems, choosing the eigenvalues to be near the origin generally results in fast responses and little oscillation. In fact, there is a pole placement strategy for discrete time systems called deat-beat control, where you compute $L$ to place all eigenvalues of $A-LC$ at the origin. In this case, the error dynamics actually converge to zero in a finite number of time steps. Aside from pole placement, you can also design state observer gains using an optimization approach. In particular, if you appeal to the duality between control and observation, you can likely compute an optimal $L$ matrix by solving an appropriate LQR problem, although there might not be a meaningful interpretation for your $R$ matrix. On the other hand, if you have a Gaussian model for the state and measurement disturbances, you can design an intuitively meaning optimal state estimator based on your model's covariance matrices. The resulting state estimator is the renowned Kalman filter.

I hope this helps!