1 ## Copyright (C) 2009, 2010, 2011 Lukas F. Reichlin
3 ## This file is part of LTI Syncope.
5 ## LTI Syncope is free software: you can redistribute it and/or modify
6 ## it under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation, either version 3 of the License, or
8 ## (at your option) any later version.
10 ## LTI Syncope is distributed in the hope that it will be useful,
11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 ## GNU General Public License for more details.
15 ## You should have received a copy of the GNU General Public License
16 ## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
19 ## @deftypefn {Function File} {@var{sys} =} sminreal (@var{sys})
20 ## @deftypefnx {Function File} {@var{sys} =} sminreal (@var{sys}, @var{tol})
21 ## Perform state-space model reduction based on structure.
22 ## Remove states which have no influence on the input-output behaviour.
23 ## The physical meaning of the states is retained.
30 ## Optional tolerance for controllability and observability.
31 ## Entries of the state-space matrices whose moduli are less or equal to @var{tol}
32 ## are assumed to be zero. Default value is 0.
38 ## Reduced state-space model.
44 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
45 ## Created: October 2009
48 function sys = sminreal (sys, tol = 0)
50 if (nargin > 2) # sminreal () not possible (inside @lti)
54 if (! isa (sys, "ss"))
55 warning ("sminreal: system not in state-space form");
56 sys = ss (sys); # needed by __sys_prune__
59 if (! (is_real_scalar (tol) && tol >= 0))
60 error ("sminreal: second argument is not a valid tolerance");
63 [a, b, c, d, e] = dssdata (sys, []);
74 co_idx = __controllable_states__ (a, b);
75 ob_idx = __controllable_states__ (a.', c.');
77 st_idx = intersect (co_idx, ob_idx);
79 sys = __sys_prune__ (sys, ":", ":", st_idx);
84 function c_idx = __controllable_states__ (a, b)
86 n = rows (a); # number of states
87 a = a & ! eye (n); # set diagonal entries to zero
89 c_vec = any (b, 2); # states directly controllable
90 c_idx = find (c_vec); # indices of directly controllable states
91 c_idx_new = 0; # any vector of length > 0 possible
93 while (all (length (c_idx) != [0, n]) && length(c_idx_new) != 0)
95 u_idx = find (! c_vec); # indices of uncontrollable states
100 repmat (c_vec.', length (u_idx), 1)
101 a(u_idx, :) & repmat (c_vec.', length (u_idx), 1)
102 any (a(u_idx, :) & repmat (c_vec.', length (u_idx), 1), 2)
103 find (any (a(u_idx, :) & repmat (c_vec.', length (u_idx), 1), 2))
106 c_idx_new = u_idx (find (any (a(u_idx, :) & repmat (c_vec.', length (u_idx), 1), 2)));
107 c_idx = union (c_idx, c_idx_new);
108 c_vec(c_idx_new) = 1;
118 %! A = ss (-2, 3, 4, 5);
120 %! C = sminreal (B); # no states should be removed
129 %! A = ss (-1, 1, 1, 0);
130 %! B = ss (-2, 3, 4, 5);
132 %! D = sminreal (C(:, 1));
133 %! E = sminreal (C(:, 2));