Hi.
I suggest the attouched patch against current svn of the "control" package. Please review it. -Added a new function named "balreal". At the moment it is tested only with continuous time systems (not with discrete-time ones), and its test should be improved, but I think it is good it is committed. I modified the INDEX file to reflect this change. -Added 2 new functions named "isct" and "isdt" and consider "is_digital" obsolete. At the moment they are only wrapper for is_digital. I tested they are compatible with Matlab documentation for simple systems (look at their tests). In the long term they should be modified to be more compatible with Matlab counterparts. I modified the INDEX file to reflect this change. -Rewritten the "gram" function to be more compatible with Matlab: added a new signature and now accept discrete-time systems. Added tests to consider all paths in the function. -Modified the "dgram" function to reflect the fact that it is left only for backward compatibility (not efficient but simple to mantain). In fact it was even removed from Matlab documentation (but it is present in recent version of Matlab, I think for compatibility reason). -Modified the INDEX file to suggest using "is_stable" instead of the missing "isstable" -Bumped package version from 1.0.8 to 1.0.9 Cheers, Luca Favatella I attouch svn diff.
Index: DESCRIPTION =================================================================== --- DESCRIPTION (revisione 5499) +++ DESCRIPTION (copia locale) @@ -1,5 +1,5 @@ Name: Control -Version: 1.0.8 +Version: 1.0.9 Date: 2008-08-23 Author: Ben Sapp Maintainer: None Index: INDEX =================================================================== --- INDEX (revisione 5499) +++ INDEX (copia locale) @@ -40,12 +40,14 @@ are dare dgram dlyap gram lyap pinv dre zgfmul zgfslv zginit zgreduce zgrownorm zgscal zgsgiv zgshsr + balreal System properties analdemo abcddim ctrb h2norm hinfnorm obsv pzmap is_abcd is_controllable is_detectable is_dgkf - is_digital is_observable is_sample is_siso is_stabilizable + is_observable is_sample is_siso is_stabilizable is_signal_list is_stable + isct isdt Time domain analysis c2d d2c dmr2d damp dcgain impulse step Frequency domain analysis @@ -65,6 +67,7 @@ Obsolete dezero dlqg minfo packsys unpacksys swaprows swapcols rotg qzval syschnames series + is_digital Internal __bodquist__ __freqresp__ __stepimp__ __abcddims__ __syschnamesl__ __syscont_disc__ __sysdefioname__ __sysdefstname__ __sysgroupn__ @@ -78,8 +81,6 @@ ssdata= use <f>sys2ss</f> tfdata=use <f>sys2tf</f> zpkdata=use <f>sys2zp</f> - isct=use !<f>is_digital</f> - isdt=use <f>is_digital</f> issiso=use <f>is_siso</f> append(sys)=use <f>sysappend</f> lft=use <f>starp</f> @@ -87,4 +88,4 @@ zero=use <f>tzero</f> or <f>tzero2</f> kalman=use <f>lqe</f> kalmd=use <f>dlqe</f> or <f>dkalman</f> - + isstable=use <f>is_stable</f> Index: inst/gram.m =================================================================== --- inst/gram.m (revisione 5499) +++ inst/gram.m (copia locale) @@ -17,30 +17,113 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {} gram (@var{a}, @var{b}) -## Return controllability gramian @var{m} of the continuous time system -## @math{dx/dt = a x + b u}. +## @deftypefn {Function File} {...@var{w} =} gram (@var{sys}, @var{mode}) +## @deftypefnx {Function File} {...@var{wc} =} gram (@var{a}, @var{b}) +## @code{gram (@var{sys}, 'c')} returns the controllability gramian of +## the (continuous- or discrete-time) system @var{sys}. +## @code{gram (@var{sys}, 'o')} returns the observability gramian of the +## (continuous- or discrete-time) system @var{sys}. +## @code{gram (@var{a}, @var{b})} returns the controllability gramian +## @var{Wc} of the continuous-time system @math{dx/dt = a x + b u}; +## i.e., @var{Wc} satisfies @math{a Wc + m Wc' + b b' = 0}. ## -## @var{m} satisfies @math{a m + m a' + b b' = 0}. ## @end deftypefn ## Author: A. S. Hodel <a.s.ho...@eng.auburn.edu> -function m = gram (a, b) + # TODO: substitute is_stable with isstable +function W = gram (argin1, argin2) + if (nargin != 2) print_usage (); - endif + else - ## Let lyap do the error checking... + if (! ischar (argin2)) ## the function was called as "gram (a, b)" + a = argin1; + b = argin2; + if (! is_stable (a)) + error ("the continuous-time system must have a stable 'a' matrix"); + else - m = lyap (a', b*b'); + ## let lyap do the error checking about dimensions + W = lyap (a', b*b'); + endif + else ## the function was called as "gram (sys, mode)" + if (! (strcmp (argin2, 'c') || strcmp (argin2, 'o'))) + print_usage (); + else + if (strcmp (argin2, 'c')) + a = argin1.a; + b = argin1.b; + elseif (strcmp (argin2, 'o')) + a = argin1.a'; + b = argin1.c'; + endif + + if (isct (argin1)) + if (! is_stable (argin1)) + error ("the continuous-time system must be stable"); + else + + ## let lyap do the error checking about dimensions + W = lyap (a', b*b'); + endif + elseif (isdt (argin1)) + if (! is_stable (argin1)) + error ("the discrete-time system must be stable"); + else + + ## let dlyap do the error checking about dimensions + W = dlyap (a, b*b'); + endif + else + error ("strange behaviour in isct/isdt: if you can reproduce \ + this, please submit a bug report"); + endif + endif + endif + + endif + endfunction %!test %! a = [-1 0 0; 1/2 -1 0; 1/2 0 -1]; %! b = [1 0; 0 -1; 0 1]; -%! m = gram (a, b); -%! assert (a * m + m * a' + b *b', zeros (size (a))) \ No newline at end of file +%! c = [0 0 1; 1 1 0]; ## it doesn't matter what the value of c is +%! Wc = gram (ss (a, b, c), 'c'); +%! assert (a * Wc + Wc * a' + b * b', zeros (size (a))) + +%!test +%! a = [-1 0 0; 1/2 -1 0; 1/2 0 -1]; +%! b = [1 0; 0 -1; 0 1]; ## it doesn't matter what the value of b is +%! c = [0 0 1; 1 1 0]; +%! Wo = gram (ss (a, b, c), 'o'); +%! assert (a' * Wo + Wo * a + c' * c, zeros (size (a))) + +%!test +%! a = [-1 0 0; 1/2 -1 0; 1/2 0 -1]; +%! b = [1 0; 0 -1; 0 1]; +%! Wc = gram (a, b); +%! assert (a * Wc + Wc * a' + b * b', zeros (size (a))) + +%!test +%! a = [-1 0 0; 1/2 1 0; 1/2 0 -1] / 2; +%! b = [1 0; 0 -1; 0 1]; +%! c = [0 0 1; 1 1 0]; ## it doesn't matter what the value of c is +%! d = zeros (rows (c), columns (b)); ## it doesn't matter what the value of d is +%! Ts = 0.1; ## Ts != 0 +%! Wc = gram (ss (a, b, c, d, Ts), 'c'); +%! assert (a * Wc * a' - Wc + b * b', zeros (size (a)), 1e-12) + +%!test +%! a = [-1 0 0; 1/2 1 0; 1/2 0 -1] / 2; +%! b = [1 0; 0 -1; 0 1]; ## it doesn't matter what the value of b is +%! c = [0 0 1; 1 1 0]; +%! d = zeros (rows (c), columns (b)); ## it doesn't matter what the value of d is +%! Ts = 0.1; ## Ts != 0 +%! Wo = gram (ss (a, b, c, d, Ts), 'o'); +%! assert (a' * Wo * a - Wo + c' * c, zeros (size (a)), 1e-12) \ No newline at end of file Index: inst/balreal.m =================================================================== --- inst/balreal.m (revisione 0) +++ inst/balreal.m (revisione 0) @@ -0,0 +1,117 @@ +## Copyright (C) 2008 Luca Favatella <slacky...@gmail.com> +## +## +## This program is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This program is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {...@var{sysb}, @var{g}] =} balreal (@var{sys}) +## @deftypefnx {Function File} {...@var{sysb}, @var{g}, @var{T}, @var{Ti}] =} balreal (@var{sys}) +## Balanced realization of the continuous-time LTI system @var{sys}. +## +## @strong{Input} +## @table @var +## @item sys +## Stable, controllable and observable continuous-time LTI system. +## @end table +## +## @strong{Outputs} +## @table @var +## @item sysb +## Balanced realization of @var{sys}. +## @item g +## Diagonal of the balanced gramians. +## @item T +## State transformation to convert @var{sys} to @var{sysb}. +## @item Ti +## Inverse of T. +## @end table +## +## @seealso{gram} +## @end deftypefn + +## Author: Luca Favatella <slacky...@gmail.com> +## Version: 0.2.4 + + # TODO + # improve continuous-time system test + # test against discrete-time systems + # clean asserts + # clean names Wc2/Wc and Wo2/Wc + # substitute is_stable with isstable + +function [sysb, g, T, Ti] = balreal (sys) + + if (nargin != 1) + print_usage (); + elseif (! is_stable (sys)) + error ("sys must be stable"); + elseif (! is_controllable (sys)) + error ("sys must be controllable"); + elseif (! is_observable (sys)) + error ("sys must be observable"); + else + ASSERT_TOL = 1e-12; ## DEBUG + + ## step 1: compute T1 + Wc2 = gram (sys, 'c'); + [Vc, Sc2, trash] = svd (Wc2); + assert (trash, Vc, ASSERT_TOL); ## DEBUG + T1 = Vc * sqrt (Sc2); + + ## step 2: apply T1 + Atilde = inv (T1) * sys.a * T1; + Btilde = inv (T1) * sys.b; + Ctilde = sys.c * T1; + Wc2tilde = gram (Atilde, Btilde); ## DEBUG + assert (Wc2tilde, eye (size (Wc2tilde)), ASSERT_TOL); ## DEBUG + + ## step 3: compute T2 + Wo2tilde = gram (Atilde', Ctilde'); + [Vo, Sigma2, trash] = svd (Wo2tilde); + assert (trash, Vo, ASSERT_TOL); ## DEBUG + T2 = Vo * (Sigma2 ^ (- 1/4)); + + ## step 4: apply T2 + Abal = inv (T2) * Atilde * T2; + Bbal = inv (T2) * Btilde; + Cbal = Ctilde * T2; + + + ## prepare return values + ## sysb + sysb = ss (Abal, Bbal, Cbal, sys.d); + + ## g + Wc2bal = gram (Abal, Bbal); + g = diag (Wc2bal); + Wo2bal = gram (Abal', Cbal'); ## DEBUG + assert (Wc2bal, Wo2bal, ASSERT_TOL); ## DEBUG + Wo2 = gram (sys, 'o'); ## DEBUG + Sigma = diag (sqrt (sort (eig (Wc2 * Wo2), 'descend'))); ## DEBUG + assert (Wc2bal, Sigma, ASSERT_TOL); ## DEBUG + assert (Wo2bal, Sigma, ASSERT_TOL); ## DEBUG + + ## T and Ti + T = inv (T1 * T2); + Ti = inv (T); + endif +endfunction + + +%!test +%! a = [-1 0 0; 1/2 -1 0; 1/2 0 -1]; +%! b = [1 0; 0 -1; 0 1]; +%! c = [0 0 1; 1 1 0]; +%! balreal (ss (a, b, c)); \ No newline at end of file Index: inst/isct.m =================================================================== --- inst/isct.m (revisione 0) +++ inst/isct.m (revisione 0) @@ -0,0 +1,49 @@ +## Copyright (C) 2008 Luca Favatella <slacky...@gmail.com> +## +## +## This program is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This program is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} isct (@var{sys}) +## Return true if the LTI system @var{sys} is continuous-time, false otherwise. +## +## @seealso{isdt} +## @end deftypefn + +## Author: Luca Favatella <slacky...@gmail.com> +## Version: 0.2.1 + +function retval = isct (sys) + if (nargin != 1) + print_usage (); + else + retval = (! is_digital (sys)); + endif +endfunction + + +%!test +%! A = [-1 0 0; 1/2 -1 0; 1/2 0 -1]; +%! B = [1 0; 0 -1; 0 1]; +%! C = [0 0 1; 1 1 0]; +%! assert (isct (ss (A, B, C))); + +%!test +%! A = [-1 0 0; 1/2 -1 0; 1/2 0 -1]; +%! B = [1 0; 0 -1; 0 1]; +%! C = [0 0 1; 1 1 0]; +%! D = zeros (rows (C), columns (B)); +%! Ts = 0.1; +%! assert (! isct (ss (A, B, C, D, Ts))); \ No newline at end of file Index: inst/dgram.m =================================================================== --- inst/dgram.m (revisione 5499) +++ inst/dgram.m (copia locale) @@ -62,10 +62,24 @@ if (nargin != 2) print_usage (); + else + + ## this is not efficient, but it is simple to mantain + ## + ## efficiency is not important because this function is here only + ## for backward compatibility + c = zeros (1, length (a)); ## it doesn't matter what the value of c is + d = zeros (rows (c), columns (b)); ## it doesn't matter what the value of d is + Ts = 0.1; ## Ts != 0 + sys = ss (a, b, c, d, Ts); + m = gram (sys, 'c'); endif - ## let dlyap do the error checking... +endfunction - m = dlyap (a, b*b'); -endfunction +%!test +%! a = [-1 0 0; 1/2 1 0; 1/2 0 -1] / 2; +%! b = [1 0; 0 -1; 0 1]; +%! Wc = dgram (a, b); +%! assert (a * Wc * a' - Wc + b * b', zeros (size (a)), 1e-12) \ No newline at end of file Index: inst/isdt.m =================================================================== --- inst/isdt.m (revisione 0) +++ inst/isdt.m (revisione 0) @@ -0,0 +1,49 @@ +## Copyright (C) 2008 Luca Favatella <slacky...@gmail.com> +## +## +## This program is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This program is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} isdt (@var{sys}) +## Return true if the LTI system @var{sys} is discrete-time, false otherwise. +## +## @seealso{isct} +## @end deftypefn + +## Author: Luca Favatella <slacky...@gmail.com> +## Version: 0.2.1 + +function retval = isdt (sys) + if (nargin != 1) + print_usage (); + else + retval = is_digital (sys); + endif +endfunction + + +%!test +%! A = [-1 0 0; 1/2 -1 0; 1/2 0 -1]; +%! B = [1 0; 0 -1; 0 1]; +%! C = [0 0 1; 1 1 0]; +%! D = zeros (rows (C), columns (B)); +%! Ts = 0.1; +%! assert (isdt (ss (A, B, C, D, Ts))); + +%!test +%! A = [-1 0 0; 1/2 -1 0; 1/2 0 -1]; +%! B = [1 0; 0 -1; 0 1]; +%! C = [0 0 1; 1 1 0]; +%! assert (! isdt (ss (A, B, C))); \ No newline at end of file
------------------------------------------------------------------------------
_______________________________________________ Octave-dev mailing list Octave-dev@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/octave-dev
