TTK4115: Linear System Theory
Realization
Requirement for realization
- The transfer function
$\hat{G}(s)$ needs to be a$\underline{\text{proper}}$ function (the degree of the denominator must be at least as high as the degree of the numerator). 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
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
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
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.
Similarity transform
The similarity transform is used to transform objects to similar object. This may be useful for different control- or analysis purposes. Similarity transform is done by having a similar state
In this course we usually use similarity transform on our well known system:
This may be transformed into:
Where
The matrix
Jordan Form
Definition
The Jordan form of a system can be derived by means of a transformation matrix
where T consist of the generalized eigenvectors of A
Diagonalization
A system can be diagonalized through the same formulas as above, just with
Discretization
Exact discretization
Ongoing...
If we want to discretize a system on the form
To the form (given a sample time T)
We need to find the matrices
$\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).$\vec{B}_d = \vec{A}^{-1}(\vec{A}_d - \vec{I})\vec{B} = \int_{0}^{T} e^{\vec{A}\tau} \vec{B} d \tau$ $\vec{C}_d = \vec{C}$ $\vec{D}_d = \vec{D}$
Euler discretization
Forward Euler Discretization:
Backward Euler Discretization:
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
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
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):
Linearizing amounts to transforming our system from
Luenberger Observer
Formula
The eigenvalues of the matrix
What does it do?
It observes things.
LQR - linear quadratic regulator
Definition
Linear quadratic regulator for finding the control input
When can we use a LQR?
We can use a LQR when we have a linear controllable system. See the sections Linearization and Controllability.
How do we use a LQR?
We need to choose the real, symmetric, positive definite
Bryson's rule
Bryson's rule can be used as a starting point to tune the regulator.
Begin the punishment
With or without using Bryson's rule, you should adjust the diagonal elements of
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:
and then finding the gain
in
Finding optional feedforward gain P
Assuming that feedback results in a stable equilibrium yields the steady-state condition
Inserting
we want
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:
The (amazing) separation principle
If a system is controllable, we can design a controller and an estimator for the system. In, fact we can assign the eigenvalues of the closed-loop system and the eigenvalues of the estimation error independently!
Stability
BIBO stability
A system is Bounded Input Bounded Output (BIBO) stable if every bounded input results in a bounded output given
-
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 (
If for any given positive definite symmetric matrix N, the Lyapunov equation has a unique
symmetric
Lyapunov equation:
What does it mean that a matrix is positive definite?
A symmetrix matrix
Checking if a matrix is positive definite
- Check if all eigenvalues are postitive or
- Check if all leading principal minors (norsk: øvre submatriser) 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.
- In a 2x2-M-matrix, this refers to
Marginal stability
Lyapunov stability where
Asymptotic stability
A subset of marginal stability.
Lyapunov stability where
Process Classification
Deterministic Process
Knowledge about the signal and it's statistical properties for
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
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:
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.
If the process is stationary, its probability density functions are invariant with time, and the autocorrelation function depends only on the time difference
Another way of writing the autocorrelation:
The last integral is difficult to evaluate, since it is difficult to know the joint density function
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
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
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
with autocorrelation function
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
Power Spectral Density
Definition
The power spectrum
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:
It can also be shown that
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:
Where
Like for the autocorrelation, the spectral density of the sum of two uncorrelated processes x and y is given by:
Finding the mean-square value of a stationary process, given its spectral function
Gauss-Markov processes
The Gauss-Markov process is a "noise" process, which fits a large number of physical processes. It is a stationary Gaussian process with an exponential autocorrelation (which completely defines the process).
Spectral density:
Mean-square value:
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:
The autocorrelation is given by:
Determination of autocorrelation and spectral density from experimental data
We don't know. Should've taken 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.
- 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}$ .
- 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 ^{-})$ .
- 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^{-}$ .
- 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$ .
- 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:
The noise and disturbance are unbiased, white and uncorrelated.
Kalman (optimal) gain:
where
This nonlinear differential equation is called the matrix Ricatti equation, and it can be solved by assuming
The stationary gain is found by setting
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?
- What are the benefits of integral effect in the estimator?
- Do several minimal realizations exist for the same transfer function? Yes, can do similarity transform, and get a A-matrix with the same dimensions.