|
3381
|
1 ## Copyright (C) 1996 Auburn University. All Rights Reserved |
|
|
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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. |
|
3213
|
18 |
|
3381
|
19 ## -*- texinfo -*- |
|
|
20 ## @deftypefn {Function File} {} tzero (@var{a}, @var{b}, @var{c}, @var{d}@{, @var{opt}@}) |
|
|
21 ## @deftypefnx {Function File} {} tzero (@var{sys}@{,@var{opt}@}) |
|
|
22 ## Compute transmission zeros of a continuous |
|
|
23 ## @example |
|
|
24 ## . |
|
|
25 ## x = Ax + Bu |
|
|
26 ## y = Cx + Du |
|
|
27 ## @end example |
|
|
28 ## or discrete |
|
|
29 ## @example |
|
|
30 ## x(k+1) = A x(k) + B u(k) |
|
|
31 ## y(k) = C x(k) + D u(k) |
|
|
32 ## @end example |
|
|
33 ## system. |
|
|
34 ## @strong{Outputs} |
|
|
35 ## @table @var |
|
|
36 ## @item zer |
|
|
37 ## transmission zeros of the system |
|
|
38 ## @item gain |
|
|
39 ## leading coefficient (pole-zero form) of SISO transfer function |
|
|
40 ## returns gain=0 if system is multivariable |
|
|
41 ## @end table |
|
|
42 ## @strong{References} |
|
|
43 ## @enumerate |
|
|
44 ## @item Emami-Naeini and Van Dooren, Automatica, 1982. |
|
|
45 ## @item Hodel, "Computation of Zeros with Balancing," 1992 Lin. Alg. Appl. |
|
|
46 ## @end enumerate |
|
|
47 ## @end deftypefn |
|
3346
|
48 |
|
|
49 |
|
3385
|
50 function [zer, gain] = tzero (A, B, C, D) |
|
245
|
51 |
|
3381
|
52 ## R. Bruce Tenison July 4, 1994 |
|
|
53 ## A. S. Hodel Aug 1995: allow for MIMO and system data structures |
|
|
54 |
|
|
55 ## get A,B,C,D and Asys variables, regardless of initial form |
|
3213
|
56 if(nargin == 4) |
|
|
57 Asys = ss2sys(A,B,C,D); |
|
|
58 elseif( (nargin == 1) && (! is_struct(A))) |
|
|
59 usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)"); |
|
|
60 elseif(nargin != 1) |
|
|
61 usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)"); |
|
|
62 else |
|
|
63 Asys = A; |
|
|
64 [A,B,C,D] = sys2ss(Asys); |
|
60
|
65 endif |
|
|
66 |
|
3213
|
67 Ao = Asys; # save for leading coefficient |
|
|
68 siso = is_siso(Asys); |
|
|
69 digital = is_digital(Asys); # check if it's mixed or not |
|
|
70 |
|
3381
|
71 ## see if it's a gain block |
|
3213
|
72 if(isempty(A)) |
|
|
73 zer = []; |
|
|
74 gain = D; |
|
|
75 return; |
|
|
76 endif |
|
60
|
77 |
|
3381
|
78 ## First, balance the system via the zero computation generalized eigenvalue |
|
|
79 ## problem balancing method (Hodel and Tiller, Linear Alg. Appl., 1992) |
|
3213
|
80 |
|
|
81 Asys = zgpbal(Asys); [A,B,C,D] = sys2ss(Asys); # balance coefficients |
|
3238
|
82 meps = 2*eps*norm([A, B; C, D],'fro'); |
|
3213
|
83 Asys = zgreduce(Asys,meps); [A, B, C, D] = sys2ss(Asys); # ENVD algorithm |
|
|
84 if(!isempty(A)) |
|
3381
|
85 ## repeat with dual system |
|
3213
|
86 Asys = ss2sys(A', C', B', D'); Asys = zgreduce(Asys,meps); |
|
|
87 |
|
3381
|
88 ## transform back |
|
3213
|
89 [A,B,C,D] = sys2ss(Asys); Asys = ss2sys(A', C', B', D'); |
|
60
|
90 endif |
|
|
91 |
|
3213
|
92 zer = []; # assume none |
|
|
93 [A,B,C,D] = sys2ss(Asys); |
|
|
94 if( !isempty(C) ) |
|
3238
|
95 [W,r,Pi] = qr([C, D]'); |
|
3213
|
96 [nonz,ztmp] = zgrownorm(r,meps); |
|
|
97 if(nonz) |
|
3381
|
98 ## We can now solve the generalized eigenvalue problem. |
|
3213
|
99 [pp,mm] = size(D); |
|
|
100 nn = rows(A); |
|
3238
|
101 Afm = [A , B ; C, D] * W'; |
|
3213
|
102 Bfm = [eye(nn), zeros(nn,mm); zeros(pp,nn+mm)]*W'; |
|
|
103 |
|
|
104 jdx = (mm+1):(mm+nn); |
|
|
105 Af = Afm(1:nn,jdx); |
|
|
106 Bf = Bfm(1:nn,jdx); |
|
|
107 zer = qz(Af,Bf); |
|
|
108 endif |
|
|
109 endif |
|
|
110 |
|
|
111 mz = length(zer); |
|
|
112 [A,B,C,D] = sys2ss(Ao); # recover original system |
|
3381
|
113 ## compute leading coefficient |
|
3213
|
114 if ( (nargout == 2) && siso) |
|
|
115 n = rows(A); |
|
|
116 if ( mz == n) |
|
|
117 gain = D; |
|
|
118 elseif ( mz < n ) |
|
|
119 gain = C*(A^(n-1-mz))*B; |
|
|
120 endif |
|
|
121 else |
|
|
122 gain = []; |
|
|
123 endif |
|
60
|
124 endfunction |
|
3213
|
125 |
