From: <sla...@us...> - 2008-11-27 22:03:04
|
Revision: 5467 http://octave.svn.sourceforge.net/octave/?rev=5467&view=rev Author: slackydeb Date: 2008-11-27 22:03:01 +0000 (Thu, 27 Nov 2008) Log Message: ----------- Bug fix: gram function not working as documented. Add a regression test. Bump version number. Modified Paths: -------------- trunk/octave-forge/main/control/DESCRIPTION trunk/octave-forge/main/control/inst/gram.m Modified: trunk/octave-forge/main/control/DESCRIPTION =================================================================== --- trunk/octave-forge/main/control/DESCRIPTION 2008-11-26 07:55:33 UTC (rev 5466) +++ trunk/octave-forge/main/control/DESCRIPTION 2008-11-27 22:03:01 UTC (rev 5467) @@ -1,5 +1,5 @@ Name: Control -Version: 1.0.7 +Version: 1.0.8 Date: 2008-08-23 Author: Ben Sapp Maintainer: None Modified: trunk/octave-forge/main/control/inst/gram.m =================================================================== --- trunk/octave-forge/main/control/inst/gram.m 2008-11-26 07:55:33 UTC (rev 5466) +++ trunk/octave-forge/main/control/inst/gram.m 2008-11-27 22:03:01 UTC (rev 5467) @@ -34,6 +34,13 @@ ## Let lyap do the error checking... - m = lyap (a, b*b'); + m = lyap (a', b*b'); 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 This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <sla...@us...> - 2008-12-28 17:47:20
|
Revision: 5501 http://octave.svn.sourceforge.net/octave/?rev=5501&view=rev Author: slackydeb Date: 2008-12-28 17:47:14 +0000 (Sun, 28 Dec 2008) Log Message: ----------- Added new functions: "balreal", "isct", "isdt". Considered the function "is_digital" as obsolete. Rewritten the function "gram". Modified the function "dgram" to be simple to mantain. Suggested using "is_stable" instead of the missing "isstable". Modified the "INDEX" file to reflect changes. Bumped package version. Modified Paths: -------------- trunk/octave-forge/main/control/DESCRIPTION trunk/octave-forge/main/control/INDEX trunk/octave-forge/main/control/inst/dgram.m trunk/octave-forge/main/control/inst/gram.m Added Paths: ----------- trunk/octave-forge/main/control/inst/balreal.m trunk/octave-forge/main/control/inst/isct.m trunk/octave-forge/main/control/inst/isdt.m Modified: trunk/octave-forge/main/control/DESCRIPTION =================================================================== --- trunk/octave-forge/main/control/DESCRIPTION 2008-12-28 10:13:22 UTC (rev 5500) +++ trunk/octave-forge/main/control/DESCRIPTION 2008-12-28 17:47:14 UTC (rev 5501) @@ -1,5 +1,5 @@ Name: Control -Version: 1.0.8 +Version: 1.0.9 Date: 2008-08-23 Author: Ben Sapp Maintainer: None Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2008-12-28 10:13:22 UTC (rev 5500) +++ trunk/octave-forge/main/control/INDEX 2008-12-28 17:47:14 UTC (rev 5501) @@ -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> Added: trunk/octave-forge/main/control/inst/balreal.m =================================================================== --- trunk/octave-forge/main/control/inst/balreal.m (rev 0) +++ trunk/octave-forge/main/control/inst/balreal.m 2008-12-28 17:47:14 UTC (rev 5501) @@ -0,0 +1,117 @@ +## Copyright (C) 2008 Luca Favatella <sla...@gm...> +## +## +## 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 <sla...@gm...> +## 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 Modified: trunk/octave-forge/main/control/inst/dgram.m =================================================================== --- trunk/octave-forge/main/control/inst/dgram.m 2008-12-28 10:13:22 UTC (rev 5500) +++ trunk/octave-forge/main/control/inst/dgram.m 2008-12-28 17:47:14 UTC (rev 5501) @@ -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 Modified: trunk/octave-forge/main/control/inst/gram.m =================================================================== --- trunk/octave-forge/main/control/inst/gram.m 2008-12-28 10:13:22 UTC (rev 5500) +++ trunk/octave-forge/main/control/inst/gram.m 2008-12-28 17:47:14 UTC (rev 5501) @@ -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...@en...> -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 Added: trunk/octave-forge/main/control/inst/isct.m =================================================================== --- trunk/octave-forge/main/control/inst/isct.m (rev 0) +++ trunk/octave-forge/main/control/inst/isct.m 2008-12-28 17:47:14 UTC (rev 5501) @@ -0,0 +1,49 @@ +## Copyright (C) 2008 Luca Favatella <sla...@gm...> +## +## +## 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 <sla...@gm...> +## 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 Added: trunk/octave-forge/main/control/inst/isdt.m =================================================================== --- trunk/octave-forge/main/control/inst/isdt.m (rev 0) +++ trunk/octave-forge/main/control/inst/isdt.m 2008-12-28 17:47:14 UTC (rev 5501) @@ -0,0 +1,49 @@ +## Copyright (C) 2008 Luca Favatella <sla...@gm...> +## +## +## 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 <sla...@gm...> +## 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 This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <sla...@us...> - 2009-02-11 21:02:30
|
Revision: 5556 http://octave.svn.sourceforge.net/octave/?rev=5556&view=rev Author: slackydeb Date: 2009-02-11 21:02:23 +0000 (Wed, 11 Feb 2009) Log Message: ----------- Little optimization fix in balreal.m. Add myself as mantainer. Bump version number and last update date in DESCRIPTION. Modified Paths: -------------- trunk/octave-forge/main/control/DESCRIPTION trunk/octave-forge/main/control/inst/balreal.m Modified: trunk/octave-forge/main/control/DESCRIPTION =================================================================== --- trunk/octave-forge/main/control/DESCRIPTION 2009-02-11 02:51:44 UTC (rev 5555) +++ trunk/octave-forge/main/control/DESCRIPTION 2009-02-11 21:02:23 UTC (rev 5556) @@ -1,8 +1,8 @@ Name: Control -Version: 1.0.9 -Date: 2008-08-23 +Version: 1.0.10 +Date: 2009-02-11 Author: Ben Sapp -Maintainer: None +Maintainer: Luca Favatella <sla...@gm...> Title: Control. Description: Additional Octave Control tools Depends: octave (>= 2.9.7) Modified: trunk/octave-forge/main/control/inst/balreal.m =================================================================== --- trunk/octave-forge/main/control/inst/balreal.m 2009-02-11 02:51:44 UTC (rev 5555) +++ trunk/octave-forge/main/control/inst/balreal.m 2009-02-11 21:02:23 UTC (rev 5556) @@ -1,4 +1,4 @@ -## Copyright (C) 2008 Luca Favatella <sla...@gm...> +## Copyright (C) 2008, 2009 Luca Favatella <sla...@gm...> ## ## ## This program is free software; you can redistribute it and/or modify it @@ -42,7 +42,7 @@ ## @end deftypefn ## Author: Luca Favatella <sla...@gm...> -## Version: 0.2.4 +## Version: 0.2.5 # TODO # improve continuous-time system test @@ -104,8 +104,8 @@ assert (Wo2bal, Sigma, ASSERT_TOL); ## DEBUG ## T and Ti - T = inv (T1 * T2); - Ti = inv (T); + Ti = T1 * T2; + T = inv (Ti); endif endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <sla...@us...> - 2009-07-06 08:26:57
|
Revision: 5986 http://octave.svn.sourceforge.net/octave/?rev=5986&view=rev Author: slackydeb Date: 2009-07-06 08:26:50 +0000 (Mon, 06 Jul 2009) Log Message: ----------- Add svplot (thanks to Lukas Reichlin <luk...@sw...>). Modified Paths: -------------- trunk/octave-forge/main/control/INDEX Added Paths: ----------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2009-07-02 12:27:35 UTC (rev 5985) +++ trunk/octave-forge/main/control/INDEX 2009-07-06 08:26:50 UTC (rev 5986) @@ -64,6 +64,7 @@ Miscellaneous axis2dlim prompt sortcom swap run_cmd strappend listidx packedform sysidx + svplot Obsolete dezero dlqg minfo packsys unpacksys swaprows swapcols rotg qzval syschnames series Added: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m (rev 0) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-07-06 08:26:50 UTC (rev 5986) @@ -0,0 +1,126 @@ +## Copyright (C) 2009 Lukas Reichlin <luk...@sw...> +## +## 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. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}) +## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{w}) +## If no output arguments are given, the singular value plot of a MIMO +## system over a range of frequencies is printed on the screen; +## otherwise, the singular values of the system data structure are +## computed and returned. +## +## @strong{Inputs} +## @table @var +## @item sys +## System data structure (must be either purely continuous or discrete; +## see @code{is_digital}). +## @item w +## Optional vector of frequency values. If @var{w} is not specified, it +## will be calculated by @code{bode_bounds}. +## @end table +## +## @strong{Outputs} +## @table @var +## @item sigma_min +## Vector of minimal singular values. +## @item sigma_max +## Vector of maximal singular values. +## @item w +## Vector of frequency values used. +## @end table +## +## @seealso{is_digital} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@sw...> +## Version: 0.1alpha + + +function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w) + + if (nargin < 1 || nargin > 2) + print_usage(); + endif + + ## Get State Space Matrices + sys = sysupdate(sys, "ss"); + [A, B, C, D] = sys2ss(sys); + + + ## Find interesting Frequency Range w if not specified + if (nargin < 2) + + j = 1; # Index Counter + + for m = 1 : size(B, 2) # For all Inputs + for p = 1 : size(C, 1) # For all Outputs + + ## sys2zp and bode_bounds don't work for MIMO Systems! + sysp = sysprune(sys, p, m); + [zer, pol, k, tsam, inname, outname] = sys2zp(sysp); + + if (tsam == 0) + DIGITAL = 0; + else + DIGITAL = 1; + endif + + [wmin, wmax] = bode_bounds(zer, pol, DIGITAL, tsam); + + w_min(j) = wmin; + w_max(j) = wmax; + j++; + endfor + endfor + + dec_min = min(w_min); # Begin Plot at 10^dec_min rad/s + dec_max = max(w_max); # End Plot at 10^dec_min rad/s + n_steps = 1000; # Number of Steps + + w = logspace(dec_min, dec_max, n_steps); # rad/s + endif + + + ## Repeat for every Frequency w + for k = 1 : size(w, 2) + + ## Frequency Response Matrix + P = C * inv(i * w(k) * eye(size(A)) - A) * B + D; + + ## Singular Value Decomposition + sigma = svd(P); + sigma_min(k) = min(sigma); + sigma_max(k) = max(sigma); + endfor + + if (nargout < 1) # Plot the Information + + ## Convert to dB for Plotting + sigma_min_db = 20 * log10(sigma_min); + sigma_max_db = 20 * log10(sigma_max); + + ## Plot Results + semilogx(w, sigma_min_db, 'b', w, sigma_max_db, 'b') + title('Singular Values') + xlabel('Frequency [rad/s]') + ylabel('Singular Values [dB]') + grid on + else # Return Values + sigma_min_r = sigma_min; + sigma_max_r = sigma_max; + w_r = w; + endif + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <sla...@us...> - 2009-07-06 09:07:01
|
Revision: 5988 http://octave.svn.sourceforge.net/octave/?rev=5988&view=rev Author: slackydeb Date: 2009-07-06 09:06:55 +0000 (Mon, 06 Jul 2009) Log Message: ----------- Add margin (thanks to Lukas Reichlin <luk...@sw...>). Modified Paths: -------------- trunk/octave-forge/main/control/INDEX Added Paths: ----------- trunk/octave-forge/main/control/inst/margin.m Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2009-07-06 08:58:15 UTC (rev 5987) +++ trunk/octave-forge/main/control/INDEX 2009-07-06 09:06:55 UTC (rev 5988) @@ -55,6 +55,7 @@ bode bode_bounds nichols freqchkw ltifr nyquist tzero tzero2 rlocus + margin Controller design dgkfdemo hinfdemo dhinfdemo dlqe dlqr lqe lqg lqr Added: trunk/octave-forge/main/control/inst/margin.m =================================================================== --- trunk/octave-forge/main/control/inst/margin.m (rev 0) +++ trunk/octave-forge/main/control/inst/margin.m 2009-07-06 09:06:55 UTC (rev 5988) @@ -0,0 +1,76 @@ +## Copyright (C) 2009 Doug Stewart and Lukas Reichlin +## +## 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 2 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; if not, write to the Free Software +## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +## -*- texinfo -*- +## @deftypefn{Function File} {[@var{gamma}, @var{phi}, @var{w_gamma}, @var{w_phi}] =} margin (@var{sys}) +## Gain and phase margins. +## +## @strong{Inputs} +## @table @var +## @item sys +## Continuous time system. +## @end table +## +## @strong{Outputs} +## @table @var +## @item gamma +## Gain margin (as gain, not dbs). +## @item phi +## Phase margin (in degrees). +## @item w_gamma +## Radian frequency for the gain margin. +## @item w_phi +## Radian frequency for the phase margin. +## @end table +## +## @seealso{bode} +## @end deftypefn + +## Version: 0.1alpha + + +function [gamma, phi, w_gamma, w_phi] = margin(sys) + + if (nargin != 1 || (! isstruct (sys))) + print_usage (); + endif + + ## Get Frequency Response + [mag, pha, w] = bode (sys); + + ## fix the phase wrap around + phw = unwrap (pha * pi/180); + pha = phw * 180/pi; + + ## find the Gain Margin and its location + [x, ix] = min (abs (pha + 180)); + + if (x > 1) + gamma = "INF"; + ix = length (pha); + else + gamma = 1 / mag(ix); + endif + + w_gamma = w(ix); + + ## find the Phase Margin and its location + [pmx, ipmx] = min (abs (mag - 1)); + + phi = 180 + pha(ipmx); + w_phi = w(ipmx); + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <sla...@us...> - 2009-07-07 18:23:34
|
Revision: 5989 http://octave.svn.sourceforge.net/octave/?rev=5989&view=rev Author: slackydeb Date: 2009-07-07 18:23:28 +0000 (Tue, 07 Jul 2009) Log Message: ----------- svplot (thanks to Lukas Reichlin <luk...@sw...>): add checks, add comments, add a test. Bump version. Modified Paths: -------------- trunk/octave-forge/main/control/DESCRIPTION trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/DESCRIPTION =================================================================== --- trunk/octave-forge/main/control/DESCRIPTION 2009-07-06 09:06:55 UTC (rev 5988) +++ trunk/octave-forge/main/control/DESCRIPTION 2009-07-07 18:23:28 UTC (rev 5989) @@ -1,6 +1,6 @@ Name: Control -Version: 1.0.11 -Date: 2009-02-11 +Version: 1.0.12 +Date: 2009-07-07 Author: Ben Sapp Maintainer: Luca Favatella <sla...@gm...> Title: Control. Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-07-06 09:06:55 UTC (rev 5988) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-07-07 18:23:28 UTC (rev 5989) @@ -45,29 +45,49 @@ ## @end deftypefn ## Author: Lukas Reichlin <luk...@sw...> -## Version: 0.1alpha +## Version: 0.1 function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w) + ## Check whether arguments are OK if (nargin < 1 || nargin > 2) print_usage (); endif + if(! isstruct (sys)) + error ("first argument sys must be a system data structure"); + endif + + if (nargin == 2) + if (! isvector (w)) + error ("second argument w must be a vector of frequencies"); + endif + endif + ## Get State Space Matrices sys = sysupdate (sys, "ss"); - [A, B, C, D] = sys2ss (sys); + [A, B, C, D] = sys2ss (sys); + I = eye (size (A)); ## Find interesting Frequency Range w if not specified if (nargin < 2) + ## Since no frequency vector w has been specified, the interesting + ## decades of the sigma plot have to be found. The already existing + ## function bode_bounds does exactly that, unfortunately for SISO + ## systems only. Therefore, a SISO system from every input m to + ## every output p is created. Then for every SISO system the + ## interesting frequency range is calculated by bode_bounds. + ## Afterwards, the min/max decade can be found by the min()/max() + ## command. + j = 1; # Index Counter for m = 1 : size (B, 2) # For all Inputs for p = 1 : size (C, 1) # For all Outputs - ## sys2zp and bode_bounds don't work for MIMO Systems! sysp = sysprune (sys, p, m); [zer, pol, k, tsam, inname, outname] = sys2zp (sysp); @@ -75,6 +95,7 @@ DIGITAL = 0; else DIGITAL = 1; + ## FIXME: Discrete systems not yet tested! endif [wmin, wmax] = bode_bounds (zer, pol, DIGITAL, tsam); @@ -87,7 +108,7 @@ dec_min = min (w_min); # Begin Plot at 10^dec_min rad/s dec_max = max (w_max); # End Plot at 10^dec_min rad/s - n_steps = 1000; # Number of Steps + n_steps = 1000; # Number of Frequencies evaluated for Plotting w = logspace (dec_min, dec_max, n_steps); # rad/s endif @@ -97,7 +118,7 @@ for k = 1 : size (w, 2) ## Frequency Response Matrix - P = C * inv (i * w(k) * eye (size(A)) - A) * B + D; + P = C * inv (i*w(k)*I - A) * B + D; ## Singular Value Decomposition sigma = svd (P); @@ -105,7 +126,7 @@ sigma_max(k) = max (sigma); endfor - if (nargout < 1) # Plot the Information + if (nargout == 0) # Plot the Information ## Convert to dB for Plotting sigma_min_db = 20 * log10 (sigma_min); @@ -124,3 +145,18 @@ endif endfunction + + +%!shared sigma_min_exp, sigma_max_exp, w_exp, sigma_min_obs, sigma_max_obs, w_obs +%! A = [1 2; 3 4]; +%! B = [5 6; 7 8]; +%! C = [4 3; 2 1]; +%! D = [8 7; 6 5]; +%! w = [2 3]; +%! sigma_min_exp = [0.698526948925716 0.608629874340667]; +%! sigma_max_exp = [7.91760889977901 8.62745836756994]; +%! w_exp = [2 3]; +%! [sigma_min_obs, sigma_max_obs, w_obs] = svplot (ss (A, B, C, D), w); +%!assert (sigma_min_obs, sigma_min_exp, 4*eps); # tolerance manually tweaked +%!assert (sigma_max_obs, sigma_max_exp, 12*eps); # tolerance manually tweaked +%!assert (w_obs, w_exp, 2*eps); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-08-15 18:44:48
|
Revision: 6097 http://octave.svn.sourceforge.net/octave/?rev=6097&view=rev Author: paramaniac Date: 2009-08-15 18:44:39 +0000 (Sat, 15 Aug 2009) Log Message: ----------- add new function gensig Modified Paths: -------------- trunk/octave-forge/main/control/INDEX Added Paths: ----------- trunk/octave-forge/main/control/inst/gensig.m Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2009-08-14 06:32:41 UTC (rev 6096) +++ trunk/octave-forge/main/control/INDEX 2009-08-15 18:44:39 UTC (rev 6097) @@ -51,7 +51,7 @@ is_signal_list is_stable isct isdt Time domain analysis - c2d d2c dmr2d damp dcgain impulse step + c2d d2c dmr2d damp dcgain impulse step gensig Frequency domain analysis frdemo bode bode_bounds nichols Added: trunk/octave-forge/main/control/inst/gensig.m =================================================================== --- trunk/octave-forge/main/control/inst/gensig.m (rev 0) +++ trunk/octave-forge/main/control/inst/gensig.m 2009-08-15 18:44:39 UTC (rev 6097) @@ -0,0 +1,90 @@ +## Copyright (C) 2009 Lukas Reichlin <luk...@sw...> +## +## 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. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} {[@var{u}, @var{t}] =} gensig (@var{sigtype}, @var{tau}) +## @deftypefnx{Function File} {[@var{u}, @var{t}] =} gensig (@var{sigtype}, @var{tau}, @var{tfinal}) +## @deftypefnx{Function File} {[@var{u}, @var{t}] =} gensig (@var{sigtype}, @var{tau}, @var{tfinal}, @var{tsam}) +## Generate periodic signal. Useful in combination with lsim. +## +## @strong{Inputs} +## @table @var +## @item sigtype = "sin" +## Sine wave. +## @item sigtype = "cos" +## Cosine wave. +## @item sigtype = "square" +## Square wave. +## @item sigtype = "pulse" +## Periodic pulse. +## @item tau +## Duration of one period in seconds. +## @item tfinal +## Optional duration of the signal in seconds. Default duration is 5 periods +## @item tsam +## Optional sampling time in seconds. Default spacing is tau/64. +## @end table +## +## @strong{Outputs} +## @table @var +## @item u +## Vector of signal values. +## @item t +## Time vector of the signal. +## @end table +## +## @seealso{lsim} +## @end deftypefn + +## Version: 0.1 + +function [u, t] = gensig (sigtype, tau, tfinal, tsam) + + if (nargin < 2 || nargin > 4) + print_usage (); + endif + + if (! ischar (sigtype)) + error ("gensig: first argument must be a string"); + endif + + if (nargin < 3) + tfinal = 5 * tau; + endif + + if (nargin < 4) + tsam = tau / 64; + endif + + t = (0 : tsam : tfinal)'; + + sigtype = lower (sigtype); + + switch (sigtype(1:2)) + case "si" + u = sin (2*pi/tau * t); + case "co" + u = cos (2*pi/tau * t); + case "sq" + u = +(rem (t, tau) >= tau/2); + case "pu" + u = +(rem (t, tau) < (1 - 1000*eps) * tsam); + otherwise + error ("gensig: invalid signal type"); + endswitch + +endfunction + +## FIXME: Add a test \ No newline at end of file This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-08-24 21:05:11
|
Revision: 6134 http://octave.svn.sourceforge.net/octave/?rev=6134&view=rev Author: paramaniac Date: 2009-08-24 21:05:03 +0000 (Mon, 24 Aug 2009) Log Message: ----------- Add new function: initial.m Modified Paths: -------------- trunk/octave-forge/main/control/INDEX Added Paths: ----------- trunk/octave-forge/main/control/inst/initial.m Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2009-08-24 19:17:20 UTC (rev 6133) +++ trunk/octave-forge/main/control/INDEX 2009-08-24 21:05:03 UTC (rev 6134) @@ -51,7 +51,7 @@ is_signal_list is_stable isct isdt Time domain analysis - c2d d2c dmr2d damp dcgain impulse step gensig + c2d d2c dmr2d damp dcgain impulse step gensig initial Frequency domain analysis frdemo bode bode_bounds nichols Added: trunk/octave-forge/main/control/inst/initial.m =================================================================== --- trunk/octave-forge/main/control/inst/initial.m (rev 0) +++ trunk/octave-forge/main/control/inst/initial.m 2009-08-24 21:05:03 UTC (rev 6134) @@ -0,0 +1,202 @@ +## Copyright (C) 1996, 1998, 2000, 2003, 2004, 2005, 2006, 2007, 2009 +## Auburn University. All rights reserved. +## +## 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{y}, @var{t}, @var{x}] =} initial (@var{sys}, @var{x0}) +## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} initial (@var{sys}, @var{x0}, @var{tfinal}) +## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} initial (@var{sys}, @var{x0}, @var{tfinal}, @var{dt}) +## Initial condition response of state-space model. +## If no output arguments are given, the response is printed on the screen; +## otherwise, the response is computed and returned. +## +## @strong{Inputs} +## @table @var +## @item sys +## System data structure. Must be either purely continuous or discrete; +## see @code{is_digital}. +## @item x0 +## Vector of initial conditions for each state. +## @item tfinal +## Optional simulation horizon. If not specified, it will be calculated by +## the poles of the system to reflect adequately the response transients. +## @item dt +## Optional sampling time. Be sure to choose it small enough to capture transient +## phenomena. If not specified, it will be calculated by the poles of the system. +## @end table +## +## @strong{Outputs} +## @table @var +## @item y +## Output response array. Has as many rows as time samples (length of t) +## and as many columns as outputs. +## @item t +## Time row vector. +## @item x +## State trajectories array. Has length(t) rows and as many columns as states. +## @end table +## +## @seealso{impulse, lsim, step} +## @example +## @group +## . +## x = A x x(0) = x0 +## +## y = C x +## @end group +## @end example +## @end deftypefn + +## Author: Lukas Reichlin <luk...@sw...> +## Created: August 16, 2009 +## based on __stepimp__.m of Kai P. Mueller and A. Scottedward Hodel +## Version: 0.1 alpha +## FIXME: Needs thorough testing for discrete systems! + +function [y_r, t_r, x_r] = initial (sys, x_0, t_final, dt) + + ## check whether arguments are OK + if (nargin < 2 || nargin > 4) + print_usage (); + endif + + if (! isstruct (sys)) + error ("initial: first argument must be a system data structure"); + endif + + if (! isvector (x_0)) + error ("initial: second argument must be a vector"); + endif + + ## get system information + sys = sysupdate (sys, "ss"); + digital = is_digital (sys, 2); + [n_c, n_d, n_in, n_out] = sysdimensions (sys); + n_st = n_c + n_d; # number of states + + if (digital == -1) + error ("initial: system must be either purely continuous or purely discrete"); + endif + + if (n_st != length (x_0)) + error ("initial: x0 must be a vector with %d elements", n_st); + endif + + ## code adapted from __stepimp__.m + if (nargin < 3) + ## we have to compute the time when the system reaches steady state + ## and the step size + eigw = eig (sys2ss (sys)); + + if (digital) + t_sam = sysgettsam (sys); + + ## perform bilinear transformation on poles in z + for k = 1 : n_d + pole = eigw(k); + if (abs (pole + 1) < 1.0e-10) + eigw(k) = 0; + else + eigw(k) = 2 / t_sam * (pole - 1) / (pole + 1); + endif + endfor + endif + + ## remove poles near zero from eigenvalue array eigw + nk = n_st; + for k = 1 : n_st + if (abs (real (eigw(k))) < 1.0e-10) + eigw(k) = 0; + nk = nk - 1; + endif + endfor + + if (nk == 0) + t_final = 10; + dt = t_final / 1000; + else + eigw = eigw(find (eigw)); + eigw_max = max (abs (eigw)); + dt = 0.2 * pi / eigw_max; + + eigw_min = min (abs (real (eigw))); + t_final = 5.0 / eigw_min; + + ## round up + yy = 10^(ceil (log10 (t_final)) - 1); + t_final = yy * ceil (t_final / yy); + + n = t_final / dt; + + if (n < 50) + dt = t_final / 50; + endif + + if (n > 2000) + dt = t_final / 2000; + endif + endif + endif + ## end of adapted code + + if (digital) + dt = sysgettsam (sys); + + if (nargin == 4) + warning ("initial: fourth argument has no effect on sampling time of digital system"); + endif + else + sys = c2d (sys, dt); + endif + + t = (0 : dt : t_final)'; + l_t = length (t); + + [F, G, C, D] = sys2ss (sys); + + ## preallocate memory + y = zeros (l_t, n_out); + x_vec = zeros (l_t, n_st); + + ## make sure that x is a row vector + x = reshape (x_0, length (x_0), 1); + + ## simulation + for k = 1 : l_t + y(k, :) = C * x; + x_vec(k, :) = x; + x = F * x; + endfor + + if (nargout == 0) # plot information + for k = 1 : n_out + subplot (n_out, 1, k) + plot (t, y(:, k)) + if (k == 1) + title ("Response to Initial Conditions") + endif + ylabel (sprintf ("Amplitude %s", sysgetsignals (sys, "out", k, 1))) + endfor + xlabel ("Time [s]") + else # return values + y_r = y; + t_r = t; + x_r = x_vec; + endif + +endfunction + +## FIXME: Add a test \ No newline at end of file This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-08-29 15:13:16
|
Revision: 6170 http://octave.svn.sourceforge.net/octave/?rev=6170&view=rev Author: paramaniac Date: 2009-08-29 15:13:04 +0000 (Sat, 29 Aug 2009) Log Message: ----------- add new function sigma.m Modified Paths: -------------- trunk/octave-forge/main/control/INDEX Added Paths: ----------- trunk/octave-forge/main/control/inst/sigma.m Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2009-08-29 10:51:40 UTC (rev 6169) +++ trunk/octave-forge/main/control/INDEX 2009-08-29 15:13:04 UTC (rev 6170) @@ -36,7 +36,7 @@ sysdup sysgroup sysmult sysprune sysreorder sysscale syssub ugain wgt1o starp series - feedback sysfeedback + feedback sysfeedback sysparallel buildsens buildretdiff Numerical functions are dare dgram dlyap gram lyap pinv dre @@ -57,7 +57,7 @@ bode bode_bounds nichols freqchkw ltifr nyquist tzero tzero2 rlocus - margin svplot + margin sigma svplot Controller design dgkfdemo hinfdemo dhinfdemo dlqe dlqr lqe lqg lqr @@ -90,4 +90,3 @@ kalman=use <f>lqe</f> kalmd=use <f>dlqe</f> or <f>dkalman</f> isstable=use <f>is_stable</f> - sigma=use <f>svplot</f> Added: trunk/octave-forge/main/control/inst/sigma.m =================================================================== --- trunk/octave-forge/main/control/inst/sigma.m (rev 0) +++ trunk/octave-forge/main/control/inst/sigma.m 2009-08-29 15:13:04 UTC (rev 6170) @@ -0,0 +1,210 @@ +## Copyright (C) 2009 Lukas Reichlin <luk...@sw...> +## +## 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. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} {[@var{sv}, @var{w}] =} sigma (@var{sys}) +## @deftypefnx{Function File} {[@var{sv}, @var{w}] =} sigma (@var{sys}, @var{w}) +## @deftypefnx{Function File} {[@var{sv}, @var{w}] =} sigma (@var{sys}, @var{[]}, @var{ptype}) +## @deftypefnx{Function File} {[@var{sv}, @var{w}] =} sigma (@var{sys}, @var{w}, @var{ptype}) +## If no output arguments are given, the singular value plot of a MIMO +## system is printed on the screen; +## otherwise, the singular values of the system data structure are +## computed and returned. +## +## @strong{Inputs} +## @table @var +## @item sys +## System data structure. Must be either purely continuous or discrete; +## see @code{is_digital}. +## @item w +## Optional vector of frequency values. If @var{w} is not specified, it +## will be calculated by @code{bode_bounds}. +## @item ptype = 0 +## Singular values of the frequency response H of system sys. Default Value. +## @item ptype = 1 +## Singular values of the frequency response inv(H); i.e. inversed system. +## @item ptype = 2 +## Singular values of the frequency response I + H; i.e. inversed sensitivity +## (or return difference) if H = P * C. +## @item ptype = 3 +## Singular values of the frequency response I + inv(H); i.e. inversed complementary +## sensitivity if H = P * C. +## @end table +## +## @strong{Outputs} +## @table @var +## @item sv +## Array of singular values. For a system with m inputs and p outputs, the array sv +## has min(m,p) rows and as many columns as frequency points (length of w). +## The singular values at the frequency w(k) are given by sv(:,k). +## @item w +## Vector of frequency values used. +## @end table +## +## @seealso{bode, svd, bode_bounds, is_digital} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@sw...> +## Version: 0.1 + +function [sv_r, w_r] = sigma (sys, w, ptype) + + ## Check whether arguments are OK + if (nargin < 1 || nargin > 3) + print_usage (); + endif + + if(! isstruct (sys)) + error ("sigma: first argument sys must be a system data structure"); + endif + + if (nargin == 1) + w = []; + elseif (! isvector (w) && ! isempty (w)) + error ("sigma: second argument w must be a vector of frequencies"); + endif + + if (is_siso (sys)) + warning ("sigma: sys is a SISO system. You might want to use bode(sys) instead."); + endif + + ## Get state space matrices + sys = sysupdate (sys, "ss"); + [A, B, C, D] = sys2ss (sys); + I = eye (size (A)); + + ## Get system information + [n_c, n_d, m, p] = sysdimensions (sys); + digital = is_digital (sys, 2); + t_sam = sysgettsam (sys); + + ## Error for mixed systems + if (digital == -1) + error ("svplot: system must be either purely continuous or purely discrete"); + endif + + ## Handle plot type + if (nargin == 3) + if (isfloat (ptype)) # Numeric constants like 2 are NOT integers in Octave! + if (ptype != 0 && ptype != 1 && ptype != 2 && ptype != 3) + error ("sigma: third argument ptype must be 0, 1, 2 or 3"); + endif + else + error ("sigma: third argument ptype must be a number"); + endif + + if (m != p) + error ("sigma: system must be square for ptype 1, 2 or 3"); + endif + J = eye (m); + else + ptype = 0; # Default value + endif + + ## Find interesting frequency range w if not specified + if (isempty (w)) + + ## Since no frequency vector w has been specified, the interesting + ## decades of the sigma plot have to be found. The already existing + ## function bode_bounds does exactly that. + + zer = tzero (sys); + pol = eig (A); + + ## Begin plot at 10^dec_min, end plot at 10^dec_max [rad/s] + [dec_min, dec_max] = bode_bounds (zer, pol, digital, t_sam); + + n_freq = 1000; # Number of frequencies evaluated for plotting + + w = logspace (dec_min, dec_max, n_freq); # [rad/s] + endif + + if (digital) # Discrete system + s = exp (i * w * t_sam); + else # Continuous system + s = i * w; + endif + + l_s = length (s); + sv = zeros (min (m, p), l_s); + + switch (ptype) + case 0 # Default system + for k = 1 : l_s # Repeat for every frequency s + H = C * inv (s(k)*I - A) * B + D; # Frequency Response Matrix + sv(:, k) = svd (H); # Singular Value Decomposition + endfor + + case 1 # Inversed system + for k = 1 : l_s + H = inv (C * inv (s(k)*I - A) * B + D); + sv(:, k) = svd (H); + endfor + + case 2 # Inversed sensitivity + for k = 1 : l_s + H = J + C * inv (s(k)*I - A) * B + D; + sv(:, k) = svd (H); + endfor + + case 3 # Inversed complementary sensitivity + for k = 1 : l_s + H = J + inv (C * inv (s(k)*I - A) * B + D); + sv(:, k) = svd (H); + endfor + endswitch + + if (nargout == 0) # Plot the information + + ## Convert to dB for plotting + sv_db = 20 * log10 (sv); + + ## Determine axes + ax_vec = axis2dlim ([[w(:), min(sv_db)(:)]; [w(:), max(sv_db)(:)]]); + ax_vec(1:2) = [min(w), max(w)]; + + ## Determine xlabel + if (digital) + xl_str = sprintf ("Frequency [rad/s] Pi / T = %g", pi/t_sam); + else + xl_str = "Frequency [rad/s]"; + endif + + ## Plot results + semilogx (w, sv_db, "b") + axis (ax_vec) + title ("Singular Values") + xlabel (xl_str) + ylabel ("Singular Values [dB]") + grid on + else # Return values + sv_r = sv; + w_r = w; + endif + +endfunction + + +%!shared sv_exp, w_exp, sv_obs, w_obs +%! A = [1 2; 3 4]; +%! B = [5 6; 7 8]; +%! C = [4 3; 2 1]; +%! D = [8 7; 6 5]; +%! w = [2 3]; +%! sv_exp = [7.91760889977901, 8.62745836756994; 0.698526948925716, 0.608629874340667]; +%! w_exp = [2 3]; +%! [sv_obs, w_obs] = sigma (ss (A, B, C, D), w); +%!assert (sv_obs, sv_exp, 16*eps); # tolerance manually tweaked +%!assert (w_obs, w_exp, 2*eps); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <sla...@us...> - 2009-08-31 09:38:49
|
Revision: 6175 http://octave.svn.sourceforge.net/octave/?rev=6175&view=rev Author: slackydeb Date: 2009-08-31 09:38:37 +0000 (Mon, 31 Aug 2009) Log Message: ----------- Remove 'svplot' keeping only 'sigma' for maintainability and compatibility. 'svplot' was different from 'sigma' in the return values and in the plots for systems with more than two inputs and/or outputs. Modified Paths: -------------- trunk/octave-forge/main/control/INDEX trunk/octave-forge/main/control/inst/buildretdiff.m trunk/octave-forge/main/control/inst/buildsens.m trunk/octave-forge/main/control/inst/sigma.m Removed Paths: ------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2009-08-30 20:41:24 UTC (rev 6174) +++ trunk/octave-forge/main/control/INDEX 2009-08-31 09:38:37 UTC (rev 6175) @@ -57,7 +57,7 @@ bode bode_bounds nichols freqchkw ltifr nyquist tzero tzero2 rlocus - margin sigma svplot + margin sigma Controller design dgkfdemo hinfdemo dhinfdemo dlqe dlqr lqe lqg lqr Modified: trunk/octave-forge/main/control/inst/buildretdiff.m =================================================================== --- trunk/octave-forge/main/control/inst/buildretdiff.m 2009-08-30 20:41:24 UTC (rev 6174) +++ trunk/octave-forge/main/control/inst/buildretdiff.m 2009-08-31 09:38:37 UTC (rev 6175) @@ -29,8 +29,6 @@ ## System data structure of the return difference. ## @end table ## -## @seealso{svplot} -## ## @example ## @group ## Q(s) = I + L(s) = I + P(s) C(s) @@ -39,7 +37,7 @@ ## ## @end deftypefn -## Version: 0.1.2 +## Version: 0.1.3 function retsys = buildretdiff (sys) Modified: trunk/octave-forge/main/control/inst/buildsens.m =================================================================== --- trunk/octave-forge/main/control/inst/buildsens.m 2009-08-30 20:41:24 UTC (rev 6174) +++ trunk/octave-forge/main/control/inst/buildsens.m 2009-08-31 09:38:37 UTC (rev 6175) @@ -29,8 +29,6 @@ ## Sensitivity system data structure. ## @end table ## -## @seealso{svplot} -## ## @example ## @group ## Y(s) 1 1 @@ -61,7 +59,7 @@ ## ## @end deftypefn -## Version: 0.1.2 +## Version: 0.1.3 function retsys = buildsens (sys) Modified: trunk/octave-forge/main/control/inst/sigma.m =================================================================== --- trunk/octave-forge/main/control/inst/sigma.m 2009-08-30 20:41:24 UTC (rev 6174) +++ trunk/octave-forge/main/control/inst/sigma.m 2009-08-31 09:38:37 UTC (rev 6175) @@ -57,7 +57,7 @@ ## @end deftypefn ## Author: Lukas Reichlin <luk...@sw...> -## Version: 0.1 +## Version: 0.1.1 function [sv_r, w_r] = sigma (sys, w, ptype) @@ -92,7 +92,7 @@ ## Error for mixed systems if (digital == -1) - error ("svplot: system must be either purely continuous or purely discrete"); + error ("sigma: system must be either purely continuous or purely discrete"); endif ## Handle plot type Deleted: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-08-30 20:41:24 UTC (rev 6174) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-08-31 09:38:37 UTC (rev 6175) @@ -1,223 +0,0 @@ -## Copyright (C) 2009 Lukas Reichlin <luk...@sw...> -## -## 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. If not, see <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}) -## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{w}) -## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{[]}, @var{ptype}) -## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{w}, @var{ptype}) -## If no output arguments are given, the singular value plot of a MIMO -## system is printed on the screen; -## otherwise, the singular values of the system data structure are -## computed and returned. -## -## @strong{Inputs} -## @table @var -## @item sys -## System data structure. Must be either purely continuous or discrete; -## see @code{is_digital}. -## @item w -## Optional vector of frequency values. If @var{w} is not specified, it -## will be calculated by @code{bode_bounds}. -## @item ptype = 0 -## Singular values of the frequency response H of system sys. Default Value. -## @item ptype = 1 -## Singular values of the frequency response inv(H); i.e. inversed system. -## @item ptype = 2 -## Singular values of the frequency response I + H; i.e. inversed sensitivity -## (or return difference) if H = P * C. -## @item ptype = 3 -## Singular values of the frequency response I + inv(H); i.e. inversed complementary -## sensitivity if H = P * C. -## @end table -## -## @strong{Outputs} -## @table @var -## @item sigma_min -## Vector of minimal singular values. -## @item sigma_max -## Vector of maximal singular values. -## @item w -## Vector of frequency values used. -## @end table -## -## @seealso{bode, svd, bode_bounds, is_digital} -## @end deftypefn - -## Author: Lukas Reichlin <luk...@sw...> -## Version: 0.3 - -function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w, ptype) - - ## Check whether arguments are OK - if (nargin < 1 || nargin > 3) - print_usage (); - endif - - if(! isstruct (sys)) - error ("svplot: first argument sys must be a system data structure"); - endif - - if (nargin == 1) - w = []; - elseif (! isvector (w) && ! isempty (w)) - error ("svplot: second argument w must be a vector of frequencies"); - endif - - if (is_siso (sys)) - warning ("svplot: sys is a SISO system. You might want to use bode(sys) instead."); - endif - - ## Get state space matrices - sys = sysupdate (sys, "ss"); - [A, B, C, D] = sys2ss (sys); - I = eye (size (A)); - - ## Get system information - digital = is_digital (sys, 2); - t_sam = sysgettsam (sys); - - ## Error for mixed systems - if (digital == -1) - error ("svplot: system must be either purely continuous or purely discrete"); - endif - - ## Handle plot type - if (nargin == 3) - if (isfloat (ptype)) # Numeric constants like 2 are NOT integers in Octave! - if (ptype != 0 && ptype != 1 && ptype != 2 && ptype != 3) - error ("svplot: third argument ptype must be 0, 1, 2 or 3"); - endif - else - error ("svplot: third argument ptype must be a number"); - endif - - [n_c, n_d, m, p] = sysdimensions (sys); - if (m != p) - error ("svplot: system must be square for ptype 1, 2 or 3"); - endif - J = eye (m); - else - ptype = 0; # Default value - endif - - ## Find interesting frequency range w if not specified - if (isempty (w)) - - ## Since no frequency vector w has been specified, the interesting - ## decades of the sigma plot have to be found. The already existing - ## function bode_bounds does exactly that. - - zer = tzero (sys); - pol = eig (A); - - ## Begin plot at 10^dec_min, end plot at 10^dec_max [rad/s] - [dec_min, dec_max] = bode_bounds (zer, pol, digital, t_sam); - - n_freq = 1000; # Number of frequencies evaluated for plotting - - w = logspace (dec_min, dec_max, n_freq); # [rad/s] - endif - - if (digital) # Discrete system - s = exp (i * w * t_sam); - else # Continuous system - s = i * w; - endif - - l_s = length (s); - sigma_min = zeros (1, l_s); - sigma_max = zeros (1, l_s); - - switch (ptype) - case 0 # Default system - for k = 1 : l_s # Repeat for every frequency s - H = C * inv (s(k)*I - A) * B + D; # Frequency Response Matrix - sigma = svd (H); # Singular Value Decomposition - sigma_min(k) = min (sigma); - sigma_max(k) = max (sigma); - endfor - - case 1 # Inversed system - for k = 1 : l_s - H = inv (C * inv (s(k)*I - A) * B + D); - sigma = svd (H); - sigma_min(k) = min (sigma); - sigma_max(k) = max (sigma); - endfor - - case 2 # Inversed sensitivity - for k = 1 : l_s - H = J + C * inv (s(k)*I - A) * B + D; - sigma = svd (H); - sigma_min(k) = min (sigma); - sigma_max(k) = max (sigma); - endfor - - case 3 # Inversed complementary sensitivity - for k = 1 : l_s - H = J + inv (C * inv (s(k)*I - A) * B + D); - sigma = svd (H); - sigma_min(k) = min (sigma); - sigma_max(k) = max (sigma); - endfor - endswitch - - if (nargout == 0) # Plot the information - - ## Convert to dB for plotting - sigma_min_db = 20 * log10 (sigma_min); - sigma_max_db = 20 * log10 (sigma_max); - - ## Determine axes - ax_vec = axis2dlim ([[w(:), sigma_min_db(:)]; [w(:), sigma_max_db(:)]]); - ax_vec(1:2) = [min(w), max(w)]; - - ## Determine xlabel - if (digital) - xl_str = sprintf ('Frequency [rad/s] Pi / T = %g', pi/t_sam); - else - xl_str = 'Frequency [rad/s]'; - endif - - ## Plot results - semilogx (w, sigma_min_db, 'b', w, sigma_max_db, 'b') - axis (ax_vec) - title ('Singular Values') - xlabel (xl_str) - ylabel ('Singular Values [dB]') - grid on - else # Return values - sigma_min_r = sigma_min; - sigma_max_r = sigma_max; - w_r = w; - endif - -endfunction - - -%!shared sigma_min_exp, sigma_max_exp, w_exp, sigma_min_obs, sigma_max_obs, w_obs -%! A = [1 2; 3 4]; -%! B = [5 6; 7 8]; -%! C = [4 3; 2 1]; -%! D = [8 7; 6 5]; -%! w = [2 3]; -%! sigma_min_exp = [0.698526948925716 0.608629874340667]; -%! sigma_max_exp = [7.91760889977901 8.62745836756994]; -%! w_exp = [2 3]; -%! [sigma_min_obs, sigma_max_obs, w_obs] = svplot (ss (A, B, C, D), w); -%!assert (sigma_min_obs, sigma_min_exp, 8*eps); # tolerance manually tweaked -%!assert (sigma_max_obs, sigma_max_exp, 16*eps); # tolerance manually tweaked -%!assert (w_obs, w_exp, 2*eps); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <sla...@us...> - 2009-11-30 00:18:43
|
Revision: 6556 http://octave.svn.sourceforge.net/octave/?rev=6556&view=rev Author: slackydeb Date: 2009-11-30 00:16:15 +0000 (Mon, 30 Nov 2009) Log Message: ----------- control: __axis2dlim__: merge from control-oo. Modified Paths: -------------- trunk/octave-forge/main/control/INDEX trunk/octave-forge/main/control/inst/bode.m trunk/octave-forge/main/control/inst/margin.m trunk/octave-forge/main/control/inst/nichols.m trunk/octave-forge/main/control/inst/nyquist.m trunk/octave-forge/main/control/inst/pzmap.m trunk/octave-forge/main/control/inst/rlocus.m trunk/octave-forge/main/control/inst/sigma.m Added Paths: ----------- trunk/octave-forge/main/control/inst/__axis2dlim__.m Removed Paths: ------------- trunk/octave-forge/main/control/inst/axis2dlim.m Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2009-11-29 19:17:14 UTC (rev 6555) +++ trunk/octave-forge/main/control/INDEX 2009-11-30 00:16:15 UTC (rev 6556) @@ -65,7 +65,7 @@ lsim place dkalman h2syn acker Miscellaneous - axis2dlim prompt sortcom swap run_cmd + prompt sortcom swap run_cmd strappend listidx packedform sysidx Obsolete dezero dlqg minfo packsys unpacksys swaprows swapcols rotg qzval @@ -76,6 +76,7 @@ __bodquist__ __freqresp__ __stepimp__ __abcddims__ __syschnamesl__ __syscont_disc__ __sysdefioname__ __sysdefstname__ __sysgroupn__ __tf2sysl__ __zp2ssg2__ __outlist__ __zgpbal__ + __axis2dlim__ Compatibility filt= use <f>tf2sys</f> zpk= use <f>zp</f> Copied: trunk/octave-forge/main/control/inst/__axis2dlim__.m (from rev 6555, trunk/octave-forge/main/control/inst/axis2dlim.m) =================================================================== --- trunk/octave-forge/main/control/inst/__axis2dlim__.m (rev 0) +++ trunk/octave-forge/main/control/inst/__axis2dlim__.m 2009-11-30 00:16:15 UTC (rev 6556) @@ -0,0 +1,71 @@ +## Copyright (C) 1998, 2000, 2004, 2005, 2007 +## Auburn University. All rights reserved. +## +## +## 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} {} __axis2dlim__ (@var{axdata}) +## Determine axis limits for 2-D data (column vectors); leaves a 10% +## margin around the plots. +## Inserts margins of +/- 0.1 if data is one-dimensional +## (or a single point). +## +## @strong{Input} +## @table @var +## @item axdata +## @var{n} by 2 matrix of data [@var{x}, @var{y}]. +## @end table +## +## @strong{Output} +## @table @var +## @item axvec +## Vector of axis limits appropriate for call to @command{axis} function. +## @end table +## @end deftypefn + +function axvec = __axis2dlim__ (axdata) + + if (nargin < 1 || isempty (axdata)) + axdata = 0; + endif + + ## compute axis limits + minv = min (axdata); + maxv = max (axdata); + delv = (maxv-minv)/2; # breadth of the plot + midv = (minv + maxv)/2; # midpoint of the plot + axmid = [midv(1), midv(1), midv(2), midv(2)]; + axdel = [-0.1, 0.1, -0.1, 0.1]; # default plot width (if less than 2-d data) + if (max (delv) == 0) + if (midv(1) != 0) + axdel(1:2) = [-0.1*midv(1), 0.1*midv(1)]; + endif + if (midv(2) != 0) + axdel(3:4) = [-0.1*midv(2), 0.1*midv(2)]; + endif + else + ## they're at least one-dimensional + tolv = max(1e-8, 1e-8*abs(midv)); + if (abs (delv(1)) >= tolv(1)) + axdel(1:2) = 1.1*[-delv(1),delv(1)]; + endif + if (abs (delv(2)) >= tolv(2)) + axdel(3:4) = 1.1*[-delv(2),delv(2)]; + endif + endif + axvec = axmid + axdel; + +endfunction Deleted: trunk/octave-forge/main/control/inst/axis2dlim.m =================================================================== --- trunk/octave-forge/main/control/inst/axis2dlim.m 2009-11-29 19:17:14 UTC (rev 6555) +++ trunk/octave-forge/main/control/inst/axis2dlim.m 2009-11-30 00:16:15 UTC (rev 6556) @@ -1,71 +0,0 @@ -## Copyright (C) 1998, 2000, 2004, 2005, 2007 -## Auburn University. All rights reserved. -## -## -## 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} {} axis2dlim (@var{axdata}) -## Determine axis limits for 2-D data (column vectors); leaves a 10% -## margin around the plots. -## Inserts margins of +/- 0.1 if data is one-dimensional -## (or a single point). -## -## @strong{Input} -## @table @var -## @item axdata -## @var{n} by 2 matrix of data [@var{x}, @var{y}]. -## @end table -## -## @strong{Output} -## @table @var -## @item axvec -## Vector of axis limits appropriate for call to @command{axis} function. -## @end table -## @end deftypefn - -function axvec = axis2dlim (axdata) - - if (nargin < 1 || isempty (axdata)) - axdata = 0; - endif - - ## compute axis limits - minv = min (axdata); - maxv = max (axdata); - delv = (maxv-minv)/2; # breadth of the plot - midv = (minv + maxv)/2; # midpoint of the plot - axmid = [midv(1), midv(1), midv(2), midv(2)]; - axdel = [-0.1, 0.1, -0.1, 0.1]; # default plot width (if less than 2-d data) - if (max (delv) == 0) - if (midv(1) != 0) - axdel(1:2) = [-0.1*midv(1), 0.1*midv(1)]; - endif - if (midv(2) != 0) - axdel(3:4) = [-0.1*midv(2), 0.1*midv(2)]; - endif - else - ## they're at least one-dimensional - tolv = max(1e-8, 1e-8*abs(midv)); - if (abs (delv(1)) >= tolv(1)) - axdel(1:2) = 1.1*[-delv(1),delv(1)]; - endif - if (abs (delv(2)) >= tolv(2)) - axdel(3:4) = 1.1*[-delv(2),delv(2)]; - endif - endif - axvec = axmid + axdel; - -endfunction Modified: trunk/octave-forge/main/control/inst/bode.m =================================================================== --- trunk/octave-forge/main/control/inst/bode.m 2009-11-29 19:17:14 UTC (rev 6555) +++ trunk/octave-forge/main/control/inst/bode.m 2009-11-30 00:16:15 UTC (rev 6556) @@ -166,7 +166,7 @@ semilogx (w, md); if (max_mag_positive) ylabel ("Gain in dB"); - axvec = axis2dlim ([w(:), md(:)]); + axvec = __axis2dlim__ ([w(:), md(:)]); axvec(1:2) = wv; axis (axvec); endif @@ -189,7 +189,7 @@ if (is_siso_sys) subplot (2, 1, 2); - axvec = axis2dlim ([w(:), phase(:)]); + axvec = __axis2dlim__ ([w(:), phase(:)]); axvec(1:2) = wv; semilogx (w, phase); axis (axvec); Modified: trunk/octave-forge/main/control/inst/margin.m =================================================================== --- trunk/octave-forge/main/control/inst/margin.m 2009-11-29 19:17:14 UTC (rev 6555) +++ trunk/octave-forge/main/control/inst/margin.m 2009-11-30 00:16:15 UTC (rev 6556) @@ -264,9 +264,9 @@ gamma_db = 20 * log10 (gamma); wv = [min(w), max(w)]; - ax_vec_mag = axis2dlim ([w(:), mag_db(:)]); + ax_vec_mag = __axis2dlim__ ([w(:), mag_db(:)]); ax_vec_mag(1:2) = wv; - ax_vec_pha = axis2dlim ([w(:), pha(:)]); + ax_vec_pha = __axis2dlim__ ([w(:), pha(:)]); ax_vec_pha(1:2) = wv; wgm = [w_gamma, w_gamma]; Modified: trunk/octave-forge/main/control/inst/nichols.m =================================================================== --- trunk/octave-forge/main/control/inst/nichols.m 2009-11-29 19:17:14 UTC (rev 6555) +++ trunk/octave-forge/main/control/inst/nichols.m 2009-11-30 00:16:15 UTC (rev 6556) @@ -139,7 +139,7 @@ __outlist__ (outname, " ")); endif - axis (axis2dlim ([phase(:), md(:)])); + axis (__axis2dlim__ ([phase(:), md(:)])); else mag2 = mag; phase2 = phase; Modified: trunk/octave-forge/main/control/inst/nyquist.m =================================================================== --- trunk/octave-forge/main/control/inst/nyquist.m 2009-11-29 19:17:14 UTC (rev 6555) +++ trunk/octave-forge/main/control/inst/nyquist.m 2009-11-30 00:16:15 UTC (rev 6556) @@ -152,7 +152,7 @@ inn{1}, outn{1}, w(1), w(end))); endif - axis (axis2dlim ([[realp(:), imagp(:)]; [realp(:), -imagp(:)]])); + axis (__axis2dlim__ ([[realp(:), imagp(:)]; [realp(:), -imagp(:)]])); ## check for interactive plots dnplot = 1; # assume done; will change later if atol is satisfied Modified: trunk/octave-forge/main/control/inst/pzmap.m =================================================================== --- trunk/octave-forge/main/control/inst/pzmap.m 2009-11-29 19:17:14 UTC (rev 6555) +++ trunk/octave-forge/main/control/inst/pzmap.m 2009-11-30 00:16:15 UTC (rev 6556) @@ -83,6 +83,6 @@ grid ("on"); ## compute axis limits - axis (axis2dlim ([zerdata; poldata])); + axis (__axis2dlim__ ([zerdata; poldata])); endfunction Modified: trunk/octave-forge/main/control/inst/rlocus.m =================================================================== --- trunk/octave-forge/main/control/inst/rlocus.m 2009-11-29 19:17:14 UTC (rev 6555) +++ trunk/octave-forge/main/control/inst/rlocus.m 2009-11-30 00:16:15 UTC (rev 6556) @@ -217,7 +217,7 @@ if (nargout == 0) rlpolv = vec(rlpol); axdata = [real(rlpolv), imag(rlpolv); real(olzer), imag(olzer)]; - axlim = axis2dlim (axdata); + axlim = __axis2dlim__ (axdata); rldata = [real(rlpolv), imag(rlpolv) ]; [stn, inname, outname] = sysgetsignals (sys); Modified: trunk/octave-forge/main/control/inst/sigma.m =================================================================== --- trunk/octave-forge/main/control/inst/sigma.m 2009-11-29 19:17:14 UTC (rev 6555) +++ trunk/octave-forge/main/control/inst/sigma.m 2009-11-30 00:16:15 UTC (rev 6556) @@ -166,7 +166,7 @@ sv_db = 20 * log10 (sv); ## Determine axes - ax_vec = axis2dlim ([[w(:), min(sv_db)(:)]; [w(:), max(sv_db)(:)]]); + ax_vec = __axis2dlim__ ([[w(:), min(sv_db)(:)]; [w(:), max(sv_db)(:)]]); ax_vec(1:2) = [min(w), max(w)]; ## Determine xlabel This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-04-18 07:33:48
|
Revision: 7208 http://octave.svn.sourceforge.net/octave/?rev=7208&view=rev Author: paramaniac Date: 2010-04-18 07:33:41 +0000 (Sun, 18 Apr 2010) Log Message: ----------- control: revamp code for oct files Modified Paths: -------------- trunk/octave-forge/main/control/DESCRIPTION trunk/octave-forge/main/control/src/slab08nd.cc trunk/octave-forge/main/control/src/slab13ad.cc trunk/octave-forge/main/control/src/slab13bd.cc trunk/octave-forge/main/control/src/slab13dd.cc trunk/octave-forge/main/control/src/slsb01bd.cc trunk/octave-forge/main/control/src/slsb02od.cc trunk/octave-forge/main/control/src/slsb03md.cc trunk/octave-forge/main/control/src/slsb04md.cc trunk/octave-forge/main/control/src/slsb04qd.cc trunk/octave-forge/main/control/src/slsb10dd.cc trunk/octave-forge/main/control/src/slsb10ed.cc trunk/octave-forge/main/control/src/slsb10fd.cc trunk/octave-forge/main/control/src/slsb10hd.cc trunk/octave-forge/main/control/src/slsg03ad.cc Added Paths: ----------- trunk/octave-forge/main/control/src/common.cc Modified: trunk/octave-forge/main/control/DESCRIPTION =================================================================== --- trunk/octave-forge/main/control/DESCRIPTION 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/DESCRIPTION 2010-04-18 07:33:41 UTC (rev 7208) @@ -1,6 +1,6 @@ Name: Control -Version: 0.3.0 -Date: 2010-03-25 +Version: 0.3.1 +Date: 2010-04-18 Author: Lukas Reichlin <luk...@gm...> Maintainer: Lukas Reichlin <luk...@gm...> Title: LTI Syncope Added: trunk/octave-forge/main/control/src/common.cc =================================================================== --- trunk/octave-forge/main/control/src/common.cc (rev 0) +++ trunk/octave-forge/main/control/src/common.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -0,0 +1,58 @@ +/* + +Copyright (C) 2010 Lukas F. Reichlin + +This file is part of LTI Syncope. + +LTI Syncope 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. + +LTI Syncope 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. If not, see <http://www.gnu.org/licenses/>. + +Common code for oct-files. + +Author: Lukas Reichlin <luk...@gm...> +Created: April 2010 +Version: 0.1 + +*/ + + +int max (int a, int b) +{ + if (a > b) + return a; + else + return b; +} + +int max (int a, int b, int c) +{ + int d = max (a, b); + + return max (c, d); +} + +int max (int a, int b, int c, int d) +{ + int e = max (a, b); + int f = max (c, d); + + return max (e, f); +} + +int min (int a, int b) +{ + if (a < b) + return a; + else + return b; +} Modified: trunk/octave-forge/main/control/src/slab08nd.cc =================================================================== --- trunk/octave-forge/main/control/src/slab08nd.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slab08nd.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: November 2009 -Version: 0.2 +Version: 0.3 */ #include <octave/oct.h> -#include "f77-fcn.h" +#include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -60,23 +61,7 @@ double* WORK, int& LWORK, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int min (int a, int b) -{ - if (a < b) - return a; - else - return b; -} - + DEFUN_DLD (slab08nd, args, nargout, "Slicot AB08ND Release 5.0") { int nargin = args.length (); @@ -91,10 +76,10 @@ // arguments in char equil = 'N'; - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray c = args(2).array_value (); - NDArray d = args(3).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); int n = a.rows (); // n: number of states int m = b.columns (); // m: number of inputs @@ -171,11 +156,9 @@ int lwork = max (1, 8*nu); OCTAVE_LOCAL_BUFFER (double, work, lwork); - dim_vector dv (1); - dv(0) = nu; - NDArray alphar (dv); - NDArray alphai (dv); - NDArray beta (dv); + ColumnVector alphar (nu); + ColumnVector alphai (nu); + ColumnVector beta (nu); int info2; Modified: trunk/octave-forge/main/control/src/slab13ad.cc =================================================================== --- trunk/octave-forge/main/control/src/slab13ad.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slab13ad.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: January 2010 -Version: 0.1 +Version: 0.2 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -44,21 +45,6 @@ double* DWORK, int& LDWORK, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int max (int a, int b, int c) -{ - int d = max (a, b); - - return max (c, d); -} DEFUN_DLD (slab13ad, args, nargout, "Slicot AB13AD Release 5.0") { @@ -75,9 +61,9 @@ char dico; char equil = 'N'; - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray c = args(2).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); int digital = args(3).int_value (); double alpha = args(4).double_value (); @@ -90,16 +76,14 @@ int m = b.columns (); // m: number of inputs int p = c.rows (); // p: number of outputs - int lda = max (1, n); - int ldb = max (1, n); - int ldc = max (1, p); + int lda = max (1, a.rows ()); + int ldb = max (1, b.rows ()); + int ldc = max (1, c.rows ()); // arguments out int ns; - dim_vector dv (1); - dv(0) = n; - NDArray hsv (dv); + ColumnVector hsv (n); // workspace int ldwork = max (1, n*(max (n, m, p) + 5) + n*(n+1)/2); Modified: trunk/octave-forge/main/control/src/slab13bd.cc =================================================================== --- trunk/octave-forge/main/control/src/slab13bd.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slab13bd.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: November 2009 -Version: 0.2 +Version: 0.3 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -45,29 +46,6 @@ int& IWARN, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int max (int a, int b, int c) -{ - int d = max (a, b); - - return max (c, d); -} - -int min (int a, int b) -{ - if (a < b) - return a; - else - return b; -} DEFUN_DLD (slab13bd, args, nargout, "Slicot AB13BD Release 5.0") { @@ -84,10 +62,10 @@ char dico; char jobn = 'H'; - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray c = args(2).array_value (); - NDArray d = args(3).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); int digital = args(4).int_value (); if (digital == 0) Modified: trunk/octave-forge/main/control/src/slab13dd.cc =================================================================== --- trunk/octave-forge/main/control/src/slab13dd.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slab13dd.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,13 +23,14 @@ Author: Lukas Reichlin <luk...@gm...> Created: November 2009 -Version: 0.2 +Version: 0.3 */ #include <octave/oct.h> #include <f77-fcn.h> #include <complex> +#include "common.cc" extern "C" { @@ -49,22 +50,6 @@ Complex* CWORK, int& LCWORK, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int min (int a, int b) -{ - if (a < b) - return a; - else - return b; -} DEFUN_DLD (slab13dd, args, nargout, "Slicot AB13DD Release 5.0") { @@ -83,10 +68,10 @@ char equil = 'N'; char jobd = 'D'; - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray c = args(2).array_value (); - NDArray d = args(3).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); double* e = 0; int digital = args(4).int_value (); double tol = args(5).double_value (); @@ -106,10 +91,8 @@ int ldd = max (1, d.rows ()); int lde = 1; - dim_vector dv (1); - dv(0) = 2; - NDArray fpeak (dv); - NDArray gpeak (dv); + ColumnVector fpeak (2); + ColumnVector gpeak (2); fpeak(0) = 0; fpeak(1) = 1; Modified: trunk/octave-forge/main/control/src/slsb01bd.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb01bd.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slsb01bd.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: November 2009 -Version: 0.2.1 +Version: 0.3 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -46,22 +47,6 @@ double* DWORK, int& LDWORK, int& IWARN, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int max (int a, int b, int c, int d) -{ - int e = max (a, b); - int f = max (c, d); - - return max (e, f); -} DEFUN_DLD (slsb01bd, args, nargout, "Slicot SB01BD Release 5.0") { @@ -77,10 +62,10 @@ // arguments in char dico; - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray wr = args(2).array_value (); - NDArray wi = args(3).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + ColumnVector wr = args(2).column_vector_value (); + ColumnVector wi = args(3).column_vector_value (); int digital = args(4).int_value (); double alpha = args(5).double_value (); double tol = args(6).double_value (); @@ -104,10 +89,7 @@ int nap; int nup; - dim_vector dv (2); - dv(0) = ldf; - dv(1) = n; - NDArray f (dv); + Matrix f (ldf, n); OCTAVE_LOCAL_BUFFER (double, z, ldz*n); Modified: trunk/octave-forge/main/control/src/slsb02od.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb02od.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slsb02od.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: February 2010 -Version: 0.1 +Version: 0.2 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -55,29 +56,6 @@ bool* BWORK, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int max (int a, int b, int c) -{ - int d = max (a, b); - - return max (c, d); -} - -int max (int a, int b, int c, int d) -{ - int e = max (a, b); - int f = max (c, d); - - return max (e, f); -} DEFUN_DLD (slsb02od, args, nargout, "Slicot SB02OD Release 5.0") { @@ -98,11 +76,11 @@ char jobl; char sort = 'S'; - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray q = args(2).array_value (); - NDArray r = args(3).array_value (); - NDArray l = args(4).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix q = args(2).matrix_value (); + Matrix r = args(3).matrix_value (); + Matrix l = args(4).matrix_value (); int digital = args(5).int_value (); int ijobl = args(6).int_value (); @@ -130,10 +108,7 @@ double rcond; int ldx = max (1, n); - dim_vector dv (2); - dv(0) = ldx; - dv(1) = n; - NDArray x (dv); + Matrix x (ldx, n); // unused output arguments OCTAVE_LOCAL_BUFFER (double, alfar, 2*n); Modified: trunk/octave-forge/main/control/src/slsb03md.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb03md.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slsb03md.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: December 2009 -Version: 0.2 +Version: 0.3 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -46,14 +47,6 @@ double* DWORK, int& LDWORK, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} DEFUN_DLD (slsb03md, args, nargout, "Slicot SB03MD Release 5.0") { @@ -72,8 +65,8 @@ char fact = 'N'; char trana = 'T'; - NDArray a = args(0).array_value (); - NDArray c = args(1).array_value (); + Matrix a = args(0).matrix_value (); + Matrix c = args(1).matrix_value (); int dt = args(2).int_value (); if (dt == 0) @@ -92,17 +85,10 @@ double sep = 0; double ferr = 0; - dim_vector dv_u (2); - dv_u(0) = ldu; - dv_u(1) = n; + Matrix u (ldu, n); + ColumnVector wr (n); + ColumnVector wi (n); - dim_vector dv (1); - dv(0) = n; - - NDArray u (dv_u); - NDArray wr (dv); - NDArray wi (dv); - // workspace int* iwork = 0; // not referenced because job = X Modified: trunk/octave-forge/main/control/src/slsb04md.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb04md.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slsb04md.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: January 2010 -Version: 0.1 +Version: 0.2 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -42,22 +43,6 @@ double* DWORK, int& LDWORK, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int max (int a, int b, int c, int d) -{ - int e = max (a, b); - int f = max (c, d); - - return max (e, f); -} DEFUN_DLD (slsb04md, args, nargout, "Slicot SB04MD Release 5.0") { @@ -71,9 +56,9 @@ else { // arguments in - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray c = args(2).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); int n = a.rows (); int m = b.rows (); @@ -84,12 +69,8 @@ int ldz = max (1, m); // arguments out - dim_vector dv (2); - dv(0) = ldz; - dv(1) = m; + Matrix z (ldz, m); - NDArray z (dv); - // workspace int ldwork = max (1, 2*n*n + 8*n, 5*m, n + m); Modified: trunk/octave-forge/main/control/src/slsb04qd.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb04qd.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slsb04qd.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: January 2010 -Version: 0.1 +Version: 0.2 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -43,22 +44,6 @@ int& INFO); } -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int max (int a, int b, int c, int d) -{ - int e = max (a, b); - int f = max (c, d); - - return max (e, f); -} - DEFUN_DLD (slsb04qd, args, nargout, "Slicot SB04QD Release 5.0") { int nargin = args.length (); @@ -71,9 +56,9 @@ else { // arguments in - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray c = args(2).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); int n = a.rows (); int m = b.rows (); @@ -84,12 +69,8 @@ int ldz = max (1, m); // arguments out - dim_vector dv (2); - dv(0) = ldz; - dv(1) = m; + Matrix z (ldz, m); - NDArray z (dv); - // workspace int ldwork = max (1, 2*n*n + 9*n, 5*m, n + m); Modified: trunk/octave-forge/main/control/src/slsb10dd.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb10dd.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slsb10dd.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: December 2009 -Version: 0.2 +Version: 0.3 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -53,29 +54,6 @@ bool* BWORK, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int max (int a, int b, int c) -{ - int d = max (a, b); - - return max (c, d); -} - -int max (int a, int b, int c, int d) -{ - int e = max (a, b); - int f = max (c, d); - - return max (e, f); -} DEFUN_DLD (slsb10dd, args, nargout, "Slicot SB10DD Release 5.0") { @@ -89,10 +67,10 @@ else { // arguments in - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray c = args(2).array_value (); - NDArray d = args(3).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); int ncon = args(4).int_value (); int nmeas = args(5).int_value (); @@ -118,41 +96,14 @@ double tol = 0; // arguments out - dim_vector dv_ak (2); - dv_ak(0) = ldak; - dv_ak(1) = n; + Matrix ak (ldak, n); + Matrix bk (ldbk, nmeas); + Matrix ck (ldck, n); + Matrix dk (lddk, nmeas); + Matrix x (ldx, n); + Matrix z (ldz, n); + ColumnVector rcond (8); - dim_vector dv_bk (2); - dv_bk(0) = ldbk; - dv_bk(1) = nmeas; - - dim_vector dv_ck (2); - dv_ck(0) = ldck; - dv_ck(1) = n; - - dim_vector dv_dk (2); - dv_dk(0) = lddk; - dv_dk(1) = nmeas; - - dim_vector dv_x (2); - dv_x(0) = ldx; - dv_x(1) = n; - - dim_vector dv_z (2); - dv_z(0) = ldz; - dv_z(1) = n; - - dim_vector dv (1); - dv(0) = 8; - - NDArray ak (dv_ak); - NDArray bk (dv_bk); - NDArray ck (dv_ck); - NDArray dk (dv_dk); - NDArray x (dv_x); - NDArray z (dv_z); - NDArray rcond (dv); - // workspace int m2 = ncon; int m1 = m - m2; Modified: trunk/octave-forge/main/control/src/slsb10ed.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb10ed.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slsb10ed.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: November 2009 -Version: 0.2 +Version: 0.3 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -51,29 +52,6 @@ int& INFO); } -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int max (int a, int b, int c) -{ - int d = max (a, b); - - return max (c, d); -} - -int max (int a, int b, int c, int d) -{ - int e = max (a, b); - int f = max (c, d); - - return max (e, f); -} - DEFUN_DLD (slsb10ed, args, nargout, "Slicot SB10ED Release 5.0") { int nargin = args.length (); @@ -86,10 +64,10 @@ else { // arguments in - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray c = args(2).array_value (); - NDArray d = args(3).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); int ncon = args(4).int_value (); int nmeas = args(5).int_value (); @@ -110,32 +88,13 @@ double tol = 0; - // arguments out - dim_vector dv_ak (2); - dv_ak(0) = ldak; - dv_ak(1) = n; + // arguments out + Matrix ak (ldak, n); + Matrix bk (ldbk, nmeas); + Matrix ck (ldck, n); + Matrix dk (lddk, nmeas); + ColumnVector rcond (7); - dim_vector dv_bk (2); - dv_bk(0) = ldbk; - dv_bk(1) = nmeas; - - dim_vector dv_ck (2); - dv_ck(0) = ldck; - dv_ck(1) = n; - - dim_vector dv_dk (2); - dv_dk(0) = lddk; - dv_dk(1) = nmeas; - - dim_vector dv (1); - dv(0) = 7; - - NDArray ak (dv_ak); - NDArray bk (dv_bk); - NDArray ck (dv_ck); - NDArray dk (dv_dk); - NDArray rcond (dv); - // workspace int m2 = ncon; int m1 = m - m2; Modified: trunk/octave-forge/main/control/src/slsb10fd.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb10fd.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slsb10fd.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: December 2009 -Version: 0.2 +Version: 0.3 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -51,29 +52,6 @@ bool* BWORK, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int max (int a, int b, int c) -{ - int d = max (a, b); - - return max (c, d); -} - -int max (int a, int b, int c, int d) -{ - int e = max (a, b); - int f = max (c, d); - - return max (e, f); -} DEFUN_DLD (slsb10fd, args, nargout, "Slicot SB10FD Release 5.0") { @@ -87,10 +65,10 @@ else { // arguments in - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray c = args(2).array_value (); - NDArray d = args(3).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); int ncon = args(4).int_value (); int nmeas = args(5).int_value (); @@ -113,31 +91,12 @@ double tol = 0; // arguments out - dim_vector dv_ak (2); - dv_ak(0) = ldak; - dv_ak(1) = n; + Matrix ak (ldak, n); + Matrix bk (ldbk, nmeas); + Matrix ck (ldck, n); + Matrix dk (lddk, nmeas); + ColumnVector rcond (4); - dim_vector dv_bk (2); - dv_bk(0) = ldbk; - dv_bk(1) = nmeas; - - dim_vector dv_ck (2); - dv_ck(0) = ldck; - dv_ck(1) = n; - - dim_vector dv_dk (2); - dv_dk(0) = lddk; - dv_dk(1) = nmeas; - - dim_vector dv (1); - dv(0) = 4; - - NDArray ak (dv_ak); - NDArray bk (dv_bk); - NDArray ck (dv_ck); - NDArray dk (dv_dk); - NDArray rcond (dv); - // workspace int m2 = ncon; Modified: trunk/octave-forge/main/control/src/slsb10hd.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb10hd.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slsb10hd.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: November 2009 -Version: 0.2 +Version: 0.3 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -50,29 +51,6 @@ bool* BWORK, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} - -int max (int a, int b, int c) -{ - int d = max (a, b); - - return max (c, d); -} - -int max (int a, int b, int c, int d) -{ - int e = max (a, b); - int f = max (c, d); - - return max (e, f); -} DEFUN_DLD (slsb10hd, args, nargout, "Slicot SB10HD Release 5.0") { @@ -86,10 +64,10 @@ else { // arguments in - NDArray a = args(0).array_value (); - NDArray b = args(1).array_value (); - NDArray c = args(2).array_value (); - NDArray d = args(3).array_value (); + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); int ncon = args(4).int_value (); int nmeas = args(5).int_value (); @@ -111,31 +89,12 @@ double tol = 0; // arguments out - dim_vector dv_ak (2); - dv_ak(0) = ldak; - dv_ak(1) = n; + Matrix ak (ldak, n); + Matrix bk (ldbk, nmeas); + Matrix ck (ldck, n); + Matrix dk (lddk, nmeas); + ColumnVector rcond (4); - dim_vector dv_bk (2); - dv_bk(0) = ldbk; - dv_bk(1) = nmeas; - - dim_vector dv_ck (2); - dv_ck(0) = ldck; - dv_ck(1) = n; - - dim_vector dv_dk (2); - dv_dk(0) = lddk; - dv_dk(1) = nmeas; - - dim_vector dv (1); - dv(0) = 4; - - NDArray ak (dv_ak); - NDArray bk (dv_bk); - NDArray ck (dv_ck); - NDArray dk (dv_dk); - NDArray rcond (dv); - // workspace int m2 = ncon; int m1 = m - m2; Modified: trunk/octave-forge/main/control/src/slsg03ad.cc =================================================================== --- trunk/octave-forge/main/control/src/slsg03ad.cc 2010-04-17 21:39:03 UTC (rev 7207) +++ trunk/octave-forge/main/control/src/slsg03ad.cc 2010-04-18 07:33:41 UTC (rev 7208) @@ -23,12 +23,13 @@ Author: Lukas Reichlin <luk...@gm...> Created: January 2010 -Version: 0.1 +Version: 0.2 */ #include <octave/oct.h> #include <f77-fcn.h> +#include "common.cc" extern "C" { @@ -50,14 +51,6 @@ double* DWORK, int& LDWORK, int& INFO); } - -int max (int a, int b) -{ - if (a > b) - return a; - else - return b; -} DEFUN_DLD (slsg03ad, args, nargout, "Slicot SG03AD Release 5.0") { @@ -77,9 +70,9 @@ char trans = 'T'; char uplo = 'U'; // ?!? - NDArray a = args(0).array_value (); - NDArray e = args(1).array_value (); - NDArray x = args(2).array_value (); + Matrix a = args(0).matrix_value (); + Matrix e = args(1).matrix_value (); + Matrix x = args(2).matrix_value (); int dt = args(3).int_value (); if (dt == 0) @@ -100,23 +93,12 @@ double sep = 0; double ferr = 0; - dim_vector dv_q (2); - dv_q(0) = ldq; - dv_q(1) = n; + Matrix q (ldq, n); + Matrix z (ldz, n); + ColumnVector alphar (n); + ColumnVector alphai (n); + ColumnVector beta (n); - dim_vector dv_z (2); - dv_z(0) = ldz; - dv_z(1) = n; - - dim_vector dv (1); - dv(0) = n; - - NDArray q (dv_q); - NDArray z (dv_z); - NDArray alphar (dv); - NDArray alphai (dv); - NDArray beta (dv); - // workspace int* iwork = 0; // not referenced because job = X This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-04-20 15:33:04
|
Revision: 7212 http://octave.svn.sourceforge.net/octave/?rev=7212&view=rev Author: paramaniac Date: 2010-04-20 15:32:57 +0000 (Tue, 20 Apr 2010) Log Message: ----------- control: style fixes Modified Paths: -------------- trunk/octave-forge/main/control/inst/@lti/norm.m trunk/octave-forge/main/control/inst/__freqbounds__.m trunk/octave-forge/main/control/inst/__timeresp__.m trunk/octave-forge/main/control/inst/hsvd.m trunk/octave-forge/main/control/inst/lsim.m trunk/octave-forge/main/control/inst/place.m trunk/octave-forge/main/control/src/slab13ad.cc trunk/octave-forge/main/control/src/slab13bd.cc trunk/octave-forge/main/control/src/slab13dd.cc trunk/octave-forge/main/control/src/slsb01bd.cc trunk/octave-forge/main/control/src/slsb02od.cc trunk/octave-forge/main/control/src/slsb03md.cc trunk/octave-forge/main/control/src/slsg03ad.cc Modified: trunk/octave-forge/main/control/inst/@lti/norm.m =================================================================== --- trunk/octave-forge/main/control/inst/@lti/norm.m 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/inst/@lti/norm.m 2010-04-20 15:32:57 UTC (rev 7212) @@ -66,12 +66,12 @@ if (isstable (sys)) [a, b, c, d] = ssdata (sys); - digital = ! isct (sys); + discrete = ! isct (sys); - if (! digital && ! all (all (d == 0))) # continuous and non-zero feedthrough + if (! discrete && ! all (all (d == 0))) # continuous and non-zero feedthrough gain = inf; else - [gain, iwarn] = slab13bd (a, b, c, d, digital); + [gain, iwarn] = slab13bd (a, b, c, d, discrete); if (iwarn) warning ("lti: norm: slab13bd: iwarn = %d", iwarn); @@ -87,13 +87,13 @@ function [gain, wpeak] = linfnorm (sys, tol = 0.01) [a, b, c, d, tsam] = ssdata (sys); - digital = ! isct (sys); + discrete = ! isct (sys); tol = max (tol, 100*eps); - [fpeak, gpeak] = slab13dd (a, b, c, d, digital, tol); + [fpeak, gpeak] = slab13dd (a, b, c, d, discrete, tol); if (fpeak(2) > 0) - if (digital) + if (discrete) wpeak = fpeak(1) / tsam; else wpeak = fpeak(1); Modified: trunk/octave-forge/main/control/inst/__freqbounds__.m =================================================================== --- trunk/octave-forge/main/control/inst/__freqbounds__.m 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/inst/__freqbounds__.m 2010-04-20 15:32:57 UTC (rev 7212) @@ -42,14 +42,14 @@ zer = zero (sys); pol = pole (sys); tsam = get (sys, "tsam"); - digital = (tsam > 0); # static gains (tsam = -1) are continuous + discrete = (tsam > 0); # static gains (tsam = -1) are continuous ## make sure zer, pol are row vectors pol = pol(:).'; zer = zer(:).'; ## check for natural frequencies away from omega = 0 - if (digital) + if (discrete) ## The 2nd conditions prevents log(0) in the next log command iiz = find (abs(zer-1) > norm(zer)*eps && abs(zer) > norm(zer)*eps); iip = find (abs(pol-1) > norm(pol)*eps && abs(pol) > norm(pol)*eps); @@ -114,8 +114,8 @@ error ("freqbounds: second argument invalid"); endswitch - ## run digital frequency all the way to pi - if (digital) + ## run discrete frequency all the way to pi + if (discrete) dec_max = log10 (pi/tsam); endif Modified: trunk/octave-forge/main/control/inst/__timeresp__.m =================================================================== --- trunk/octave-forge/main/control/inst/__timeresp__.m 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/inst/__timeresp__.m 2010-04-20 15:32:57 UTC (rev 7212) @@ -35,19 +35,19 @@ [A, B, C, D, tsam] = ssdata (sys); - digital = ! isct (sys); # static gains are treated as analog systems + discrete = ! isct (sys); # static gains are treated as analog systems - if (digital) + if (discrete) if (! isempty (dt)) - warning ("timeresp: argument dt has no effect on sampling time of digital system"); + warning ("timeresp: argument dt has no effect on sampling time of discrete system"); endif dt = tsam; endif - [tfinal, dt] = __simhorizon__ (A, digital, tfinal, dt); + [tfinal, dt] = __simhorizon__ (A, discrete, tfinal, dt); - if (! digital) + if (! discrete) sys = c2d (sys, dt, "zoh"); endif @@ -119,7 +119,7 @@ u = zeros (m, 1); u(j) = 1; - if (digital) + if (discrete) x = zeros (n, 1); # zero by definition y(1, :, j) = D * u / dt; x_arr(1, :, j) = x; @@ -139,7 +139,7 @@ endfor endfor - if (digital) + if (discrete) y *= dt; x_arr *= dt; endif @@ -170,7 +170,7 @@ cols = m; endif - if (digital) # discrete system + if (discrete) # discrete system for k = 1 : p for j = 1 : cols @@ -230,7 +230,7 @@ endfunction -function [tfinal, dt] = __simhorizon__ (A, digital, tfinal, Ts) +function [tfinal, dt] = __simhorizon__ (A, discrete, tfinal, Ts) ## code based on __stepimp__.m of Kai P. Mueller and A. Scottedward Hodel @@ -243,7 +243,7 @@ n = rows (A); eigw = eig (A); - if (digital) + if (discrete) ## perform bilinear transformation on poles in z for k = 1 : n pol = eigw(k); @@ -269,14 +269,14 @@ tfinal = T_DEF; endif - if (! digital) + if (! discrete) dt = tfinal / N_DEF; endif else eigw = eigw(find (eigw)); eigw_max = max (abs (eigw)); - if (! digital) + if (! discrete) dt = 0.2 * pi / eigw_max; endif @@ -289,7 +289,7 @@ tfinal = yy * ceil (tfinal / yy); endif - if (! digital) + if (! discrete) N = tfinal / dt; if (N < N_MIN) Modified: trunk/octave-forge/main/control/inst/hsvd.m =================================================================== --- trunk/octave-forge/main/control/inst/hsvd.m 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/inst/hsvd.m 2010-04-20 15:32:57 UTC (rev 7212) @@ -44,15 +44,15 @@ [a, b, c] = ssdata (sys); - digital = ! isct (sys); + discrete = ! isct (sys); - if (digital) + if (discrete) alpha = 1 - val; else alpha = - val; endif - [hsv, ns] = slab13ad (a, b, c, digital, alpha); + [hsv, ns] = slab13ad (a, b, c, discrete, alpha); idx = 1 : ns; hsv = hsv (idx); Modified: trunk/octave-forge/main/control/inst/lsim.m =================================================================== --- trunk/octave-forge/main/control/inst/lsim.m 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/inst/lsim.m 2010-04-20 15:32:57 UTC (rev 7212) @@ -39,12 +39,12 @@ [A, B, C, D, tsam] = ssdata (sys); - digital = ! isct (sys); # static gains are treated as analog systems + discrete = ! isct (sys); # static gains are treated as analog systems urows = rows (u); ucols = columns (u); - if (digital) # discrete system + if (discrete) # discrete system if (isempty (t)) dt = tsam; tfinal = tsam * (urows - 1); @@ -52,7 +52,7 @@ dt = tsam; tfinal = t; else - warning ("lsim: spacing of time vector has no effect on sampling time of digital system"); + warning ("lsim: spacing of time vector has no effect on sampling time of discrete system"); dt = tsam; tfinal = t(end); endif @@ -116,7 +116,7 @@ outname = __markemptynames__ (outname); endif - if (digital) # discrete system + if (discrete) # discrete system for k = 1 : p subplot (p, 1, k); stairs (t, y(:, k)); Modified: trunk/octave-forge/main/control/inst/place.m =================================================================== --- trunk/octave-forge/main/control/inst/place.m 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/inst/place.m 2010-04-20 15:32:57 UTC (rev 7212) @@ -84,7 +84,7 @@ p = b; sys = a; [a, b] = ssdata (sys); - digital = ! isct (sys); # treat tsam = -1 as continuous system + discrete = ! isct (sys); # treat tsam = -1 as continuous system endif else # place (a, b, p), place (a, b, p, alpha), place (a, b, p, alpha, tol) if (nargin < 3) # nargin > 5 already tested @@ -93,7 +93,7 @@ if (! isnumeric (a) || ! isnumeric (b) || ! issquare (a) || rows (a) != rows (b)) error ("place: matrices a and b not conformal"); endif - digital = 0; # assume continuous system + discrete = 0; # assume continuous system endif endif @@ -114,7 +114,7 @@ endif if (isempty (alpha)) - if (digital) + if (discrete) alpha = 0; else alpha = - norm (a, inf); @@ -125,7 +125,7 @@ tol = 0; endif - [f, iwarn, nfp, nap, nup] = slsb01bd (a, b, wr, wi, digital, alpha, tol); + [f, iwarn, nfp, nap, nup] = slsb01bd (a, b, wr, wi, discrete, alpha, tol); f = -f; # A + B*F --> A - B*F if (iwarn) Modified: trunk/octave-forge/main/control/src/slab13ad.cc =================================================================== --- trunk/octave-forge/main/control/src/slab13ad.cc 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/src/slab13ad.cc 2010-04-20 15:32:57 UTC (rev 7212) @@ -64,10 +64,10 @@ Matrix a = args(0).matrix_value (); Matrix b = args(1).matrix_value (); Matrix c = args(2).matrix_value (); - int digital = args(3).int_value (); + int discrete = args(3).int_value (); double alpha = args(4).double_value (); - if (digital == 0) + if (discrete == 0) dico = 'C'; else dico = 'D'; Modified: trunk/octave-forge/main/control/src/slab13bd.cc =================================================================== --- trunk/octave-forge/main/control/src/slab13bd.cc 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/src/slab13bd.cc 2010-04-20 15:32:57 UTC (rev 7212) @@ -66,9 +66,9 @@ Matrix b = args(1).matrix_value (); Matrix c = args(2).matrix_value (); Matrix d = args(3).matrix_value (); - int digital = args(4).int_value (); + int discrete = args(4).int_value (); - if (digital == 0) + if (discrete == 0) dico = 'C'; else dico = 'D'; Modified: trunk/octave-forge/main/control/src/slab13dd.cc =================================================================== --- trunk/octave-forge/main/control/src/slab13dd.cc 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/src/slab13dd.cc 2010-04-20 15:32:57 UTC (rev 7212) @@ -73,10 +73,10 @@ Matrix c = args(2).matrix_value (); Matrix d = args(3).matrix_value (); double* e = 0; - int digital = args(4).int_value (); + int discrete = args(4).int_value (); double tol = args(5).double_value (); - if (digital == 0) + if (discrete == 0) dico = 'C'; else dico = 'D'; Modified: trunk/octave-forge/main/control/src/slsb01bd.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb01bd.cc 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/src/slsb01bd.cc 2010-04-20 15:32:57 UTC (rev 7212) @@ -66,11 +66,11 @@ Matrix b = args(1).matrix_value (); ColumnVector wr = args(2).column_vector_value (); ColumnVector wi = args(3).column_vector_value (); - int digital = args(4).int_value (); + int discrete = args(4).int_value (); double alpha = args(5).double_value (); double tol = args(6).double_value (); - if (digital == 1) + if (discrete == 1) dico = 'D'; else dico = 'C'; Modified: trunk/octave-forge/main/control/src/slsb02od.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb02od.cc 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/src/slsb02od.cc 2010-04-20 15:32:57 UTC (rev 7212) @@ -81,10 +81,10 @@ Matrix q = args(2).matrix_value (); Matrix r = args(3).matrix_value (); Matrix l = args(4).matrix_value (); - int digital = args(5).int_value (); + int discrete = args(5).int_value (); int ijobl = args(6).int_value (); - if (digital == 0) + if (discrete == 0) dico = 'C'; else dico = 'D'; Modified: trunk/octave-forge/main/control/src/slsb03md.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb03md.cc 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/src/slsb03md.cc 2010-04-20 15:32:57 UTC (rev 7212) @@ -67,9 +67,9 @@ Matrix a = args(0).matrix_value (); Matrix c = args(1).matrix_value (); - int dt = args(2).int_value (); + int discrete = args(2).int_value (); - if (dt == 0) + if (discrete == 0) dico = 'C'; else dico = 'D'; Modified: trunk/octave-forge/main/control/src/slsg03ad.cc =================================================================== --- trunk/octave-forge/main/control/src/slsg03ad.cc 2010-04-20 03:57:28 UTC (rev 7211) +++ trunk/octave-forge/main/control/src/slsg03ad.cc 2010-04-20 15:32:57 UTC (rev 7212) @@ -73,9 +73,9 @@ Matrix a = args(0).matrix_value (); Matrix e = args(1).matrix_value (); Matrix x = args(2).matrix_value (); - int dt = args(3).int_value (); + int discrete = args(3).int_value (); - if (dt == 0) + if (discrete == 0) dico = 'C'; else dico = 'D'; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-08-18 12:27:21
|
Revision: 7552 http://octave.svn.sourceforge.net/octave/?rev=7552&view=rev Author: paramaniac Date: 2010-08-18 12:27:15 +0000 (Wed, 18 Aug 2010) Log Message: ----------- control: doc test_control.m Modified Paths: -------------- trunk/octave-forge/main/control/DESCRIPTION trunk/octave-forge/main/control/INDEX trunk/octave-forge/main/control/inst/test_control.m Modified: trunk/octave-forge/main/control/DESCRIPTION =================================================================== --- trunk/octave-forge/main/control/DESCRIPTION 2010-08-18 12:06:26 UTC (rev 7551) +++ trunk/octave-forge/main/control/DESCRIPTION 2010-08-18 12:27:15 UTC (rev 7552) @@ -1,6 +1,6 @@ Name: Control -Version: 0.3.1 -Date: 2010-04-18 +Version: 0.3.2 +Date: 2010-08-18 Author: Lukas Reichlin <luk...@gm...> Maintainer: Lukas Reichlin <luk...@gm...> Title: LTI Syncope Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2010-08-18 12:06:26 UTC (rev 7551) +++ trunk/octave-forge/main/control/INDEX 2010-08-18 12:27:15 UTC (rev 7552) @@ -96,4 +96,6 @@ Octave-only isctrb isobsv - @lti/mconnect \ No newline at end of file + @lti/mconnect +Tests + test_control \ No newline at end of file Modified: trunk/octave-forge/main/control/inst/test_control.m =================================================================== --- trunk/octave-forge/main/control/inst/test_control.m 2010-08-18 12:06:26 UTC (rev 7551) +++ trunk/octave-forge/main/control/inst/test_control.m 2010-08-18 12:27:15 UTC (rev 7552) @@ -1,5 +1,29 @@ -## execute all available tests at once +## Copyright (C) 2010 Lukas F. Reichlin +## +## This file is part of LTI Syncope. +## +## LTI Syncope 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. +## +## LTI Syncope 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. If not, see <http://www.gnu.org/licenses/>. +## -*- texinfo -*- +## @deftypefn {Script File} test_control +## Execute all available tests at once. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: May 2010 +## Version: 0.1 + test ltimodels test hinfsyn test h2syn This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-08-29 01:23:36
|
Revision: 7596 http://octave.svn.sourceforge.net/octave/?rev=7596&view=rev Author: paramaniac Date: 2010-08-29 01:23:30 +0000 (Sun, 29 Aug 2010) Log Message: ----------- control: clarifications Modified Paths: -------------- trunk/octave-forge/main/control/DESCRIPTION trunk/octave-forge/main/control/inst/gensig.m Modified: trunk/octave-forge/main/control/DESCRIPTION =================================================================== --- trunk/octave-forge/main/control/DESCRIPTION 2010-08-28 10:01:43 UTC (rev 7595) +++ trunk/octave-forge/main/control/DESCRIPTION 2010-08-29 01:23:30 UTC (rev 7596) @@ -3,8 +3,8 @@ Date: 2010-08-25 Author: Lukas Reichlin <luk...@gm...> Maintainer: Lukas Reichlin <luk...@gm...> -Title: Computer Aided Control Systems Design -Description: Object-oriented control tools for GNU Octave +Title: Control Systems +Description: Computer Aided Control Systems Design package for GNU Octave Depends: octave (>= 3.2.3) Autoload: yes License: GPL version 3 or later Modified: trunk/octave-forge/main/control/inst/gensig.m =================================================================== --- trunk/octave-forge/main/control/inst/gensig.m 2010-08-28 10:01:43 UTC (rev 7595) +++ trunk/octave-forge/main/control/inst/gensig.m 2010-08-29 01:23:30 UTC (rev 7596) @@ -74,17 +74,15 @@ t = (0 : tsam : tfinal)'; - sigtype = lower (sigtype); - - switch (sigtype(1:2)) + switch (lower (sigtype(1:2))) case "si" u = sin (2*pi/tau * t); case "co" u = cos (2*pi/tau * t); case "sq" - u = (rem (t, tau) >= tau/2); + u = rem (t, tau) >= tau/2; case "pu" - u = (rem (t, tau) < (1 - 1000*eps) * tsam); + u = rem (t, tau) < (1 - 1000*eps) * tsam; otherwise error ("gensig: invalid signal type"); endswitch This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-09-07 08:41:04
|
Revision: 7673 http://octave.svn.sourceforge.net/octave/?rev=7673&view=rev Author: paramaniac Date: 2010-09-07 08:40:56 +0000 (Tue, 07 Sep 2010) Log Message: ----------- control: use SLICOT for ss:minreal Modified Paths: -------------- trunk/octave-forge/main/control/DESCRIPTION trunk/octave-forge/main/control/inst/@lti/minreal.m trunk/octave-forge/main/control/inst/@ss/__minreal__.m trunk/octave-forge/main/control/inst/ltimodels.m trunk/octave-forge/main/control/src/Makefile Added Paths: ----------- trunk/octave-forge/main/control/src/AB07MD.f trunk/octave-forge/main/control/src/TB01PD.f trunk/octave-forge/main/control/src/TB01UD.f trunk/octave-forge/main/control/src/TB01XD.f trunk/octave-forge/main/control/src/sltb01pd.cc Modified: trunk/octave-forge/main/control/DESCRIPTION =================================================================== --- trunk/octave-forge/main/control/DESCRIPTION 2010-09-07 06:25:09 UTC (rev 7672) +++ trunk/octave-forge/main/control/DESCRIPTION 2010-09-07 08:40:56 UTC (rev 7673) @@ -1,6 +1,6 @@ Name: Control Version: 2.0.0 -Date: 2010-08-25 +Date: 2010-09-07 Author: Lukas Reichlin <luk...@gm...> Maintainer: Lukas Reichlin <luk...@gm...> Title: Control Systems Modified: trunk/octave-forge/main/control/inst/@lti/minreal.m =================================================================== --- trunk/octave-forge/main/control/inst/@lti/minreal.m 2010-09-07 06:25:09 UTC (rev 7672) +++ trunk/octave-forge/main/control/inst/@lti/minreal.m 2010-09-07 08:40:56 UTC (rev 7673) @@ -1,4 +1,4 @@ -## Copyright (C) 2009 Lukas F. Reichlin +## Copyright (C) 2009 - 2010 Lukas F. Reichlin ## ## This file is part of LTI Syncope. ## @@ -23,10 +23,18 @@ ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 -## Version: 0.1 +## Version: 0.2 function sys = minreal (sys, tol = "def") + if (nargin > 2) # nargin == 0 not possible because minreal is inside @lti + print_usage (); + endif + + if (isa (tol, "lti") || (tol != "def" && ! isreal (tol) && ! isscalar (tol))) + error ("minreal: second argument must be a real scalar"); + endif + sys = __minreal__ (sys, tol); endfunction \ No newline at end of file Modified: trunk/octave-forge/main/control/inst/@ss/__minreal__.m =================================================================== --- trunk/octave-forge/main/control/inst/@ss/__minreal__.m 2010-09-07 06:25:09 UTC (rev 7672) +++ trunk/octave-forge/main/control/inst/@ss/__minreal__.m 2010-09-07 08:40:56 UTC (rev 7673) @@ -1,4 +1,4 @@ -## Copyright (C) 2009 Lukas F. Reichlin +## Copyright (C) 2009 - 2010 Lukas F. Reichlin ## ## This file is part of LTI Syncope. ## @@ -17,55 +17,27 @@ ## -*- texinfo -*- ## Minimal realization of SS models. The physical meaning of states is lost. +## Uses SLICOT TB01PD by courtesy of NICONET e.V. <http://www.slicot.org> -## Algorithm based on sysmin by A. Scottedward Hodel ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 -## Version: 0.1 +## Version: 0.2 function retsys = __minreal__ (sys, tol) - A = sys.a; - B = sys.b; - C = sys.c; - - if (! isempty (A)) - if (tol == "def") - [cflg, Uc] = isctrb (A, B); - else - [cflg, Uc] = isctrb (A, B, tol); - endif - - if (! cflg) - if (! isempty (Uc)) - A = Uc.' * A * Uc; - B = Uc.' * B; - C = C * Uc; - else - A = B = C = []; - endif - endif + if (tol == "def") + tol = 0; + elseif (tol > 1) + error ("ss: minreal: require tol <= 1"); endif - if (! isempty (A)) - if (tol == "def") - [oflg, Uo] = isobsv (A, C); - else - [oflg, Uo] = isobsv (A, C, tol); - endif + [a, b, c, nr] = sltb01pd (sys.a, sys.b, sys.c, tol); - if (! oflg) - if (! isempty (Uo)) - A = Uo.' * A * Uo; - B = Uo.' * B; - C = C * Uo; - else - A = B = C = []; - endif - endif - endif + a = a(1:nr, 1:nr); + b = b(1:nr, :); + c = c(:, 1:nr); - retsys = ss (A, B, C, sys.d); + retsys = ss (a, b, c, sys.d); retsys.lti = sys.lti; # retain i/o names and tsam endfunction Modified: trunk/octave-forge/main/control/inst/ltimodels.m =================================================================== --- trunk/octave-forge/main/control/inst/ltimodels.m 2010-09-07 06:25:09 UTC (rev 7672) +++ trunk/octave-forge/main/control/inst/ltimodels.m 2010-09-07 08:40:56 UTC (rev 7673) @@ -251,13 +251,13 @@ %! 23.0 5.0 7.0 14.0 16.0 %! 4.0 6.0 13.0 20.0 22.0 %! 10.0 12.0 19.0 21.0 3.0 -%! 11.0 18.0 25.0 2.0 9.0]; +%! 11.0 18.0 25.0 2.0 9.0 ]; %! %! B = [ -1.0 -4.0 %! 4.0 9.0 %! -9.0 -16.0 %! 16.0 25.0 -%! -25.0 -36.0]; +%! -25.0 -36.0 ]; %! %! tol = 0; %! @@ -270,13 +270,50 @@ %! 4.4741 -12.5544 5.3509 5.9403 1.4360 %! 14.4576 7.6855 23.1452 26.3872 -29.9557 %! 0.0000 1.4805 27.4668 22.6564 -0.0072 -%! 0.0000 0.0000 -30.4822 0.6745 18.8680]; +%! 0.0000 0.0000 -30.4822 0.6745 18.8680 ]; %! %! Bce = [ 31.1199 47.6865 %! 3.2480 0.0000 %! 0.0000 0.0000 %! 0.0000 0.0000 -%! 0.0000 0.0000]; +%! 0.0000 0.0000 ]; %! %!assert (Ac, Ace, 1e-4); %!assert (Bc, Bce, 1e-4); + +## ss: minreal (SLICOT TB01PD) +%!shared M, Me +%! A = [ 1.0 2.0 0.0 +%! 4.0 -1.0 0.0 +%! 0.0 0.0 1.0 ]; +%! +%! B = [ 1.0 +%! 0.0 +%! 1.0 ]; +%! +%! C = [ 0.0 1.0 -1.0 +%! 0.0 0.0 1.0 ]; +%! +%! D = zeros (2, 1); +%! +%! sys = ss (A, B, C, D); +%! sysmin = minreal (sys, 0.0); +%! [Ar, Br, Cr, Dr] = ssdata (sysmin); +%! M = [Ar, Br; Cr, Dr]; +%! +%! Ae = [ 1.0000 -1.4142 1.4142 +%! -2.8284 0.0000 1.0000 +%! 2.8284 1.0000 0.0000 ]; +%! +%! Be = [-1.0000 +%! 0.7071 +%! 0.7071 ]; +%! +%! Ce = [ 0.0000 0.0000 -1.4142 +%! 0.0000 0.7071 0.7071 ]; +%! +%! De = zeros (2, 1); +%! +%! Me = [Ae, Be; Ce, De]; +%! +%!assert (M, Me, 1e-4); Added: trunk/octave-forge/main/control/src/AB07MD.f =================================================================== --- trunk/octave-forge/main/control/src/AB07MD.f (rev 0) +++ trunk/octave-forge/main/control/src/AB07MD.f 2010-09-07 08:40:56 UTC (rev 7673) @@ -0,0 +1,224 @@ + SUBROUTINE AB07MD( JOBD, N, M, P, A, LDA, B, LDB, C, LDC, D, LDD, + $ INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2009 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To find the dual of a given state-space representation. +C +C ARGUMENTS +C +C Mode Parameters +C +C JOBD CHARACTER*1 +C Specifies whether or not a non-zero matrix D appears in +C the given state space model: +C = 'D': D is present; +C = 'Z': D is assumed a zero matrix. +C +C Input/Output Parameters +C +C N (input) INTEGER +C The order of the state-space representation. N >= 0. +C +C M (input) INTEGER +C The number of system inputs. M >= 0. +C +C P (input) INTEGER +C The number of system outputs. P >= 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C On entry, the leading N-by-N part of this array must +C contain the original state dynamics matrix A. +C On exit, the leading N-by-N part of this array contains +C the dual state dynamics matrix A'. +C +C LDA INTEGER +C The leading dimension of array A. LDA >= MAX(1,N). +C +C B (input/output) DOUBLE PRECISION array, dimension +C (LDB,MAX(M,P)) +C On entry, the leading N-by-M part of this array must +C contain the original input/state matrix B. +C On exit, the leading N-by-P part of this array contains +C the dual input/state matrix C'. +C +C LDB INTEGER +C The leading dimension of array B. LDB >= MAX(1,N). +C +C C (input/output) DOUBLE PRECISION array, dimension (LDC,N) +C On entry, the leading P-by-N part of this array must +C contain the original state/output matrix C. +C On exit, the leading M-by-N part of this array contains +C the dual state/output matrix B'. +C +C LDC INTEGER +C The leading dimension of array C. +C LDC >= MAX(1,M,P) if N > 0. +C LDC >= 1 if N = 0. +C +C D (input/output) DOUBLE PRECISION array, dimension +C (LDD,MAX(M,P)) +C On entry, if JOBD = 'D', the leading P-by-M part of this +C array must contain the original direct transmission +C matrix D. +C On exit, if JOBD = 'D', the leading M-by-P part of this +C array contains the dual direct transmission matrix D'. +C The array D is not referenced if JOBD = 'Z'. +C +C LDD INTEGER +C The leading dimension of array D. +C LDD >= MAX(1,M,P) if JOBD = 'D'. +C LDD >= 1 if JOBD = 'Z'. +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C < 0: if INFO = -i, the i-th argument had an illegal +C value. +C +C METHOD +C +C If the given state-space representation is the M-input/P-output +C (A,B,C,D), its dual is simply the P-input/M-output (A',C',B',D'). +C +C REFERENCES +C +C None +C +C NUMERICAL ASPECTS +C +C None +C +C CONTRIBUTOR +C +C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Dec. 1996. +C Supersedes Release 2.0 routine AB07AD by T.W.C.Williams, Kingston +C Polytechnic, United Kingdom, March 1982. +C +C REVISIONS +C +C V. Sima, Research Institute for Informatics, Bucharest, Feb. 2004. +C +C KEYWORDS +C +C Dual system, state-space model, state-space representation. +C +C ****************************************************************** +C +C .. Scalar Arguments .. + CHARACTER JOBD + INTEGER INFO, LDA, LDB, LDC, LDD, M, N, P +C .. Array Arguments .. + DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*) +C .. Local Scalars .. + LOGICAL LJOBD + INTEGER J, MINMP, MPLIM +C .. External functions .. + LOGICAL LSAME + EXTERNAL LSAME +C .. External subroutines .. + EXTERNAL DCOPY, DSWAP, XERBLA +C .. Intrinsic Functions .. + INTRINSIC MAX, MIN +C .. Executable Statements .. +C + INFO = 0 + LJOBD = LSAME( JOBD, 'D' ) + MPLIM = MAX( M, P ) + MINMP = MIN( M, P ) +C +C Test the input scalar arguments. +C + IF( .NOT.LJOBD .AND. .NOT.LSAME( JOBD, 'Z' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( M.LT.0 ) THEN + INFO = -3 + ELSE IF( P.LT.0 ) THEN + INFO = -4 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -6 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( ( N.GT.0 .AND. LDC.LT.MAX( 1, MPLIM ) ) .OR. + $ ( N.EQ.0 .AND. LDC.LT.1 ) ) THEN + INFO = -10 + ELSE IF( ( LJOBD .AND. LDD.LT.MAX( 1, MPLIM ) ) .OR. + $ ( .NOT.LJOBD .AND. LDD.LT.1 ) ) THEN + INFO = -12 + END IF +C + IF ( INFO.NE.0 ) THEN +C +C Error return. +C + CALL XERBLA( 'AB07MD', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + IF ( MAX( N, MINMP ).EQ.0 ) + $ RETURN +C + IF ( N.GT.0 ) THEN +C +C Transpose A, if non-scalar. +C + DO 10 J = 1, N - 1 + CALL DSWAP( N-J, A(J+1,J), 1, A(J,J+1), LDA ) + 10 CONTINUE +C +C Replace B by C' and C by B'. +C + DO 20 J = 1, MPLIM + IF ( J.LE.MINMP ) THEN + CALL DSWAP( N, B(1,J), 1, C(J,1), LDC ) + ELSE IF ( J.GT.P ) THEN + CALL DCOPY( N, B(1,J), 1, C(J,1), LDC ) + ELSE + CALL DCOPY( N, C(J,1), LDC, B(1,J), 1 ) + END IF + 20 CONTINUE +C + END IF +C + IF ( LJOBD .AND. MINMP.GT.0 ) THEN +C +C Transpose D, if non-scalar. +C + DO 30 J = 1, MPLIM + IF ( J.LT.MINMP ) THEN + CALL DSWAP( MINMP-J, D(J+1,J), 1, D(J,J+1), LDD ) + ELSE IF ( J.GT.P ) THEN + CALL DCOPY( P, D(1,J), 1, D(J,1), LDD ) + ELSE IF ( J.GT.M ) THEN + CALL DCOPY( M, D(J,1), LDD, D(1,J), 1 ) + END IF + 30 CONTINUE +C + END IF +C + RETURN +C *** Last line of AB07MD *** + END Modified: trunk/octave-forge/main/control/src/Makefile =================================================================== --- trunk/octave-forge/main/control/src/Makefile 2010-09-07 06:25:09 UTC (rev 7672) +++ trunk/octave-forge/main/control/src/Makefile 2010-09-07 08:40:56 UTC (rev 7673) @@ -1,6 +1,7 @@ all: slab08nd.oct slab13dd.oct slsb10hd.oct slsb10ed.oct slab13bd.oct \ slsb01bd.oct slsb10fd.oct slsb10dd.oct slsb03md.oct slsb04md.oct \ - slsb04qd.oct slsg03ad.oct slsb02od.oct slab13ad.oct slab01od.oct + slsb04qd.oct slsg03ad.oct slsb02od.oct slab13ad.oct slab01od.oct \ + sltb01pd.oct # transmission zeros of state-space models slab08nd.oct: slab08nd.cc @@ -111,5 +112,11 @@ mkoctfile slab01od.cc \ AB01OD.f AB01ND.f MB03OY.f MB01PD.f MB01QD.f +# Minimal realization of state-space models +sltb01pd.oct: sltb01pd.cc + mkoctfile sltb01pd.cc \ + TB01PD.f TB01XD.f TB01ID.f AB07MD.f TB01UD.f \ + MB03OY.f MB01PD.f MB01QD.f + clean: rm *.o core octave-core *.oct *~ Added: trunk/octave-forge/main/control/src/TB01PD.f =================================================================== --- trunk/octave-forge/main/control/src/TB01PD.f (rev 0) +++ trunk/octave-forge/main/control/src/TB01PD.f 2010-09-07 08:40:56 UTC (rev 7673) @@ -0,0 +1,352 @@ + SUBROUTINE TB01PD( JOB, EQUIL, N, M, P, A, LDA, B, LDB, C, LDC, + $ NR, TOL, IWORK, DWORK, LDWORK, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2009 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To find a reduced (controllable, observable, or minimal) state- +C space representation (Ar,Br,Cr) for any original state-space +C representation (A,B,C). The matrix Ar is in upper block +C Hessenberg form. +C +C ARGUMENTS +C +C Mode Parameters +C +C JOB CHARACTER*1 +C Indicates whether the user wishes to remove the +C uncontrollable and/or unobservable parts as follows: +C = 'M': Remove both the uncontrollable and unobservable +C parts to get a minimal state-space representation; +C = 'C': Remove the uncontrollable part only to get a +C controllable state-space representation; +C = 'O': Remove the unobservable part only to get an +C observable state-space representation. +C +C EQUIL CHARACTER*1 +C Specifies whether the user wishes to preliminarily balance +C the triplet (A,B,C) as follows: +C = 'S': Perform balancing (scaling); +C = 'N': Do not perform balancing. +C +C Input/Output Parameters +C +C N (input) INTEGER +C The order of the original state-space representation, i.e. +C the order of the matrix A. N >= 0. +C +C M (input) INTEGER +C The number of system inputs. M >= 0. +C +C P (input) INTEGER +C The number of system outputs. P >= 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C On entry, the leading N-by-N part of this array must +C contain the original state dynamics matrix A. +C On exit, the leading NR-by-NR part of this array contains +C the upper block Hessenberg state dynamics matrix Ar of a +C minimal, controllable, or observable realization for the +C original system, depending on the value of JOB, JOB = 'M', +C JOB = 'C', or JOB = 'O', respectively. +C +C LDA INTEGER +C The leading dimension of array A. LDA >= MAX(1,N). +C +C B (input/output) DOUBLE PRECISION array, dimension (LDB,M), +C if JOB = 'C', or (LDB,MAX(M,P)), otherwise. +C On entry, the leading N-by-M part of this array must +C contain the original input/state matrix B; if JOB = 'M', +C or JOB = 'O', the remainder of the leading N-by-MAX(M,P) +C part is used as internal workspace. +C On exit, the leading NR-by-M part of this array contains +C the transformed input/state matrix Br of a minimal, +C controllable, or observable realization for the original +C system, depending on the value of JOB, JOB = 'M', +C JOB = 'C', or JOB = 'O', respectively. +C If JOB = 'C', only the first IWORK(1) rows of B are +C nonzero. +C +C LDB INTEGER +C The leading dimension of array B. LDB >= MAX(1,N). +C +C C (input/output) DOUBLE PRECISION array, dimension (LDC,N) +C On entry, the leading P-by-N part of this array must +C contain the original state/output matrix C; if JOB = 'M', +C or JOB = 'O', the remainder of the leading MAX(M,P)-by-N +C part is used as internal workspace. +C On exit, the leading P-by-NR part of this array contains +C the transformed state/output matrix Cr of a minimal, +C controllable, or observable realization for the original +C system, depending on the value of JOB, JOB = 'M', +C JOB = 'C', or JOB = 'O', respectively. +C If JOB = 'M', or JOB = 'O', only the last IWORK(1) columns +C (in the first NR columns) of C are nonzero. +C +C LDC INTEGER +C The leading dimension of array C. +C LDC >= MAX(1,M,P) if N > 0. +C LDC >= 1 if N = 0. +C +C NR (output) INTEGER +C The order of the reduced state-space representation +C (Ar,Br,Cr) of a minimal, controllable, or observable +C realization for the original system, depending on +C JOB = 'M', JOB = 'C', or JOB = 'O'. +C +C Tolerances +C +C TOL DOUBLE PRECISION +C The tolerance to be used in rank determination when +C transforming (A, B, C). If the user sets TOL > 0, then +C the given value of TOL is used as a lower bound for the +C reciprocal condition number (see the description of the +C argument RCOND in the SLICOT routine MB03OD); a +C (sub)matrix whose estimated condition number is less than +C 1/TOL is considered to be of full rank. If the user sets +C TOL <= 0, then an implicitly computed, default tolerance +C (determined by the SLICOT routine TB01UD) is used instead. +C +C Workspace +C +C IWORK INTEGER array, dimension (N+MAX(M,P)) +C On exit, if INFO = 0, the first nonzero elements of +C IWORK(1:N) return the orders of the diagonal blocks of A. +C +C DWORK DOUBLE PRECISION array, dimension (LDWORK) +C On exit, if INFO = 0, DWORK(1) returns the optimal value +C of LDWORK. +C +C LDWORK INTEGER +C The length of the array DWORK. +C LDWORK >= MAX(1, N + MAX(N, 3*M, 3*P)). +C For optimum performance LDWORK should be larger. +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C < 0: if INFO = -i, the i-th argument had an illegal +C value. +C +C METHOD +C +C If JOB = 'M', the matrices A and B are operated on by orthogonal +C similarity transformations (made up of products of Householder +C transformations) so as to produce an upper block Hessenberg matrix +C A1 and a matrix B1 with all but its first rank(B) rows zero; this +C separates out the controllable part of the original system. +C Applying the same algorithm to the dual of this subsystem, +C therefore separates out the controllable and observable (i.e. +C minimal) part of the original system representation, with the +C final Ar upper block Hessenberg (after using pertransposition). +C If JOB = 'C', or JOB = 'O', only the corresponding part of the +C above procedure is applied. +C +C REFERENCES +C +C [1] Van Dooren, P. +C The Generalized Eigenstructure Problem in Linear System +C Theory. (Algorithm 1) +C IEEE Trans. Auto. Contr., AC-26, pp. 111-129, 1981. +C +C NUMERICAL ASPECTS +C 3 +C The algorithm requires 0(N ) operations and is backward stable. +C +C CONTRIBUTOR +C +C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1998. +C +C REVISIONS +C +C A. Varga, DLR Oberpfaffenhofen, July 1998. +C A. Varga, DLR Oberpfaffenhofen, April 28, 1999. +C V. Sima, Research Institute for Informatics, Bucharest, Mar. 2004. +C +C KEYWORDS +C +C Hessenberg form, minimal realization, multivariable system, +C orthogonal transformation, state-space model, state-space +C representation. +C +C ****************************************************************** +C +C .. Parameters .. + INTEGER LDIZ + PARAMETER ( LDIZ = 1 ) + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. Scalar Arguments .. + CHARACTER EQUIL, JOB + INTEGER INFO, LDA, LDB, LDC, LDWORK, M, N, NR, P + DOUBLE PRECISION TOL +C .. Array Arguments .. + INTEGER IWORK(*) + DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*) +C .. Local Scalars .. + LOGICAL LEQUIL, LNJOBC, LNJOBO + INTEGER I, INDCON, ITAU, IZ, JWORK, KL, MAXMP, NCONT, + $ WRKOPT + DOUBLE PRECISION MAXRED +C .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +C .. External Subroutines .. + EXTERNAL AB07MD, TB01ID, TB01UD, TB01XD, XERBLA +C .. Intrinsic Functions .. + INTRINSIC INT, MAX, MIN +C .. Executable Statements .. +C + INFO = 0 + MAXMP = MAX( M, P ) + LNJOBC = .NOT.LSAME( JOB, 'C' ) + LNJOBO = .NOT.LSAME( JOB, 'O' ) + LEQUIL = LSAME( EQUIL, 'S' ) +C +C Test the input scalar arguments. +C + IF( LNJOBC .AND. LNJOBO .AND. .NOT.LSAME( JOB, 'M' ) ) THEN + INFO = -1 + ELSE IF( .NOT.LEQUIL .AND. .NOT.LSAME( EQUIL, 'N' ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( M.LT.0 ) THEN + INFO = -4 + ELSE IF( P.LT.0 ) THEN + INFO = -5 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -9 + ELSE IF( LDC.LT.1 .OR. ( N.GT.0 .AND. LDC.LT.MAXMP ) ) THEN + INFO = -11 + ELSE IF( LDWORK.LT.MAX( 1, N + MAX( N, 3*MAXMP ) ) ) THEN + INFO = -16 + END IF +C + IF ( INFO.NE.0 ) THEN +C +C Error return. +C + CALL XERBLA( 'TB01PD', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + IF ( N.EQ.0 .OR. ( LNJOBC .AND. MIN( N, P ).EQ.0 ) .OR. + $ ( LNJOBO .AND. MIN( N, M ).EQ.0 ) ) THEN + NR = 0 +C + DO 5 I = 1, N + IWORK(I) = 0 + 5 CONTINUE +C + DWORK(1) = ONE + RETURN + END IF +C +C If required, balance the triplet (A,B,C) (default MAXRED). +C Workspace: need N. +C +C (Note: Comments in the code beginning "Workspace:" describe the +C minimal amount of real workspace needed at that point in the code, +C as well as the preferred amount for good performance.) +C + IF ( LEQUIL ) THEN + MAXRED = ZERO + CALL TB01ID( 'A', N, M, P, MAXRED, A, LDA, B, LDB, C, LDC, + $ DWORK, INFO ) + WRKOPT = N + ELSE + WRKOPT = 1 + END IF +C + IZ = 1 + ITAU = 1 + JWORK = ITAU + N + IF ( LNJOBO ) THEN +C +C Separate out controllable subsystem (of order NCONT): +C A <-- Z'*A*Z, B <-- Z'*B, C <-- C*Z. +C +C Workspace: need N + MAX(N, 3*M, P). +C prefer larger. +C + CALL TB01UD( 'No Z', N, M, P, A, LDA, B, LDB, C, LDC, NCONT, + $ INDCON, IWORK, DWORK(IZ), LDIZ, DWORK(ITAU), TOL, + $ IWORK(N+1), DWORK(JWORK), LDWORK-JWORK+1, INFO ) +C + WRKOPT = INT( DWORK(JWORK) ) + JWORK - 1 + ELSE + NCONT = N + END IF +C + IF ( LNJOBC ) THEN +C +C Separate out the observable subsystem (of order NR): +C Form the dual of the subsystem of order NCONT (which is +C controllable, if JOB = 'M'), leaving rest as it is. +C + CALL AB07MD( 'Z', NCONT, M, P, A, LDA, B, LDB, C, LDC, DWORK, + $ 1, INFO ) +C +C And separate out the controllable part of this dual subsystem. +C +C Workspace: need NCONT + MAX(NCONT, 3*P, M). +C prefer larger. +C + CALL TB01UD( 'No Z', NCONT, P, M, A, LDA, B, LDB, C, LDC, NR, + $ INDCON, IWORK, DWORK(IZ), LDIZ, DWORK(ITAU), TOL, + $ IWORK(N+1), DWORK(JWORK), LDWORK-JWORK+1, INFO ) +C + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 ) +C +C Transpose and reorder (to get a block upper Hessenberg +C matrix A), giving, for JOB = 'M', the controllable and +C observable (i.e., minimal) part of original system. +C + IF( INDCON.GT.0 ) THEN + KL = IWORK(1) - 1 + IF ( INDCON.GE.2 ) + $ KL = KL + IWORK(2) + ELSE + KL = 0 + END IF + CALL TB01XD( 'Zero D', NR, P, M, KL, MAX( 0, NR-1 ), A, LDA, + $ B, LDB, C, LDC, DWORK, 1, INFO ) + ELSE + NR = NCONT + END IF +C +C Annihilate the trailing components of IWORK(1:N). +C + DO 10 I = INDCON + 1, N + IWORK(I) = 0 + 10 CONTINUE +C +C Set optimal workspace dimension. +C + DWORK(1) = WRKOPT + RETURN +C *** Last line of TB01PD *** + END Added: trunk/octave-forge/main/control/src/TB01UD.f =================================================================== --- trunk/octave-forge/main/control/src/TB01UD.f (rev 0) +++ trunk/octave-forge/main/control/src/TB01UD.f 2010-09-07 08:40:56 UTC (rev 7673) @@ -0,0 +1,491 @@ + SUBROUTINE TB01UD( JOBZ, N, M, P, A, LDA, B, LDB, C, LDC, NCONT, + $ INDCON, NBLK, Z, LDZ, TAU, TOL, IWORK, DWORK, + $ LDWORK, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2009 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To find a controllable realization for the linear time-invariant +C multi-input system +C +C dX/dt = A * X + B * U, +C Y = C * X, +C +C where A, B, and C are N-by-N, N-by-M, and P-by-N matrices, +C respectively, and A and B are reduced by this routine to +C orthogonal canonical form using (and optionally accumulating) +C orthogonal similarity transformations, which are also applied +C to C. Specifically, the system (A, B, C) is reduced to the +C triplet (Ac, Bc, Cc), where Ac = Z' * A * Z, Bc = Z' * B, +C Cc = C * Z, with +C +C [ Acont * ] [ Bcont ] +C Ac = [ ], Bc = [ ], +C [ 0 Auncont ] [ 0 ] +C +C and +C +C [ A11 A12 . . . A1,p-1 A1p ] [ B1 ] +C [ A21 A22 . . . A2,p-1 A2p ] [ 0 ] +C [ 0 A32 . . . A3,p-1 A3p ] [ 0 ] +C Acont = [ . . . . . . . ], Bc = [ . ], +C [ . . . . . . ] [ . ] +C [ . . . . . ] [ . ] +C [ 0 0 . . . Ap,p-1 App ] [ 0 ] +C +C where the blocks B1, A21, ..., Ap,p-1 have full row ranks and +C p is the controllability index of the pair. The size of the +C block Auncont is equal to the dimension of the uncontrollable +C subspace of the pair (A, B). +C +C ARGUMENTS +C +C Mode Parameters +C +C JOBZ CHARACTER*1 +C Indicates whether the user wishes to accumulate in a +C matrix Z the orthogonal similarity transformations for +C reducing the system, as follows: +C = 'N': Do not form Z and do not store the orthogonal +C transformations; +C = 'F': Do not form Z, but store the orthogonal +C transformations in the factored form; +C = 'I': Z is initialized to the unit matrix and the +C orthogonal transformation matrix Z is returned. +C +C Input/Output Parameters +C +C N (input) INTEGER +C The order of the original state-space representation, +C i.e. the order of the matrix A. N >= 0. +C +C M (input) INTEGER +C The number of system inputs, or of columns of B. M >= 0. +C +C P (input) INTEGER +C The number of system outputs, or of rows of C. P >= 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C On entry, the leading N-by-N part of this array must +C contain the original state dynamics matrix A. +C On exit, the leading NCONT-by-NCONT part contains the +C upper block Hessenberg state dynamics matrix Acont in Ac, +C given by Z' * A * Z, of a controllable realization for +C the original system. The elements below the first block- +C subdiagonal are set to zero. The leading N-by-N part +C contains the matrix Ac. +C +C LDA INTEGER +C The leading dimension of array A. LDA >= MAX(1,N). +C +C B (input/output) DOUBLE PRECISION array, dimension (LDB,M) +C On entry, the leading N-by-M part of this array must +C contain the input matrix B. +C On exit, the leading NCONT-by-M part of this array +C contains the transformed input matrix Bcont in Bc, given +C by Z' * B, with all elements but the first block set to +C zero. The leading N-by-M part contains the matrix Bc. +C +C LDB INTEGER +C The leading dimension of array B. LDB >= MAX(1,N). +C +C C (input/output) DOUBLE PRECISION array, dimension (LDC,N) +C On entry, the leading P-by-N part of this array must +C contain the output matrix C. +C On exit, the leading P-by-N part of this array contains +C the transformed output matrix Cc, given by C * Z. +C +C LDC INTEGER +C The leading dimension of array C. LDC >= MAX(1,P). +C +C NCONT (output) INTEGER +C The order of the controllable state-space representation. +C +C INDCON (output) INTEGER +C The controllability index of the controllable part of the +C system representation. +C +C NBLK (output) INTEGER array, dimension (N) +C The leading INDCON elements of this array contain the +C the orders of the diagonal blocks of Acont. +C +C Z (output) DOUBLE PRECISION array, dimension (LDZ,N) +C If JOBZ = 'I', then the leading N-by-N part of this +C array contains the matrix of accumulated orthogonal +C similarity transformations which reduces the given system +C to orthogonal canonical form. +C If JOBZ = 'F', the elements below the diagonal, with the +C array TAU, represent the orthogonal transformation matrix +C as a product of elementary reflectors. The transformation +C matrix can then be obtained by calling the LAPACK Library +C routine DORGQR. +C If JOBZ = 'N', the array Z is not referenced and can be +C supplied as a dummy array (i.e. set parameter LDZ = 1 and +C declare this array to be Z(1,1) in the calling program). +C +C LDZ INTEGER +C The leading dimension of array Z. If JOBZ = 'I' or +C JOBZ = 'F', LDZ >= MAX(1,N); if JOBZ = 'N', LDZ >= 1. +C +C TAU (output) DOUBLE PRECISION array, dimension (N) +C The elements of TAU contain the scalar factors of the +C elementary reflectors used in the reduction of B and A. +C +C Tolerances +C +C TOL DOUBLE PRECISION +C The tolerance to be used in rank determination when +C transforming (A, B). If the user sets TOL > 0, then +C the given value of TOL is used as a lower bound for the +C reciprocal condition number (see the description of the +C argument RCOND in the SLICOT routine MB03OD); a +C (sub)matrix whose estimated condition number is less than +C 1/TOL is considered to be of full rank. If the user sets +C TOL <= 0, then an implicitly computed, default tolerance, +C defined by TOLDEF = N*N*EPS, is used instead, where EPS +C is the machine precision (see LAPACK Library routine +C DLAMCH). +C +C Workspace +C +C IWORK INTEGER array, dimension (M) +C +C DWORK DOUBLE PRECISION array, dimension (LDWORK) +C On exit, if INFO = 0, DWORK(1) returns the optimal value +C of LDWORK. +C +C LDWORK INTEGER +C The length of the array DWORK. +C LDWORK >= MAX(1, N, 3*M, P). +C For optimum performance LDWORK should be larger. +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C < 0: if INFO = -i, the i-th argument had an illegal +C value. +C +C METHOD +C +C Matrix B is first QR-decomposed and the appropriate orthogonal +C similarity transformation applied to the matrix A. Leaving the +C first rank(B) states unchanged, the remaining lower left block +C of A is then QR-decomposed and the new orthogonal matrix, Q1, +C is also applied to the right of A to complete the similarity +C transformation. By continuing in this manner, a completely +C controllable state-space pair (Acont, Bcont) is found for the +C given (A, B), where Acont is upper block Hessenberg with each +C subdiagonal block of full row rank, and Bcont is zero apart from +C its (independent) first rank(B) rows. +C All orthogonal transformations determined in this process are also +C applied to the matrix C, from the right. +C NOTE that the system controllability indices are easily +C calculated from the dimensions of the blocks of Acont. +C +C REFERENCES +C +C [1] Konstantinov, M.M., Petkov, P.Hr. and Christov, N.D. +C Orthogonal Invariants and Canonical Forms for Linear +C Controllable Systems. +C Proc. 8th IFAC World Congress, Kyoto, 1, pp. 49-54, 1981. +C +C [2] Paige, C.C. +C Properties of numerical algorithms related to computing +C controllablity. +C IEEE Trans. Auto. Contr., AC-26, pp. 130-138, 1981. +C +C [3] Petkov, P.Hr., Konstantinov, M.M., Gu, D.W. and +C Postlethwaite, I. +C Optimal Pole Assignment Design of Linear Multi-Input Systems. +C Leicester University, Report 99-11, May 1996. +C +C NUMERICAL ASPECTS +C 3 +C The algorithm requires 0(N ) operations and is backward stable. +C +C FURTHER COMMENTS +C +C If the system matrices A and B are badly scaled, it would be +C useful to scale them with SLICOT routine TB01ID, before calling +C the routine. +C +C CONTRIBUTOR +C +C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1998. +C +C REVISIONS +C +C V. Sima, Katholieke Univ. Leuven, Belgium, May 1999, Nov. 2003. +C A. Varga, DLR Oberpfaffenhofen, March 2002, Nov. 2003. +C +C KEYWORDS +C +C Controllability, minimal realization, orthogonal canonical form, +C orthogonal transformation. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. Scalar Arguments .. + CHARACTER JOBZ + INTEGER INDCON, INFO, LDA, LDB, LDC, LDWORK, LDZ, M, N, + $ NCONT, P + DOUBLE PRECISION TOL +C .. Array Arguments .. + DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), TAU(*), + $ Z(LDZ,*) + INTEGER IWORK(*), NBLK(*) +C .. Local Scalars .. + LOGICAL LJOBF, LJOBI, LJOBZ + INTEGER IQR, ITAU, J, MCRT, NBL, NCRT, NI, NJ, RANK, + $ WRKOPT + DOUBLE PRECISION ANORM, BNORM, FNRM, TOLDEF +C .. Local Arrays .. + DOUBLE PRECISION SVAL(3) +C .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2 + EXTERNAL DLAMCH, DLANGE, DLAPY2, LSAME +C .. External Subroutines .. + EXTERNAL DCOPY, DLACPY, DLAPMT, DLASET, DORGQR, DORMQR, + $ MB01PD, MB03OY, XERBLA +C .. Intrinsic Functions .. + INTRINSIC DBLE, INT, MAX, MIN +C .. +C .. Executable Statements .. +C + INFO = 0 + LJOBF = LSAME( JOBZ, 'F' ) + LJOBI = LSAME( JOBZ, 'I' ) + LJOBZ = LJOBF.OR.LJOBI +C +C Test the input scalar arguments. +C + IF( .NOT.LJOBZ .AND. .NOT.LSAME( JOBZ, 'N' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( M.LT.0 ) THEN + INFO = -3 + ELSE IF( P.LT.0 ) THEN + INFO = -4 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -6 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( LDC.LT.MAX( 1, P ) ) THEN + INFO = -10 + ELSE IF( .NOT.LJOBZ .AND. LDZ.LT.1 .OR. + $ LJOBZ .AND. LDZ.LT.MAX( 1, N ) ) THEN + INFO = -15 + ELSE IF( LDWORK.LT.MAX( 1, N, 3*M, P ) ) THEN + INFO = -20 + END IF +C + IF ( INFO.NE.0 ) THEN +C +C Error return. +C + CALL XERBLA( 'TB01UD', -INFO ) + RETURN + END IF +C + NCONT = 0 + INDCON = 0 +C +C Calculate the absolute norms of A and B (used for scaling). +C + ANORM = DLANGE( 'M', N, N, A, LDA, DWORK ) + BNORM = DLANGE( 'M', N, M, B, LDB, DWORK ) +C +C Quick return if possible. +C + IF ( MIN( N, M ).EQ.0 .OR. BNORM.EQ.ZERO ) THEN + IF( N.GT.0 ) THEN + IF ( LJOBI ) THEN + CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) + ELSE IF ( LJOBF ) THEN + CALL DLASET( 'Full', N, N, ZERO, ZERO, Z, LDZ ) + CALL DLASET( 'Full', N, 1, ZERO, ZERO, TAU, N ) + END IF + END IF + DWORK(1) = ONE + RETURN + END IF +C +C Scale (if needed) the matrices A and B. +C + CALL MB01PD( 'S', 'G', N, N, 0, 0, ANORM, 0, NBLK, A, LDA, INFO ) + CALL MB01PD( 'S', 'G', N, M, 0, 0, BNORM, 0, NBLK, B, LDB, INFO ) +C +C Compute the Frobenius norm of [ B A ] (used for rank estimation). +C + FNRM = DLAPY2( DLANGE( 'F', N, M, B, LDB, DWORK ), + $ DLANGE( 'F', N, N, A, LDA, DWORK ) ) +C + TOLDEF = TOL + IF ( TOLDEF.LE.ZERO ) THEN +C +C Use the default tolerance in controllability determination. +C + TOLDEF = DBLE( N*N )*DLAMCH( 'EPSILON' ) + END IF +C + IF ( FNRM.LT.TOLDEF ) + $ FNRM = ONE +C + WRKOPT = 1 + NI = 0 + ITAU = 1 + NCRT = N + MCRT = M + IQR = 1 +C +C (Note: Comments in the code beginning "Workspace:" describe the +C minimal amount of real workspace needed at that point in the +C code, as well as the preferred amount for good performance. +C NB refers to the optimal block size for the immediately +C following subroutine, as returned by ILAENV.) +C + 10 CONTINUE +C +C Rank-revealing QR decomposition with column pivoting. +C The calculation is performed in NCRT rows of B starting from +C the row IQR (initialized to 1 and then set to rank(B)+1). +C Workspace: 3*MCRT. +C + CALL MB03OY( NCRT, MCRT, B(IQR,1), LDB, TOLDEF, FNRM, RANK, + $ SVAL, IWORK, TAU(ITAU), DWORK, INFO ) +C + IF ( RANK.NE.0 ) THEN + NJ = NI + NI = NCONT + NCONT = NCONT + RANK + INDCON = INDCON + 1 + NBLK(INDCON) = RANK +C +C Premultiply and postmultiply the appropriate block row +C and block column of A by Q' and Q, respectively. +C Workspace: need NCRT; +C prefer NCRT*NB. +C + CALL DORMQR( 'Left', 'Transpose', NCRT, NCRT, RANK, + $ B(IQR,1), LDB, TAU(ITAU), A(NI+1,NI+1), LDA, + $ DWORK, LDWORK, INFO ) + WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) ) +C +C Workspace: need N; +C prefer N*NB. +C + CALL DORMQR( 'Right', 'No transpose', N, NCRT, RANK, + $ B(IQR,1), LDB, TAU(ITAU), A(1,NI+1), LDA, + $ DWORK, LDWORK, INFO ) + WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) ) +C +C Postmultiply the appropriate block column of C by Q. +C Workspace: need P; +C prefer P*NB. +C + CALL DORMQR( 'Right', 'No transpose', P, NCRT, RANK, + $ B(IQR,1), LDB, TAU(ITAU), C(1,NI+1), LDC, + $ DWORK, LDWORK, INFO ) + WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) ) +C +C If required, save transformations. +C + IF ( LJOBZ.AND.NCRT.GT.1 ) THEN + CALL DLACPY( 'L', NCRT-1, MIN( RANK, NCRT-1 ), + $ B(IQR+1,1), LDB, Z(NI+2,ITAU), LDZ ) + END IF +C +C Zero the subdiagonal elements of the current matrix. +C + IF ( RANK.GT.1 ) + $ CALL DLASET( 'L', RANK-1, RANK-1, ZERO, ZERO, B(IQR+1,1), + $ LDB ) +C +C Backward permutation of the columns of B or A. +C + IF ( INDCON.EQ.1 ) THEN + CALL DLAPMT( .FALSE., RANK, M, B(IQR,1), LDB, IWORK ) + IQR = RANK + 1 + ELSE + DO 20 J = 1, MCRT + CALL DCOPY( RANK, B(IQR,J), 1, A(NI+1,NJ+IWORK(J)), + $ 1 ) + 20 CONTINUE + END IF +C + ITAU = ITAU + RANK + IF ( RANK.NE.NCRT ) THEN + MCRT = RANK + NCRT = NCRT - RANK + CALL DLACPY( 'G', NCRT, MCRT, A(NCONT+1,NI+1), LDA, + $ B(IQR,1), LDB ) + CALL DLASET( 'G', NCRT, MCRT, ZERO, ZERO, + $ A(NCONT+1,NI+1), LDA ) + GO TO 10 + END IF + END IF +C +C If required, accumulate transformations. +C Workspace: need N; prefer N*NB. +C + IF ( LJOBI ) THEN + CALL DORGQR( N, N, ITAU-1, Z, LDZ, TAU, DWORK, + $ LDWORK, INFO ) + WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) ) + END IF +C +C Annihilate the trailing blocks of B. +C + IF( IQR.LE.N ) + $ CALL DLASET( 'G', N-IQR+1, M, ZERO, ZERO, B(IQR,1), LDB ) +C +C Annihilate the trailing elements of TAU, if JOBZ = 'F'. +C + IF ( LJOBF ) THEN + DO 30 J = ITAU, N + TAU(J) = ZERO + 30 CONTINUE + END IF +C +C Undo scaling of A and B. +C + IF ( INDCON.LT.N ) THEN + NBL = INDCON + 1 + NBLK(NBL) = N - NCONT + ELSE + NBL = 0 + END IF + CALL MB01PD( 'U', 'H', N, N, 0, 0, ANORM, NBL, NBLK, A, LDA, + $ INFO ) + CALL MB01PD( 'U', 'G', NBLK(1), M, 0, 0, BNORM, 0, NBLK, B, LDB, + $ INFO ) +C +C Set optimal workspace dimension. +C + DWORK(1) = WRKOPT + RETURN +C *** Last line of TB01UD *** + END Added: trunk/octave-forge/main/control/src/TB01XD.f =================================================================== --- trunk/octave-forge/main/control/src/TB01XD.f (rev 0) +++ trunk/octave-forge/main/control/src/TB01XD.f 2010-09-07 08:40:56 UTC (rev 7673) @@ -0,0 +1,284 @@ + SUBROUTINE TB01XD( JOBD, N, M, P, KL, KU, A, LDA, B, LDB, C, LDC, + $ D, LDD, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2009 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To apply a special transformation to a system given as a triple +C (A,B,C), +C +C A <-- P * A' * P, B <-- P * C', C <-- B' * P, +C +C where P is a matrix with 1 on the secondary diagonal, and with 0 +C in the other entries. Matrix A can be specified as a band matrix. +C Optionally, matrix D of the system can be transposed. This +C transformation is actually a special similarity transformation of +C the dual system. +C +C ARGUMENTS +C +C Mode Parameters +C +C JOBD CHARACTER*1 +C Specifies whether or not a non-zero matrix D appears in +C the given state space model: +C = 'D': D is present; +C = 'Z': D is assumed a zero matrix. +C +C Input/Output Parameters +C +C N (input) INTEGER +C The order of the matrix A, the number of rows of matrix B +C and the number of columns of matrix C. +C N represents the dimension of the state vector. N >= 0. +C +C M (input) INTEGER. +C The number of columns of matrix B. +C M represents the dimension of input vector. M >= 0. +C +C P (input) INTEGER. +C The number of rows of matrix C. +C P represents the dimension of output vector. P >= 0. +C +C KL (input) INTEGER +C The number of subdiagonals of A to be transformed. +C MAX( 0, N-1 ) >= KL >= 0. +C +C KU (input) INTEGER +C The number of superdiagonals of A to be transformed. +C MAX( 0, N-1 ) >= KU >= 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C On entry, the leading N-by-N part of this array must +C contain the system state matrix A. +C On exit, the leading N-by-N part of this array contains +C the transformed (pertransposed) matrix P*A'*P. +C +C LDA INTEGER +C The leading dimension of the array A. LDA >= MAX(1,N). +C +C B (input/output) DOUBLE PRECISION array, dimension +C (LDB,MAX(M,P)) +C On entry, the leading N-by-M part of this array must +C contain the original input/state matrix B. +C On exit, the leading N-by-P part of this array contains +C the dual input/state matrix P*C'. +C +C LDB INTEGER +C The leading dimension of the array B. +C LDB >= MAX(1,N) if M > 0 or P > 0. +C LDB >= 1 if M = 0 and P = 0. +C +C C (input/output) DOUBLE PRECISION array, dimension (LDC,N) +C On entry, the leading P-by-N part of this array must +C contain the original state/output matrix C. +C On exit, the leading M-by-N part of this array contains +C the dual state/output matrix B'*P. +C +C LDC INTEGER +C The leading dimension of array C. +C LDC >= MAX(1,M,P) if N > 0. +C LDC >= 1 if N = 0. +C +C D (input/output) DOUBLE PRECISION array, dimension +C (LDD,MAX(M,P)) +C On entry, if JOBD = 'D', the leading P-by-M part of this +C array must contain the original direct transmission +C matrix D. +C On exit, if JOBD = 'D', the leading M-by-P part of this +C array contains the transposed direct transmission matrix +C D'. The array D is not referenced if JOBD = 'Z'. +C +C LDD INTEGER +C The leading dimension of array D. +C LDD >= MAX(1,M,P) if JOBD = 'D'. +C LDD >= 1 if JOBD = 'Z'. +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit. +C < 0: if INFO = -i, the i-th argument had an illegal +C value. +C +C METHOD +C +C The rows and/or columns of the matrices of the triplet (A,B,C) +C and, optionally, of the matrix D are swapped in a special way. +C +C NUMERICAL ASPECTS +C +C None. +C +C CONTRIBUTOR +C +C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1998. +C Partly based on routine DMPTR (A. Varga, German Aerospace +C Research Establishment, DLR, Aug. 1992). +C +C +C REVISIONS +C +C 07-31-1998, 04-25-1999, A. Varga. +C 03-16-2004, V. Sima. +C +C KEYWORDS +C +C Matrix algebra, matrix operations, similarity transformation. +C +C ********************************************************************* +C +C .. +C .. Scalar Arguments .. + CHARACTER JOBD + INTEGER INFO, KL, KU, LDA, LDB, LDC, LDD, M, N, P +C .. +C .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ), + $ D( LDD, * ) +C .. +C .. Local Scalars .. + LOGICAL LJOBD + INTEGER J, J1, LDA1, MAXMP, MINMP, NM1 +C .. +C .. External functions .. + LOGICAL LSAME + EXTERNAL LSAME +C .. +C .. External Subroutines .. + EXTERNAL DCOPY, DSWAP, XERBLA +C .. +C .. Intrinsic Functions .. + INTRINSIC MAX, MIN +C .. +C .. Executable Statements .. +C +C Test the scalar input arguments. +C + INFO = 0 + LJOBD = LSAME( JOBD, 'D' ) + MAXMP = MAX( M, P ) + MINMP = MIN( M, P ) + NM1 = N - 1 +C + IF( .NOT.LJOBD .AND. .NOT.LSAME( JOBD, 'Z' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( M.LT.0 ) THEN + INFO = -3 + ELSE IF( P.LT.0 ) THEN + INFO = -4 + ELSE IF( KL.LT.0 .OR. KL.GT.MAX( 0, NM1 ) ) THEN + INFO = -5 + ELSE IF( KU.LT.0 .OR. KU.GT.MAX( 0, NM1 ) ) THEN + INFO = -6 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( ( MAXMP.GT.0 .AND. LDB.LT.MAX( 1, N ) ) .OR. + $ ( MINMP.EQ.0 .AND. LDB.LT.1 ) ) THEN + INFO = -10 + ELSE IF( LDC.LT.1 .OR. ( N.GT.0 .AND. LDC.LT.MAXMP ) ) THEN + INFO = -12 + ELSE IF( LDD.LT.1 .OR. ( LJOBD .AND. LDD.LT.MAXMP ) ) THEN + INFO = -14 + END IF +C + IF( INFO.NE.0 ) THEN +C +C Error return. +C + CALL XERBLA( 'TB01XD', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + IF ( LJOBD ) THEN +C +C Replace D by D', if non-scalar. +C + DO 5 J = 1, MAXMP + IF ( J.LT.MINMP ) THEN + CALL DSWAP( MINMP-J, D(J+1,J), 1, D(J,J+1), LDD ) + ELSE IF ( J.GT.P ) THEN + CALL DCOPY( P, D(1,J), 1, D(J,1), LDD ) + ELSE IF ( J.GT.M ) THEN + CALL DCOPY( M, D(J,1), LDD, D(1,J), 1 ) + END IF + 5 CONTINUE +C + END IF +C + IF( N.EQ.0 ) + $ RETURN +C +C Replace matrix A by P*A'*P. +C + IF ( KL.EQ.NM1 .AND. KU.EQ.NM1 ) THEN +C +C Full matrix A. +C + DO 10 J = 1, NM1 + CALL DSWAP( N-J, A( 1, J ), 1, A( N-J+1, J+1 ), -LDA ) + 10 CONTINUE +C + ELSE +C +C Band matrix A. +C + LDA1 = LDA + 1 +C +C Pertranspose the KL subdiagonals. +C + DO 20 J = 1, MIN( KL, N-2 ) + J1 = ( N - J )/2 + CALL DSWAP( J1, A(J+1,1), LDA1, A(N-J1+1,N-J1+1-J), -LDA1 ) + 20 CONTINUE +C +C Pertranspose the KU superdiagonals. +C + DO 30 J = 1, MIN( KU, N-2 ) + J1 = ( N - J )/2 + CALL DSWAP( J1, A(1,J+1), LDA1, A(N-J1+1-J,N-J1+1), -LDA1 ) + 30 CONTINUE +C +C Pertranspose the diagonal. +C + J1 = N/2 + CALL DSWAP( J1, A(1,1), LDA1, A(N-J1+1,N-J1+1), -LDA1 ) +C + END IF +C +C Replace matrix B by P*C' and matrix C by B'*P. +C + DO 40 J = 1, MAXMP + IF ( J.LE.MINMP ) THEN + CALL DSWAP( N, B(1,J), 1, C(J,1), -LDC ) + ELSE IF ( J.GT.P ) THEN + CALL DCOPY( N, B(1,J), 1, C(J,1), -LDC ) + ELSE + CALL DCOPY( N, C(J,1), -LDC, B(1,J), 1 ) + END IF + 40 CONTINUE +C + RETURN +C *** Last line of TB01XD *** + END Added: trunk/octave-forge/main/control/src/sltb01pd.cc =================================================================== --- trunk/octave-forge/main/control/src/sltb01pd.cc (rev 0) +++ trunk/octave-forge/main/control/src/sltb01pd.cc 2010-09-07 08:40:56 UTC (rev 7673) @@ -0,0 +1,123 @@ +/* + +Copyright (C) 2010 Lukas F. Reichlin + +This file is part of LTI Syncope. + +LTI Syncope 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. + +LTI Syncope 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. If not, see <http://www.gnu.org/licenses/>. + +Minimal realization of state-space models. +Uses SLICOT TB01PD by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: September 2010 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include "common.cc" + +extern "C" +{ + int F77_FUNC (tb01pd, TB01PD) + (char& JOB, char& EQUIL, + int& N, int& M, int& P, + double* A, int& LDA, + double* B, int& LDB, + double* C, int& LDC, + int& NR, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& INFO); +} + +DEFUN_DLD (sltb01pd, args, nargout, "Slicot TB01PD Release 5.0") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin != 4) + { + print_usage (); + } + else + { + // arguments in + char job = 'M'; + char equil = 'N'; + + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + double tol = args(3).double_value (); + + int n = a.rows (); // n: number of states + int m = b.columns (); // m: number of inputs + int p = c.rows (); // p: number of outputs + + int lda = max (1, n); + int ldb = max (1, n); + int ldc; + + if (n == 0) + ldc = 1; + else + ldc = max (1, m, p); + + // arguments out + int nr; + + // workspace + int liwork = n + max (m, p); + int ldwork = max (1, n + max (n, 3*m, 3*p)); + + OCTAVE_LOCAL_BUFFER (int, iwork, liwork); + OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); + + // error indicators + int info; + + + // SLICOT routine TB01PD + F77_XFCN (tb01pd, TB01PD, + (job, equil, + n, m, p, + a.fortran_vec (), lda, + b.fortran_vec (), ldb, + c.fortran_vec (), ldc, + nr, + tol, + iwork, + dwork, ldwork, + info)); + + if (f77_exception_encountered) + error ("ss: minreal: sltb01pd: exception in SLICOT subroutine TB01PD"); + + if (info != 0) + error ("ss: minreal: sltb01pd: TB01PD returned info = %d", info); + + // return values + retval(0) = a; + retval(1) = b; + retval(2) = c; + retval(3) = octave_value (nr); + } + + return retval; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-09-07 23:48:27
|
Revision: 7680 http://octave.svn.sourceforge.net/octave/?rev=7680&view=rev Author: paramaniac Date: 2010-09-07 23:48:21 +0000 (Tue, 07 Sep 2010) Log Message: ----------- control: resize arrays inside oct-files Modified Paths: -------------- trunk/octave-forge/main/control/inst/@ss/__minreal__.m trunk/octave-forge/main/control/inst/hsvd.m trunk/octave-forge/main/control/src/slab13ad.cc trunk/octave-forge/main/control/src/sltb01pd.cc Modified: trunk/octave-forge/main/control/inst/@ss/__minreal__.m =================================================================== --- trunk/octave-forge/main/control/inst/@ss/__minreal__.m 2010-09-07 19:30:23 UTC (rev 7679) +++ trunk/octave-forge/main/control/inst/@ss/__minreal__.m 2010-09-07 23:48:21 UTC (rev 7680) @@ -33,10 +33,6 @@ [a, b, c, nr] = sltb01pd (sys.a, sys.b, sys.c, tol); - a = a(1:nr, 1:nr); - b = b(1:nr, :); - c = c(:, 1:nr); - retsys = ss (a, b, c, sys.d); retsys.lti = sys.lti; # retain i/o names and tsam Modified: trunk/octave-forge/main/control/inst/hsvd.m =================================================================== --- trunk/octave-forge/main/control/inst/hsvd.m 2010-09-07 19:30:23 UTC (rev 7679) +++ trunk/octave-forge/main/control/inst/hsvd.m 2010-09-07 23:48:21 UTC (rev 7680) @@ -26,7 +26,7 @@ ## Author: Lukas Reichlin <luk...@gm...> ## Created: January 2010 -## Version: 0.1 +## Version: 0.2 function hsv_r = hsvd (sys, prop = "offset", val = 1e-8) @@ -54,13 +54,10 @@ [hsv, ns] = slab13ad (a, b, c, discrete, alpha); - idx = 1 : ns; - hsv = hsv (idx); - if (nargout) hsv_r = hsv; else - bar (idx + (rows (a) - ns), hsv); + bar ((1:ns) + (rows (a) - ns), hsv); title ("Hankel Singular Values of Stable Part"); xlabel ("State"); ylabel ("State Energy"); Modified: trunk/octave-forge/main/control/src/slab13ad.cc =================================================================== --- trunk/octave-forge/main/control/src/slab13ad.cc 2010-09-07 19:30:23 UTC (rev 7679) +++ trunk/octave-forge/main/control/src/slab13ad.cc 2010-09-07 23:48:21 UTC (rev 7680) @@ -112,6 +112,9 @@ if (info != 0) error ("hsvd: slab13ad: AB13AD returned info = %d", info); + + // resize + hsv.resize (ns); // return values retval(0) = hsv; Modified: trunk/octave-forge/main/control/src/sltb01pd.cc =================================================================== --- trunk/octave-forge/main/control/src/sltb01pd.cc 2010-09-07 19:30:23 UTC (rev 7679) +++ trunk/octave-forge/main/control/src/sltb01pd.cc 2010-09-07 23:48:21 UTC (rev 7680) @@ -111,7 +111,12 @@ if (info != 0) error ("ss: minreal: sltb01pd: TB01PD returned info = %d", info); - + + // resize + a.resize (nr, nr); + b.resize (nr, m); + c.resize (p, nr); + // return values retval(0) = a; retval(1) = b; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-09-08 16:14:46
|
Revision: 7686 http://octave.svn.sourceforge.net/octave/?rev=7686&view=rev Author: paramaniac Date: 2010-09-08 16:14:40 +0000 (Wed, 08 Sep 2010) Log Message: ----------- control: assemble complex vector inside oct-file Modified Paths: -------------- trunk/octave-forge/main/control/inst/@ss/__zero__.m trunk/octave-forge/main/control/src/slab08nd.cc Modified: trunk/octave-forge/main/control/inst/@ss/__zero__.m =================================================================== --- trunk/octave-forge/main/control/inst/@ss/__zero__.m 2010-09-08 10:05:43 UTC (rev 7685) +++ trunk/octave-forge/main/control/inst/@ss/__zero__.m 2010-09-08 16:14:40 UTC (rev 7686) @@ -1,4 +1,4 @@ -## Copyright (C) 2009 Lukas F. Reichlin +## Copyright (C) 2009 - 2010 Lukas F. Reichlin ## ## This file is part of LTI Syncope. ## @@ -22,16 +22,11 @@ ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 -## Version: 0.1 +## Version: 0.2 function [zer, gain] = __zero__ (sys) - if (isempty (sys.a)) - zer = zeros (0, 1); - else - [alphar, alphai, beta] = slab08nd (sys.a, sys.b, sys.c, sys.d); - zer = (alphar + i*alphai) ./ beta; - endif + zer = slab08nd (sys.a, sys.b, sys.c, sys.d); lzer = length (zer); n = rows (sys.a); Modified: trunk/octave-forge/main/control/src/slab08nd.cc =================================================================== --- trunk/octave-forge/main/control/src/slab08nd.cc 2010-09-08 10:05:43 UTC (rev 7685) +++ trunk/octave-forge/main/control/src/slab08nd.cc 2010-09-08 16:14:40 UTC (rev 7686) @@ -23,13 +23,14 @@ Author: Lukas Reichlin <luk...@gm...> Created: November 2009 -Version: 0.3 +Version: 0.4 */ #include <octave/oct.h> #include <f77-fcn.h> #include "common.cc" +#include <complex> extern "C" { @@ -118,7 +119,7 @@ int info; // tolerance - double tol = 2.2204e-16; // TODO: use LAPACK dlamch + double tol = 0; // AB08ND uses DLAMCH for default tolerance // SLICOT routine AB08ND F77_XFCN (ab08nd, AB08ND, @@ -179,11 +180,21 @@ if (info2 != 0) error ("ss: zero: slab08nd: DGGEV returned info = %d", info2); - + + // assemble complex vector - adapted from DEFUN complex in data.cc + ColumnVector zeror (nu); + ColumnVector zeroi (nu); + + zeror = quotient (alphar, beta); + zeroi = quotient (alphai, beta); + + ComplexColumnVector zero (nu, Complex ()); + + for (octave_idx_type i = 0; i < nu; i++) + zero.xelem (i) = Complex (zeror(i), zeroi(i)); + // return values - retval(0) = alphar; - retval(1) = alphai; - retval(2) = beta; + retval(0) = zero; } return retval; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-09-09 20:00:29
|
Revision: 7695 http://octave.svn.sourceforge.net/octave/?rev=7695&view=rev Author: paramaniac Date: 2010-09-09 20:00:23 +0000 (Thu, 09 Sep 2010) Log Message: ----------- control: calculate gain of state-space models inside oct-file Modified Paths: -------------- trunk/octave-forge/main/control/inst/@ss/__zero__.m trunk/octave-forge/main/control/src/slab08nd.cc Modified: trunk/octave-forge/main/control/inst/@ss/__zero__.m =================================================================== --- trunk/octave-forge/main/control/inst/@ss/__zero__.m 2010-09-09 19:01:41 UTC (rev 7694) +++ trunk/octave-forge/main/control/inst/@ss/__zero__.m 2010-09-09 20:00:23 UTC (rev 7695) @@ -26,21 +26,6 @@ function [zer, gain] = __zero__ (sys) - zer = slab08nd (sys.a, sys.b, sys.c, sys.d); + [zer, gain] = slab08nd (sys.a, sys.b, sys.c, sys.d); - lzer = length (zer); - n = rows (sys.a); - m = columns (sys.b); - p = rows (sys.c); - - if (m == 1 && p == 1) - if (lzer == n) - gain = sys.d; - else - gain = sys.c * (sys.a^(n-1-lzer)) * sys.b; - endif - else - gain = []; - endif - endfunction \ No newline at end of file Modified: trunk/octave-forge/main/control/src/slab08nd.cc =================================================================== --- trunk/octave-forge/main/control/src/slab08nd.cc 2010-09-09 19:01:41 UTC (rev 7694) +++ trunk/octave-forge/main/control/src/slab08nd.cc 2010-09-09 20:00:23 UTC (rev 7695) @@ -31,6 +31,7 @@ #include <f77-fcn.h> #include "common.cc" #include <complex> +#include <xpow.h> extern "C" { @@ -181,6 +182,21 @@ if (info2 != 0) error ("ss: zero: slab08nd: DGGEV returned info = %d", info2); + // calculate gain + octave_value gain; + + if (m == 1 && p == 1) + { + if (nu == n) + gain = d; + else + gain = c * xpow (a, double (n-1-nu)) * b; + } + else + { + gain = Matrix (0, 0); + } + // assemble complex vector - adapted from DEFUN complex in data.cc ColumnVector zeror (nu); ColumnVector zeroi (nu); @@ -195,6 +211,7 @@ // return values retval(0) = zero; + retval(1) = gain; } return retval; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-09-13 11:08:45
|
Revision: 7707 http://octave.svn.sourceforge.net/octave/?rev=7707&view=rev Author: paramaniac Date: 2010-09-13 11:08:39 +0000 (Mon, 13 Sep 2010) Log Message: ----------- control: improve argument checking Modified Paths: -------------- trunk/octave-forge/main/control/inst/isctrb.m trunk/octave-forge/main/control/inst/isdetectable.m trunk/octave-forge/main/control/inst/isobsv.m trunk/octave-forge/main/control/inst/isstabilizable.m trunk/octave-forge/main/control/src/slab01od.cc Modified: trunk/octave-forge/main/control/inst/isctrb.m =================================================================== --- trunk/octave-forge/main/control/inst/isctrb.m 2010-09-13 10:33:45 UTC (rev 7706) +++ trunk/octave-forge/main/control/inst/isctrb.m 2010-09-13 11:08:39 UTC (rev 7707) @@ -51,9 +51,9 @@ ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 -## Version: 0.2 +## Version: 0.2.1 -function [bool, u] = isctrb (a, b = 0, tol = 0) +function [bool, u] = isctrb (a, b = [], tol = []) if (nargin < 1 || nargin > 3) print_usage (); @@ -70,10 +70,10 @@ rows (a), columns (a), rows (b), columns (b)); endif - ## check tol dimensions - if (! isscalar (tol)) - error ("isctrb: tol(%dx%d) must be a scalar", - rows (tol), columns (tol)); + if (! isreal (tol) && ! isscalar (tol)) + error ("isctrb: tol must be a real scalar"); + elseif (isempty (tol)) + tol = 0; # default tolerance endif [ac, bc, u, ncont] = slab01od (a, b, tol); Modified: trunk/octave-forge/main/control/inst/isdetectable.m =================================================================== --- trunk/octave-forge/main/control/inst/isdetectable.m 2010-09-13 10:33:45 UTC (rev 7706) +++ trunk/octave-forge/main/control/inst/isdetectable.m 2010-09-13 11:08:39 UTC (rev 7707) @@ -56,7 +56,7 @@ ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 -## Version: 0.2 +## Version: 0.2.1 function bool = isdetectable (a, c = [], tol = [], dflg = 0) @@ -67,7 +67,7 @@ print_usage (); endif tol = c; - dflg = isdt (a); + dflg = ! isct (a); [a, b, c] = ssdata (a); elseif (nargin == 1) # isdetectable (a, c, ...) print_usage (); Modified: trunk/octave-forge/main/control/inst/isobsv.m =================================================================== --- trunk/octave-forge/main/control/inst/isobsv.m 2010-09-13 10:33:45 UTC (rev 7706) +++ trunk/octave-forge/main/control/inst/isobsv.m 2010-09-13 11:08:39 UTC (rev 7707) @@ -51,7 +51,7 @@ ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 -## Version: 0.2 +## Version: 0.2.1 function [bool, u] = isobsv (a, c = [], tol = []) @@ -67,11 +67,7 @@ print_usage (); endif - if (isempty (tol)) - [bool, u] = isctrb (a.', c.'); - else - [bool, u] = isctrb (a.', c.', tol); - endif + [bool, u] = isctrb (a.', c.', tol); endfunction Modified: trunk/octave-forge/main/control/inst/isstabilizable.m =================================================================== --- trunk/octave-forge/main/control/inst/isstabilizable.m 2010-09-13 10:33:45 UTC (rev 7706) +++ trunk/octave-forge/main/control/inst/isstabilizable.m 2010-09-13 11:08:39 UTC (rev 7707) @@ -66,7 +66,7 @@ ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 -## Version: 0.2 +## Version: 0.2.1 function bool = isstabilizable (a, b = [], tol = [], dflg = 0) @@ -77,7 +77,7 @@ print_usage (); endif tol = b; - dflg = isdt (a); + dflg = ! isct (a); [a, b] = ssdata (a); elseif (nargin == 1) # isstabilizable (a, b, ...) print_usage (); @@ -85,8 +85,10 @@ error ("isstabilizable: a must be square and conformal to b") endif - if (isempty (tol)) - tol = 0; + if (! isreal (tol) && ! isscalar (tol)) + error ("isstabilizable: tol must be a real scalar"); + elseif (isempty (tol)) + tol = 0; # default tolerance endif ## controllability staircase form Modified: trunk/octave-forge/main/control/src/slab01od.cc =================================================================== --- trunk/octave-forge/main/control/src/slab01od.cc 2010-09-13 10:33:45 UTC (rev 7706) +++ trunk/octave-forge/main/control/src/slab01od.cc 2010-09-13 11:08:39 UTC (rev 7707) @@ -117,7 +117,12 @@ if (info != 0) error ("slab01od: AB01OD returned info = %d", info); - + + // resize + a.resize (n, n); + b.resize (n, m); + u.resize (n, n); + // return values retval(0) = a; retval(1) = b; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-09-13 16:31:29
|
Revision: 7708 http://octave.svn.sourceforge.net/octave/?rev=7708&view=rev Author: paramaniac Date: 2010-09-13 16:31:22 +0000 (Mon, 13 Sep 2010) Log Message: ----------- control: fix segmentation fault in minreal for state-space models Modified Paths: -------------- trunk/octave-forge/main/control/inst/@ss/__minreal__.m trunk/octave-forge/main/control/inst/ltimodels.m trunk/octave-forge/main/control/src/sltb01pd.cc Removed Paths: ------------- trunk/octave-forge/main/control/inst/test_minreal.m Modified: trunk/octave-forge/main/control/inst/@ss/__minreal__.m =================================================================== --- trunk/octave-forge/main/control/inst/@ss/__minreal__.m 2010-09-13 11:08:39 UTC (rev 7707) +++ trunk/octave-forge/main/control/inst/@ss/__minreal__.m 2010-09-13 16:31:22 UTC (rev 7708) @@ -22,7 +22,7 @@ ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 ## Version: 0.2 -%{ + function retsys = __minreal__ (sys, tol) if (tol == "def") @@ -37,62 +37,3 @@ retsys.lti = sys.lti; # retain i/o names and tsam endfunction -%} - -## ============================================================================== -## TODO: Fix frequent segfaults in sltb01pd.oct, afterwards delete the code below -## ============================================================================== - -## Algorithm based on sysmin by A. Scottedward Hodel -## Author: Lukas Reichlin <luk...@gm...> -## Created: October 2009 -## Version: 0.1 - -function retsys = __minreal__ (sys, tol) - - A = sys.a; - B = sys.b; - C = sys.c; - - if (! isempty (A)) - if (tol == "def") - [cflg, Uc] = isctrb (A, B); - else - [cflg, Uc] = isctrb (A, B, tol); - endif - - if (! cflg) - if (! isempty (Uc)) - A = Uc.' * A * Uc; - B = Uc.' * B; - C = C * Uc; - else - A = B = C = []; - endif - endif - endif - - if (! isempty (A)) - if (tol == "def") - [oflg, Uo] = isobsv (A, C); - else - [oflg, Uo] = isobsv (A, C, tol); - endif - - if (! oflg) - if (! isempty (Uo)) - A = Uo.' * A * Uo; - B = Uo.' * B; - C = C * Uo; - else - A = B = C = []; - endif - endif - endif - - retsys = ss (A, B, C, sys.d); - retsys.lti = sys.lti; # retain i/o names and tsam - -endfunction - -## ============================================================================== Modified: trunk/octave-forge/main/control/inst/ltimodels.m =================================================================== --- trunk/octave-forge/main/control/inst/ltimodels.m 2010-09-13 11:08:39 UTC (rev 7707) +++ trunk/octave-forge/main/control/inst/ltimodels.m 2010-09-13 16:31:22 UTC (rev 7708) @@ -23,9 +23,9 @@ ## Test suite and help for LTI models. ## @end deftypefn -## Author: Luca Favatella <sla...@gm...> +## Author: Lukas Reichlin <luk...@gm...> ## Created: November 2009 -## Version: 0.1 +## Version: 0.2 function ltimodels (systype = "general") @@ -197,7 +197,7 @@ %!assert (z, z_exp, 1e-4); -## ss: minreal +## ss: minreal (SLICOT TB01PD) %!shared C, D %! %! A = ss (-2, 3, 4, 5); @@ -210,45 +210,43 @@ %!assert (C.c, D.c); %!assert (C.d, D.d); +%!shared M, Me +%! A = [ 1.0 2.0 0.0 +%! 4.0 -1.0 0.0 +%! 0.0 0.0 1.0 ]; +%! +%! B = [ 1.0 +%! 0.0 +%! 1.0 ]; +%! +%! C = [ 0.0 1.0 -1.0 +%! 0.0 0.0 1.0 ]; +%! +%! D = zeros (2, 1); +%! +%! sys = ss (A, B, C, D); +%! sysmin = minreal (sys, 0.0); +%! [Ar, Br, Cr, Dr] = ssdata (sysmin); +%! M = [Ar, Br; Cr, Dr]; +%! +%! Ae = [ 1.0000 -1.4142 1.4142 +%! -2.8284 0.0000 1.0000 +%! 2.8284 1.0000 0.0000 ]; +%! +%! Be = [-1.0000 +%! 0.7071 +%! 0.7071 ]; +%! +%! Ce = [ 0.0000 0.0000 -1.4142 +%! 0.0000 0.7071 0.7071 ]; +%! +%! De = zeros (2, 1); +%! +%! Me = [Ae, Be; Ce, De]; +%! +%!assert (M, Me, 1e-4); -## ss: minreal (SLICOT TB01PD) -#%!shared M, Me -#%! A = [ 1.0 2.0 0.0 -#%! 4.0 -1.0 0.0 -#%! 0.0 0.0 1.0 ]; -#%! -#%! B = [ 1.0 -#%! 0.0 -#%! 1.0 ]; -#%! -#%! C = [ 0.0 1.0 -1.0 -#%! 0.0 0.0 1.0 ]; -#%! -#%! D = zeros (2, 1); -#%! -#%! sys = ss (A, B, C, D); -#%! sysmin = minreal (sys, 0.0); -#%! [Ar, Br, Cr, Dr] = ssdata (sysmin); -#%! M = [Ar, Br; Cr, Dr]; -#%! -#%! Ae = [ 1.0000 -1.4142 1.4142 -#%! -2.8284 0.0000 1.0000 -#%! 2.8284 1.0000 0.0000 ]; -#%! -#%! Be = [-1.0000 -#%! 0.7071 -#%! 0.7071 ]; -#%! -#%! Ce = [ 0.0000 0.0000 -1.4142 -#%! 0.0000 0.7071 0.7071 ]; -#%! -#%! De = zeros (2, 1); -#%! -#%! Me = [Ae, Be; Ce, De]; -#%! -#%!assert (M, Me, 1e-4); - ## ss: sminreal %!shared B, C %! Deleted: trunk/octave-forge/main/control/inst/test_minreal.m =================================================================== --- trunk/octave-forge/main/control/inst/test_minreal.m 2010-09-13 11:08:39 UTC (rev 7707) +++ trunk/octave-forge/main/control/inst/test_minreal.m 2010-09-13 16:31:22 UTC (rev 7708) @@ -1,58 +0,0 @@ -## ============================================================================== -## The tests below usually segfault after consecutive calls -## ============================================================================== - -test test_minreal - -## ss: minreal -%!shared Ar, Br, Cr -%! -%! A = [ -2 -2.4 -%! 0 -4.4 ]; -%! -%! B = [ 0.6 -%! 0.6 ]; -%! -%! C = [ 4 -4 ]; -%! -%! [Ar, Br, Cr, nr] = sltb01pd (A, B, C, 0); -%! -%!assert (isempty (Ar), true); -%!assert (isempty (Br), true); -%!assert (isempty (Cr), true); - - -## ss: minreal (SLICOT TB01PD) -%!shared M, Me -%! A = [ 1.0 2.0 0.0 -%! 4.0 -1.0 0.0 -%! 0.0 0.0 1.0 ]; -%! -%! B = [ 1.0 -%! 0.0 -%! 1.0 ]; -%! -%! C = [ 0.0 1.0 -1.0 -%! 0.0 0.0 1.0 ]; -%! -%! D = zeros (2, 1); -%! -%! [Ar, Br, Cr, nr] = sltb01pd (A, B, C, 0.0); -%! M = [Ar, Br; Cr, D]; -%! -%! Ae = [ 1.0000 -1.4142 1.4142 -%! -2.8284 0.0000 1.0000 -%! 2.8284 1.0000 0.0000 ]; -%! -%! Be = [-1.0000 -%! 0.7071 -%! 0.7071 ]; -%! -%! Ce = [ 0.0000 0.0000 -1.4142 -%! 0.0000 0.7071 0.7071 ]; -%! -%! De = zeros (2, 1); -%! -%! Me = [Ae, Be; Ce, De]; -%! -%!assert (M, Me, 1e-4); \ No newline at end of file Modified: trunk/octave-forge/main/control/src/sltb01pd.cc =================================================================== --- trunk/octave-forge/main/control/src/sltb01pd.cc 2010-09-13 11:08:39 UTC (rev 7707) +++ trunk/octave-forge/main/control/src/sltb01pd.cc 2010-09-13 16:31:22 UTC (rev 7708) @@ -23,7 +23,7 @@ Author: Lukas Reichlin <luk...@gm...> Created: September 2010 -Version: 0.1 +Version: 0.2 */ @@ -79,6 +79,8 @@ else ldc = max (1, m, p); + b.resize (ldb, max (m, p)); + // arguments out int nr = 0; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-09-13 22:59:34
|
Revision: 7710 http://octave.svn.sourceforge.net/octave/?rev=7710&view=rev Author: paramaniac Date: 2010-09-13 22:59:27 +0000 (Mon, 13 Sep 2010) Log Message: ----------- control: various small enhancements to ARE solvers Modified Paths: -------------- trunk/octave-forge/main/control/inst/care.m trunk/octave-forge/main/control/inst/dare.m trunk/octave-forge/main/control/src/slsb02od.cc Modified: trunk/octave-forge/main/control/inst/care.m =================================================================== --- trunk/octave-forge/main/control/inst/care.m 2010-09-13 19:32:31 UTC (rev 7709) +++ trunk/octave-forge/main/control/inst/care.m 2010-09-13 22:59:27 UTC (rev 7710) @@ -25,25 +25,25 @@ ## @strong{Inputs} ## @table @var ## @item a -## Real matrix. +## Real matrix (n-by-n). ## @item b -## Real matrix. +## Real matrix (n-by-m). ## @item q -## Real matrix. +## Real matrix (n-by-n). ## @item r -## Real matrix. +## Real matrix (m-by-m). ## @item s -## Optional real matrix. If @var{s} is not specified, a zero matrix is assumed. +## Optional real matrix (n-by-m). If @var{s} is not specified, a zero matrix is assumed. ## @end table ## ## @strong{Outputs} ## @table @var ## @item x -## Unique stabilizing solution of the continuous-time Riccati equation. +## Unique stabilizing solution of the continuous-time Riccati equation (n-by-n). ## @item l -## Closed-loop poles. +## Closed-loop poles (n-by-1). ## @item g -## Corresponding gain matrix. +## Corresponding gain matrix (m-by-n). ## @end table ## ## @example @@ -68,39 +68,41 @@ ## Author: Lukas Reichlin <luk...@gm...> ## Created: November 2009 -## Version: 0.3 +## Version: 0.4 function [x, l, g] = care (a, b, q, r, s = []) ## TODO: Add SLICOT SG02AD (Solution of continuous- or discrete-time ## algebraic Riccati equations for descriptor systems) + ## TODO: extract feedback matrix g from SB02OD (and SG02AD) + if (nargin < 4 || nargin > 5) print_usage (); endif - if (! issquare (a)) - error ("care: a is not square"); + if (! isreal (a) || ! issquare (a)) + error ("care: a must be real and square"); endif + + if (! isreal (b) || rows (a) != rows (b)) + error ("care: b must be real and conformal to a"); + endif - if (! issquare (q)) - error ("care: q is not square"); + if (! isreal (q) || ! issquare (q)) + error ("care: q must be real and square"); endif - if (! issquare (r)) - error ("care: r is not square"); + if (! isreal (r) || ! issquare (r)) + error ("care: r must be real and square"); endif - if (rows (a) != rows (b)) - error ("care: (a, b) not conformable"); - endif - if (columns (r) != columns (b)) error ("care: (b, r) not conformable"); endif - if (! isempty (s) && any (size (s) != size (b))) - error ("care: s (%dx%d) must be identically dimensioned with b (%dx%d)", + if (! isempty (s) && (! isreal (s) || any (size (s) != size (b)))) + error ("care: s(%dx%d) must be real and identically dimensioned with b(%dx%d)", rows (s), columns (s), rows (b), columns (b)); endif @@ -124,24 +126,15 @@ ## solve the riccati equation if (isempty (s)) - ## unique stabilizing solution - x = slsb02od (a, b, q, r, b, false, false); - - ## corresponding gain matrix - g = r \ (b.'*x); + [x, l] = slsb02od (a, b, q, r, b, false, false); + + g = r \ (b.'*x); # gain matrix else - ## unique stabilizing solution - x = slsb02od (a, b, q, r, s, false, true); - - ## corresponding gain matrix - g = r \ (b.'*x + s.'); + [x, l] = slsb02od (a, b, q, r, s, false, true); + + g = r \ (b.'*x + s.'); # gain matrix endif - ## closed-loop poles - l = eig (a - b*g); - - ## TODO: use alphar, alphai and beta from SB02OD - endfunction Modified: trunk/octave-forge/main/control/inst/dare.m =================================================================== --- trunk/octave-forge/main/control/inst/dare.m 2010-09-13 19:32:31 UTC (rev 7709) +++ trunk/octave-forge/main/control/inst/dare.m 2010-09-13 22:59:27 UTC (rev 7710) @@ -25,25 +25,25 @@ ## @strong{Inputs} ## @table @var ## @item a -## Real matrix. +## Real matrix (n-by-n). ## @item b -## Real matrix. +## Real matrix (n-by-m). ## @item q -## Real matrix. +## Real matrix (n-by-n). ## @item r -## Real matrix. +## Real matrix (m-by-m). ## @item s -## Optional real matrix. If @var{s} is not specified, a zero matrix is assumed. +## Optional real matrix (n-by-m). If @var{s} is not specified, a zero matrix is assumed. ## @end table ## ## @strong{Outputs} ## @table @var ## @item x -## Unique stabilizing solution of the discrete-time Riccati equation. +## Unique stabilizing solution of the discrete-time Riccati equation (n-by-n). ## @item l -## Closed-loop poles. +## Closed-loop poles (n-by-1). ## @item g -## Corresponding gain matrix. +## Corresponding gain matrix (m-by-n). ## @end table ## ## @example @@ -68,39 +68,41 @@ ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 -## Version: 0.3 +## Version: 0.4 function [x, l, g] = dare (a, b, q, r, s = []) ## TODO: Add SLICOT SG02AD (Solution of continuous- or discrete-time ## algebraic Riccati equations for descriptor systems) + ## TODO: extract feedback matrix g from SB02OD (and SG02AD) + if (nargin < 4 || nargin > 5) print_usage (); endif - if (! issquare (a)) - error ("dare: a is not square"); + if (! isreal (a) || ! issquare (a)) + error ("dare: a must be real and square"); endif + + if (! isreal (b) || rows (a) != rows (b)) + error ("dare: b must be real and conformal to a"); + endif - if (! issquare (q)) - error ("dare: q is not square"); + if (! isreal (q) || ! issquare (q)) + error ("dare: q must be real and square"); endif - if (! issquare (r)) - error ("dare: r is not square"); + if (! isreal (r) || ! issquare (r)) + error ("dare: r must be real and square"); endif - if (rows (a) != rows (b)) - error ("dare: (a, b) not conformable"); - endif - if (columns (r) != columns (b)) error ("dare: (b, r) not conformable"); endif - if (! isempty (s) && any (size (s) != size (b))) - error ("dare: s (%dx%d) must be identically dimensioned with b (%dx%d)", + if (! isempty (s) && (! isreal (s) || any (size (s) != size (b)))) + error ("dare: s(%dx%d) must be real and identically dimensioned with b(%dx%d)", rows (s), columns (s), rows (b), columns (b)); endif @@ -124,24 +126,15 @@ ## solve the riccati equation if (isempty (s)) - ## unique stabilizing solution - x = slsb02od (a, b, q, r, b, true, false); + [x, l] = slsb02od (a, b, q, r, b, true, false); - ## corresponding gain matrix - g = (r + b.'*x*b) \ (b.'*x*a); + g = (r + b.'*x*b) \ (b.'*x*a); # gain matrix else - ## unique stabilizing solution - x = slsb02od (a, b, q, r, s, true, true); - - ## corresponding gain matrix - g = (r + b.'*x*b) \ (b.'*x*a + s.'); + [x, l] = slsb02od (a, b, q, r, s, true, true); + + g = (r + b.'*x*b) \ (b.'*x*a + s.'); # gain matrix endif - ## closed-loop poles - l = eig (a - b*g); - - ## TODO: use alphar, alphai and beta from SB02OD - endfunction @@ -167,5 +160,5 @@ %! ge = [ 0.4092 1.7283]; %! %!assert (x, xe, 1e-4); -%!assert (l, le, 1e-4); +%!assert (sort (l), sort (le), 1e-4); %!assert (g, ge, 1e-4); \ No newline at end of file Modified: trunk/octave-forge/main/control/src/slsb02od.cc =================================================================== --- trunk/octave-forge/main/control/src/slsb02od.cc 2010-09-13 19:32:31 UTC (rev 7709) +++ trunk/octave-forge/main/control/src/slsb02od.cc 2010-09-13 22:59:27 UTC (rev 7710) @@ -23,13 +23,14 @@ Author: Lukas Reichlin <luk...@gm...> Created: February 2010 -Version: 0.2 +Version: 0.2.1 */ #include <octave/oct.h> #include <f77-fcn.h> #include "common.cc" +#include <complex> extern "C" { @@ -56,12 +57,12 @@ bool* BWORK, int& INFO); } - + DEFUN_DLD (slsb02od, args, nargout, "Slicot SB02OD Release 5.0") { int nargin = args.length (); octave_value_list retval; - + if (nargin != 7) { print_usage (); @@ -75,7 +76,7 @@ char uplo = 'U'; char jobl; char sort = 'S'; - + Matrix a = args(0).matrix_value (); Matrix b = args(1).matrix_value (); Matrix q = args(2).matrix_value (); @@ -88,54 +89,55 @@ dico = 'C'; else dico = 'D'; - + if (ijobl == 0) jobl = 'Z'; else jobl = 'N'; - + int n = a.rows (); // n: number of states int m = b.columns (); // m: number of inputs int p = 0; // p: number of outputs, not used because FACT = 'N' - + int lda = max (1, n); int ldb = max (1, n); int ldq = max (1, n); int ldr = max (1, m); int ldl = max (1, n); - + // arguments out double rcond; - + int ldx = max (1, n); Matrix x (ldx, n); - + + int nu = 2*n; + ColumnVector alfar (nu); + ColumnVector alfai (nu); + ColumnVector beta (nu); + + int lds = max (1, 2*n + m); + Matrix s (lds, lds); + // unused output arguments - OCTAVE_LOCAL_BUFFER (double, alfar, 2*n); - OCTAVE_LOCAL_BUFFER (double, alfai, 2*n); - OCTAVE_LOCAL_BUFFER (double, beta, 2*n); - - int lds = max (1, 2*n + m); - OCTAVE_LOCAL_BUFFER (double, s, lds*lds); - int ldt = max (1, 2*n + m); OCTAVE_LOCAL_BUFFER (double, t, ldt * 2*n); - + int ldu = max (1, 2*n); OCTAVE_LOCAL_BUFFER (double, u, ldu * 2*n); - + // tolerance double tol = 0; // use default value - + // workspace int liwork = max (1, m, 2*n); OCTAVE_LOCAL_BUFFER (int, iwork, liwork); - + int ldwork = max (7*(2*n + 1) + 16, 16*n, 2*n + m, 3*m); OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); - + OCTAVE_LOCAL_BUFFER (bool, bwork, 2*n); - + // error indicator int info; @@ -153,9 +155,9 @@ l.fortran_vec (), ldl, rcond, x.fortran_vec (), ldx, - alfar, alfai, - beta, - s, lds, + alfar.fortran_vec (), alfai.fortran_vec (), + beta.fortran_vec (), + s.fortran_vec (), lds, t, ldt, u, ldu, tol, @@ -166,13 +168,30 @@ if (f77_exception_encountered) error ("are: slsb02od: exception in SLICOT subroutine SB02OD"); - + if (info != 0) error ("are: slsb02od: SB02OD returned info = %d", info); + + // assemble complex vector - adapted from DEFUN complex in data.cc + alfar.resize (n); + alfai.resize (n); + beta.resize (n); + + ColumnVector zeror (n); + ColumnVector zeroi (n); + + zeror = quotient (alfar, beta); + zeroi = quotient (alfai, beta); + ComplexColumnVector zero (n, Complex ()); + + for (octave_idx_type i = 0; i < n; i++) + zero.xelem (i) = Complex (zeror(i), zeroi(i)); + // return value retval(0) = x; + retval(1) = zero; } - + return retval; } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-09-14 10:59:45
|
Revision: 7713 http://octave.svn.sourceforge.net/octave/?rev=7713&view=rev Author: paramaniac Date: 2010-09-14 10:59:36 +0000 (Tue, 14 Sep 2010) Log Message: ----------- control: add lyapchol and dlyapchol solvers Modified Paths: -------------- trunk/octave-forge/main/control/INDEX trunk/octave-forge/main/control/inst/dlyap.m trunk/octave-forge/main/control/inst/lyap.m trunk/octave-forge/main/control/inst/test_control.m trunk/octave-forge/main/control/src/Makefile Added Paths: ----------- trunk/octave-forge/main/control/inst/dlyapchol.m trunk/octave-forge/main/control/inst/lyapchol.m trunk/octave-forge/main/control/src/SB03OD.f trunk/octave-forge/main/control/src/SG03BD.f trunk/octave-forge/main/control/src/SG03BU.f trunk/octave-forge/main/control/src/SG03BV.f trunk/octave-forge/main/control/src/SG03BW.f trunk/octave-forge/main/control/src/SG03BX.f trunk/octave-forge/main/control/src/SG03BY.f trunk/octave-forge/main/control/src/slsb03od.cc trunk/octave-forge/main/control/src/slsg03bd.cc Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2010-09-14 09:11:45 UTC (rev 7712) +++ trunk/octave-forge/main/control/INDEX 2010-09-14 10:59:36 UTC (rev 7713) @@ -92,7 +92,9 @@ care dare dlyap + dlyapchol lyap + lyapchol Octave-only isctrb isobsv Modified: trunk/octave-forge/main/control/inst/dlyap.m =================================================================== --- trunk/octave-forge/main/control/inst/dlyap.m 2010-09-14 09:11:45 UTC (rev 7712) +++ trunk/octave-forge/main/control/inst/dlyap.m 2010-09-14 10:59:36 UTC (rev 7713) @@ -32,6 +32,8 @@ ## AXA' - EXE' + B = 0 (Generalized Lyapunov Equation) ## @end group ## @end example +## +## @seealso{dlyapchol, lyap, lyapchol} ## @end deftypefn ## Author: Lukas Reichlin <luk...@gm...> Added: trunk/octave-forge/main/control/inst/dlyapchol.m =================================================================== --- trunk/octave-forge/main/control/inst/dlyapchol.m (rev 0) +++ trunk/octave-forge/main/control/inst/dlyapchol.m 2010-09-14 10:59:36 UTC (rev 7713) @@ -0,0 +1,95 @@ +## Copyright (C) 2010 Lukas F. Reichlin +## +## This file is part of LTI Syncope. +## +## LTI Syncope 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. +## +## LTI Syncope 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. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} {@var{u} =} dlyapchol (@var{a}, @var{b}) +## @deftypefnx{Function File} {@var{u} =} dlyapchol (@var{a}, @var{b}, @var{e}) +## Compute Cholesky factor of discrete-time Lyapunov equations. +## Uses SLICOT SB03OD and SG03BD by courtesy of NICONET e.V. +## <http://www.slicot.org> +## +## @example +## @group +## A U' U A' - X + B B' = 0 (Lyapunov Equation) +## +## A U' U A' - E U' U E' + B B' = 0 (Generalized Lyapunov Equation) +## @end group +## @end example +## +## @seealso{dlyap, lyap, lyapchol} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: January 2010 +## Version: 0.2 + +function [u, scale] = dlyapchol (a, b, e) + + switch (nargin) + case 2 + + if (! isreal (a) || isempty (a) || ! issquare (a)) + error ("dlyapchol: a must be real and square"); + endif + + if (! isreal (b) || isempty (b)) + error ("dlyapchol: b must be real and square") + endif + + if (rows (a) != rows (b)) + error ("dlyapchol: a and b must have the same number of rows"); + endif + + [u, scale] = slsb03od (a.', b.', true); + + ## NOTE: TRANS = 'T' not suitable because we need U' U, not U U' + + case 3 + + if (! isreal (a) || isempty (a) || ! issquare (a)) + error ("dlyapchol: a must be real and square"); + endif + + if (! isreal (e) || isempty (e) || ! issquare (e)) + error ("dlyapchol: e must be real and square"); + endif + + if (! isreal (b) || isempty (b)) + error ("dlyapchol: b must be real"); + endif + + if (rows (b) != rows (a) || rows (e) != rows (a)) + error ("dlyapchol: a, b, e not conformal"); + endif + + [u, scale] = slsg03bd (a.', e.', b.', true); + + ## NOTE: TRANS = 'T' not suitable because we need U' U, not U U' + + otherwise + print_usage (); + + endswitch + + if (scale < 1) + warning ("dlyapchol: solution scaled by %g to prevent overflow", scale); + endif + +endfunction + + +## TODO: add tests Modified: trunk/octave-forge/main/control/inst/lyap.m =================================================================== --- trunk/octave-forge/main/control/inst/lyap.m 2010-09-14 09:11:45 UTC (rev 7712) +++ trunk/octave-forge/main/control/inst/lyap.m 2010-09-14 10:59:36 UTC (rev 7713) @@ -32,6 +32,8 @@ ## AXE' + EXA' + B = 0 (Generalized Lyapunov Equation) ## @end group ## @end example +## +## @seealso{lyapchol, dlyap, dlyapchol} ## @end deftypefn ## Author: Lukas Reichlin <luk...@gm...> Added: trunk/octave-forge/main/control/inst/lyapchol.m =================================================================== --- trunk/octave-forge/main/control/inst/lyapchol.m (rev 0) +++ trunk/octave-forge/main/control/inst/lyapchol.m 2010-09-14 10:59:36 UTC (rev 7713) @@ -0,0 +1,143 @@ +## Copyright (C) 2010 Lukas F. Reichlin +## +## This file is part of LTI Syncope. +## +## LTI Syncope 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. +## +## LTI Syncope 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. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} {@var{u} =} lyapchol (@var{a}, @var{b}) +## @deftypefnx{Function File} {@var{u} =} lyapchol (@var{a}, @var{b}, @var{e}) +## Compute Cholesky factor of continuous-time Lyapunov equations. +## Uses SLICOT SB03OD and SG03BD by courtesy of NICONET e.V. +## <http://www.slicot.org> +## +## @example +## @group +## A U' U + U' U A' + B B' = 0 (Lyapunov Equation) +## +## A U' U E' + E U' U A' + B B' = 0 (Generalized Lyapunov Equation) +## @end group +## @end example +## +## @seealso{lyap, dlyap, dlyapchol} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: January 2010 +## Version: 0.2 + +function [u, scale] = lyapchol (a, b, e) + + switch (nargin) + case 2 + + if (! isreal (a) || isempty (a) || ! issquare (a)) + error ("lyapchol: a must be real and square"); + endif + + if (! isreal (b) || isempty (b)) + error ("lyapchol: b must be real and square") + endif + + if (rows (a) != rows (b)) + error ("lyapchol: a and b must have the same number of rows"); + endif + + [u, scale] = slsb03od (a.', b.', false); + + ## NOTE: TRANS = 'T' not suitable because we need U' U, not U U' + + case 3 + + if (! isreal (a) || isempty (a) || ! issquare (a)) + error ("lyapchol: a must be real and square"); + endif + + if (! isreal (e) || isempty (e) || ! issquare (e)) + error ("lyapchol: e must be real and square"); + endif + + if (! isreal (b) || isempty (b)) + error ("lyapchol: b must be real"); + endif + + if (rows (b) != rows (a) || rows (e) != rows (a)) + error ("lyapchol: a, b, e not conformal"); + endif + + [u, scale] = slsg03bd (a.', e.', b.', false); + + ## NOTE: TRANS = 'T' not suitable because we need U' U, not U U' + + otherwise + print_usage (); + + endswitch + + if (scale < 1) + warning ("lyapchol: solution scaled by %g to prevent overflow", scale); + endif + +endfunction + + +%!shared U, U_exp, X, X_exp +%! +%! A = [ -1.0 37.0 -12.0 -12.0 +%! -1.0 -10.0 0.0 4.0 +%! 2.0 -4.0 7.0 -6.0 +%! 2.0 2.0 7.0 -9.0 ].'; +%! +%! B = [ 1.0 2.5 1.0 3.5 +%! 0.0 1.0 0.0 1.0 +%! -1.0 -2.5 -1.0 -1.5 +%! 1.0 2.5 4.0 -5.5 +%! -1.0 -2.5 -4.0 3.5 ].'; +%! +%! U = lyapchol (A, B); +%! +%! X = U.' * U; # use lyap at home! +%! +%! U_exp = [ 1.0000 0.0000 0.0000 0.0000 +%! 3.0000 1.0000 0.0000 0.0000 +%! 2.0000 -1.0000 1.0000 0.0000 +%! -1.0000 1.0000 -2.0000 1.0000 ].'; +%! +%! X_exp = [ 1.0000 3.0000 2.0000 -1.0000 +%! 3.0000 10.0000 5.0000 -2.0000 +%! 2.0000 5.0000 6.0000 -5.0000 +%! -1.0000 -2.0000 -5.0000 7.0000 ]; +%! +%!assert (U, U_exp, 1e-4); +%!assert (X, X_exp, 1e-4); + +%!shared U, U_exp, X, X_exp +%! +%! A = [ -1.0 3.0 -4.0 +%! 0.0 5.0 -2.0 +%! -4.0 4.0 1.0 ].'; +%! +%! E = [ 2.0 1.0 3.0 +%! 2.0 0.0 1.0 +%! 4.0 5.0 1.0 ].'; +%! +%! B = [ 2.0 -1.0 7.0 ].'; +%! +%! U = lyapchol (A, B, E); +%! +%! U_exp = [ 1.6003 -0.4418 -0.1523 +%! 0.0000 0.6795 -0.2499 +%! 0.0000 0.0000 0.2041 ]; +%! +%!assert (U, U_exp, 1e-4); Modified: trunk/octave-forge/main/control/inst/test_control.m =================================================================== --- trunk/octave-forge/main/control/inst/test_control.m 2010-09-14 09:11:45 UTC (rev 7712) +++ trunk/octave-forge/main/control/inst/test_control.m 2010-09-14 10:59:36 UTC (rev 7713) @@ -41,6 +41,8 @@ test dlyap test gram test covar +test lyapchol +## test dlyapchol # TODO: add tests ## various oct-files test place Modified: trunk/octave-forge/main/control/src/Makefile =================================================================== --- trunk/octave-forge/main/control/src/Makefile 2010-09-14 09:11:45 UTC (rev 7712) +++ trunk/octave-forge/main/control/src/Makefile 2010-09-14 10:59:36 UTC (rev 7713) @@ -1,7 +1,7 @@ all: slab08nd.oct slab13dd.oct slsb10hd.oct slsb10ed.oct slab13bd.oct \ slsb01bd.oct slsb10fd.oct slsb10dd.oct slsb03md.oct slsb04md.oct \ slsb04qd.oct slsg03ad.oct slsb02od.oct slab13ad.oct slab01od.oct \ - sltb01pd.oct + sltb01pd.oct slsb03od.oct slsg03bd.oct # transmission zeros of state-space models slab08nd.oct: slab08nd.cc @@ -107,16 +107,29 @@ SB03OR.f SB03OY.f SB04PX.f MB04NY.f MB04OY.f \ SB03OV.f -# Staircase form using orthogonal transformations +# staircase form using orthogonal transformations slab01od.oct: slab01od.cc mkoctfile slab01od.cc \ AB01OD.f AB01ND.f MB03OY.f MB01PD.f MB01QD.f -# Minimal realization of state-space models +# minimal realization of state-space models sltb01pd.oct: sltb01pd.cc mkoctfile sltb01pd.cc \ TB01PD.f TB01XD.f TB01ID.f AB07MD.f TB01UD.f \ MB03OY.f MB01PD.f MB01QD.f +# Cholesky factor of Lyapunov equations +slsb03od.oct: slsb03od.cc + mkoctfile slsb03od.cc \ + SB03OD.f select.f SB03OU.f SB03OT.f MB04ND.f \ + MB04OD.f SB03OR.f SB03OY.f SB04PX.f MB04NY.f \ + MB04OY.f SB03OV.f + +# Cholesky factor of generalized Lyapunov equations +slsg03bd.oct: slsg03bd.cc + mkoctfile slsg03bd.cc \ + SG03BD.f SG03BV.f SG03BU.f SG03BW.f SG03BX.f \ + SG03BY.f MB02UU.f MB02UV.f + clean: rm *.o core octave-core *.oct *~ Added: trunk/octave-forge/main/control/src/SB03OD.f =================================================================== --- trunk/octave-forge/main/control/src/SB03OD.f (rev 0) +++ trunk/octave-forge/main/control/src/SB03OD.f 2010-09-14 10:59:36 UTC (rev 7713) @@ -0,0 +1,662 @@ + SUBROUTINE SB03OD( DICO, FACT, TRANS, N, M, A, LDA, Q, LDQ, B, + $ LDB, SCALE, WR, WI, DWORK, LDWORK, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2009 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To solve for X = op(U)'*op(U) either the stable non-negative +C definite continuous-time Lyapunov equation +C 2 +C op(A)'*X + X*op(A) = -scale *op(B)'*op(B) (1) +C +C or the convergent non-negative definite discrete-time Lyapunov +C equation +C 2 +C op(A)'*X*op(A) - X = -scale *op(B)'*op(B) (2) +C +C where op(K) = K or K' (i.e., the transpose of the matrix K), A is +C an N-by-N matrix, op(B) is an M-by-N matrix, U is an upper +C triangular matrix containing the Cholesky factor of the solution +C matrix X, X = op(U)'*op(U), and scale is an output scale factor, +C set less than or equal to 1 to avoid overflow in X. If matrix B +C has full rank then the solution matrix X will be positive-definite +C and hence the Cholesky factor U will be nonsingular, but if B is +C rank deficient then X may be only positive semi-definite and U +C will be singular. +C +C In the case of equation (1) the matrix A must be stable (that +C is, all the eigenvalues of A must have negative real parts), +C and for equation (2) the matrix A must be convergent (that is, +C all the eigenvalues of A must lie inside the unit circle). +C +C ARGUMENTS +C +C Mode Parameters +C +C DICO CHARACTER*1 +C Specifies the type of Lyapunov equation to be solved as +C follows: +C = 'C': Equation (1), continuous-time case; +C = 'D': Equation (2), discrete-time case. +C +C FACT CHARACTER*1 +C Specifies whether or not the real Schur factorization +C of the matrix A is supplied on entry, as follows: +C = 'F': On entry, A and Q contain the factors from the +C real Schur factorization of the matrix A; +C = 'N': The Schur factorization of A will be computed +C and the factors will be stored in A and Q. +C +C TRANS CHARACTER*1 +C Specifies the form of op(K) to be used, as follows: +C = 'N': op(K) = K (No transpose); +C = 'T': op(K) = K**T (Transpose). +C +C Input/Output Parameters +C +C N (input) INTEGER +C The order of the matrix A and the number of columns in +C matrix op(B). N >= 0. +C +C M (input) INTEGER +C The number of rows in matrix op(B). M >= 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C On entry, the leading N-by-N part of this array must +C contain the matrix A. If FACT = 'F', then A contains +C an upper quasi-triangular matrix S in Schur canonical +C form; the elements below the upper Hessenberg part of the +C array A are not referenced. +C On exit, the leading N-by-N upper Hessenberg part of this +C array contains the upper quasi-triangular matrix S in +C Schur canonical form from the Shur factorization of A. +C The contents of array A is not modified if FACT = 'F'. +C +C LDA INTEGER +C The leading dimension of array A. LDA >= MAX(1,N). +C +C Q (input or output) DOUBLE PRECISION array, dimension +C (LDQ,N) +C On entry, if FACT = 'F', then the leading N-by-N part of +C this array must contain the orthogonal matrix Q of the +C Schur factorization of A. +C Otherwise, Q need not be set on entry. +C On exit, the leading N-by-N part of this array contains +C the orthogonal matrix Q of the Schur factorization of A. +C The contents of array Q is not modified if FACT = 'F'. +C +C LDQ INTEGER +C The leading dimension of array Q. LDQ >= MAX(1,N). +C +C B (input/output) DOUBLE PRECISION array, dimension (LDB,N) +C if TRANS = 'N', and dimension (LDB,max(M,N)), if +C TRANS = 'T'. +C On entry, if TRANS = 'N', the leading M-by-N part of this +C array must contain the coefficient matrix B of the +C equation. +C On entry, if TRANS = 'T', the leading N-by-M part of this +C array must contain the coefficient matrix B of the +C equation. +C On exit, the leading N-by-N part of this array contains +C the upper triangular Cholesky factor U of the solution +C matrix X of the problem, X = op(U)'*op(U). +C If M = 0 and N > 0, then U is set to zero. +C +C LDB INTEGER +C The leading dimension of array B. +C LDB >= MAX(1,N,M), if TRANS = 'N'; +C LDB >= MAX(1,N), if TRANS = 'T'. +C +C SCALE (output) DOUBLE PRECISION +C The scale factor, scale, set less than or equal to 1 to +C prevent the solution overflowing. +C +C WR (output) DOUBLE PRECISION array, dimension (N) +C WI (output) DOUBLE PRECISION array, dimension (N) +C If FACT = 'N', and INFO >= 0 and INFO <= 2, WR and WI +C contain the real and imaginary parts, respectively, of +C the eigenvalues of A. +C If FACT = 'F', WR and WI are not referenced. +C +C Workspace +C +C DWORK DOUBLE PRECISION array, dimension (LDWORK) +C On exit, if INFO = 0, or INFO = 1, DWORK(1) returns the +C optimal value of LDWORK. +C +C LDWORK INTEGER +C The length of the array DWORK. +C If M > 0, LDWORK >= MAX(1,4*N + MIN(M,N)); +C If M = 0, LDWORK >= 1. +C For optimum performance LDWORK should sometimes be larger. +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C < 0: if INFO = -i, the i-th argument had an illegal +C value; +C = 1: if the Lyapunov equation is (nearly) singular +C (warning indicator); +C if DICO = 'C' this means that while the matrix A +C (or the factor S) has computed eigenvalues with +C negative real parts, it is only just stable in the +C sense that small perturbations in A can make one or +C more of the eigenvalues have a non-negative real +C part; +C if DICO = 'D' this means that while the matrix A +C (or the factor S) has computed eigenvalues inside +C the unit circle, it is nevertheless only just +C convergent, in the sense that small perturbations +C in A can make one or more of the eigenvalues lie +C outside the unit circle; +C perturbed values were used to solve the equation; +C = 2: if FACT = 'N' and DICO = 'C', but the matrix A is +C not stable (that is, one or more of the eigenvalues +C of A has a non-negative real part), or DICO = 'D', +C but the matrix A is not convergent (that is, one or +C more of the eigenvalues of A lies outside the unit +C circle); however, A will still have been factored +C and the eigenvalues of A returned in WR and WI. +C = 3: if FACT = 'F' and DICO = 'C', but the Schur factor S +C supplied in the array A is not stable (that is, one +C or more of the eigenvalues of S has a non-negative +C real part), or DICO = 'D', but the Schur factor S +C supplied in the array A is not convergent (that is, +C one or more of the eigenvalues of S lies outside the +C unit circle); +C = 4: if FACT = 'F' and the Schur factor S supplied in +C the array A has two or more consecutive non-zero +C elements on the first sub-diagonal, so that there is +C a block larger than 2-by-2 on the diagonal; +C = 5: if FACT = 'F' and the Schur factor S supplied in +C the array A has a 2-by-2 diagonal block with real +C eigenvalues instead of a complex conjugate pair; +C = 6: if FACT = 'N' and the LAPACK Library routine DGEES +C has failed to converge. This failure is not likely +C to occur. The matrix B will be unaltered but A will +C be destroyed. +C +C METHOD +C +C The method used by the routine is based on the Bartels and Stewart +C method [1], except that it finds the upper triangular matrix U +C directly without first finding X and without the need to form the +C normal matrix op(B)'*op(B). +C +C The Schur factorization of a square matrix A is given by +C +C A = QSQ', +C +C where Q is orthogonal and S is an N-by-N block upper triangular +C matrix with 1-by-1 and 2-by-2 blocks on its diagonal (which +C correspond to the eigenvalues of A). If A has already been +C factored prior to calling the routine however, then the factors +C Q and S may be supplied and the initial factorization omitted. +C +C If TRANS = 'N', the matrix B is factored as (QR factorization) +C _ _ _ _ _ +C B = P ( R ), M >= N, B = P ( R Z ), M < N, +C ( 0 ) +C _ _ +C where P is an M-by-M orthogonal matrix and R is a square upper +C _ _ _ _ _ +C triangular matrix. Then, the matrix B = RQ, or B = ( R Z )Q (if +C M < N) is factored as +C _ _ +C B = P ( R ), M >= N, B = P ( R Z ), M < N. +C +C If TRANS = 'T', the matrix B is factored as (RQ factorization) +C _ +C _ _ ( Z ) _ +C B = ( 0 R ) P, M >= N, B = ( _ ) P, M < N, +C ( R ) +C _ _ +C where P is an M-by-M orthogonal matrix and R is a square upper +C _ _ _ _ _ +C triangular matrix. Then, the matrix B = Q'R, or B = Q'( Z' R' )' +C (if M < N) is factored as +C _ _ +C B = ( R ) P, M >= N, B = ( Z ) P, M < N. +C ( R ) +C +C These factorizations are utilised to either transform the +C continuous-time Lyapunov equation to the canonical form +C 2 +C op(S)'*op(V)'*op(V) + op(V)'*op(V)*op(S) = -scale *op(F)'*op(F), +C +C or the discrete-time Lyapunov equation to the canonical form +C 2 +C op(S)'*op(V)'*op(V)*op(S) - op(V)'*op(V) = -scale *op(F)'*op(F), +C +C where V and F are upper triangular, and +C +C F = R, M >= N, F = ( R Z ), M < N, if TRANS = 'N'; +C ( 0 0 ) +C +C F = R, M >= N, F = ( 0 Z ), M < N, if TRANS = 'T'. +C ( 0 R ) +C +C The transformed equation is then solved for V, from which U is +C obtained via the QR factorization of V*Q', if TRANS = 'N', or +C via the RQ factorization of Q*V, if TRANS = 'T'. +C +C REFERENCES +C +C [1] Bartels, R.H. and Stewart, G.W. +C Solution of the matrix equation A'X + XB = C. +C Comm. A.C.M., 15, pp. 820-826, 1972. +C +C [2] Hammarling, S.J. +C Numerical solution of the stable, non-negative definite +C Lyapunov equation. +C IMA J. Num. Anal., 2, pp. 303-325, 1982. +C +C NUMERICAL ASPECTS +C 3 +C The algorithm requires 0(N ) operations and is backward stable. +C +C FURTHER COMMENTS +C +C The Lyapunov equation may be very ill-conditioned. In particular, +C if A is only just stable (or convergent) then the Lyapunov +C equation will be ill-conditioned. A symptom of ill-conditioning +C is "large" elements in U relative to those of A and B, or a +C "small" value for scale. A condition estimate can be computed +C using SLICOT Library routine SB03MD. +C +C SB03OD routine can be also used for solving "unstable" Lyapunov +C equations, i.e., when matrix A has all eigenvalues with positive +C real parts, if DICO = 'C', or with moduli greater than one, +C if DICO = 'D'. Specifically, one may solve for X = op(U)'*op(U) +C either the continuous-time Lyapunov equation +C 2 +C op(A)'*X + X*op(A) = scale *op(B)'*op(B), (3) +C +C or the discrete-time Lyapunov equation +C 2 +C op(A)'*X*op(A) - X = scale *op(B)'*op(B), (4) +C +C provided, for equation (3), the given matrix A is replaced by -A, +C or, for equation (4), the given matrices A and B are replaced by +C inv(A) and B*inv(A), if TRANS = 'N' (or inv(A)*B, if TRANS = 'T'), +C respectively. Although the inversion generally can rise numerical +C problems, in case of equation (4) it is expected that the matrix A +C is enough well-conditioned, having only eigenvalues with moduli +C greater than 1. However, if A is ill-conditioned, it could be +C preferable to use the more general SLICOT Lyapunov solver SB03MD. +C +C CONTRIBUTOR +C +C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Aug. 1997. +C Supersedes Release 2.0 routine SB03CD by Sven Hammarling, +C NAG Ltd, United Kingdom. +C +C REVISIONS +C +C Dec. 1997, April 1998, May 1998, May 1999, Oct. 2001 (V. Sima). +C March 2002 (A. Varga). +C +C KEYWORDS +C +C Lyapunov equation, orthogonal transformation, real Schur form, +C Sylvester equation. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. Scalar Arguments .. + CHARACTER DICO, FACT, TRANS + INTEGER INFO, LDA, LDB, LDQ, LDWORK, M, N + DOUBLE PRECISION SCALE +C .. Array Arguments .. + DOUBLE PRECISION A(LDA,*), B(LDB,*), DWORK(*), Q(LDQ,*), WI(*), + $ WR(*) +C .. Local Scalars .. + LOGICAL CONT, LTRANS, NOFACT + INTEGER I, IFAIL, INFORM, ITAU, J, JWORK, K, L, MINMN, + $ NE, SDIM, WRKOPT + DOUBLE PRECISION EMAX, TEMP +C .. Local Arrays .. + LOGICAL BWORK(1) +C .. External Functions .. + LOGICAL LSAME, SELECT + DOUBLE PRECISION DLAPY2 + EXTERNAL DLAPY2, LSAME, SELECT +C .. External Subroutines .. + EXTERNAL DCOPY, DGEES, DGEMM, DGEMV, DGEQRF, DGERQF, + $ DLACPY, DLASET, DTRMM, SB03OU, XERBLA +C .. Intrinsic Functions .. + INTRINSIC INT, MAX, MIN +C .. Executable Statements .. +C +C Test the input scalar arguments. +C + CONT = LSAME( DICO, 'C' ) + NOFACT = LSAME( FACT, 'N' ) + LTRANS = LSAME( TRANS, 'T' ) + MINMN = MIN( M, N ) +C + INFO = 0 + IF( .NOT.CONT .AND. .NOT.LSAME( DICO, 'D' ) ) THEN + INFO = -1 + ELSE IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN + INFO = -2 + ELSE IF( .NOT.LTRANS .AND. .NOT.LSAME( TRANS, 'N' ) ) THEN + INFO = -3 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( M.LT.0 ) THEN + INFO = -5 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN + INFO = -9 + ELSE IF( ( LDB.LT.MAX( 1, N ) ) .OR. + $ ( LDB.LT.MAX( 1, N, M ) .AND. .NOT.LTRANS ) ) THEN + INFO = -11 + ELSE IF( LDWORK.LT.1 .OR. ( M.GT.0 .AND. LDWORK.LT.4*N + MINMN ) ) + $ THEN + INFO = -16 + END IF +C + IF ( INFO.NE.0 ) THEN +C +C Error return. +C + CALL XERBLA( 'SB03OD', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + IF( MINMN.EQ.0 ) THEN + IF( M.EQ.0 ) + $ CALL DLASET( 'Full', N, N, ZERO, ZERO, B, LDB ) + SCALE = ONE + DWORK(1) = ONE + RETURN + END IF +C +C Start the solution. +C +C (Note: Comments in the code beginning "Workspace:" describe the +C minimal amount of real workspace needed at that point in the +C code, as well as the preferred amount for good performance. +C NB refers to the optimal block size for the immediately +C following subroutine, as returned by ILAENV.) +C + IF ( NOFACT ) THEN +C +C Find the Schur factorization of A, A = Q*S*Q'. +C Workspace: need 3*N; +C prefer larger. +C + CALL DGEES( 'Vectors', 'Not ordered', SELECT, N, A, LDA, SDIM, + $ WR, WI, Q, LDQ, DWORK, LDWORK, BWORK, INFORM ) + IF ( INFORM.NE.0 ) THEN + INFO = 6 + RETURN + END IF + WRKOPT = DWORK(1) +C +C Check the eigenvalues for stability. +C + IF ( CONT ) THEN + EMAX = WR(1) +C + DO 20 J = 2, N + IF ( WR(J).GT.EMAX ) + $ EMAX = WR(J) + 20 CONTINUE +C + ELSE + EMAX = DLAPY2( WR(1), WI(1) ) +C + DO 40 J = 2, N + TEMP = DLAPY2( WR(J), WI(J) ) + IF ( TEMP.GT.EMAX ) + $ EMAX = TEMP + 40 CONTINUE +C + END IF +C + IF ( ( CONT ) .AND. ( EMAX.GE.ZERO ) .OR. + $ ( .NOT.CONT ) .AND. ( EMAX.GE.ONE ) ) THEN + INFO = 2 + RETURN + END IF + ELSE + WRKOPT = 0 + END IF +C +C Perform the QR or RQ factorization of B, +C _ _ _ _ _ +C B = P ( R ), or B = P ( R Z ), if TRANS = 'N', or +C ( 0 ) +C _ +C _ _ ( Z ) _ +C B = ( 0 R ) P, or B = ( _ ) P, if TRANS = 'T'. +C ( R ) +C Workspace: need MIN(M,N) + N; +C prefer MIN(M,N) + N*NB. +C + ITAU = 1 + JWORK = ITAU + MINMN + IF ( LTRANS ) THEN + CALL DGERQF( N, M, B, LDB, DWORK(ITAU), DWORK(JWORK), + $ LDWORK-JWORK+1, IFAIL ) + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1, MINMN*N) + JWORK = ITAU +C +C Form in B +C _ _ _ _ _ _ +C B := Q'R, m >= n, B := Q'*( Z' R' )', m < n, with B an +C n-by-min(m,n) matrix. +C Use a BLAS 3 operation if enough workspace, and BLAS 2, +C _ +C otherwise: B is formed column by column. +C + IF ( LDWORK.GE.JWORK+MINMN*N-1 ) THEN + K = JWORK +C + DO 60 I = 1, MINMN + CALL DCOPY( N, Q(N-MINMN+I,1), LDQ, DWORK(K), 1 ) + K = K + N + 60 CONTINUE +C + CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Non-unit', + $ N, MINMN, ONE, B(N-MINMN+1,M-MINMN+1), LDB, + $ DWORK(JWORK), N ) + IF ( M.LT.N ) + $ CALL DGEMM( 'Transpose', 'No transpose', N, M, N-M, + $ ONE, Q, LDQ, B, LDB, ONE, DWORK(JWORK), N ) + CALL DLACPY( 'Full', N, MINMN, DWORK(JWORK), N, B, LDB ) + ELSE + NE = N - MINMN +C + DO 80 J = 1, MINMN + NE = NE + 1 + CALL DCOPY( NE, B(1,M-MINMN+J), 1, DWORK(JWORK), 1 ) + CALL DGEMV( 'Transpose', NE, N, ONE, Q, LDQ, + $ DWORK(JWORK), 1, ZERO, B(1,J), 1 ) + 80 CONTINUE +C + END IF + ELSE + CALL DGEQRF( M, N, B, LDB, DWORK(ITAU), DWORK(JWORK), + $ LDWORK-JWORK+1, IFAIL ) + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1, MINMN*N) + JWORK = ITAU +C +C Form in B +C _ _ _ _ _ _ +C B := RQ, m >= n, B := ( R Z )*Q, m < n, with B an +C min(m,n)-by-n matrix. +C Use a BLAS 3 operation if enough workspace, and BLAS 2, +C _ +C otherwise: B is formed row by row. +C + IF ( LDWORK.GE.JWORK+MINMN*N-1 ) THEN + CALL DLACPY( 'Full', MINMN, N, Q, LDQ, DWORK(JWORK), MINMN ) + CALL DTRMM( 'Left', 'Upper', 'No transpose', 'Non-unit', + $ MINMN, N, ONE, B, LDB, DWORK(JWORK), MINMN ) + IF ( M.LT.N ) + $ CALL DGEMM( 'No transpose', 'No transpose', M, N, N-M, + $ ONE, B(1,M+1), LDB, Q(M+1,1), LDQ, ONE, + $ DWORK(JWORK), MINMN ) + CALL DLACPY( 'Full', MINMN, N, DWORK(JWORK), MINMN, B, LDB ) + ELSE + NE = MINMN + MAX( 0, N-M ) +C + DO 100 J = 1, MINMN + CALL DCOPY( NE, B(J,J), LDB, DWORK(JWORK), 1 ) + CALL DGEMV( 'Transpose', NE, N, ONE, Q(J,1), LDQ, + $ DWORK(JWORK), 1, ZERO, B(J,1), LDB ) + NE = NE - 1 + 100 CONTINUE +C + END IF + END IF + JWORK = ITAU + MINMN +C +C Solve for U the transformed Lyapunov equation +C 2 _ _ +C op(S)'*op(U)'*op(U) + op(U)'*op(U)*op(S) = -scale *op(B)'*op(B), +C +C or +C 2 _ _ +C op(S)'*op(U)'*op(U)*op(S) - op(U)'*op(U) = -scale *op(B)'*op(B) +C +C Workspace: need MIN(M,N) + 4*N; +C prefer larger. +C + CALL SB03OU( .NOT.CONT, LTRANS, N, MINMN, A, LDA, B, LDB, + $ DWORK(ITAU), B, LDB, SCALE, DWORK(JWORK), + $ LDWORK-JWORK+1, INFO ) + IF ( INFO.GT.1 ) THEN + INFO = INFO + 1 + RETURN + END IF + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 ) + JWORK = ITAU +C +C Form U := U*Q' or U := Q*U in the array B. +C Use a BLAS 3 operation if enough workspace, and BLAS 2, otherwise. +C Workspace: need N; +C prefer N*N; +C + IF ( LDWORK.GE.JWORK+N*N-1 ) THEN + IF ( LTRANS ) THEN + CALL DLACPY( 'Full', N, N, Q, LDQ, DWORK(JWORK), N ) + CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Non-unit', N, + $ N, ONE, B, LDB, DWORK(JWORK), N ) + ELSE + K = JWORK +C + DO 120 I = 1, N + CALL DCOPY( N, Q(1,I), 1, DWORK(K), N ) + K = K + 1 + 120 CONTINUE +C + CALL DTRMM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, + $ N, ONE, B, LDB, DWORK(JWORK), N ) + END IF + CALL DLACPY( 'Full', N, N, DWORK(JWORK), N, B, LDB ) + WRKOPT = MAX( WRKOPT, JWORK + N*N - 1 ) + ELSE + IF ( LTRANS ) THEN +C +C U is formed column by column ( U := Q*U ). +C + DO 140 I = 1, N + CALL DCOPY( I, B(1,I), 1, DWORK(JWORK), 1 ) + CALL DGEMV( 'No transpose', N, I, ONE, Q, LDQ, + $ DWORK(JWORK), 1, ZERO, B(1,I), 1 ) + 140 CONTINUE + ELSE +C +C U is formed row by row ( U' := Q*U' ). +C + DO 160 I = 1, N + CALL DCOPY( N-I+1, B(I,I), LDB, DWORK(JWORK), 1 ) + CALL DGEMV( 'No transpose', N, N-I+1, ONE, Q(1,I), LDQ, + $ DWORK(JWORK), 1, ZERO, B(I,1), LDB ) + 160 CONTINUE + END IF + END IF +C +C Lastly find the QR or RQ factorization of U, overwriting on B, +C to give the required Cholesky factor. +C Workspace: need 2*N; +C prefer N + N*NB; +C + JWORK = ITAU + N + IF ( LTRANS ) THEN + CALL DGERQF( N, N, B, LDB, DWORK(ITAU), DWORK(JWORK), + $ LDWORK-JWORK+1, IFAIL ) + ELSE + CALL DGEQRF( N, N, B, LDB, DWORK(ITAU), DWORK(JWORK), + $ LDWORK-JWORK+1, IFAIL ) + END IF + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 ) +C +C Make the diagonal elements of U non-negative. +C + IF ( LTRANS ) THEN +C + DO 200 J = 1, N + IF ( B(J,J).LT.ZERO ) THEN +C + DO 180 I = 1, J + B(I,J) = -B(I,J) + 180 CONTINUE +C + END IF + 200 CONTINUE +C + ELSE + K = JWORK +C + DO 240 J = 1, N + DWORK(K) = B(J,J) + L = JWORK +C + DO 220 I = 1, J + IF ( DWORK(L).LT.ZERO ) B(I,J) = -B(I,J) + L = L + 1 + 220 CONTINUE +C + K = K + 1 + 240 CONTINUE + END IF +C + IF( N.GT.1 ) + $ CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO, B(2,1), LDB ) +C +C Set the optimal workspace. +C + DWORK(1) = WRKOPT +C + RETURN +C *** Last line of SB03OD *** + END Added: trunk/octave-forge/main/control/src/SG03BD.f =================================================================== --- trunk/octave-forge/main/control/src/SG03BD.f (rev 0) +++ trunk/octave-forge/main/control/src/SG03BD.f 2010-09-14 10:59:36 UTC (rev 7713) @@ -0,0 +1,814 @@ + SUBROUTINE SG03BD( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q, + $ LDQ, Z, LDZ, B, LDB, SCALE, ALPHAR, ALPHAI, + $ BETA, DWORK, LDWORK, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2009 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To compute the Cholesky factor U of the matrix X, +C +C T +C X = op(U) * op(U), +C +C which is the solution of either the generalized +C c-stable continuous-time Lyapunov equation +C +C T T +C op(A) * X * op(E) + op(E) * X * op(A) +C +C 2 T +C = - SCALE * op(B) * op(B), (1) +C +C or the generalized d-stable discrete-time Lyapunov equation +C +C T T +C op(A) * X * op(A) - op(E) * X * op(E) +C +C 2 T +C = - SCALE * op(B) * op(B), (2) +C +C without first finding X and without the need to form the matrix +C op(B)**T * op(B). +C +C op(K) is either K or K**T for K = A, B, E, U. A and E are N-by-N +C matrices, op(B) is an M-by-N matrix. The resulting matrix U is an +C N-by-N upper triangular matrix with non-negative entries on its +C main diagonal. SCALE is an output scale factor set to avoid +C overflow in U. +C +C In the continuous-time case (1) the pencil A - lambda * E must be +C c-stable (that is, all eigenvalues must have negative real parts). +C In the discrete-time case (2) the pencil A - lambda * E must be +C d-stable (that is, the moduli of all eigenvalues must be smaller +C than one). +C +C ARGUMENTS +C +C Mode Parameters +C +C DICO CHARACTER*1 +C Specifies which type of the equation is considered: +C = 'C': Continuous-time equation (1); +C = 'D': Discrete-time equation (2). +C +C FACT CHARACTER*1 +C Specifies whether the generalized real Schur +C factorization of the pencil A - lambda * E is supplied +C on entry or not: +C = 'N': Factorization is not supplied; +C = 'F': Factorization is supplied. +C +C TRANS CHARACTER*1 +C Specifies whether the transposed equation is to be solved +C or not: +C = 'N': op(A) = A, op(E) = E; +C = 'T': op(A) = A**T, op(E) = E**T. +C +C Input/Output Parameters +C +C N (input) INTEGER +C The order of the matrix A. N >= 0. +C +C M (input) INTEGER +C The number of rows in the matrix op(B). M >= 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C On entry, if FACT = 'F', then the leading N-by-N upper +C Hessenberg part of this array must contain the +C generalized Schur factor A_s of the matrix A (see +C definition (3) in section METHOD). A_s must be an upper +C quasitriangular matrix. The elements below the upper +C Hessenberg part of the array A are not referenced. +C If FACT = 'N', then the leading N-by-N part of this +C array must contain the matrix A. +C On exit, the leading N-by-N part of this array contains +C the generalized Schur factor A_s of the matrix A. (A_s is +C an upper quasitriangular matrix.) +C +C LDA INTEGER +C The leading dimension of the array A. LDA >= MAX(1,N). +C +C E (input/output) DOUBLE PRECISION array, dimension (LDE,N) +C On entry, if FACT = 'F', then the leading N-by-N upper +C triangular part of this array must contain the +C generalized Schur factor E_s of the matrix E (see +C definition (4) in section METHOD). The elements below the +C upper triangular part of the array E are not referenced. +C If FACT = 'N', then the leading N-by-N part of this +C array must contain the coefficient matrix E of the +C equation. +C On exit, the leading N-by-N part of this array contains +C the generalized Schur factor E_s of the matrix E. (E_s is +C an upper triangular matrix.) +C +C LDE INTEGER +C The leading dimension of the array E. LDE >= MAX(1,N). +C +C Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) +C On entry, if FACT = 'F', then the leading N-by-N part of +C this array must contain the orthogonal matrix Q from +C the generalized Schur factorization (see definitions (3) +C and (4) in section METHOD). +C If FACT = 'N', Q need not be set on entry. +C On exit, the leading N-by-N part of this array contains +C the orthogonal matrix Q from the generalized Schur +C factorization. +C +C LDQ INTEGER +C The leading dimension of the array Q. LDQ >= MAX(1,N). +C +C Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) +C On entry, if FACT = 'F', then the leading N-by-N part of +C this array must contain the orthogonal matrix Z from +C the generalized Schur factorization (see definitions (3) +C and (4) in section METHOD). +C If FACT = 'N', Z need not be set on entry. +C On exit, the leading N-by-N part of this array contains +C the orthogonal matrix Z from the generalized Schur +C factorization. +C +C LDZ INTEGER +C The leading dimension of the array Z. LDZ >= MAX(1,N). +C +C B (input/output) DOUBLE PRECISION array, dimension (LDB,N1) +C On entry, if TRANS = 'T', the leading N-by-M part of this +C array must contain the matrix B and N1 >= MAX(M,N). +C If TRANS = 'N', the leading M-by-N part of this array +C must contain the matrix B and N1 >= N. +C On exit, the leading N-by-N part of this array contains +C the Cholesky factor U of the solution matrix X of the +C problem, X = op(U)**T * op(U). +C If M = 0 and N > 0, then U is set to zero. +C +C LDB INTEGER +C The leading dimension of the array B. +C If TRANS = 'T', LDB >= MAX(1,N). +C If TRANS = 'N', LDB >= MAX(1,M,N). +C +C SCALE (output) DOUBLE PRECISION +C The scale factor set to avoid overflow in U. +C 0 < SCALE <= 1. +C +C ALPHAR (output) DOUBLE PRECISION array, dimension (N) +C ALPHAI (output) DOUBLE PRECISION array, dimension (N) +C BETA (output) DOUBLE PRECISION array, dimension (N) +C If INFO = 0, 3, 5, 6, or 7, then +C (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, are the +C eigenvalues of the matrix pencil A - lambda * E. +C +C Workspace +C +C DWORK DOUBLE PRECISION array, dimension (LDWORK) +C On exit, if INFO = 0, DWORK(1) returns the optimal value +C of LDWORK. +C +C LDWORK INTEGER +C The dimension of the array DWORK. +C LDWORK >= MAX(1,4*N,6*N-6), if FACT = 'N'; +C LDWORK >= MAX(1,2*N,6*N-6), if FACT = 'F'. +C For good performance, LDWORK should be larger. +C +C Error indicator +C +C INFO INTEGER +C = 0: successful exit; +C < 0: if INFO = -i, the i-th argument had an illegal +C value; +C = 1: the pencil A - lambda * E is (nearly) singular; +C perturbed values were used to solve the equation +C (but the reduced (quasi)triangular matrices A and E +C are unchanged); +C = 2: FACT = 'F' and the matrix contained in the upper +C Hessenberg part of the array A is not in upper +C quasitriangular form; +C = 3: FACT = 'F' and there is a 2-by-2 block on the main +C diagonal of the pencil A_s - lambda * E_s whose +C eigenvalues are not conjugate complex; +C = 4: FACT = 'N' and the pencil A - lambda * E cannot be +C reduced to generalized Schur form: LAPACK routine +C DGEGS has failed to converge; +C = 5: DICO = 'C' and the pencil A - lambda * E is not +C c-stable; +C = 6: DICO = 'D' and the pencil A - lambda * E is not +C d-stable; +C = 7: the LAPACK routine DSYEVX utilized to factorize M3 +C failed to converge in the discrete-time case (see +C section METHOD for SLICOT Library routine SG03BU). +C This error is unlikely to occur. +C +C METHOD +C +C An extension [2] of Hammarling's method [1] to generalized +C Lyapunov equations is utilized to solve (1) or (2). +C +C First the pencil A - lambda * E is reduced to real generalized +C Schur form A_s - lambda * E_s by means of orthogonal +C transformations (QZ-algorithm): +C +C A_s = Q**T * A * Z (upper quasitriangular) (3) +C +C E_s = Q**T * E * Z (upper triangular). (4) +C +C If the pencil A - lambda * E has already been factorized prior to +C calling the routine however, then the factors A_s, E_s, Q and Z +C may be supplied and the initial factorization omitted. +C +C Depending on the parameters TRANS and M the N-by-N upper +C triangular matrix B_s is defined as follows. In any case Q_B is +C an M-by-M orthogonal matrix, which need not be accumulated. +C +C 1. If TRANS = 'N' and M < N, B_s is the upper triangular matrix +C from the QR-factorization +C +C ( Q_B O ) ( B * Z ) +C ( ) * B_s = ( ), +C ( O I ) ( O ) +C +C where the O's are zero matrices of proper size and I is the +C identity matrix of order N-M. +C +C 2. If TRANS = 'N' and M >= N, B_s is the upper triangular matrix +C from the (rectangular) QR-factorization +C +C ( B_s ) +C Q_B * ( ) = B * Z, +C ( O ) +C +C where O is the (M-N)-by-N zero matrix. +C +C 3. If TRANS = 'T' and M < N, B_s is the upper triangular matrix +C from the RQ-factorization +C +C ( Q_B O ) +C (B_s O ) * ( ) = ( Q**T * B O ). +C ( O I ) +C +C 4. If TRANS = 'T' and M >= N, B_s is the upper triangular matrix +C from the (rectangular) RQ-factorization +C +C ( B_s O ) * Q_B = Q**T * B, +C +C where O is the N-by-(M-N) zero matrix. +C +C Assuming SCALE = 1, the transformation of A, E and B described +C above leads to the reduced continuous-time equation +C +C T T +C op(A_s) op(U_s) op(U_s) op(E_s) +C +C T T +C + op(E_s) op(U_s) op(U_s) op(A_s) +C +C T +C = - op(B_s) op(B_s) (5) +C +C or to the reduced discrete-time equation +C +C T T +C op(A_s) op(U_s) op(U_s) op(A_s) +C +C T T +C - op(E_s) op(U_s) op(U_s) op(E_s) +C +C T +C = - op(B_s) op(B_s). (6) +C +C For brevity we restrict ourself to equation (5) and the case +C TRANS = 'N'. The other three cases can be treated in a similar +C fashion. +C +C We use the following partitioning for the matrices A_s, E_s, B_s +C and U_s +C +C ( A11 A12 ) ( E11 E12 ) +C A_s = ( ), E_s = ( ), +C ( 0 A22 ) ( 0 E22 ) +C +C ( B11 B12 ) ( U11 U12 ) +C B_s = ( ), U_s = ( ). (7) +C ( 0 B22 ) ( 0 U22 ) +C +C The size of the (1,1)-blocks is 1-by-1 (iff A_s(2,1) = 0.0) or +C 2-by-2. +C +C We compute U11 and U12**T in three steps. +C +C Step I: +C +C From (5) and (7) we get the 1-by-1 or 2-by-2 equation +C +C T T T T +C A11 * U11 * U11 * E11 + E11 * U11 * U11 * A11 +C +C T +C = - B11 * B11. +C +C For brevity, details are omitted here. See [2]. The technique +C for computing U11 is similar to those applied to standard +C Lyapunov equations in Hammarling's algorithm ([1], section 6). +C +C Furthermore, the auxiliary matrices M1 and M2 defined as +C follows +C +C -1 -1 +C M1 = U11 * A11 * E11 * U11 +C +C -1 -1 +C M2 = B11 * E11 * U11 +C +C are computed in a numerically reliable way. +C +C Step II: +C +C The generalized Sylvester equation +C +C T T T T +C A22 * U12 + E22 * U12 * M1 = +C +C T T T T T +C - B12 * M2 - A12 * U11 - E12 * U11 * M1 +C +C is solved for U12**T. +C +C Step III: +C +C It can be shown that +C +C T T T T +C A22 * U22 * U22 * E22 + E22 * U22 * U22 * A22 = +C +C T T +C - B22 * B22 - y * y (8) +C +C holds, where y is defined as +C +C T T T T T T +C y = B12 - ( E12 * U11 + E22 * U12 ) * M2 . +C +C If B22_tilde is the square triangular matrix arising from the +C (rectangular) QR-factorization +C +C ( B22_tilde ) ( B22 ) +C Q_B_tilde * ( ) = ( ), +C ( O ) ( y**T ) +C +C where Q_B_tilde is an orthogonal matrix of order N, then +C +C T T T +C - B22 * B22 - y * y = - B22_tilde * B22_tilde. +C +C Replacing the right hand side in (8) by the term +C - B22_tilde**T * B22_tilde leads to a reduced generalized +C Lyapunov equation of lower dimension compared to (5). +C +C The recursive application of the steps I to III yields the +C solution U_s of the equation (5). +C +C It remains to compute the solution matrix U of the original +C problem (1) or (2) from the matrix U_s. To this end we transform +C the solution back (with respect to the transformation that led +C from (1) to (5) (from (2) to (6)) and apply the QR-factorization +C (RQ-factorization). The upper triangular solution matrix U is +C obtained by +C +C Q_U * U = U_s * Q**T (if TRANS = 'N') +C +C or +C +C U * Q_U = Z * U_s (if TRANS = 'T') +C +C where Q_U is an N-by-N orthogonal matrix. Again, the orthogonal +C matrix Q_U need not be accumulated. +C +C REFERENCES +C +C [1] Hammarling, S.J. +C Numerical solution of the stable, non-negative definite +C Lyapunov equation. +C IMA J. Num. Anal., 2, pp. 303-323, 1982. +C +C [2] Penzl, T. +C Numerical solution of generalized Lyapunov equations. +C Advances in Comp. Math., vol. 8, pp. 33-48, 1998. +C +C NUMERICAL ASPECTS +C +C The number of flops required by the routine is given by the +C following table. Note that we count a single floating point +C arithmetic operation as one flop. +C +C | FACT = 'F' FACT = 'N' +C ---------+-------------------------------------------------- +C M <= N | (13*N**3+6*M*N**2 (211*N**3+6*M*N**2 +C | +6*M**2*N-2*M**3)/3 +6*M**2*N-2*M**3)/3 +C | +C M > N | (11*N**3+12*M*N**2)/3 (209*N**3+12*M*N**2)/3 +C +C FURTHER COMMENTS +C +C The Lyapunov equation may be very ill-conditioned. In particular, +C if DICO = 'D' and the pencil A - lambda * E has a pair of almost +C reciprocal eigenvalues, or DICO = 'C' and the pencil has an almost +C degenerate pair of eigenvalues, then the Lyapunov equation will be +C ill-conditioned. Perturbed values were used to solve the equation. +C A condition estimate can be obtained from the routine SG03AD. +C When setting the error indicator INFO, the routine does not test +C for near instability in the equation but only for exact +C instability. +C +C CONTRIBUTOR +C +C T. Penzl, Technical University Chemnitz, Germany, Aug. 1998. +C +C REVISIONS +C +C Sep. 1998 (V. Sima). +C May 1999 (V. Sima). +C March 2002 (A. Varga). +C Feb. 2004 (V. Sima). +C +C KEYWORDS +C +C Lyapunov equation +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION MONE, ONE, TWO, ZERO + PARAMETER ( MONE = -1.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, + $ ZERO = 0.0D+0 ) +C .. Scalar Arguments .. + DOUBLE PRECISION SCALE + INTEGER INFO, LDA, LDB, LDE, LDQ, LDWORK, LDZ, M, N + CHARACTER DICO, FACT, TRANS +C .. Array Arguments .. + DOUBLE PRECISION A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*), + $ BETA(*), DWORK(*), E(LDE,*), Q(LDQ,*), Z(LDZ,*) +C .. Local Scalars .. + DOUBLE PRECISION S1, S2, SAFMIN, WI, WR1, WR2 + INTEGER I, INFO1, MINMN, MINWRK, OPTWRK + LOGICAL ISDISC, ISFACT, ISTRAN +C .. Local Arrays .. + DOUBLE PRECISION E1(2,2) +C .. External Functions .. + DOUBLE PRECISION DLAMCH, DLAPY2 + LOGICAL LSAME + EXTERNAL DLAMCH, DLAPY2, LSAME +C .. External Subroutines .. + EXTERNAL DCOPY, DGEGS, DGEMM, DGEMV, DGEQRF, DGERQF, + $ DLACPY, DLAG2, DLASET, DSCAL, DTRMM, SG03BU, + $ SG03BV, XERBLA +C .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, INT, MAX, MIN, SIGN +C .. Executable Statements .. +C +C Decode input parameters. +C + ISDISC = LSAME( DICO, 'D' ) + ISFACT = LSAME( FACT, 'F' ) + ISTRAN = LSAME( TRANS, 'T' ) +C +C Compute minimal workspace. +C + IF (ISFACT ) THEN + MINWRK = MAX( 1, 2*N, 6*N-6 ) + ELSE + MINWRK = MAX( 1, 4*N, 6*N-6 ) + END IF +C +C Check the scalar input parameters. +C + IF ( .NOT.( ISDISC .OR. LSAME( DICO, 'C' ) ) ) THEN + INFO = -1 + ELSEIF ( .NOT.( ISFACT .OR. LSAME( FACT, 'N' ) ) ) THEN + INFO = -2 + ELSEIF ( .NOT.( ISTRAN .OR. LSAME( TRANS, 'N' ) ) ) THEN + INFO = -3 + ELSEIF ( N .LT. 0 ) THEN + INFO = -4 + ELSEIF ( M .LT. 0 ) THEN + INFO = -5 + ELSEIF ( LDA .LT. MAX( 1, N ) ) THEN + INFO = -7 + ELSEIF ( LDE .LT. MAX( 1, N ) ) THEN + INFO = -9 + ELSEIF ( LDQ .LT. MAX( 1, N ) ) THEN + INFO = -11 + ELSEIF ( LDZ .LT. MAX( 1, N ) ) THEN + INFO = -13 + ELSEIF ( ( ISTRAN .AND. ( LDB .LT. MAX( 1, N ) ) ) .OR. + $ ( .NOT.ISTRAN .AND. ( LDB .LT. MAX( 1, M, N ) ) ) ) THEN + INFO = -15 + ELSEIF ( LDWORK .LT. MINWRK ) THEN + INFO = -21 + ELSE + INFO = 0 + END IF + IF ( INFO .NE. 0 ) THEN + CALL XERBLA( 'SG03BD', -INFO ) + RETURN + END IF +C + SCALE = ONE +C +C Quick return if possible. +C + MINMN = MIN( M, N ) + IF ( MINMN .EQ. 0 ) THEN + IF ( N.GT.0 ) + $ CALL DLASET( 'Full', N, N, ZERO, ZERO, B, LDB ) + DWORK(1) = ONE + RETURN + ENDIF +C + IF ( ISFACT ) THEN +C +C Make sure the upper Hessenberg part of A is quasitriangular. +C + DO 20 I = 1, N-2 + IF ( A(I+1,I).NE.ZERO .AND. A(I+2,I+1).NE.ZERO ) THEN + INFO = 2 + RETURN + END IF + 20 CONTINUE + END IF +C + IF ( .NOT.ISFACT ) THEN +C +C Reduce the pencil A - lambda * E to generalized Schur form. +C +C A := Q**T * A * Z (upper quasitriangular) +C E := Q**T * E * Z (upper triangular) +C +C ( Workspace: >= MAX(1,4*N) ) +C + CALL DGEGS( 'Vectors', 'Vectors', N, A, LDA, E, LDE, ALPHAR, + $ ALPHAI, BETA, Q, LDQ, Z, LDZ, DWORK, LDWORK, + $ INFO1 ) + IF... [truncated message content] |
From: <par...@us...> - 2010-09-14 14:52:02
|
Revision: 7718 http://octave.svn.sourceforge.net/octave/?rev=7718&view=rev Author: paramaniac Date: 2010-09-14 14:51:54 +0000 (Tue, 14 Sep 2010) Log Message: ----------- control: add some remarks Modified Paths: -------------- trunk/octave-forge/main/control/DESCRIPTION trunk/octave-forge/main/control/doc/slicot_readme.txt Modified: trunk/octave-forge/main/control/DESCRIPTION =================================================================== --- trunk/octave-forge/main/control/DESCRIPTION 2010-09-14 12:55:14 UTC (rev 7717) +++ trunk/octave-forge/main/control/DESCRIPTION 2010-09-14 14:51:54 UTC (rev 7718) @@ -1,10 +1,10 @@ Name: Control Version: 2.0.0 -Date: 2010-09-07 +Date: 2010-09-14 Author: Lukas Reichlin <luk...@gm...> Maintainer: Lukas Reichlin <luk...@gm...> Title: Control Systems -Description: Computer Aided Control Systems Design package for GNU Octave +Description: control systems package increasingly based on the SLICOT Fortran library Depends: octave (>= 3.2.3) Autoload: yes License: GPL version 3 or later Modified: trunk/octave-forge/main/control/doc/slicot_readme.txt =================================================================== --- trunk/octave-forge/main/control/doc/slicot_readme.txt 2010-09-14 12:55:14 UTC (rev 7717) +++ trunk/octave-forge/main/control/doc/slicot_readme.txt 2010-09-14 14:51:54 UTC (rev 7718) @@ -1,7 +1,19 @@ +Important +********* + +If you encounter any bugs in the octave control package, please +DONT bother the great people behind SLICOT. Most likely it's a +problem of the control package, not the proven SLICOT library. +Therefore, please send any reports to the octave mailing list +and the maintainer of the control package. +See the file packinfo/DESCRIPTION for the current maintainer +of the control package and his mail address. + + Acknowledgments *************** -LTI SYNCOPE - LTI SYstems aNd COntrol PackagE - is partially based +LTI SYNCOPE - LTI SYstems aNd COntrol PackagE - is ever-increasingly based on SLICOT subroutines by courtesy of the Numerics In COntrol NETwork NICONET e.V. <http://www.icm.tu-bs.de/NICONET/>. This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-09-15 11:11:02
|
Revision: 7728 http://octave.svn.sourceforge.net/octave/?rev=7728&view=rev Author: paramaniac Date: 2010-09-15 11:10:55 +0000 (Wed, 15 Sep 2010) Log Message: ----------- control: relocate functions, improve argument checks Modified Paths: -------------- trunk/octave-forge/main/control/INDEX Added Paths: ----------- trunk/octave-forge/main/control/inst/@lti/dcgain.m trunk/octave-forge/main/control/inst/@lti/freqresp.m Removed Paths: ------------- trunk/octave-forge/main/control/inst/dcgain.m trunk/octave-forge/main/control/inst/freqresp.m Modified: trunk/octave-forge/main/control/INDEX =================================================================== --- trunk/octave-forge/main/control/INDEX 2010-09-15 09:59:03 UTC (rev 7727) +++ trunk/octave-forge/main/control/INDEX 2010-09-15 11:10:55 UTC (rev 7728) @@ -26,7 +26,7 @@ @lti/series strseq System Gain and Dynamics - dcgain + @lti/dcgain @lti/norm @lti/pole pzmap @@ -41,7 +41,7 @@ Frequency Domain Analysis bode bodemag - freqresp + @lti/freqresp margin nichols nyquist Added: trunk/octave-forge/main/control/inst/@lti/dcgain.m =================================================================== --- trunk/octave-forge/main/control/inst/@lti/dcgain.m (rev 0) +++ trunk/octave-forge/main/control/inst/@lti/dcgain.m 2010-09-15 11:10:55 UTC (rev 7728) @@ -0,0 +1,50 @@ +## Copyright (C) 2009 Lukas F. Reichlin +## +## This file is part of LTI Syncope. +## +## LTI Syncope 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. +## +## LTI Syncope 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. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{k} =} dcgain (@var{sys}) +## DC gain of LTI model. +## +## @strong{Inputs} +## @table @var +## @item sys +## LTI system. +## @end table +## +## @strong{Outputs} +## @table @var +## @item k +## DC gain matrice. For a system with m inputs and p outputs, the array @var{k} +## has dimensions [p, m]. +## @end table +## +## @seealso{freqresp} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function gain = dcgain (sys) + + if (nargin != 1) # sys is always an LTI model + print_usage (); + endif + + gain = __freqresp__ (sys, 0); + +endfunction \ No newline at end of file Added: trunk/octave-forge/main/control/inst/@lti/freqresp.m =================================================================== --- trunk/octave-forge/main/control/inst/@lti/freqresp.m (rev 0) +++ trunk/octave-forge/main/control/inst/@lti/freqresp.m 2010-09-15 11:10:55 UTC (rev 7728) @@ -0,0 +1,57 @@ +## Copyright (C) 2009 - 2010 Lukas F. Reichlin +## +## This file is part of LTI Syncope. +## +## LTI Syncope 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. +## +## LTI Syncope 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. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} {@var{H} =} freqresp (@var{sys}, @var{w}) +## Evaluate frequency response at given frequencies. +## +## @strong{Inputs} +## @table @var +## @item sys +## LTI system. +## @item w +## Vector of frequency values. +## @end table +## +## @strong{Outputs} +## @table @var +## @item H +## Array of frequency response. For a system with m inputs and p outputs, the array @var{H} +## has dimensions [p, m, length (w)]. +## The frequency response at the frequency w(k) is given by H(:,:,k). +## @end table +## +## @seealso{dcgain} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.2 + +function H = freqresp (sys, w) + + if (nargin != 2) # case freqresp () not possible + print_usage (); + endif + + if (! isreal (w) || ! isvector (w)) # catches freqresp (sys, sys) and freqresp (w, sys) as well + error ("freqresp: second argument must be a real vector"); + endif + + H = __freqresp__ (sys, w); + +endfunction \ No newline at end of file Deleted: trunk/octave-forge/main/control/inst/dcgain.m =================================================================== --- trunk/octave-forge/main/control/inst/dcgain.m 2010-09-15 09:59:03 UTC (rev 7727) +++ trunk/octave-forge/main/control/inst/dcgain.m 2010-09-15 11:10:55 UTC (rev 7728) @@ -1,50 +0,0 @@ -## Copyright (C) 2009 Lukas F. Reichlin -## -## This file is part of LTI Syncope. -## -## LTI Syncope 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. -## -## LTI Syncope 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. If not, see <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{k} =} dcgain (@var{sys}) -## DC gain of LTI model. -## -## @strong{Inputs} -## @table @var -## @item sys -## LTI system. -## @end table -## -## @strong{Outputs} -## @table @var -## @item k -## DC gain matrice. For a system with m inputs and p outputs, the array @var{k} -## has dimensions [p, m]. -## @end table -## -## @seealso{freqresp} -## @end deftypefn - -## Author: Lukas Reichlin <luk...@gm...> -## Created: October 2009 -## Version: 0.1 - -function gain = dcgain (sys) - - if (nargin != 1) - print_usage (); - endif - - gain = __freqresp__ (sys, 0); - -endfunction \ No newline at end of file Deleted: trunk/octave-forge/main/control/inst/freqresp.m =================================================================== --- trunk/octave-forge/main/control/inst/freqresp.m 2010-09-15 09:59:03 UTC (rev 7727) +++ trunk/octave-forge/main/control/inst/freqresp.m 2010-09-15 11:10:55 UTC (rev 7728) @@ -1,53 +0,0 @@ -## Copyright (C) 2009 Lukas F. Reichlin -## -## This file is part of LTI Syncope. -## -## LTI Syncope 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. -## -## LTI Syncope 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. If not, see <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn{Function File} {@var{H} =} freqresp (@var{sys}, @var{w}) -## Evaluate frequency response at given frequencies. -## -## @strong{Inputs} -## @table @var -## @item sys -## LTI system. -## @item w -## Vector of frequency values. -## @end table -## -## @strong{Outputs} -## @table @var -## @item H -## Array of frequency response. For a system with m inputs and p outputs, the array @var{H} -## has dimensions [p, m, length (w)]. -## The frequency response at the frequency w(k) is given by H(:,:,k). -## @end table -## -## @seealso{dcgain} -## @end deftypefn - -## Author: Lukas Reichlin <luk...@gm...> -## Created: October 2009 -## Version: 0.1 - -function H = freqresp (sys, w) - - if (nargin != 2) - print_usage (); - endif - - H = __freqresp__ (sys, w); - -endfunction \ No newline at end of file This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |