|
3213
|
1 # Copyright (C) 1996 A. Scottedward Hodel |
|
|
2 # |
|
|
3 # This file is part of Octave. |
|
|
4 # |
|
|
5 # Octave is free software; you can redistribute it and/or modify it |
|
|
6 # under the terms of the GNU General Public License as published by the |
|
|
7 # Free Software Foundation; either version 2, or (at your option) any |
|
|
8 # later version. |
|
|
9 # |
|
|
10 # Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
12 # FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
13 # for more details. |
|
|
14 # |
|
|
15 # You should have received a copy of the GNU General Public License |
|
|
16 # along with Octave; see the file COPYING. If not, write to the Free |
|
|
17 # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. |
|
|
18 |
|
|
19 function [zer, gain] = tzero(A,B,C,D) |
|
|
20 # [zer{,gain}] = tzero(A,B,C,D) -or- |
|
|
21 # [zer{,gain}] = tzero(Asys) |
|
|
22 # Compute transmission zeros of a continuous |
|
|
23 # . |
|
|
24 # x = Ax + Bu |
|
|
25 # y = Cx + Du |
|
|
26 # |
|
|
27 # or discrete |
|
|
28 # x(k+1) = A x(k) + B u(k) |
|
|
29 # y(k) = C x(k) + D u(k) |
|
|
30 # |
|
|
31 # system. |
|
|
32 # |
|
|
33 # outputs: |
|
|
34 # zer: transmission zeros of the system |
|
|
35 # gain: leading coefficient (pole-zero form) of SISO transfer function |
|
|
36 # returns gain=0 if system is multivariable |
|
|
37 # References: |
|
|
38 # Hodel, "Computation of Zeros with Balancing," 1992 Lin. Alg. Appl. |
|
|
39 |
|
|
40 # R. Bruce Tenison July 4, 1994 |
|
|
41 # A. S. Hodel Aug 1995: allow for MIMO and system data structures |
|
|
42 # $Revision: 1.16 $ |
|
|
43 # $Log: tzero.m,v $ |
|
|
44 # Revision 1.16 1998-11-06 16:15:37 jwe |
|
|
45 # *** empty log message *** |
|
|
46 # |
|
|
47 # Revision 1.7 1998/08/24 15:50:30 hodelas |
|
|
48 # updated documentation |
|
|
49 # |
|
|
50 # Revision 1.4 1998/08/12 20:34:36 hodelas |
|
|
51 # Updated to use system access calls instead of direct structure access |
|
|
52 # |
|
|
53 # Revision 1.3 1998/07/21 14:53:11 hodelas |
|
|
54 # use isempty instead of size tests; use sys calls to reduce direct |
|
|
55 # access to system structure elements |
|
|
56 # |
|
|
57 # Revision 1.2 1997/02/13 11:58:05 hodel |
|
|
58 # tracked down error in zgfslv; added Log message |
|
245
|
59 |
|
3213
|
60 # get A,B,C,D and Asys variables, regardless of initial form |
|
|
61 if(nargin == 4) |
|
|
62 Asys = ss2sys(A,B,C,D); |
|
|
63 elseif( (nargin == 1) && (! is_struct(A))) |
|
|
64 usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)"); |
|
|
65 elseif(nargin != 1) |
|
|
66 usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)"); |
|
|
67 else |
|
|
68 Asys = A; |
|
|
69 [A,B,C,D] = sys2ss(Asys); |
|
60
|
70 endif |
|
|
71 |
|
3213
|
72 Ao = Asys; # save for leading coefficient |
|
|
73 siso = is_siso(Asys); |
|
|
74 digital = is_digital(Asys); # check if it's mixed or not |
|
|
75 |
|
|
76 # see if it's a gain block |
|
|
77 if(isempty(A)) |
|
|
78 zer = []; |
|
|
79 gain = D; |
|
|
80 return; |
|
|
81 endif |
|
60
|
82 |
|
3213
|
83 # First, balance the system via the zero computation generalized eigenvalue |
|
|
84 # problem balancing method (Hodel and Tiller, Linear Alg. Appl., 1992) |
|
|
85 |
|
|
86 Asys = zgpbal(Asys); [A,B,C,D] = sys2ss(Asys); # balance coefficients |
|
|
87 meps = 2*eps*norm([A B; C D],'fro'); |
|
|
88 Asys = zgreduce(Asys,meps); [A, B, C, D] = sys2ss(Asys); # ENVD algorithm |
|
|
89 if(!isempty(A)) |
|
|
90 # repeat with dual system |
|
|
91 Asys = ss2sys(A', C', B', D'); Asys = zgreduce(Asys,meps); |
|
|
92 |
|
|
93 # transform back |
|
|
94 [A,B,C,D] = sys2ss(Asys); Asys = ss2sys(A', C', B', D'); |
|
60
|
95 endif |
|
|
96 |
|
3213
|
97 zer = []; # assume none |
|
|
98 [A,B,C,D] = sys2ss(Asys); |
|
|
99 if( !isempty(C) ) |
|
|
100 [W,r,Pi] = qr([C D]'); |
|
|
101 [nonz,ztmp] = zgrownorm(r,meps); |
|
|
102 if(nonz) |
|
|
103 # We can now solve the generalized eigenvalue problem. |
|
|
104 [pp,mm] = size(D); |
|
|
105 nn = rows(A); |
|
|
106 Afm = [A , B ; C D] * W'; |
|
|
107 Bfm = [eye(nn), zeros(nn,mm); zeros(pp,nn+mm)]*W'; |
|
|
108 |
|
|
109 jdx = (mm+1):(mm+nn); |
|
|
110 Af = Afm(1:nn,jdx); |
|
|
111 Bf = Bfm(1:nn,jdx); |
|
|
112 zer = qz(Af,Bf); |
|
|
113 endif |
|
|
114 endif |
|
|
115 |
|
|
116 mz = length(zer); |
|
|
117 [A,B,C,D] = sys2ss(Ao); # recover original system |
|
|
118 #compute leading coefficient |
|
|
119 if ( (nargout == 2) && siso) |
|
|
120 n = rows(A); |
|
|
121 if ( mz == n) |
|
|
122 gain = D; |
|
|
123 elseif ( mz < n ) |
|
|
124 gain = C*(A^(n-1-mz))*B; |
|
|
125 endif |
|
|
126 else |
|
|
127 gain = []; |
|
|
128 endif |
|
60
|
129 endfunction |
|
3213
|
130 |
