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 '20 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 '21 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 '21 at 16:07

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 '17 at 12:42