7
$\begingroup$

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!

--- Edit with answer.

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} $$

Improve this question
$\endgroup$
7
  • $\begingroup$ 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? $\endgroup$ Commented Sep 14, 2020 at 3:35
  • $\begingroup$ 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? $\endgroup$ Commented Aug 10, 2021 at 15:35
  • $\begingroup$ 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. $\endgroup$ Commented Nov 15, 2021 at 16:07
  • $\begingroup$ @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. $\endgroup$ Commented Jul 27, 2022 at 9:35
  • $\begingroup$ @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. $\endgroup$ Commented Jul 27, 2022 at 9:43

2 Answers 2

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Commented Nov 10, 2017 at 12:42
4
$\begingroup$

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 $$

Also, the output $$ \begin{aligned} y &= Cx + Du \\ &= C\left[ \left(I - \frac{1}{\alpha}A\right)^{-1} \left(x_\text{d}+\frac{1}{\alpha}Bu\right) \right] + Du \\ &= C \left(I - \frac{1}{\alpha}A\right)^{-1}x_\text{d} + \left[ \frac{1}{\alpha}C\left(I-\frac{1}{\alpha}A\right)^{-1}B + D \right] u \end{aligned} $$

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 \\ C_\text{d} &= C \left(I - \frac{1}{\alpha}A\right)^{-1} \\ D_\text{d} &= D + \frac{1}{\alpha}C B \end{aligned} $$

Then, $$ \begin{aligned} x_\text{d}^{k+1} &= A_\text{d}x_\text{d}^k + B_\text{d}u^k \\ y^k &= C_\text{d}x_\text{d}^k + D_\text{d}u^k \end{aligned} $$

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 = 0; 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; Ad = (eye(2)+A/alph)/(eye(2)-A/alph); Bd = (1/alph)*(Ad + eye(2))*B; Cd = C/(eye(2)-A/alph); Dd = D + Cd*B/alph; % Compare fprintf('Compare Ad:\n') disp([Ad sysd.A]) fprintf('Compare Bd:\n') disp([Bd sysd.B]) fprintf('Compare Cd:\n') disp([Cd sysd.C]) fprintf('Compare Bd:\n') disp([Dd sysd.D])

This gives $$A_\text{d} = \begin{bmatrix}0.9526 & 0.0096 \\-9.4865 & 0.9224 \end{bmatrix}$$ $$B_\text{d} = \begin{bmatrix}0.0000 \\ 0.0096 \end{bmatrix}$$ $$C_\text{d} = \begin{bmatrix}0.9763 & 0.0048\end{bmatrix}$$ and $$D_\text{d} = 2.403\times 10^{-5}$$

Improve this answer
$\endgroup$
4
  • $\begingroup$ This is a good answer, but I had already upvoted it. It's actually a really good answer. $\endgroup$ Commented Mar 20, 2024 at 15:56
  • $\begingroup$ The coordinate transformation used also changes the C and D matrices. Could you for completeness sake add those to your answer as well? $\endgroup$ Commented Mar 12 at 10:36
  • $\begingroup$ Echoing @fibonatic comment that C and D matrices are important to be defined in this method for an accurate match to C2D results. Assuming Cd=C and Dd=D does not yield the expected results $\endgroup$ Commented May 13 at 18:20
  • $\begingroup$ In response to comments by fibonatic and @AIS, I have added the matrices $C_\text{d}$ and $D_\text{d}$ for output. Sorry I wasn't able to get to this sooner. $\endgroup$ Commented May 14 at 14:40

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.