1

I have a state space system with matrices A,B,C and D.

I can either create a state space system, sys1 = ss(A,B,C,D), of it or compute the transfer function matrix, sys2 = C*inv(z*I - A)*B + D

However when I draw the bode plot of both systems, they are different while they should be the same.

What is going wrong here? Does anyone have a clue? I know btw that the bodeplot generated by sys1 is correct.

The system can be downloaded here: https://dl.dropboxusercontent.com/u/20782274/system.mat

clear all;
close all;
clc;

Ts = 0.01;
z = tf('z',Ts);

% Discrete system
A = [0 1 0; 0 0 1; 0.41 -1.21 1.8];
B = [0; 0; 0.01];
C = [7 -73 170];
D = 1;

% Set as state space
sys1 = ss(A,B,C,D,Ts);

% Compute transfer function
sys2 = C*inv(z*eye(3) - A)*B + D;

% Compute the actual transfer function
[num,den] = ss2tf(A,B,C,D);
sys3 = tf(num,den,Ts);

% Show bode
bode(sys1,'b',sys2,'r--',sys3,'g--');

Edit: I made a small mistake, the transfer function matrix is sys2 = C*inv(z*I - A)*B + D, instead of sys2 = C*inv(z*I - A)*B - D which I did wrote done before. The problem still holds.

Edit 2: I have noticted that when I compute the denominator, it is correct.

syms z;
collect(det(z*eye(3) - A),z)
WG-
  • 1,046
  • 2
  • 14
  • 27

4 Answers4

1

Your assumption that sys2 = C*inv(z*I- A)*B + D is incorrect. The correct equivalent to your state-space system (A,B,C,D) is sys2 = C*inv(s*I- A)*B + D. If you want to express it in terms of z, you'll need to invert the relationship z = exp(s*T). sys1 is the correct representation of your state-space system. What I would suggest for sys2 is to do as follows:

sys1 = ss(mjlsCE.A,mjlsCE.B,mjlsCE.C,mjlsCE.D,Ts);
sys1_c = d2c(sys1);
s = tf('s');
sys2_c = sys1_c.C*inv(s*eye(length(sys1_c.A)) - sys1_c.A)*sys1_c.B + sys1_c.D;
sys2_d = c2d(sys2_c,Ts);

That should give you the correct result.

am304
  • 13,758
  • 2
  • 22
  • 40
  • I believe that is not correct. See: http://www.dcsc.tudelft.nl/~sc4060/transp/discreteNOTES.pdf page 47. Computing the transfer function matrix for a continuous- or discrete-time systems is the same only `s` is interchanged with `z`, i.e., `cont_sys = C*inv(s*I - A)*B + D`, `disc_sys = C*inv(z*I - A)*B + D`. Note also the small edit I made to my openingspost. Further also note the A,B,C,D matrices are already discrete! – WG- Oct 28 '14 at 12:11
  • Actually, my mistake `z` and `s` are related by `z = exp(s*T)`. Your assumption remains incorrect, but I will update my answer. – am304 Oct 28 '14 at 12:16
  • the bode plots are still not the same, please check it for yourself using my script and your corrections to convince yourself. I believe this is because I used 'zero-matched pole' (http://www.mathworks.nl/help/control/ref/c2d.html) discritization. Furthermore, I am still not convinced really :p. Indeed the relation between z = exp(s*T), however I am not dealing with any continuous-time system here, I am dealing with only a discrete system, e.g., the system matrices A,B,C,D are already discrete. None of the matrices are changing after using `sys1 = ss(A,B,C,D,Ts)` – WG- Oct 28 '14 at 12:36
  • I only have access to Octave, and it doesn't allow multiplying a transfer function and a matrix `s*eye(length(mjlsCE.A))` or `z*eye(length(mjlsCE.A))`, so I can' try it unfortunately. However, if you are using a discrete system, then stick with the `sys1` approach, it is the correct way to go about it. BTW, a matrix is not "discrete", it's just a collection of numbers. It can be used to describe a discrete or a continuous system, or both. – am304 Oct 28 '14 at 12:40
  • That is a shame. And indeed you are right a matrix is not discrete, but what I meant what that the matrices for the continuous-time system have already been converted to the matrices for the discrete-time system, i.e., (A,B,C,D) -> (Ad,Bd,Cd,Dd). And unfortunatly sticking with sys1 does not work for me, since I need to compute also different transfer function matrix with a different method but involves the same techniques as in `C*inv(z*I - A)*B + D` – WG- Oct 28 '14 at 12:46
  • I see what you mean about the matrices being defined for the discrete system only. I will edit my answer. – am304 Oct 28 '14 at 13:08
  • `Error using DynamicSystem/d2c (line 101) The "zoh" and "foh" methods cannot be used for discrete models with poles near z=0.` ;-) further even it would work it is still a work-around it is not an answer to my actual problem of why my transfer function matrix is different from my state space system, while in fact they should be the same. Thank you very much for your efforts in helping me :-) – WG- Oct 28 '14 at 13:19
  • I simplified the problem. Since I noticted the same problem occures when using the example from dcsc.tudelft.nl/~sc4060/transp/discreteNOTES.pdf on page 47 – WG- Oct 28 '14 at 13:50
1

Due to inacurracy of the inverse function extra unobservable poles and zeros are added to the system. For this reason you need to compute the minimal realization of your transfer function matrix.

Meaning 

% Compute transfer function
sys2 = minreal(C*inv(z*eye(3) - A)*B + D);
WG-
  • 1,046
  • 2
  • 14
  • 27
1

What you are noticing is actually a numerical instability regarding pole-zero pair cancellations. If you run the following code:

A = [0, 1, 0; 0, 0, 1; 0.41, -1.21, 1.8] ;
B = [0; 0; 0.01] ;
C = [7, -73, 170] ;
D = 1 ;

sys_ss = ss(A, B, C, D) ;

sys_tf_simp = tf(sys_ss) ;

s = tf('s') ;
sys_tf_full = tf(C*inv(s*eye(3) - A)*B + D) ;

zero(sys_tf_simp)
zero(sys_tf_full)

pole(sys_tf_simp)
pole(sys_tf_full)

you will see that the transfer function formulated by matrices directly has a lot more poles and zeros than the one formulated by MatLab's tf function. You will also notice that every single pair of these "extra" poles and zeros are equal- meaning that they cancel with each other if you were to simply the rational expression. MatLab's tf presents the simplified form, with equal pole-zero pairs cancelled out. This is algebraically equivalent to the unsimplified form, but not numerically.

When you call bode on the unsimplified transfer function, MatLab begins its numerical plotting routine with the pole-zero pairs not cancelled algebraically. If the computer was perfect, the result would be the same as in the simplified case. However, numerical error when evaluating the numerator and denominators effectively leaves some of the pole-zero pairs "uncancelled" and as many of these poles are in the far right side of the s plane, they drastically influence the output behavior.

Check out this link for info on this same problem but from the perspective of design: http://ctms.engin.umich.edu/CTMS/index.php?aux=Extras_PZ

In your original code, you can think of the output drawn in green as what the naive designer wanted to see when he cancelled all his unstable poles with zeros, but the output drawn in red is what he actually got because in practice, finite-precision and real-world tolerances prevent the poles and zeros from cancelling perfectly.

jnez71
  • 138
  • 6
0

Why is an unobservable / uncontrollable pole? I think this issue comes only because the inverse of a transfer function matrix is inaccurate in Matlab.

Note:

  1. A is 3x3 and the minimal realization has also order 3.
  2. What you did is the inverse of a transfer function matrix, not a symbolic or numeric matrix.
# Discrete system
Ts = 0.01;
A = [0 1 0; 0 0 1; 0.41 -1.21 1.8];
B = [0; 0; 0.01];
C = [7 -73 170];
D = 1;

z = tf('z', Ts)) # z is a discrete tf
A1 = z*eye(3) - A # a tf matrix with a direct feedthrough matrix A

# inverse it, multiply with C and B from left and right, and plus D
G = D + C*inv(A1)*B

G is now a scalar (SISO) transfer function.

Without "minreal", G has order 9 (funny, I don't know how Matlab computes it, perhaps the "Adj(.)/det(.)" method). Matlab cannot cancel the common factors in the numerator and the denominator, because z is of class 'tf' rather than a symbolic variable.

Do you agree or do I have misunderstanding?

gundamlh
  • 183
  • 1
  • 9