Simulating linear systems with non-zero initial conditions in state space

Matlab’s lsim function for simulating linear systems will give you the option to provide an initial condition if your system is in state-space but not for transfer-functions.

One problem is that you either start in state-space, in which case you already know your initial conditions in terms of the initial state, or you have somehow to figure out how to convert initial conditions stated in terms of the system output as an initial state. This is a common source of confusion to readers that are not familiar with state-space techniques.

Before proceeding, one might want to take a look at Chapter 5 for an introductory overview of state-space systems.

Let’s start by looking at how you can use Matlab to transform a transfer-function, say

$$G(s) = \frac{1}{s^2 + s + 1},$$

into one of its many possible state-space realizations. For that one can use ss or tf2ss. For example, after running

Gtf = tf([1], [1 1 1])
Gss = ss(Gtf)

the variable Gss will contain a state-space realization of the transfer-function Gtf. You can check that by using tf to transform it back as in

tf(Gss)

Now recall that in state-space Gss is described by the matrices $(A,B,C,D)$, accessible as Gss.a, Gss.b, Gss.c, and Gss.d, or using

[a,b,c,d] = ssdata(Gss)

Once in state-space, the function lsim will accept an initial state option, as in lsim(Gss, u, t, x0), where x0 is a vector with the initial condition. If you have a vector of initial conditions $y_0$, how can you transform it to the corresponding initial state $x_0$? That is the question addressed next.

We will use linearity and look only at the part of the response due to the initial condition, that is we shall consider a zero input in the next paragraphs. As seen in Chapter 5, the response to a zero input and nonzero initial state $x(0) = x_0$ is given by

$$y(t)=C e^{At} x_0$$

Say that your system is of order $n$ you have initial conditions $(y(0), \dot{y}(0), \cdots, y^{(n-1)}(0))$ and you would like to know how to calculate a corresponding initial state, $x_0$. Because

$$\dot{y}(t)=C A e^{A t} x_0, \quad \cdots \quad, y^{(n-1)}(t) = C A^{n-1} e^{A t} x_0$$

one can do that by solving a system of linear equations. Here’s how it works. Simply rearrange

$$y(0) = C x_0, \quad \dot{y}(0) = C A x_0, \quad \cdots \quad, y^{(n-1)}=CA^{n-1} x_0$$

in the form of a set of linear equations

$$\mathbf{O} x_0 = \mathbf{y}_0, \quad \mathbf{O}=\begin{bmatrix}C\\CA\\\vdots\\CA^{n-1}\end{bmatrix}, \quad \mathbf{y}_0 =\begin{pmatrix}y(0)\\\dot{y}(0)\\\vdots\\y^{(n-1)}(0)\end{pmatrix}$$

and solve for $x_0 = \mathbf{O}^{-1} \mathbf{y}_0$ to calculate a suitable initial state $x_0$. The above matrix $\mathbf{O}$ is known as the observability matrix and it is so important that Matlab has a function to calculate it: obsv. The following excerpt puts it all together to simulate a step response from a non-zero initial condition.

Gtf = tf([1],[1 1 1]);
Gss = ss(Gtf);
y0 = [-1; 0]; % initial conditions y(0) = -1, y'(0) = 0
O = obsv(Gss.a, Gss.c);
x0 = O \ y0;
T = 2;
t = linspace(0, T, 100);
u = ones(size(t));
y = lsim(Gss, u, t, x0);

P.S.: You might be asking when is the above matrix $\mathbf{O}$ invertible? The answer is it will be invertible if the corresponding state-space realization is observable (See Section 5.3). For SISO systems this will correspond to not having pole-zero cancellations.