How do I simulate a linearized state space model with MATLAB's lsim command?

Question

Consider that you got these equations:

$$2.5 \dot{x_1} = 0.0050 - u_1*0.01*\sqrt{x_1}$$ $$2.5 \dot{x_2} = u_1*0.01*\sqrt{x_1} - u_2*0.01*\sqrt{x_2}$$

That will result a nonlinear state space model, due to $\sqrt{}$. Linearize that model in $$x_1 = 1.5, x_2 = 0.5 , u_1 = 0.40825, u_2 = 0.70711$$

And you will get this state space model: $$ \begin{bmatrix} \dot{x_1} - 0\\ \dot{x_2} - 0 \end{bmatrix} = \begin{bmatrix} -0.00167 & 0 \\ 0.00167 & -0.005 \end{bmatrix} \begin{bmatrix} x_1 - 1.5\\ x_2 - 0.5 \end{bmatrix} + \begin{bmatrix} -0.0122 & 0 \\ 0.0122 & -0.00707 \end{bmatrix} \begin{bmatrix} u_1 - 0.40825\\ u_2 - 0.70711 \end{bmatrix}\ $$

Question: $$ $$ How do I simulate this linearized state space equation by using MATLAb's / Octave's command lsim ? https://se.mathworks.com/help/control/ref/lsim.html

Answer 1

Well, you're asking specifically how to use the lsim command and not about the linearization of your model, so I'm working from the assumption that your model is correctly linearized. SO, starting there, you would use the lsim command as follows:

lsim(sys,u,t,x0)

From the documentation:

The vector t specifies the time samples for the simulation ... and consists of regularly spaced time samples:

t = 0:dt:Tfinal;

Then,

The input u is an array having as many rows as time samples (length(t)) and as many columns as system inputs. For instance, if sys is a SISO system, then u is a t-by-1 vector. If sys has three inputs, then u is a t-by-3 array. Each row u(i,:) specifies the input value(s) at the time sample t(i).

And finally,

x0 is an initial condition for the system states. This syntax applies only when sys is a state-space model.

So, first define your system. Your question isn't formatted in a state space form that Matlab can accept because your state and input vectors $\bf{x}$ and $\bf{u}$ have constants in them. I'm going to make this a little easier on me by just masking everything behind constants. Multiply everything out and you get:

$$ \left[\begin{array}{cc} a & b \\ c & d \\ \end{array}\right]\left[\begin{array}{c} x_1 + e \\ x_2 + f \\ \end{array}\right] = \left[\begin{array}{cc} a(x_1 + e) + b(x_2 + f) \\ c(x_1 + e) + d(x_2 + f) \\ \end{array}\right] $$

$$ \left[\begin{array}{cc} a(x_1 + e) + b(x_2 + f) \\ c(x_1 + e) + d(x_2 + f) \\ \end{array}\right] = \left[\begin{array}{cc} a x_1 + ae + b x_2 + bf \\ c x_1 + ce + d x_2 + df \\ \end{array}\right] = \left[\begin{array}{cc} a x_1 + b x_2 + ae + bf \\ c x_1 + d x_2 + ce + df \\ \end{array}\right] $$

So, the tricky bit here might be to recognize that those last terms, ae+bf and ce+df are multiplied by one, so really you wind up with something like:

$$ \left[\begin{array}{cc} (a)x_1 + (b)x_2 + (ae + bf)1 \\ (c)x_1 + (d)x_2 + (ce + df)1 \\ \end{array}\right] $$

Now, if you were to have a variable that didn't change over time, you could set its initial value to one and then use that. You do this by adding a state to your state matrix $\bf{A}$ and leaving the corresponding row for that state equal to zeros. Now the new state ($x_3$) doesn't change over time, and as long as you set its initial value to one, you are left with the form you need in proper state space form:

$$ \left[\begin{array}{cc} (a)x_1 + (b)x_2 + (ae + bf)x_3 \\ (c)x_1 + (d)x_2 + (ce + df)x_3 \\ \end{array}\right] $$

You can easily separate that out: $$ \left[\begin{array}{ccc} a & b & ae + bf \\ c & d & ce + df \\ 0 & 0 & 0 \\ \end{array}\right]\left[\begin{array}{c} x_1 \\ x_2 \\ x_3 \\ \end{array}\right] $$

Now you need to be careful to include the constants you also have in your input vector:

$$ \left[\begin{array}{cc} g & h \\ j & k \\ \end{array}\right]\left[\begin{array}{c} u_1 + m \\ u_2 + n \\ \end{array}\right] = \left[\begin{array}{cc} g(u_1 + m) + h(u_2 + n) \\ j(u_1 + m) + k(u_2 + n) \\ \end{array}\right] $$

Same steps as before:

$$ \left[\begin{array}{cc} g(u_1 + m) + h(u_2 + n) \\ j(u_1 + m) + k(u_2 + n) \\ \end{array}\right] = \left[\begin{array}{cc} g(u_1) + gm + h(u_2) + hn \\ j(u_1) + jm + k(u_2) + kn \\ \end{array}\right] $$ $$ = \left[\begin{array}{cc} g(u_1) + h(u_2) + (gm + hn)x_3 \\ j(u_1) + k(u_2) + (jm + kn)x_3 \\ \end{array}\right] $$

Put it all together:

$$ \left[\begin{array}{c} \dot{x_1} \\ \dot{x_2} \\ \dot{x_3} \\ \end{array}\right] = \left[\begin{array}{ccc} a & b & ae + bf + gm + hn \\ c & d & ce + df + jm + kn \\ 0 & 0 & 0 \\ \end{array}\right]\left[\begin{array}{c} x_1 \\ x_2 \\ x_3 \\ \end{array}\right] + \left[\begin{array}{cc} g & h \\ j & k \\ 0 & 0 \\ \end{array}\right]\left[\begin{array}{c} u_1 \\ u_2 \\ \end{array}\right] $$

This was probably the tricky part in getting your equations to the point of being able to use lsim. From here, everything else is self-explanatory.

1. Determine your sample interval dt or your sample rate fs and calculate the sample interval as dt = 1/fs;.
2. Define your time vector t = 0:dt:tFinal;.
3. Define your initial conditions. If you want to start with $x_1$ and $x_2$ at zero, then you would use x0 = [0;0;1]; because again, $x_3$ must be one.
4. Define your inputs. If you want to step only $u_1$, then you would use u = [ones(length(t),1),zeros(length(t),1)];. Swap those to step $u_2$, or whatever other combination you want.
5. Define your matrices and then your system; A = [...]; B = [...]; etc., then sys = ss(A,B,C,D);, where the matrices are as I've defined above.
6. Invoke lsim - lsim(sys,u,t,x0);.

Done!