|
7017
|
1 ## Copyright (C) 1996, 2000, 2002, 2004, 2005, 2006, 2007 |
|
|
2 ## Auburn University. All rights reserved. |
|
3432
|
3 ## |
|
|
4 ## This file is part of Octave. |
|
|
5 ## |
|
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
7016
|
7 ## under the terms of the GNU General Public License as published by |
|
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
9 ## your option) any later version. |
|
3432
|
10 ## |
|
7016
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
14 ## General Public License for more details. |
|
3432
|
15 ## |
|
|
16 ## You should have received a copy of the GNU General Public License |
|
7016
|
17 ## along with Octave; see the file COPYING. If not, see |
|
|
18 ## <http://www.gnu.org/licenses/>. |
|
3432
|
19 |
|
3451
|
20 ## -*- texinfo -*- |
|
5016
|
21 ## @deftypefn {Function File} {[@var{zer}, @var{gain}] =} tzero (@var{a}, @var{b}, @var{c}, @var{d}, @var{opt}) |
|
|
22 ## @deftypefnx {Function File} {[@var{zer}, @var{gain}] =} tzero (@var{sys}, @var{opt}) |
|
|
23 ## Compute transmission zeros of a continuous system: |
|
|
24 ## @iftex |
|
|
25 ## @tex |
|
|
26 ## $$ \dot x = Ax + Bu $$ |
|
|
27 ## $$ y = Cx + Du $$ |
|
|
28 ## @end tex |
|
|
29 ## @end iftex |
|
|
30 ## @ifinfo |
|
3432
|
31 ## @example |
|
|
32 ## . |
|
|
33 ## x = Ax + Bu |
|
|
34 ## y = Cx + Du |
|
|
35 ## @end example |
|
5016
|
36 ## @end ifinfo |
|
|
37 ## or of a discrete one: |
|
|
38 ## @iftex |
|
|
39 ## @tex |
|
|
40 ## $$ x_{k+1} = Ax_k + Bu_k $$ |
|
|
41 ## $$ y_k = Cx_k + Du_k $$ |
|
|
42 ## @end tex |
|
|
43 ## @end iftex |
|
|
44 ## @ifinfo |
|
3432
|
45 ## @example |
|
|
46 ## x(k+1) = A x(k) + B u(k) |
|
|
47 ## y(k) = C x(k) + D u(k) |
|
|
48 ## @end example |
|
5016
|
49 ## @end ifinfo |
|
|
50 ## |
|
3432
|
51 ## @strong{Outputs} |
|
|
52 ## @table @var |
|
|
53 ## @item zer |
|
|
54 ## transmission zeros of the system |
|
|
55 ## @item gain |
|
5016
|
56 ## leading coefficient (pole-zero form) of @acronym{SISO} transfer function |
|
3432
|
57 ## returns gain=0 if system is multivariable |
|
|
58 ## @end table |
|
|
59 ## @strong{References} |
|
|
60 ## @enumerate |
|
|
61 ## @item Emami-Naeini and Van Dooren, Automatica, 1982. |
|
5016
|
62 ## @item Hodel, @cite{Computation of Zeros with Balancing}, 1992 Lin. Alg. Appl. |
|
3432
|
63 ## @end enumerate |
|
|
64 ## @end deftypefn |
|
|
65 |
|
|
66 ## Author: R. Bruce Tenison <btenison@eng.auburn.edu> |
|
|
67 ## Created: July 4, 1994 |
|
|
68 ## A. S. Hodel Aug 1995: allow for MIMO and system data structures |
|
|
69 |
|
|
70 function [zer, gain] = tzero (A, B, C, D) |
|
|
71 |
|
|
72 ## get A,B,C,D and Asys variables, regardless of initial form |
|
7126
|
73 if (nargin == 4) |
|
|
74 Asys = ss (A, B, C, D); |
|
|
75 elseif (nargin == 1 && ! isstruct (A)) |
|
|
76 error ("tzero: expecting argument to be system structure"); |
|
|
77 elseif (nargin != 1) |
|
6046
|
78 print_usage (); |
|
3432
|
79 else |
|
|
80 Asys = A; |
|
7126
|
81 [A, B, C, D] = sys2ss (Asys); |
|
3432
|
82 endif |
|
|
83 |
|
|
84 Ao = Asys; # save for leading coefficient |
|
7126
|
85 siso = is_siso (Asys); |
|
|
86 digital = is_digital (Asys); # check if it's mixed or not |
|
3432
|
87 |
|
|
88 ## see if it's a gain block |
|
7126
|
89 if (isempty (A)) |
|
3432
|
90 zer = []; |
|
|
91 gain = D; |
|
|
92 return; |
|
|
93 endif |
|
|
94 |
|
|
95 ## First, balance the system via the zero computation generalized eigenvalue |
|
|
96 ## problem balancing method (Hodel and Tiller, Linear Alg. Appl., 1992) |
|
|
97 |
|
7126
|
98 ## balance coefficients |
|
|
99 Asys = __zgpbal__ (Asys); |
|
|
100 [A, B, C, D] = sys2ss (Asys); |
|
3432
|
101 meps = 2*eps*norm ([A, B; C, D], "fro"); |
|
7126
|
102 ## ENVD algorithm |
|
|
103 Asys = zgreduce (Asys, meps); |
|
|
104 [A, B, C, D] = sys2ss (Asys); |
|
|
105 if (! isempty (A)) |
|
3432
|
106 ## repeat with dual system |
|
7126
|
107 Asys = ss (A', C', B', D'); |
|
|
108 Asys = zgreduce (Asys, meps); |
|
3432
|
109 |
|
|
110 ## transform back |
|
7126
|
111 [A, B, C, D] = sys2ss (Asys); |
|
|
112 Asys = ss (A', C', B', D'); |
|
3432
|
113 endif |
|
|
114 |
|
|
115 zer = []; # assume none |
|
7126
|
116 [A, B, C, D] = sys2ss (Asys); |
|
|
117 if (! isempty (C)) |
|
|
118 [W, r, Pi] = qr ([C, D]'); |
|
|
119 [nonz, ztmp] = zgrownorm (r, meps); |
|
|
120 if (nonz) |
|
3432
|
121 ## We can now solve the generalized eigenvalue problem. |
|
7126
|
122 [pp, mm] = size (D); |
|
|
123 nn = rows (A); |
|
3432
|
124 Afm = [A , B ; C, D] * W'; |
|
|
125 Bfm = [eye(nn), zeros(nn,mm); zeros(pp,nn+mm)]*W'; |
|
|
126 |
|
|
127 jdx = (mm+1):(mm+nn); |
|
|
128 Af = Afm(1:nn,jdx); |
|
|
129 Bf = Bfm(1:nn,jdx); |
|
7126
|
130 zer = qz (Af, Bf); |
|
3432
|
131 endif |
|
|
132 endif |
|
|
133 |
|
7126
|
134 mz = length (zer); |
|
|
135 [A, B, C, D] = sys2ss (Ao); # recover original system |
|
3432
|
136 ## compute leading coefficient |
|
7126
|
137 if (nargout == 2 && siso) |
|
|
138 n = rows (A); |
|
|
139 if (mz == n) |
|
3432
|
140 gain = D; |
|
7126
|
141 elseif (mz < n) |
|
3432
|
142 gain = C*(A^(n-1-mz))*B; |
|
|
143 endif |
|
|
144 else |
|
|
145 gain = []; |
|
|
146 endif |
|
|
147 endfunction |
|
|
148 |
