TTK4115: Linear System Theory
$$
\newcommand{\dt}{\,\mathrm{d}t}
\newcommand{\dx}{\,\mathrm{d}x}
\newcommand{\dy}{\,\mathrm{d}y}
\newcommand{\dh}{\,\mathrm{d}h}
\newcommand{\pt}{\partial t}
\newcommand{\px}{\partial x}
\newcommand{\py}{\partial y}
\newcommand{\QEDA}{\hfill\ensuremath{\blacksquare}}
\newcommand{\QEDB}{\hfill\ensuremath{\square}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\bmat}[1]{\begin{bmatrix}#1\end{bmatrix}}
\renewcommand{\vec}[1]{\mathbf{#1}}
$$
#Realization
##Requirement for realization
- The transfer function $\hat{G}(s)$ needs to be a $\underline{\text{proper}}$ function (the degree of the numerator must be at least a high as the degree of the denominator). This is because no real device can amplify signals more and more at infinite frequencies!
- It also needs to be $\underline{\text{rational}}$.
##Definition
Given a rational proper matrix function $\hat{G}(s)$, a realization is any state-space
model $(\vec{A};\vec{B};\vec{C};\vec{D})$ such that $\hat{G}(s)$ is the corresponding transfer matrix, i.e.
$$\hat{G}(s) = \vec{C}(s\vec{I} - \vec{A})^{-1}\vec{B} + \vec{D}$$
A realization is said to be minimal if its state has the least achievable dimension; in
particular, this corresponds to the requirement that state-space model $(\vec{A};\vec{B};\vec{C};\vec{D})$
is both $\underline{\text{controllable}}$ and $\underline{\text{observable}}$.
Furthermore, this means that in a minimal realization, the dimension of A is equal to the degree of the transfer function (unobservable and uncontrollable states will disappear during Laplace transformation).
####Other:
We got
$$ G(s) = G_{sp}(s) + G_\infty(s) $$
where $G_{sp}(s)$ is the strictly proper par of the transfer function.
## Minimal Realization
A realization is minimal if the system on a canonical form is observable and controllable. It is also a minimal realization if the dimension of the nxn-matrix A is the same order as the denominator in the transfer matrix.
#Jordan Form
## Definition
The Jordan form of a system can be derived by means of a transformation matrix $ \vec{T} $ such that
$$
\vec{J} = \vec{T}^{-1} \vec{A} \vec{T}, \qquad \vec{\hat{B}} = \vec{T}^{-1} \vec{B}, \qquad \vec{\hat{C}} = \vec{C} \vec{T}
$$
where T consist of the generalized eigenvectors of A
$$
\vec{T}= [ \textbf{v}_1 \quad \textbf{v}_2 \quad ... \quad \textbf{v}_n ]
$$
#Discretization
## Exact discretization
Ongoing...
If we want to discretize a system on the form
$$ \vec{\dot{x}}(t) = \vec{A}\vec{x}(t) +\vec{B}\vec{u}(t) $$
$$ y(t) =\vec{C}\vec{x}(t) + \vec{D}\vec{u}(t) $$
To the form (given a sample time T)
$$ \vec{x}[k +1] = \vec{A}_d\vec{x}[k] +\vec{B}_d\vec{u}[k] $$
$$ y[k] =\vec{C}_d\vec{x}[k] + \vec{D}_d\vec{u}[k] $$
We need to find the matrices $\vec{A}_d, \vec{B}_d, \vec{C}_d, \vec{D}_d $:
1. $\vec{A}_d = e^{\vec{A}T} = Qe^{\hat{\vec{A}}T}Q^{-1} = \mathcal{L}^{-1}((\vec{sI} - \vec{A})^{-1})_{t=T}$ (Choose your liking, a third option is to use The Cayley Hamilton Theorem in calculating square matrix functions).
2. $\vec{B}_d = \vec{A}^{-1}(\vec{A}_d - \vec{I})\vec{B} = \int_{0}^{T} e^{\vec{A}\tau} \vec{B} d \tau$
3. $\vec{C}_d = \vec{C}$
4. $\vec{D}_d = \vec{D}$
## Euler discretization
Forward Euler Discretization:
$$ \dot{x} \approx \frac{x[k+1]-x[k]}{T_s} $$
Backward Euler Discretization:
$$ \dot{x} \approx \frac{x[k]- x[k-1]}{T_s} $$
None of these last discretizations are commonly used in this subject.
#Observability
We check if we are able to observe our states. Essential if we want to make a estimator, e.g. a Kalman filter. It is often the case that we can observe states even though we can't observe them directly, often through derivation. For example, you can find the velocity of an object by differentiating the position measurements.
A simple way of checking observability, is checking if the observability matrix $ \mathcal{O}=
\begin{bmatrix}
\mathbf{C}\\
\mathbf{C}\mathbf{A}\\
\vdots\\
\mathbf{C}\mathbf{A^{n-1}}
\end{bmatrix} $ has full rank.
Finding the observability matrix in Matlab:
Ob = obsv(A,C);
#Controllability
We check that our actuations are able to affect our states.
A simple way of checking this, is to check if the controllability matrix
$\mathcal{C}=
\begin{bmatrix}
\mathbf{B} & \mathbf{A}\mathbf{B} & \quad ... \quad \mathbf{A^{n-1}}\mathbf{B}
\end{bmatrix}$ has full rank.
Finding the controllability matrix in Matlab:
Co = ctrb(A,B);
#Linearization
If we have a non-linear system, we can use linearize the system around an operating point (often chosen to be equilibrium):
$ \mathbf{x} = \mathbf{x_p} \quad \mathbf{u} = \mathbf{u_p} $.
Linearizing amounts to transforming our system from $\mathbf{\dot{x}} = \mathbf{A} \mathbf{x} + \mathbf{B} \mathbf{u}$ to $\mathbf{\dot{\tilde{x}}} = \mathbf{\tilde{A}} \mathbf{\tilde{x}} + \mathbf{\tilde{B}} \mathbf{\tilde{u}}$ where
$
\mathbf{\tilde{A}} =
\begin{bmatrix}
\frac{\partial h_{1}}{\partial x_{1}} & \dots & \frac{\partial h_{1}}{\partial x_{n}}\\
\vdots & \ddots & \vdots \\
\frac{\partial h_{n}}{\partial x_{1}} & \dots & \frac{\partial h_{n}}{\partial x_{n}}
\end{bmatrix}\Bigg|_{\mathbf{x} = \mathbf{x_p},\mathbf{u} = \mathbf{u_p}}
$
$
\mathbf{\tilde{B}} =
\begin{bmatrix}
\frac{\partial h_{1}}{\partial u_{1}} & \frac{\partial h_{1}}{\partial u_{p}}\\
\vdots & \vdots \\
\frac{\partial h_{n}}{\partial u_{1}} & \frac{\partial h_{n}}{\partial u_{p}}
\end{bmatrix}\Bigg|_{\mathbf{x} = \mathbf{x_p},\mathbf{u} = \mathbf{u_p}}
$.
#LQR - linear quadratic regulator
## Definition
Linear quadratic regulator for finding the control input $\mathbf{u} = \mathbf{P}\mathbf{r}- \mathbf{K} \mathbf{x}$, based on optimizing the cost function
$J = \int_{0}^{\infty} \left(\mathbf{x^T}(t)\mathbf{Q}\mathbf{x}(t) + \mathbf{u^T}(t)\mathbf{R}\mathbf{u}(t) \right) dt$, under the constraints of the system. If the upper bound of the integral in the cost function is a finite time, the resulting control will be timevarying.
## When can we use a LQR?
We can use a LQR when we have a linear controllable system. See the sections [Linearization](#Linearization) and [Controllability](#Controllability).
## How do we use a LQR?
We need to choose the real, symmetric, positive definite $\mathbf{Q}$ and $\mathbf{R}$ in order to tune the regulator.
###Bryson's rule
Bryson's rule can be used as a starting point to tune the regulator.
$
Q_{ii} = \frac{1}{\text{maximum accepted value of }x_i^2}
$
$
R_{jj} = \frac{1}{\text{maximum accepted value of }u_j^2}
$
###Begin the punishment
With or without using Bryson's rule, you should adjust the diagonal elements of $\mathbf{Q}$ and $\mathbf{R}$ based on the response of your system. Increasing the diagonal elements of $\mathbf{Q}$ amounts to punishing the state errors of the corresponding states in $\mathbf{x}$ (increasing $Q_{22}$ punishes errors in state $x_{2}$). Decreasing the diagonal elements of $\mathbf{R}$ amounts to punishing the errors in requires actuation in $\mathbf{u}$. Increasing a diagonal element of $\mathbf{R}$ in many cases corresponds to decreasing a diagonal element of $\mathbf{Q}$. Scaling both matrices by a constant will not affect the cost function, as the constant can be moved outside the integral.
###Finding K
The desired feedback gain K can be found by using the MATLAB command K = lqr(A,B,Q,R).
Or you could solve the algebraic Ricatti equation for finding the symmetric positive-definite P:
$$\mathbf{A}^T\mathbf{P} + \mathbf{P}\mathbf{A} + \mathbf{Q} - \mathbf{P}\mathbf{B}\mathbf{R}^{-1}\mathbf{B}^T\mathbf{P} = 0$$
and then finding the gain
$$\mathbf{K} = \mathbf{R}^{-1}\mathbf{B}^T\mathbf{P}$$
in
$$\mathbf{u}(t) = -\mathbf{K}\mathbf{x}(t) $$.
###Finding optional feedforward gain P
Assuming that feedback results in a stable equilibrium yields the steady-state condition
$$
\mathbf{0} = \mathbf{\dot{x}} = \mathbf{A}\mathbf{x}_\infty + \mathbf{B}\mathbf{u} = \mathbf{A}\mathbf{x}_\infty + \mathbf{B}(\mathbf{P}\mathbf{r_0} - \mathbf{K}\mathbf{x}_\infty)
$$
$$
(\mathbf{A}-\mathbf{B}\mathbf{K})\mathbf{x}_\infty = -\mathbf{B}\mathbf{P}\mathbf{r_0}
$$
Inserting $\mathbf{y} = \mathbf{C}\mathbf{x}$ yields
$$
\mathbf{y}_\infty = [\mathbf{C}(\mathbf{B}\mathbf{K} - \mathbf{A})^{-1}\mathbf{B}]\mathbf{P}\mathbf{r_0}
$$
we want $\mathbf{y}_\infty = \mathbf{r_0}$, hence we choose
$$
\mathbf{P} = [\mathbf{C}(\mathbf{B}\mathbf{K} - \mathbf{A})^{-1}\mathbf{B}]^{-1}
$$
##Integral effect in LQR
Integral effect in a LQR will compensate for process disturbances and (asymptotically) set the correct bias if reference feed forward is omitted.
#State estimation
Example of simulink state estimation implementation:

Example of state estimation pole placement:

#Stability
##BIBO stability
A system is Bounded Input Bounded Output (BIBO) stable if every bounded input results in a bounded output given $\mathbf{x_0}=\mathbf{0}$ (zero-state response).
Asymptotic stability implies BIBO stability (but not necessarily the other way).
- A $\underline{\text{continuous}}$ system is BIBO stable if every pole of the transfer function has a negative real part. The poles must be in the $\underline{\text{strict left half}}$ of the s-plane
- If the system is $\underline{\text{discrete}}$, the poles have to be $\underline{\text{inside the unit circle}}$ in the z-plane for BIBO stability.
- Given a $\underline{\text{impulse response}}$ of a system, the system i BIBO stable if, and only if the impulse response $g(t)$ is absolutely integrable on the domain $[0,\infty)$. That means there has to exist a $M > 0$ such that:
$$ \int_0^{\infty} |g(t)| dt \leq M < \infty $$
Simply integrate your impulse response from zero to infniity. If the result is less than infinity, your system is BIBO stable.
##Lyapunov stability
A system is Lyapunov stable if every finite initial state results in a bounded response (but not necessarily goes to zero = marginal stability). This is therefore the stability of the zero-input response, with no forcing input ($\mathbf{u}=\mathbf{0}$).
If for any given positive definite symmetric matrix N, the Lyapunov equation has a unique
symmetric $\underline{\text{positive definite}}$ solution M (Sylvester criterion), the system is asymptotically stable.
Lyapunov equation:
$$ {\vec{A}}^{T}\vec{M} + \vec{M}\vec{A} = -\vec{N}$$
####What does it mean that a matrix is positive definite?
A symmetrix matrix $\mathbf{M}$ is said to be positive definite if the scalar $\mathbf{x}^T\mathbf{M}\mathbf{x}$ is positive for all nontrivial $\mathbf{x}$.
####Checking if a matrix is positive definite
- Check if all eigenvalues are postitive or
- Check if all leading principal minors are postive
- In a 2x2-M-matrix, this refers to $M_{11}$ and $det(M)$.
- In a 3x3-M-matrix, this refers to $M_{11}$ , $det(M_{11}, M_{12}, M_{21}, M_{22})$ and $det(M)$
- Etc.
##Marginal stability
Lyapunov stability where $\mathbf{x}$ doesn't go to 0 as time goes to $\infty$ (it's only bounded). All eigenvalues of A have zero or negative real parts. And the Jordan blocks of the eigenvalues with zero real parts must be scalar (one-by-one).
##Asymptotic stability
Lyapunov stability where $\mathbf{x}$ goes to 0 as time goes to $\infty$. All eigenvalues of A have negative real parts.
#Process Classification
## Deterministic Process
Knowledge about the signal and it's statistical properties for $t \leq t_0$ makes identification of the parametres such that we can know the rest of the process for $ t \geq t_0 $.
## Wide Sense-Stationary Process
A process is stationary if all its density functions are time-invariant. In the case of WSS, this only applies to the mean, variance and autocorrelation. The mean and the variance have to be constant and the autocorrelation function have to be dependent only on the time difference $ t_2 - t_1 $.
## Ergodicity
An ergodic process is where time averages are equal to ensample averages. This is to show that no statistical properties changes due to change in time and space. This is normally found after a long realization of the process.
# Random Signals
##Some general properties of normal random variables
- In general, two statistically independent variables are uncorrelated. The opposite is not generally true. However, two normal random variables that are uncorrelated are statistically independent.
- The probability density function describing a vector random variable $\mathbf{X}$ is completely defined by specifying the mean and the covariance matrix.
- The covariance matrix of $\mathbf{X}$ is positive definite. The magnitudes of all correlation coefficients are less than unity.
- Normality is perserved in a linear transformation.
## Autocovariance Function
Autocovariance tells us how much a signal correlates with itself across time. It includes a possibly time varying mean of X:
$A_{x}(t_1,t_2) = cov(X(t_1,t_2)) = E[(X(t_1) - m_{x}(t_1)) (X(t_2) - m_{x}(t_2))]$.
## Autocorrelation Function
Autocorrelation also tells us how well the process is correlated with itself at two different times. However, it assumes zero mean, and is more frequently used the autocovariance.
$$ R_x (t_1, t_2) = E[X(t_1) X(t_2)] $$
If the process is stationary, its probability density functions are invariant with time, and the autocorrelation function depends only on the time difference $ \tau = t_2 - t_1 $.
$$ R_x (\tau) = E[X(t) X(t+\tau)] $$
Qualitatively, if the autocorrelation decreases rapidly with $\tau$, the process changes rapidly with time. Conversly, a slowly changing process will have an autocorrelation function that decreases slowly with $\tau$. Therefore, it contains information about the frequency content, as we will see in [the section about power spectral density](#power-spectral-density).
Another way of writing the autocorrelation:
$$ R_x (\tau) =\int_{-\infty}^\infty \int_{-\infty}^\infty x(t_1) x(t_2) f_{X(t_1) X(t_2)}(x(t_1),x(t_2))dx(t_1)dx(t_2) $$
The last integral is difficult to evaluate, since it is difficult to know the joint density function $f_{X(t_1) X(t_2)}(x(t_1),x(t_2))$ explicitly.
If the [ergodicity](#ergodicity) hypothesis applies, it is easier to compute the autocorrelation as a time average, rather than an ensamble average:
$$\mathbf{\Re_{x}} = \lim_{T \to \infty} \frac{1}{T} \int_{0}^T X_{A}(t) X_{A}(t+\tau)dt $$ where $X_{A}(t)$ denotes a sample realization of the $X(t)$ process. This _time autocorrelation function_ is of course only equal to the usual autocorrelation function when the process is ergodic.
### Properties of the autocorrelation function (for stationary processes)
- $ R_x (0) = E[(X(t))^2] $ = mean-square value (for a zero-mean process, this is equal to the variance)
- $ R_x (\tau) $ even function of $\tau$ (dir of time shift is immaterial)
- $\left | R_x(\tau) \right | \leq R_x(0)$
- $\lim_{\tau \to \infty} R_x(\tau) = 0$ if $X(t)$ does not contain any periodic components (this also implies zero mean for the process, as a constant is a special case of a periodic function)
- $\mathfrak{F}[R_x(\tau)]$ is real, symmetric and nonnegative
## Crosscorrelation Function
Just as the autocorrelation function tells us something about how process is correlated with itself, the crosscorrelation function provides information about the mutual correlation between two processes $ X(t) $ and $ Y(t) $. In the stationary case:
$$ R_{XY} (\tau) = E[X(t) Y(t+\tau)] $$
It should be noted that interchanging the order of the subscripts of the crosscorrelation function has the effect of changing the sign of the argument, and that the crosscorrelation lacks the symmetry possessed by the autocorrelation function.
For zero-mean _uncorrelated_ stationary random processes $X$ and $Y$ the additative combination has this property:
$$R_{X+Y} (\tau) = R_{X} (\tau) + R_{Y} (\tau)$$
## White Noise
White noise is deined to be a stationary random process having a constant spectral density function. The term "white" is a carryover from optics, where white light is light containing all visible frequencies. Denoting the white noise spectral amplitude as A, we have
$$ S_{wn} (j \omega) = A $$
with autocorrelation function
$$ R_{wn} (\tau) = A \delta (\tau) $$
This is a process with infinite variance, which is more a useful mathematical abstraction than a physical process. All physical systems are bandlimited, therefore bandlimited white noise should be used. However, the mathematical forms are more complicated in the bandlimited case.
##Shaping filter
The filter $G(s)$ shaping unity white noise into a given spectral density funtion $S_x(s)$.
#Power Spectral Density
## Definition
The power spectrum $ S_x(j \omega) $ of a time series $ x(t) $ describes the distribution of power into frequency components composing that specific signal. By Fourier analysis, any physical signal can be decomposed into a number of discrete frequencies or a spectrum over a continous range.
In other words, the power spectral density function shows the strength of the variations (energy) as a function of frequency. It shows at which frequencies variations are strong and opposite. The unit of the Power Spectral Density is energy per frequency. You can compute the PSD is by Fast Fourier Transform (FFT). In this subject, it is done by computing the autocorrelation function and then transform it. This is called the Wiener-Khinchine relation:
$$
S_X (j \omega) = \mathfrak{F}\{R_X(\tau)\} = \int_{-\infty}^\infty R_X(\tau) e^{-j\omega\tau} d\tau
$$
It can also be shown that
$$
S_X (j \omega) = \lim_{T\to \infty} E\left [ \frac{1}{T}|\mathfrak{F}\{X_T(t)\}|^{2} \right ]
$$
where $X_T(t)$ is the truncated version of the process (truncated to zero outside span of time T), because a stationary process is not absolutely integrable wandering of to infinity, which would cause our Fourier transform to diverge.
Since a white noise is a stochastic process and with infinite energy, it has a flat PSD.
The power of a signal is defined as:
$$
P = \lim_{T\to \infty} \frac{1}{2T} \int\limits_{-T}^T x(t)^2 dt
$$
$$
S_x (j \omega) = {|G(j \omega)|}^2 \cdot S_h(j \omega) = G(j \omega) G(-j \omega) \cdot S_h(j \omega)
$$
Where $ S_h(j \omega) $ is the spectral density of the input noise.
Like for the autocorrelation, the spectral density of the sum of two uncorrelated processes x and y is given by:
$$S_{X+Y}(j\omega) = S_{X}(j\omega) + S_{Y}(j\omega)$$
##Finding the mean-square value of a stationary process, given its spectral function
$$
R_X (0) = E[X(t) X(t)] = \frac{1}{2\pi} \int_{-\infty}^\infty S_X (j \omega) d\omega
$$
#Gauss-Markov processes
The Gauss-Markov process is a "noise" process, which fits a large nummber of physical processes. It is a stationary Gaussian process with an exponential autocorrelation (which completely defines the process).
$$
R_X (\tau) = \sigma^2e^{-\beta|t|}
$$
Spectral density:
$$
S_X (j \omega) =\frac{2\sigma^2\beta}{\omega^2+\beta^2}
$$
Mean-square value:
$$
E[(X(t))^2] = \frac{1}{2\pi} \int_{-\infty}^\infty S_X (j \omega) d\omega = \sigma^2
$$
#Narrowband Gaussian process
When a narrowband system is excited by wideband Gaussian noise, the output is narrowband Gaussian noise.
#Wiener (or Brownian-motion) process
The rigorous way of describing white noise!
Nonstationary Gaussian (random-walk) continuous process, the ouput of an integrator driven by Gaussian white noise.
The mean of the ouput is 0, as the mean of white noise is 0. The variance can be calculated as the mean-square value using the sifting property. The calculations for unity white noise shows that the variance increases linearly with time:
$$E[(X(t))^2] = t$$
The autocorrelation is given by:
$$
R_X(t_1,t_2) = E[X(t_1) X(t_2)] = min\{t_1,t_2\}
$$
#Determination of autocorrelation and spectral density from experimental data
Take a DSP course.
#Kalman filtering
Kalman filtering applies a recursive method for estimation of a random process. The optimization criterion used is minimization of the mean-squared estimation error of the random variable x.
#Algorithm for discrete Kalman filter
In Kalman filtering, we assume that the key statistical parameters are determined, and the remaining problem is that of estimating the random process as it evolves with time.
Given: a priori estimate error covariance, initial a priori state estimate, process noise covariance and measurement noise covariance.
1. Compute the Kalman gain using the estimate error covariance and the measurement noise:
- $\mathbf{L}_k = \mathbf{P}_k^{-}\mathbf{C}^T(\mathbf{C}\mathbf{P}_k^{-}\mathbf{C}^T+\mathbf{{\bar{R}_v}})^{-1}$.
2. Update the estimate with the measurement:
- $\mathbf{\hat{x}}_k =\mathbf{\hat{x}}_k ^{-}+\mathbf{L}_k(\mathbf{y}_k-\mathbf{C}\mathbf{\hat{x}}_k ^{-})$.
3. Compute error covariance for updated estimate:
- $\mathbf{P}_k = (\mathbb{I}-\mathbf{L}_k\mathbf{C})\mathbf{P}_k^{-}(\mathbb{I}-\mathbf{L}_k\mathbf{C})^{T}+\mathbf{L}_k\mathbf{\bar{R}_v}\mathbf{L}_k^{T}$).
(Alternatively for optimal gain: $ \mathbf{P}_k = (\mathbb{I}-\mathbf{L}_k\mathbf{C})\mathbf{P}_k^{-}$.
4. Project ahead:
- $\mathbf{P}_{k+1}^{-} = \mathbf{A}\mathbf{P}_k^{-}\mathbf{A}^{T}+\mathbf{E}\mathbf{Q_w}\mathbf{E}^{T}$
- $\mathbf{\hat{x}}_{k+1}^{-} = \mathbf{A}\mathbf{\hat{x}}_k+\mathbf{B}\mathbf{u}_k$.
5. Repeat with $k = k +1$.
#Continuous Kalman filter
This filter can actually be derived from the discrete form. It is not as useful, though.
The process:
$$ \mathbf{\dot{x}} = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u} + \mathbf{G} \mathbf{w} $$
$$ \mathbf{y} = \mathbf{C}\mathbf{x} + \mathbf{v} $$
$$E[\mathbf{v}(t)\mathbf{v}(\tau)^T] = \delta (t-\tau)\mathbf{R}$$
$$E[\mathbf{w}(t)\mathbf{w}(\tau)^T] = \delta (t-\tau)\mathbf{Q}$$
The noise and disturbance are unbiased, white and uncorrelated.
Kalman (optimal) gain:
$$ \mathbf{L}(t) = \mathbf{P}(t)\mathbf{C}^T\mathbf{R}^{-1}$$
where
$$\mathbf{\dot{P}} = \mathbf{A}\mathbf{P} +\mathbf{P}\mathbf{A}^T + \mathbf{G}\mathbf{Q}\mathbf{G}^T-\mathbf{P}\mathbf{C}^T\mathbf{R}^{-1}\mathbf{C}\mathbf{P}$$.
This nonlinear differential equation is called the matrix Ricatti equation, and it can be solved by assuming $\mathbf{P} = \mathbf{X}\mathbf{Z}^{-1}$, $\mathbf{Z}(0) = \mathbf{I}$, and splitting into a pair of linear differential equations.
The stationary gain is found by setting $$\mathbf{\dot{P}}=\mathbf{0}$$.
#Useful tips (Timesaving etc..)
- The rank of a matrix will not change after pre- or postmultiplication by an invertable matrix.
- Sifting property for dirac delta function (for calculating mean-square values of white noise sent through filter transfer function etc.)
- Rules for using the transposed of a matrix:
- $(A^T)^T = A$
- $(A + B)^T = A^T + B^T$
- $(AB)^T = B^TA^T$
- $det(A^T) = det(A)$
- It is sometimes the case that you don't have to calculate the entire controllability/observability matrix to determine if it has full rank.
- For any block-triangular matrix
$$
\mathbf{G} =
\begin{bmatrix}
G_{11} & G_{12}\\
0 & G_{22}\\
\end{bmatrix}
$$
or
$$
\mathbf{G} =
\begin{bmatrix}
G_{11} & 0 & 0\\
G_{21} & G_{22}& 0\\
G_{31} & G_{32}& G_{33}\\
\end{bmatrix}
$$
it holds that $det(G) = det(G_{11})det(G_{22})...det(G_{nn})$. It is the union of the determinants of its diagonal elements. Leading to that the eigenvalues of $G$; $\lambda_1, \lambda_2...\lambda_n$, is equal the the eigenvalues of $G_{11}, G_{22} ... G_{nn}$
- The inverse Fourier transform of a constant is:
$$ \mathcal{F}^{-1}(a) = a \delta (t) $$
- There are risks of numerical problems when $\mathbf{R_k}$ is zero in the Kalman filter. This can be avioded by letting $\mathbf{R_k}$ be small compared to the $(\mathbf{C}\mathbf{P}_k^{-}\mathbf{C}^T)$ term.
# Contributors
Mainly written by Kristian Fjelde Pedersen, Anders Vatland, Christian Peter Bech Aschehoug and Rune Nordmo
#Questions for Q&A
- Are there any situations where having values in other places than the diagonal in Q and R (in LQR) is useful?
- Why is the covarince matrix for unity white noise in the Kalman filter equal to 1?
- Why is the LQR timevarying when the upper limit of the cost function integral is finite?