# Bilinear transformation of continuous time state space system

I'm trying to understand the derivation of the bilinear transform for a set of continuous time state-space matrices. I've found plenty of websites which list steps to perform the conversion (here 1 or here 2 or here 3) - but haven't found any derivation or insight into how they came about. When I attempt to do it myself, I keep getting terms in z which I can't seem to get rid of.

To make the problem clear, if I define the bilinear transform as:

$$s=\frac{\alpha\left(z-1\right)}{z+1}$$

Where alpha is usually given as 2/T.

The transfer function of a continuous-time state-space system can be given as:

$$H(s) = \mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B} + \mathbf{D}$$

How do I derive the expressions for the matrices of the discrete time version bilinear transformed model:

$$H(z) = \mathbf{C_d}(z\mathbf{I}-\mathbf{A_d})^{-1}\mathbf{B_d} + \mathbf{D_d}$$

When I make the variable change, I keep getting extra terms in z which I don't know what to do with. I'm kinda hoping this is just me failing to have a strong grasp of matrix operations. Any insight (including if I am going about this completely the wrong way) would be really appreciated!

Thanks to Klaz for posting an answer to this first. I managed to work through the problem today and thought I would post my own alternate working here.

Given the typical state-space equations:

$$s\mathbf{Q}(s) = \mathbf{A}\mathbf{Q}(s) + \mathbf{B}\mathbf{X}(s)$$ $$\mathbf{Y}(s) = \mathbf{C}\mathbf{Q}(s) + \mathbf{D}\mathbf{X}(s)$$

Focus on the state update equation first and make the bilinear substitution:

$$\frac{\alpha\left(z-1\right)}{z+1} \mathbf{Q}(z) = \mathbf{A}\mathbf{Q}(z) + \mathbf{B}\mathbf{X}(z)$$

We manipulate this expression to make the LHS $z\mathbf{Q}(z)$:

$$\alpha\left(z-1\right) \mathbf{Q}(z) = \mathbf{A}\left(z+1\right)\mathbf{Q}(z) + \mathbf{B}\left(z+1\right)\mathbf{X}(z)$$

$$z\left(\alpha\mathbf{I}-A\right)\mathbf{Q}(z) = \left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \mathbf{B}\left(z+1\right)\mathbf{X}(z)$$

$$z\mathbf{Q}(z) = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B}\left(z+1\right)\mathbf{X}(z)$$

This instantly gives us $\mathbf{A_d}$ as:

$$\mathbf{A_d} = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right)$$

We know that the linear $(z+1)$ term of $\mathbf{X}(z)$ can be applied wherever $\mathbf{Q}(z)$ is used. We apply this in the output expression to give the new output discrete equation as:

$$z\mathbf{Q}(z) = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B}\mathbf{X}(z)$$ $$\mathbf{Y}(z) = \mathbf{C}\left(z+1\right)\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$ $$\mathbf{Y}(z) = \mathbf{C}z\mathbf{Q}(z) + \mathbf{C}\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$

We can now pull out $\mathbf{B_d}$ as:

$$\mathbf{B_d} = \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B}$$

Because we know an expression for $z\mathbf{Q}(z)$, we can formulate $\mathbf{Y}(z)$ as:

$$\mathbf{Y}(z) = \mathbf{C}\left(\mathbf{A_d}\mathbf{Q}(z) + \mathbf{B_d} \mathbf{X}(z) \right) + \mathbf{C}\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$

$$\mathbf{Y}(z) = \mathbf{C}\left(\mathbf{A_d}+\mathbf{I}\right)\mathbf{Q}(z) + \left(\mathbf{C}\mathbf{B_d} + \mathbf{D}\right)\mathbf{X}(z)$$

Giving:

$$\mathbf{C_d} = \mathbf{C}\left(\mathbf{A_d}+\mathbf{I}\right)$$ $$\mathbf{D_d} = \mathbf{C}\mathbf{B_d} + \mathbf{D}$$

• I noticed that your solution is slight different from the result in 1, I also find different solution here. Obviously, moving the linear term $(z+1)$ from $B$ matrix to $C$ matrix would not change the total transfer function. But would inner state matrix $Q(z)$ lose its structure? It is no longer discretization of original $Q(s)$ but something else? Sep 14, 2020 at 3:35
• I don't understand the part: We know that the linear $(z+1)$ term of $\mathbf{X}(z)$ can be applied wherever $\mathbf{Q}(z)$ is used. Why? Aug 10, 2021 at 15:35
• First, if you're doing control systems then your knowledge of the parameters is inexact, so you may as well use something that makes the math easy. The forward Euler method works. I.e. $x_n = (\mathbf I + T_s \mathbf A)x_{n-1} + T_s \mathbf B$. Second, if you feel that exactitude is important, and since it's 2021 and we're surrounded by a sea of computers, use the exact method. Note that it's easy to get over the problem with $\mathbf A$ being singular. Nov 15, 2021 at 16:07
• @JuanGonzalezBurgos yes, this was very poor wording on my part. If 𝐘(𝑧)=𝐀(𝑧) and $\mathbf{A}(z)=(z+1)\mathbf{X}(z)$, because the expressions are linear, we can move the (𝑧+1) from the second expression into the first e.g. 𝐘(𝑧)=(𝑧+1)𝐀(𝑧) and $\mathbf{A}(z)=\mathbf{X}(z)$. This is equivalent to a change of variable - but I was lazy and didn't change the variable. Jul 27, 2022 at 9:35
• @TimWescott thanks for the additional comments. The particular use-case surrounding my question didn't arise from control systems, but from attempting to build algorithms that would efficiently realise parallel-form (as opposed to cascade) discrete implementations of high-order filters. State-space representation is IMO widely under-appreciated for it's uses outside of control systems. Even regular digital IIR filters can have their performance dramatically improved by designing the filter in state space. Jul 27, 2022 at 9:43

I've had the same question last week, but I've managed to find how to derive it (getting rid of those $$z$$ terms is indeed tricky). I will give here detailed demonstration of how to arrive to the result given in 1 (with, in your notation, $$\alpha = 2 \lambda$$).

So we define our new discrete-time function transfer as

$$\begin{array}{rcl} H_d(z) &=& H(\frac{\alpha\left(z-1\right)}{z+1}) \\ &=& \mathbf{C} \left[ \frac{\alpha\left(z-1\right)}{z+1}\mathbf{I}-\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\ &=& \frac{z+1}{\alpha} \mathbf{C} \left[ \left(z-1\right)\mathbf{I}- \frac{z+1}{\alpha}\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\ &=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \left(\mathbf{I} - \frac{1}{\alpha} \mathbf{A} \right) - \left(\mathbf{I} + \frac{1}{\alpha} \mathbf{A}\right) \right]^{-1}\mathbf{B} + \mathbf{D} \\ \end{array}$$

For sake of notation, let $$\mathbf{P} = \mathbf{I} - \frac{1}{\alpha} \mathbf{A}$$ and $$\mathbf{Q} = \mathbf{I} + \frac{1}{\alpha} \mathbf{A}$$ (important remark useful later on: $$\mathbf{P} + \mathbf{Q} = 2 \mathbf{I}$$). Then

$$\begin{array}{rcl} H_d(z) &=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{P} - \mathbf{Q} \right]^{-1}\mathbf{B} + \mathbf{D} \\ &=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \mathbf{P}^{-1} \mathbf{B} + \mathbf{D} \\ &=& \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \left( \sqrt{\frac{2}{\alpha}} \mathbf{P}^{-1} \mathbf{B} \right) + \mathbf{D} \\ &=& \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \end{array}$$

where we defined

$$\mathbf{A_d} = \mathbf{P}^{-1} \mathbf{Q}$$ and $$\mathbf{B_d} = \sqrt{\frac{2}{\alpha}} \mathbf{P}^{-1} \mathbf{B}$$. Furthermore, we have

$$\begin{array}{rcl} \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z + 1 \right) \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z \mathbf{I} + \mathbf{I} + \mathbf{A_d} - \mathbf{A_d} \right) \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \left( \mathbf{I} + \mathbf{A_d} \right) \right] \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \mathbf{P}^{-1} \left( \mathbf{P} + \mathbf{Q} \right) \right] \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + 2 \mathbf{P}^{-1} \right] \end{array}$$

which gives us

$$\begin{array}{rcl} H_d(z) &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + 2 \mathbf{P}^{-1} \right] \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d} + \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1} \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \\ &=& \mathbf{C_d}(z\mathbf{I}-\mathbf{A_d})^{-1}\mathbf{B_d} + \mathbf{D_d} \end{array}$$ with $$\mathbf{C_d} = \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1}$$ and $$\mathbf{D_d} = \mathbf{D} + \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d}$$. Replacing $$\mathbf{P}$$ and $$\alpha = 2 \lambda$$ gives the results in 1.

PS: the same idea can be used to prove the discrete-time state-space representation found using Generalized Bilinear Transform ($$s \leftarrow \alpha \frac{z - 1}{\beta z + \left(1-\beta\right)}$$) or First-order Holder methods.

• Thanks Klaz for taking the time to post a solution! I managed to work through this one today after almost giving up. I think I went about it in a slightly different way and have added my own working to the original question. Nov 10, 2017 at 12:42

I wanted to add to this some specifics on the MATLAB implementation.

Continuous time system: $$\dot{x} = A x + B u$$

Laplace transform: $$sx = Ax + Bu$$

Tustin discretization: $$\frac{\alpha(z-1)}{z+1}x = Ax + Bu$$ where $$\alpha = 2f_\text{s}$$, $$f_\text{s}$$ being the sampling frequency.

When $$z \ne -1$$, multiplying through by $$(z+1)$$, and rearranging, $$z\left[x - \frac{1}{\alpha}(Ax+Bu)\right] = x + \frac{1}{\alpha}(Ax + Bu)$$

Motivated by this, do the following coordinate transformation (this is what is done in MATLAB, see here): \begin{aligned} x_\text{d} &= x - \frac{1}{\alpha}(Ax+Bu) \\ &= \left(I - \frac{1}{\alpha}A\right)x - \frac{1}{\alpha}Bu \end{aligned}

Then the inverse coordinate transformation, $$x = \left(I - \frac{1}{\alpha}A\right)^{-1} \left(x_\text{d}+\frac{1}{\alpha}Bu\right)$$

Substituting back, $$z x_\text{d} = \left(I + \frac{1}{\alpha}A\right) \left(I - \frac{1}{\alpha}A\right)^{-1}x_\text{d} + \left[\left(I + \frac{1}{\alpha}A\right) \left(I - \frac{1}{\alpha}A\right)^{-1} + I\right] \frac{1}{\alpha}Bu$$

Define \begin{aligned} A_\text{d} &= \left(I + \frac{1}{\alpha}A\right) \left(I - \frac{1}{\alpha}A\right)^{-1} \\ B_\text{d} &= \frac{1}{\alpha}(A_\text{d} + I)B \end{aligned}

Then, $$x_\text{d}^{k+1} = A_\text{d}x_\text{d}^k + B_\text{d}u^k$$

To check if this is indeed what $$\texttt{c2d}$$ in MATLAB is doing, try the following:

% Construct a second order system to test
f = 5; % frequency
zet = 0.05; % damping ratio
w = 2*pi*f; % circular frequency

% state space representation
A = [0 1; -w^2 -2*zet*w];
B = [0; 1];
C = [1 0];
D = [];
sys = ss(A, B, C, D);

% Discrete time system using Tustin
Ts = 0.01; % sampling time (100Hz)
sysd = c2d(sys, 0.01, 'Tustin');

% Construct discrete time matrices using the formulas derived
alph = 2/Ts;

This gives $$A_\text{d} = \begin{bmatrix}0.9526 & 0.0096 \\-9.4865 & 0.9224 \end{bmatrix}$$ and $$B_\text{d} = \begin{bmatrix}0.0000 \\ 0.0096 \end{bmatrix}$$