Control of Linear Systems

State-space representation is a mathematical model used to describe a linear time-invariant system using a set of first-order difference (or differential) equations. It expresses the system’s behavior in terms of state variables, inputs, and outputs.

State-space representation of a linear system

Linear time-varying systems

State-space representation of a continuous-time linear time-invariant system

\[ \dot{x}(t) = A(t) x(t) + B(t) u(t) \] \[ y(t) = C(t) x(t) + D(t) u(t) \]

Here, \( x(t) \in \mathbb{R}^n \), \( u(t) \in \mathbb{R}^p \), \( y(t) \in \mathbb{R}^m \) are the state, input, and output vectors at time \(t\); \( A(t) \in \mathbb{R}^{n\times n}\), \(B(t) \in \mathbb{R}^{n \times p}\), \(C(t) \in \mathbb{R}^{m \times n}\), \(D(t) \in \mathbb{R}^{m\times p}\) are time-varying matrices; \(\dot{x}(t)\) denotes differentiation of \(x\) with respect to time \(t\).

State-space representation of a discrete-time time linear time-invariant system

\[ x[k+1] = A[k] x[k] + B[k] u[k] \] \[ y[k] = C[k] x[k] + D[k] u[k] \]

Here, \( x[k] \in \mathbb{R}^n \), \( u[k] \in \mathbb{R}^p \), \( y[k] \in \mathbb{R}^m \) are the state, input, and output vectors at time \(t\); \( A[k] \in \mathbb{R}^{n\times n}\), \(B[k] \in \mathbb{R}^{n \times p}\), \(C \in \mathbb{R}^{m \times n}\), \(D[k] \in \mathbb{R}^{m\times p}\) are time-varying matrices.

Linear time-invariant systems

State-space representation of a continuous-time linear time-invariant system

\[ \dot{x}(t) = A x(t) + B u(t) \] \[ y(t) = C x(t) + D u(t) \]

Here, \( x(t) \in \mathbb{R}^n \), \( u(t) \in \mathbb{R}^p \), \( y(t) \in \mathbb{R}^m \) are the state, input, and output vectors at time \(t\); \( A \in \mathbb{R}^{n\times n}\), \(B \in \mathbb{R}^{n \times p}\), \(C \in \mathbb{R}^{m \times n}\), \(D \in \mathbb{R}^{m\times p}\) are constant matrices; \(\dot{x}(t)\) denotes differentiation of \(x\) with respect to time \(t\).

State-space representation of a discrete-time time linear time-invariant system

\[ {x}[k+1] = A x[k] + B u[k] \] \[ y[k] = C x[k] + D u[k] \]

Here, \( x[k] \in \mathbb{R}^n \), \( u[k] \in \mathbb{R}^p \), \( y[k] \in \mathbb{R}^m \) are the state, input, and output vectors at time \(t\); \( A \in \mathbb{R}^{n\times n}\), \(B \in \mathbb{R}^{n \times p}\), \(C \in \mathbb{R}^{m \times n}\), \(D \in \mathbb{R}^{m\times p}\) are constant matrices.

For a detailed discussion on , refer notes by MIT on state-space systems.

Reduced-order polynomial of a matrix using Cayley-Hamilton theorem

Cayley-Hamilton theorem

If the characteristic polynomial of \( A \in \mathbb{R}^{n \times n} \) is

\[ \text{det}(\lambda I - A) = \lambda^n + a_{n-1} \lambda^{n-1} + \cdots + a_1 \lambda + a_0 \]

Then

\[ A^n + a_{n-1} A^{n-1} + \cdots + a_1 A + a_0 I = 0 \]

As a consequence of , we get the following result on reducing the order of a polynomial in matrix \(A\).

Reducing the order of a polynomial in matrix using Cayley-Hamilton theorem

For any polynomial in matrix \(A\), say \(P(A)\), we can find some \(a_i\in\mathbb{R}\), \(i=0,1,\ldots,n-1\) such that

\[ P(A)=a_{n-1} A^{n-1} + \cdots + a_1 A + a_0 I \]

The unknowns \(a_i\in\mathbb{R}\), \(i=0,1,\ldots,n-1\) can be found using coefficient matching. will be used to calculate exponential and power of a matrix in and , respectively.

For a formal discussion on , refer notes by MIT on computing the matrix exponential using the Cayley-Hamilton method

Computing exponential of a matrix

Exponential of a matrix, \(e^{At}\), appears in the solution of a state-space equation of a continuous-time, linear time-invariant system (). In this section, we will explain how to calculate exponential of a matrix using Cayley-Hamilton theorem ().

Example: Computing exponential of a matrix using Cayley-Hamilton theorem

Given matrix \(A=\begin{bmatrix} 0 & 1 \\ -2 & -3 \end{bmatrix}\), let's compute \(\exp{(A)}\). From , we know

\[ \exp{(A)} = a_1 A + a_0 I \]

The unknowns \(a_0, a_1\) can be found using the eigenvalues of \(A\). This is because eigenvectors of a matrix \(A\) are also the eigenvectors of \(\exp{(A)}\). The characteristic polynomial of matrix \(A\) is given by

\[ \det(\lambda I - A) = \begin{vmatrix} \lambda & -1 \\ 2 & \lambda + 3 \end{vmatrix} = \lambda(\lambda + 3) + 2 = \lambda^2 + 3\lambda + 2 \]

This implies that the eigenvalues of \(A\) are \(1,2\). Let \(v_1,v_2\) be the corresponding eigenvectors. Multiply , on the right, with \(v_1,v_2\) yields

\[ e^2 v_2 = 2 a_1 v_2 + a_0 v_2 \] \[ e v_1 = 2 a_1 v_1 + a_0 v_1 \]

Using coefficients matching, we obtain the following linear equations

\[ \begin{bmatrix}e \\ e^2\end{bmatrix} = \begin{bmatrix}1 & 1\\ 2 & 1\end{bmatrix} \begin{bmatrix}a_1 \\ a_0\end{bmatrix} \] \[ \begin{bmatrix}a_1 \\ a_0\end{bmatrix} = \begin{bmatrix}1 & 1\\ 2 & 1\end{bmatrix}^{-1} \begin{bmatrix}e \\ e^2\end{bmatrix} = \begin{bmatrix}1 & -1\\ -2 & 1\end{bmatrix}^{-1} \begin{bmatrix}e \\ e^2\end{bmatrix} \]

The solution of is

\[ a_1 = e-e^2 \] \[ a_0 = -2e + e^2 \]

Exponential of a matrix with repeated eigenvalues

For repeated eigenvalues, differentiate with respect to the matrix \(A\) and multiply again with the eigenvector to obtain another equation. Differentiate as many times as the multiplicity of the eigenvalue to obtain equations for coefficients matching.

Python code: Exponential of a matrix

import numpy as np
from scipy.linalg import expm  # Import matrix exponential function

# Define a 2x2 matrix A
A = np.array([[1, 1],
              [2, 1]])

# Compute the matrix exponential e^A
expA = expm(A)

# Display the result
print("Matrix exponential of A:\n", expA)

Computing power of a matrix

In the solution of a state-space equation of a discrete-time, linear time-invariant system (), we would be required to compute \(A^q, q\in\mathbb{I}\). We can use Cayley-Hamilton theorem () to calculate a reduced order expression of \(A^q\), similar to .

Example: Computing polynomial of a matrix using Cayley-Hamilton theorem

Given matrix \(A=\begin{bmatrix} 0 & 1 \\ -2 & -3 \end{bmatrix}\), let's compute \(\exp{(A)}\). From , we know

\[ A^q = a_1 A + a_0 I \]

The unknowns \(a_0, a_1\) can be found using the eigenvalues of \(A\). This is because eigenvectors of a matrix \(A\) are also the eigenvectors of \(A^q\). The characteristic polynomial of matrix \(A\) is given by

\[ \det(\lambda I - A) = \begin{vmatrix} \lambda & -1 \\ 2 & \lambda + 3 \end{vmatrix} = \lambda(\lambda + 3) + 2 = \lambda^2 + 3\lambda + 2 \]

This implies that the eigenvalues of \(A\) are \(1,2\). Let \(v_1,v_2\) be the corresponding eigenvectors. Multiply , on the right, with \(v_1,v_2\) yields

\[ 2^q v_2 = 2 a_1 v_2 + a_0 v_2 \] \[ v_1 = 2 a_1 v_1 + a_0 v_1 \]

Using coefficients matching, we obtain the following linear equations

The solution of is

\[ a_1 = 1-2^q \] \[ a_0 = -2 + 2^q \]

Power of a matrix with repeated eigenvalues

Python code: Power of a matrix

import numpy as np

# Define the matrix A
A = np.array([[1, 1],
              [2, 1]])

# Define the power n
n = 3

# Compute A raised to the power n
A_power_n = np.linalg.matrix_power(A, n)

# Print the result
print(f"A^{n}:\n", A_power_n)

Common control paradigms

We will discuss following control techniques.

Proportional integral derivative (PID) control (): This is suitable for single input, single output systems. PID is computationally inexpensive to implement.
Linear quadratic regulator (LQR) (): This is suitable for multiple input, multiple output systems. LQR can also be used to implement soft system constraints. LQR is computationally efficient to implement due to closed-form expression for optimal control.
Model predictive control (MPC) (): This is suitable for multiple input, multiple output systems with hard system constraints. MPC is computationally expensive as it requires solving a quadratic program at each time step.
H-∞ control (): This is a robust controller that ensures that performance of the system in the presence of an exogenous input is bounded. H-∞ control requires a solution of linear matrix inequality.

PID (proportional, integeral and derivative) control

PID control

Given a reference input \(r[k]\), PID control for a discrete-time linear time-invariant system () is given by

\[ u[k] = K_P \, e[k] + K_I \sum_{i=0}^{k} e[i] + K_D \left( e[k] - e[k-1] \right) \]

Here \( e[k] = r[k] - y[k] \) is the error at step \( k \); \( K_P, K_I, K_D \) are the proportional, integral, and derivative gains, respectively.

Proportional, integral and derivation gains

Proportional gain \(K_P\) reacts to present error. Large \( K_P \) increases response speed but may cause overshoot.
Integral gain \(K_I\) reacts to accumulated past error. It improves steady-state accuracy but may introduce oscillation or instability.
Derivative gain \(K_D\) reacts to the rate of error. It adds damping, reduces overshoot, but is sensitive to noise.

Tuning a PID controller

Proper tuning is essential to avoid instability or slow response. Tuning can be achieved either by trial-and-error with step response observation, or through optimization-based or model-based tuning involves pole placement. Stability of the closed-loop system depends on poles of the resulting transfer function.

Linear quadratic regulator (LQR)

Consider a discrete-time linear time-invariant system () with initial condition \(x[0]\).

Linear quadratic regulator (LQR) with finite horizon

A linear quadratic regulator with finite horizon of length \(N\) minimizes the following quadratic cost function.

\[ J = \sum_{k=0}^{N} \left( x[k]^T Q x[k] + u[k]^T R u[k] \right) \]

Here, \(Q \geq 0\) is the state penalty matrix, \(R > 0\) is the control penalty matrix.

Solution of the linear quadratic regulator for finite horizon

A linear quadratic regulator with finite horizon of length \(N\) finds a state-feedback control law of the form \(u[k] = -K x[k]\) that minimizes the quaratic cost \(J\) in . The optimal solution is given in .

For the linear quadratic regulator (LQR) with a finite horizon of length \( N \), the optimal control law at each time step \( k \) is

\[ u[k] = -K[k] \, x[k] \]

The feedback gain matrix \( K[k] \) is given by

\[ K[k] = \left(R + B^T P[k+1] B\right)^{-1} B^T P[k+1] A \]

The sequence \( P[k] \) is computed recursively backward in time with terminal condition \( P[N] = Q \), using:

\[ P[k] = Q + A^T P[k+1] A - A^T P[k+1] B \left(R + B^T P[k+1] B\right)^{-1} B^T P[k+1] A \]

Linear quadratic regulator (LQR) with infinite horizon

A linear quadratic regulator with infinite horizon minimizes the following quadratic cost function.

\[ J = \sum_{k=0}^{\infty} \left( x[k]^T Q x[k] + u[k]^T R u[k] \right) \]

Here, \(Q \geq 0\) is the state penalty matrix, \(R > 0\) is the control penalty matrix.

Solution of the linear quadratic regulator for infinite horizon

A linear quadratic regulator with infinite horizon finds a state-feedback control law of the form \(u[k] = -K x[k]\) that minimizes the quaratic cost \(J\) in . The optimal solution is given in .

For the linear quadratic regulator (LQR) with infinite horizon,the optimal control law at each time step \( k \) is

\[ u[k] = -K x[k] \]

The optimal gain matrix \(K\) is given by

\[ K = (R + B^T P B)^{-1} B^T P A \]

\(P\) satisfies the Ricatti equation

\[ P = A^T P A - A^T P B (R + B^T P B)^{-1} B^T P A + Q \]

The closed-loop system, given by \(x[k+1] = (A - B K) x[k]\), is asymptotically stable, i.e., \(x[k]\) converge to zero as \( k \to \infty \). To understand when the solution of Ricatti equation exists and is unique, we define the concept of stabilizability, and detectability.

Stabilizablity

A matrix pair \( (A, B) \) is stabilizable if all eigenvalues of \( A \) with \( |\lambda| \geq 1 \) can be controlled via \( B \), i.e.,

\[ \text{rank} \begin{bmatrix} \lambda I - A & B \end{bmatrix} = n \quad \text{for all } |\lambda| \geq 1 \]

Detectability

A matrix pair \( (Q, A) \) is detectable if all eigenvalues of \( A \) with \( |\lambda| \geq 1 \) are observable through \( Q \), i.e.,

\[ \text{rank} \begin{bmatrix} \lambda I - A \\ Q^{1/2} \end{bmatrix} = n \quad \text{for all } |\lambda| \geq 1 \]

If \((A, B)\) is stabilizable and \((Q, A)\) is detectable then solution \(P\) of the Ricatti equation in exists and is unique.

Model predictive control (MPC)

MPC solves an optimization problem that predicts future system behavior over a finite horizon using a mathematical model. Based on the prediction, it computes an optimal sequence of control inputs to minimize a cost function while satisfying constraints. Only the first input in the sequence is applied, and the process repeats at the next time step with updated state measurements. This is called the receding horizon principle.

Implementation of model predictive control

At each time step \( k \), MPC solves an optimization problem to minimize a cost function over a finite prediction horizon \( N \):

\[ \min_{\{u[k], \ldots, u[k+N-1]\}} \sum_{i=k}^{k+N-1} \left( x[i]^T Q x[i] + u[i]^T R u[i] \right) + x[N]^T Q_f x[N] \]

such that \[ x_{\min} \leq x[k+i] \leq x_{\max}, \quad u_{\min} \leq u[k+i] \leq u_{\max}, \quad \text{for } i = 0, \dots, N-1 \]

Here \( Q \geq 0 \) is a state weighting matrix; \( R > 0 \) is an input weighting matrix; \( Q_f \geq 0 \) is a terminal state weight; \(x[k]\) is the current state of the system. MPC can incorporate additional system constraints; this would require solving a quadratic program at each time step. After optimization, MPC applies only the first control input \( u[0] \), and then repeats the optimization at the next time step with new state of the sytem as the initial condition.

The advantage of MPC is its ability to incorporate system constraints. Moreover, under model uncertainity, online optimization at each time step provides a feedback to the controller.

H-∞ control

Consider the following discrete-time linear system:

\[ x[k+1] = A x[k] + B_1 w[k] + B_2 u[k] \] \[ z[k] = C_1 x[k] + D_{12} u[k] \]

Here, \( x[k] \in \mathbb{R}^n \) is the system state; \( u[k] \in \mathbb{R}^m \) is the control input; \( w[k] \in \mathbb{R}^p \) is an exogenous disturbance input; \(z[k]\) is the performance output.

An H-∞ control law \(u[k]=K x[k]\) ensures that the closed-loop transfer function from \( w[k] \) to \( z[k] \) satisfies

\[ \| T_{zw}(z) \|_\infty = \min_{u}\max_{\lVert w\rVert\neq 0} \frac{\lVert z\rVert}{\lVert w\rVert} < \gamma \]

Here, \(\gamma>0\) is the desired performance bound.

H-∞ controller minimizes the influence of exogenous input \(w[k]\) on the performance output \(z[k]\). This ensures that the system is robust to disturbances, keeping the effect of \( w[k] \) on \( z[k] \) below a specified bound \( \gamma > 0 \). The min-max form expresses the robust performance objective: the controller tries to minimize the worst-case effect of disturbances.

provides a condition for existence of an H-∞ control law.

Existence of an H-∞ control law

Step 2: An H-∞ control law exists if there exists a matrix \( P \succ 0 \) such that

\[ (A + B_2 K)^T P (A + B_2 K) - P + (C_1 + D_{12} K)^T (C_1 + D_{12} K) + \gamma^{-2} P B_1 B_1^T P \preceq 0 \]

An equivalent condition for existence of an H-&infin is if there exists a solution- pair \((X\succ 0, F)\) such that the following Linear matrix inequality (LMI) hold

\[ X - A X A^T - A F^T B_2^T - B_2 F A^T - B_2 F F^T B_2^T - \gamma^{-2} B_1 B_1^T - (C_1 X + D_{12} F)^T (C_1 X + D_{12} F) \succeq 0 \]

If a solution-pair \((X\succ 0, F)\) exists, the H-&infin controller is obtained as

\[ K = F X^{-1} \]

Author

Anurag Gupta is an M.S. graduate in Electrical and Computer Engineering from Cornell University. He also holds an M.Tech degree in Systems and Control Engineering and a B.Tech degree in Electrical Engineering from the Indian Institute of Technology, Bombay.

Comment

Past Comments

No comments yet. Be the first!