4 Time response of state space models
4.1 Solution to vector differential equations
Recall that for a scalar differential equation (i.e., \(x(t) \in \reals\)): \[ \dot x(t) = a x(t), \quad \hbox{initial condition } x(0)=x_0 \] has the solution \[ x(t) = e^{at} x_0. \]
Can we say the same for vector systems? That is, can we define matrix exponential \(e^{At}\) in such a way that the vector differential equation (i.e., \(x(t) \in \reals^n\)): \[ \dot x(t) = A x(t), \quad \hbox{initial condition } x(0) = x_0 \] where \(A \in \reals^{n × n}\) has the solution \[ x(t) = e^{At} x_0? \]
Mathematically, it turns out that this is straight forward. Recall that for a scalar \(a\), we have \[e^{at} = 1 + at + \frac{(at)^2}{2!} + \frac{(at)^3}{3!} + \cdots\] So, a natural choice is to define matrix exponential as \[\begin{align*} e^{At} &= I + At + \frac{(At)^2}{2!} + \frac{(At)^3}{3!} + \cdots\\ &= I + At + \frac{A^2t^2}{2!} + \frac{A^3t^3}{3!} + \cdots \end{align*}\] where the second equation uses the fact that \((At)^2 = At At = A^2 t^2\) because \(t\) is a scalar.
With this definition, we have \[\begin{align*} \frac{d}{dt} e^{At} &= A + A^2 t + \frac{A^3 t^2}{2!} + \cdots \\ &= A[I + At + \frac{A^2t^2}{2!} + \frac{A^3t^3}{3!} + \cdots] \\ &= A e^{At}. \end{align*}\]
Thus, if we take the candidate solution \(x(t) = e^{At} x_0\), we have that \[ \frac{d}{dt} x(t) = \left[ \frac{d}{dt} e^{At} \right] x_0 = A e^{At} x_0 = A x(t). \] Thus, our candidate solution satisfies the vector differential equation!
If we define \(e^{At}\) as \[\begin{equation}\label{eq:matrix-exponential} e^{At} = I + At + \frac{A^2t^2}{2!} + \frac{A^3t^3}{3!} + \cdots \end{equation}\]
Then, we can write the solution of the vector differential equation (i.e., \(x(t) \in \reals^n\)) \[ \dot x(t) = A x(t), \quad \hbox{initial condition } x(0) = x_0 \] where \(A \in \reals^{n × n}\) as \[ x(t) = e^{At} x_0. \]
We now present some examples where the matrix exponential can be computed easily.
Example 4.1 Compute \(e^{At}\) for \[A = \MATRIX{0 & 1 \\ 0 & 0}. \]
Based on the solution, solve the vector differential equation \[ \dot x(t) = A x(t), \quad x_0 = \MATRIX{α_1 \\ α_2}. \]
We compute the terms of \(\eqref{eq:matrix-exponential}\):
- \(A^2 = \MATRIX{0 & 0 \\ 0 & 0}\).
which means that \(A^k = 0\) for \(k \ge 2\). Thus, the right hand side of \(\eqref{eq:matrix-exponential}\) contains only a finite number of nonzero terms: \[\begin{align*} e^{At} &= I + At \\ &= \MATRIX{1 & 0 \\ 0 & 1} + \MATRIX{0 & t \\ 0 & 0 } \\ &=\MATRIX{1 & t \\ 0 & 1}. \end{align*}\]
Thus, the solution of the vector differential equation \(\dot x(t) = A x(t)\) is given by \[ x(t) = e^{At} x_0 = \MATRIX{1 & t \\ 0 & 1} \MATRIX{α_1 \\ α_2} = \MATRIX{α_1 + α_2 t \\ α_2}. \]
Example 4.2 Compute \(e^{At}\) for \[A = \MATRIX{0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0}. \]
Based on the solution, solve the vector differential equation \[ \dot x(t) = A x(t), \quad x_0 = \MATRIX{α_1 \\ α_2 \\ α_3}. \]
We compute the terms of \(\eqref{eq:matrix-exponential}\):
- \(A^2 = \MATRIX{0 & 0 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 }\).
- \(A^3 = \MATRIX{0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 }\).
which means that \(A^k = 0\) for \(k \ge 3\). Thus, the right hand side of \(\eqref{eq:matrix-exponential}\) contains only a finite number of nonzero terms: \[\begin{align*} e^{At} &= I + At + \frac 12 A^2 t^2 \\ &= \MATRIX{1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 } + \MATRIX{0 & t & 0 \\ 0 & 0 & t \\ 0 & 0 & 0} + \frac 12 \MATRIX{0 & 0 & t^2 \\ 0 & 0 & 0 \\ 0 & 0 & 0 } \\ &=\MATRIX{1 & t & \frac 12 t^2 \\ 0 & 1 & t \\ 0 & 0 & 1}. \end{align*}\]
Thus, the solution of the vector differential equation \(\dot x(t) = A x(t)\) is given by \[ x(t) = e^{At} x_0 = \MATRIX{1 & t & \frac 12 t^2 \\ 0 & 1 & t \\ 0 & 0 & 1} \MATRIX{α_1 \\ α_2 \\ α_3} = \MATRIX{α_1 + α_2 t + \frac 12 α_3^2 \\ α_2 + α_3 t \\ α_3}. \]
Example 4.3 Compute \(e^{At}\) for \[A = \MATRIX{λ_1 & 0 & 0 \\ 0 & λ_2 & 0 \\ 0 & 0 & λ_3}. \]
Based on the solution, solve the vector differential equation \[ \dot x(t) = A x(t), \quad x_0 = \MATRIX{α_1 \\ α_2 \\ α_3}. \]
We compute the terms of \(\eqref{eq:matrix-exponential}\):
- \(A^2 = \MATRIX{λ_1^2 & 0 & 0 \\ 0 & λ_2^2 & 0 \\ 0 & 0 & λ_3^2}.\)
- \(A^3 = \MATRIX{λ_1^3 & 0 & 0 \\ 0 & λ_2^3 & 0 \\ 0 & 0 & λ_3^3}.\)
- and so on.
Thus, we have \[\begin{align*} e^{At} &= I + At + \frac{A^2t^2}{2!} + \frac{A^3t^3}{3!} + \cdots \\ &= \MATRIX{1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1} + \MATRIX{λ_1 t & 0 & 0 \\ 0 & λ_2 t & 0 \\ 0 & 0 & λ_3 t} + \frac 12 \MATRIX{λ_1^2 t^2 & 0 & 0 \\ 0 & λ_2^2 t^2 & 0 \\ 0 & 0 & λ_3^2 t^2} + \cdots \\[10pt] &= \MATRIX{1 + λ_1 t + \frac 12 λ_1^2 t^2 + \cdots & 0 & 0 \\ 0 & 1 + λ_2 t + \frac 12 λ_2^2 t^2 + \cdots & 0 \\ 0 & 0 & 1 + λ_3 t + \frac 12 λ_3^2 t^2 + \cdots} \\[10pt] &= \MATRIX{e^{λ_1t} & 0 & 0 \\ 0 & e^{λ_2 t} & 0 \\ 0 & 0 & e^{λ_3 t}} \end{align*}\]
Thus, the solution of the vector differential equation \(\dot x(t) = A x(t)\) is given by \[ x(t) = e^{At} x_0 = \MATRIX{e^{λ_1t} & 0 & 0 \\ 0 & e^{λ_2 t} & 0 \\ 0 & 0 & e^{λ_3 t}} \MATRIX{α_1 \\ α_2 \\ α_3} = \MATRIX{α_1 e^{λ_1 t} \\ α_2 e^{λ_2 t} \\ α_3 e^{λ_3 t}}. \]
But outside of such few special cases, computing \(e^{At}\) via definition \(\eqref{eq:matrix-exponential}\) is not computationally feasible. We now present computationally efficient methods to compute the matrix exponential.
4.2 Computing matrix exponential
4.2.1 Method 1: Eigenvalue diagonalization method
As illustrated by Example 4.3, compute matrix exponential is easy for diagonal matrix. So, if the matrix \(A\) is diagonalizable (i.e., has distinct eigenvalues) we can do a change of coordinates and compute the matrix exponential in the eigen-coordinates.
In particular, suppose \(A\) has distinct eivenvalues \(λ_1, \dots, λ_n\) and \(v_1, \dots, v_n\) are the corresponding eigvenvectors. Thus, $ A v_k = λ_k v_k$ for all \(k \in \{1,\dots, n\}\). Writing this in matrix form, we have \[\begin{align*} A \MATRIX{v_1 & v_2 & \cdots & v_n} &= \MATRIX{A v_1 & A v_2 & \cdots & A v_n} \\ &= \MATRIX{λ_1 v_1 & λ_2 v_2 & \cdots & λ_n v_n} \\ &= \MATRIX{v_1 & v_2 & \cdots & v_n} \MATRIX{ λ_1 & 0 & \cdots & 0 \\ 0 & λ_2 & \ddots & 0 \\ 0 & \cdots & \ddots & 0 \\ 0 & \cdots & \cdots & λ_n} \end{align*}\] Now define \[T = \MATRIX{v_1 & v_2 & \cdots & v_n} \quad\text{and}\quad Λ = \MATRIX{ λ_1 & 0 & \cdots & 0 \\ 0 & λ_2 & \ddots & 0 \\ 0 & \cdots & \ddots & 0 \\ 0 & \cdots & \cdots & λ_n}. \] So, the above equation can be writen as \(AT = T Λ\) or \[ \bbox[5pt,border: 1px solid] {T^{-1} A T = Λ} \] Observe that
- \(A = T Λ T^{-1}\)
- \(A^2 = T Λ^2 T^{-1}\)
- \(A^3 = T Λ^3 T^{-1}\)
- and so on.
Therefore, \[\begin{align*} e^{At} &= I + At + \frac{A^2 t^2}{2!} + \frac{A^3 t^3}{3!} + \cdots \\ &= I + T Λ T^{-1} t + \frac{T Λ^2 T^{-1} t^2}{2!} + \frac{T Λ^3 T^{-1} t^3}{3!} + \cdots \\ &= T\bigl( I + Λt + \frac{Λ^2 t^2}{2!} + \frac{Λ^3 t^3}{3!} + \cdots\bigr) T^{-1} \\ &= T e^{Λ t} T^{-1} \end{align*}\]
Thus, once we know the eigenvalues and eigenvectors of \(A\) (and if all eigenvalues are distinct), then \[ \bbox[5pt,border: 1px solid] {e^{At} = T e^{Λt} T^{-1}} \]
Exercise 4.1 Use the eigenvalue diagonalizable method to compute \(e^{At}\) for \(A = \MATRIX{-6 & 4 \\ -5 & 3 }\).
We start by computing the eigenvalues and eigenvectors of \(A\).
To compute the eigenvalues: \[sI - A = \MATRIX{s + 6 & -4 \\ 5 & s-3 }\] Therefore, the characteristic equation is \[\det(sI - A) = s^2 + 3s + 2 = (s+1)(s+2) \] Hence, the eigenvalues of \(A\) are \(λ_1 = -1\) and \(λ_2 = -2\).
We now compute the eigenvectors. Recall that for any eigenvalue \(λ\), the eigen-vector satisfies \((λI - A) v = 0\). We start wtih \(λ_1 = -1\). Then, \[ \MATRIX{5 & -4 \\ 5 & -4 } \MATRIX{ v_{11} \\ v_{12} } = 0 \] We set \(v_{11} = 4\). Then, \(v_{12} = 5\). Thus, the eigenvector \[ v_1 = \MATRIX{ 4 \\ 5}. \]
Similarly, for \(λ_2 = -2\), we have \[ \MATRIX{4 & -4 \\ 5 & -5 } \MATRIX{ v_{21} \\ v_{22} } = 0 \] We set \(v_{21} = 1\). Then, \(v_{22} = 1\). Thus, the eigenvector \[ v_2 = \MATRIX{ 1 \\ 1}. \]
Thus, \(T = \MATRIX{ 4 & 1 \\ 5 & 1 }\) and therefore \(T^{-1} = \MATRIX{-1 & 1 \\ 5 & -4 }\). Hence, \[ e^{At} = \MATRIX{ 4 & 1 \\ 5 & 1 } \MATRIX{e^{-t} & 0 \\ 0 & e^{-2t} } \MATRIX{-1 & 1 \\ 5 & -4 } = \MATRIX{ -4 e^{-t} + 5 e^{-2t} & 4 e^{-t} - 4 e^{-2t} \\ -5 e^{-t} + 5 e^{-2t} & 5 e^{-t} - 4 e^{-2t} } \]
Exercise 4.2 Use the eigenvalue diagonalizable method to compute \(e^{At}\) for \(A = \MATRIX{1 & 2 \\ -1 & 4}\).
We start by computing the eigenvalues and eigenvectors of \(A\).
To compute the eigenvalues: \[sI - A = \MATRIX{s - 1 & -2 \\ 1 & s-4 }\] Therefore, the characteristic equation is \[\det(sI - A) = s^2 - 5s + 6 = (s-2)(s-3). \] Hence, the eigenvalues of \(A\) are \(λ_1 = 2\) and \(λ_2 = 3\).
We now compute the eigenvectors. Recall that for any eigenvalue \(λ\), the eigen-vector satisfies \((λI - A) v = 0\). We start wtih \(λ_1 = 2\). Then, \[ \MATRIX{1 & -2 \\ 1 & -2 } \MATRIX{ v_{11} \\ v_{12} } = 0 \] We set \(v_{11} = 2\). Then, \(v_{12} = 1\). Thus, the eigenvector \[ v_1 = \MATRIX{ 2 \\ 1}. \]
Similarly, for \(λ_2 = 3\), we have \[ \MATRIX{2 & -2 \\ 1 & -1 } \MATRIX{ v_{21} \\ v_{22} } = 0 \] We set \(v_{21} = 1\). Then, \(v_{22} = 1\). Thus, the eigenvector \[ v_2 = \MATRIX{ 1 \\ 1}. \]
Thus, \(T = \MATRIX{ 2 & 1 \\ 1 & 1 }\) and therefore \(T^{-1} = \MATRIX{1 & -1 \\ -1 & 2 }\). Hence, \[ e^{At} = \MATRIX{ 2 & 1 \\ 1 & 1 } \MATRIX{e^{2t} & 0 \\ 0 & e^{3t} } \MATRIX{1 & -1 \\ -1 & 2 } = \MATRIX{ 2 e^{2t} - e^{3t} & -2 e^{2t} + 2 e^{3t} \\ e^{2t} - e^{3t} & - e^{2t} + 2 e^{3t} } \]
4.2.2 Method 2: Laplace transform method
Recall that we have \[ \dot x(t) = A x(t). \] Taking (unilateral) Laplace transforms of both sides gives \[ s X(s) - x_0 = A X(s). \] Rearranging terms, we get \[ (sI - A)X(s) = x_0 \implies X(s) = (sI-A)^{-1} x_0. \] Taking inverse Laplace transforms, we get \[x(t) = \mathcal L^{-1}( (sI - A)^{-1} ) x_0. \] Comparing it with the solution obtained earlier, we get \[ \bbox[5pt,border: 1px solid] { e^{At} = \mathcal L^{-1}( (sI - A)^{-1} ).} \] In the above, we interpret the inverse Laplace transform of a matrix to mean the inverse Laplace transform of each entry.
Exercise 4.3 Use the Laplace transform method to compute \(e^{At}\) for \(A = \MATRIX{-6 & 4 \\ -5 & 3 }\).
We first compute \[sI - A = \MATRIX{s + 6 & -4 \\ 5 & s-3 }\] Therefore, \[\det(sI - A) = s^2 + 3s + 2 = (s+1)(s+2) \] and \[ (sI - A)^{-1} = \MATRIX{ \dfrac{s-3}{(s+1)(s+2)} & \dfrac{4}{(s+1)(s+2)} \\ \dfrac{-5}{(s+1)(s+2)} & \dfrac{s+6}{(s+1)(s+2)} } \] We now do partial fraction expansion of each term: \[ (sI - A)^{-1} = \MATRIX{ \dfrac{-4}{s+1} + \dfrac{5}{s+2} & \dfrac{ 4}{s+1} - \dfrac{4}{s+2} \\ \dfrac{-5}{s+1} + \dfrac{5}{s+2} & \dfrac{5}{s+1} - \dfrac{4}{s+2} } \] Taking the inverse Laplace transform of each term, we get \[ e^{At} = \mathcal L^{-1}(sI - A)^{-1} = \MATRIX{ -4 e^{-t} + 5 e^{-2t} & 4 e^{-t} - 4 e^{-2t} \\ -5 e^{-t} + 5 e^{-2t} & 5 e^{-t} - 4 e^{-2t}} \]
4.3 Internal stability
As stated in the beginning of the lecture, matrix exponential allows us to solve the vector differential equation \[ \dot x(t) = A x(t), \quad \text{initial condition } x(0) = x_0 \] in the same manner as a scalar linear differential equation. The solution is given by \[ x(t) = e^{At} x_0. \]
Suppose \(A\) has distinct eigenvalues \((λ_1, \dots, λ_n)\) with corresponding (linearly independent) eigenvalues \((v_1, \dots, v_n)\).
Recall that for any eigenvector \(v_j\), \(A^k v_j = λ_j^k v_j\). Thus, if the system starts from the intial condition \(x_0 = v_j\), then \[\begin{align*} x(t) &= e^{At} x_0 = e^{At} v_j \\ &= \biggl[ I + At + \frac{A^2 t^2}{2!} + \frac{A^3 t^3}{3!} + \cdots \biggr] v_j \\ &= \biggl[ v_j + λ_j t v_j + \frac{λ_j^2 t^2}{2!} v_j + \frac{λ_j^3 t^3 }{3!} v_j + \cdots \biggr] \\ &= e^{λ_jt} v_j \end{align*}\] Therefore, if we start along an eigenvector of \(A\), the system trajectory \(x(t)\) remains along that direction, with the length scaling as \(e^{λ_jt}\).
In general, since the eigenvectors \((v_1, \dots, v_n)\) are linearly independent, an arbitrary initial condition \(x_0\) can be written as \[ x_0 = α_1 v_1 + \cdots + α_n v_n. \] Therefore, \[\begin{align*} x(t) &= e^{At} x_0 = e^{At} \bigl[ α_1 v_1 + \cdots + α_n v_n \bigr] \\ &= α_1 e^{At} v_1 + \cdots + α_n e^{At} v_n \\ &= α_1 e^{λ_1t} v_1 + \cdots + α_n e^{λ_n t} v_n \end{align*}\]
Thus, the response of the dynamical system \(\dot x(t) = A x(t)\) is a combination of motions along the eigenvectors of the matrix \(A\). Each eigendirection is called a mode of the system. A particular mode is excited by choosing an initial condition to have a component \(α_i\) along that eigendirection.
An implication of the above is the following: the state \(x(t) \to 0\) as \(t \to ∞\) for all initial states \(x_0\) if and only if all eigenvalues of \(A\) lie in the open left-hand plane.
A SSM where all eigenvalues of the \(A\) matrix lie in the open-left hand plane is called internally stable.
The TF of a SSM is given by \[ G(s) = C(sI - A)^{-1}B = \frac{C \text{adj}(sI - A) B}{\det(sI - A)}. \] The elements of adjoint of \((sI - A)\) are all polynomials in \(s\). Thus, the numerator is a polynomial in \(s\). So, is the denominator. Moreover, the denominator equals the characteristic equation of \(A\); thus, the roots of the denominator are the eigenvalues of \(A\).
The numerator and denominator may have common roots that cancel each other. So, in general, \[ \{\text{eigenvalues of } A \} \supseteq \{\text{poles of } G(s) \}. \] Hence, if \(A\) is internally stable (i.e., all its eigenvalues are in the open left-hand plane) then \(G(s)\) is BIBO stable (i.e., all its poles are in the open left-hand plane).
However, the converse is not true there may be pole-zero cancellations. For example, consider \(A = \MATRIX{-1 & 0 \\ 0 & 2 }\), \(B = \MATRIX{1 \\ 0}\), \(C = \MATRIX{1 & 0}\). Then, \[ G(s) = \frac{s-2}{(s+1)(s-2)} = \frac{1}{s+1}. \] Notice that the TF is BIBO stable but the SSM is not internally stable! Any initial condition that activates the mode corresponding to eigenvalue \(2\) will cause \(x(t)\) to diverge to infinity.
4.4 Time response of state space models
Now consider a SSM given by \[\begin{align*} \dot x(t) &= A x(t) + B u(t) \\ y(t) &= C x(t) \end{align*}\]
Suppose the system starts at \(t=0\) with an initial state \(x(0) = x_0\) and we apply the input \(u(t)\). How do we find the output?
Taking Laplace transform of the SSM, we get: \[\begin{align*} sX(s) - x_0 &= A X(s) + B U(s) \\ Y(s) &= C X(s) \end{align*}\] Solving for \(X(s)\), we get \[X(s) = (sI - A)^{-1} x_0 + (sI - A)^{-1} B U(s). \] Substituting in \(Y(s) = C X(s)\), we get \[ \bbox[5pt,border: 1px solid] {Y(s) = \underbrace{C (sI - A)^{-1} x_0}_{\text{zero-input response}} + \underbrace{C (sI - A)^{-1} B U(s)}_{\text{zero-state response}}} \] We can the use the above expression compute \(y(t)\) by taking the inverse Laplace transform.
It is sometimes useful to write the expression in time domain (but we will not use this expression for computations). To do so, recall that \(C(sI - A)^{-1}B\) is the transfer function \(G(s)\) of the system. Therefore, its inverse Laplace transform is the impulse response: \[ g(t) = C e^{At} B. \] Then, we can compute the inverse Laplace transform of \(G(s) U(s)\) using the convolution formula and write \[ \bbox[5pt,border: 1px solid] {y(t) = \underbrace{C e^{At} x_0}_{\text{zero-input response}} + \underbrace{\int_{0}^t C e^{A (t-τ) } B u(τ) d τ}_{\text{zero-state response}}} \]