From: <par...@us...> - 2009-10-18 09:54:35
|
Revision: 6325 http://octave.svn.sourceforge.net/octave/?rev=6325&view=rev Author: paramaniac Date: 2009-10-18 09:54:21 +0000 (Sun, 18 Oct 2009) Log Message: ----------- add files Added Paths: ----------- trunk/octave-forge/extra/control-oo/inst/@lti/ trunk/octave-forge/extra/control-oo/inst/@lti/__ltigroup__.m trunk/octave-forge/extra/control-oo/inst/@lti/__ltiprune__.m trunk/octave-forge/extra/control-oo/inst/@lti/__propnames__.m trunk/octave-forge/extra/control-oo/inst/@lti/append.m trunk/octave-forge/extra/control-oo/inst/@lti/blkdiag.m trunk/octave-forge/extra/control-oo/inst/@lti/c2d.m trunk/octave-forge/extra/control-oo/inst/@lti/connect.m trunk/octave-forge/extra/control-oo/inst/@lti/display.m trunk/octave-forge/extra/control-oo/inst/@lti/eig.m trunk/octave-forge/extra/control-oo/inst/@lti/feedback.m trunk/octave-forge/extra/control-oo/inst/@lti/get.m trunk/octave-forge/extra/control-oo/inst/@lti/horzcat.m trunk/octave-forge/extra/control-oo/inst/@lti/inv.m trunk/octave-forge/extra/control-oo/inst/@lti/isct.m trunk/octave-forge/extra/control-oo/inst/@lti/isdt.m trunk/octave-forge/extra/control-oo/inst/@lti/issiso.m trunk/octave-forge/extra/control-oo/inst/@lti/isstable.m trunk/octave-forge/extra/control-oo/inst/@lti/lti.m trunk/octave-forge/extra/control-oo/inst/@lti/mconnect.m trunk/octave-forge/extra/control-oo/inst/@lti/minus.m trunk/octave-forge/extra/control-oo/inst/@lti/mldivide.m trunk/octave-forge/extra/control-oo/inst/@lti/mpower.m trunk/octave-forge/extra/control-oo/inst/@lti/mrdivide.m trunk/octave-forge/extra/control-oo/inst/@lti/mtimes.m trunk/octave-forge/extra/control-oo/inst/@lti/parallel.m trunk/octave-forge/extra/control-oo/inst/@lti/plus.m trunk/octave-forge/extra/control-oo/inst/@lti/pole.m trunk/octave-forge/extra/control-oo/inst/@lti/series.m trunk/octave-forge/extra/control-oo/inst/@lti/set.m trunk/octave-forge/extra/control-oo/inst/@lti/size.m trunk/octave-forge/extra/control-oo/inst/@lti/ssdata.m trunk/octave-forge/extra/control-oo/inst/@lti/subsasgn.m trunk/octave-forge/extra/control-oo/inst/@lti/subsref.m trunk/octave-forge/extra/control-oo/inst/@lti/tfdata.m trunk/octave-forge/extra/control-oo/inst/@lti/uminus.m trunk/octave-forge/extra/control-oo/inst/@lti/vertcat.m trunk/octave-forge/extra/control-oo/inst/@lti/zero.m trunk/octave-forge/extra/control-oo/inst/@ss/ trunk/octave-forge/extra/control-oo/inst/@ss/__c2d__.m trunk/octave-forge/extra/control-oo/inst/@ss/__get__.m trunk/octave-forge/extra/control-oo/inst/@ss/__getsysdata__.m trunk/octave-forge/extra/control-oo/inst/@ss/__pole__.m trunk/octave-forge/extra/control-oo/inst/@ss/__propnames__.m trunk/octave-forge/extra/control-oo/inst/@ss/__set__.m trunk/octave-forge/extra/control-oo/inst/@ss/__sys2tf__.m trunk/octave-forge/extra/control-oo/inst/@ss/__sysconnect__.m trunk/octave-forge/extra/control-oo/inst/@ss/__sysgroup__.m trunk/octave-forge/extra/control-oo/inst/@ss/__sysinv__.m trunk/octave-forge/extra/control-oo/inst/@ss/__sysprune__.m trunk/octave-forge/extra/control-oo/inst/@ss/__zero__.m trunk/octave-forge/extra/control-oo/inst/@ss/display.m trunk/octave-forge/extra/control-oo/inst/@ss/ss.m trunk/octave-forge/extra/control-oo/inst/@tf/ trunk/octave-forge/extra/control-oo/inst/@tf/__c2d__.m trunk/octave-forge/extra/control-oo/inst/@tf/__get__.m trunk/octave-forge/extra/control-oo/inst/@tf/__getsysdata__.m trunk/octave-forge/extra/control-oo/inst/@tf/__pole__.m trunk/octave-forge/extra/control-oo/inst/@tf/__propnames__.m trunk/octave-forge/extra/control-oo/inst/@tf/__set__.m trunk/octave-forge/extra/control-oo/inst/@tf/__sys2ss__.m trunk/octave-forge/extra/control-oo/inst/@tf/__sysconnect__.m trunk/octave-forge/extra/control-oo/inst/@tf/__sysgroup__.m trunk/octave-forge/extra/control-oo/inst/@tf/__sysinv__.m trunk/octave-forge/extra/control-oo/inst/@tf/__sysprune__.m trunk/octave-forge/extra/control-oo/inst/@tf/__zero__.m trunk/octave-forge/extra/control-oo/inst/@tf/display.m trunk/octave-forge/extra/control-oo/inst/@tf/tf.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/ trunk/octave-forge/extra/control-oo/inst/@tfpoly/__equalizer__.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/__remleadzer__.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/display.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/eq.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/get.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/length.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/minus.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/mpower.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/mtimes.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/ne.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/numel.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/plus.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/roots.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/subsref.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/tfpoly.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/tfpoly2str.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/uminus.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/uplus.m Added: trunk/octave-forge/extra/control-oo/inst/@lti/__ltigroup__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/__ltigroup__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/__ltigroup__.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,47 @@ +## 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 -*- +## Block diagonal concatenation of two LTI models. +## This file is part of the Model Abstraction Layer. +## For internal use only. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function retlti = __ltigroup__ (lti1, lti2) + + retlti = lti (); + + retlti.inname = [lti1.inname; + lti2.inname]; + + retlti.outname = [lti1.outname; + lti2.outname]; + + if (lti1.tsam == lti2.tsam) + retlti.tsam = lti1.tsam; + elseif (lti1.tsam == -1) + retlti.tsam = lti2.tsam; + elseif (lti2.tsam == -1) + retlti.tsam = lti1.tsam; + else + error ("ltigroup: systems must have identical sampling times"); + endif + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/__ltiprune__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/__ltiprune__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/__ltiprune__.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,32 @@ +## 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 -*- +## Submodel extraction and reordering for LTI objects. +## This file is part of the Model Abstraction Layer. +## For internal use only. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function lti = __ltiprune__ (lti, out_idx, in_idx) + + lti.inname = lti.inname(in_idx); + lti.outname = lti.outname(out_idx); + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/__propnames__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/__propnames__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/__propnames__.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,37 @@ +## Copyright (C) 2009 Lukas F. 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 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{propv}, @var{asgnvalv}] =} __propnames__ (@var{sys}) +## Return the list of properties as well as the assignable values for a lti object sys. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function [propv, asgnvalv] = __propnames__ (sys) + + ## cell vector of lti-specific properties + propv = {"tsam"; + "inname"; + "outname"}; + + ## cell vector of lti-specific assignable values + asgnvalv = {"scalar (sample time in seconds)"; + "m-by-1 cell vector of strings"; + "p-by-1 cell vector of strings"}; + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/append.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/append.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/append.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,37 @@ +## 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{sys} =} append (@var{sys1}, @var{sys2}) +## Group LTI models by appending their inputs and outputs. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function sys = append (varargin) + + sys = varargin{1}; + + if (nargin > 1) + for k = 2 : nargin + sys = __sysgroup__ (sys, varargin{k}); + endfor + endif + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/blkdiag.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/blkdiag.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/blkdiag.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,31 @@ +## 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{sys} =} blkdiag (@var{sys1}, @var{sys2}) +## Block-diagonal concatenation of LTI models. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function sys = blkdiag (varargin) + + sys = append (varargin{:}); + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/c2d.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/c2d.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/c2d.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,49 @@ +## 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{sys} =} c2d (@var{sys}, @var{tsam}) +## @deftypefnx {Function File} {@var{sys} =} c2d (@var{sys}, @var{tsam}, @var{method}) +## Convert the continuous lti model into its discrete time equivalent. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function sys = c2d (sys, tsam, method = "std") + + if (nargin < 2 || nargin > 3) + print_usage (); + endif + + if (! isa (sys, "lti")) + error ("c2d: first argument is not an lti model"); + endif + + if (! issample (tsam)) + error ("c2d: second argument is not a valid sample time"); + endif + + if (! ischar (method)) + error ("c2d: third argument is not a string"); + endif + + sys = __c2d__ (sys, tsam, method); + sys.tsam = tsam; + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/connect.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/connect.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/connect.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,59 @@ +## 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{sys} =} connect (@var{sys}, @var{cm}, @var{inputs}, @var{outputs}) +## Arbitrary interconnections between the inputs and outputs of an LTI model. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function sys = connect (sys, cm, in_idx, out_idx) + + if (nargin != 4) + print_usage (); + endif + + [p, m] = size (sys); + [cmrows, cmcols] = size (cm); + + ## TODO: proper argument checking + ## TODO: name-based interconnections + + if (! isreal (cm)) + error ("connect: second argument must be a matrix with real coefficients"); + endif + + cm = round (cm); + M = zeros (m, p); + in = cm(:, 1); + out = cm(:, 2 : cmcols); + + for a = 1 : cmrows + for b = 1 : cmcols-1 + if (out(a, b) != 0) + M(in(a, 1), abs (out(a, b))) = sign (out(a, b)); + endif + endfor + endfor + + sys = __sysconnect__ (sys, M); + sys = __sysprune__ (sys, out_idx, in_idx); + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/display.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/display.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/display.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,33 @@ +## 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 -*- +## Display routine for LTI objects. Called by its child classes. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function display (ltisys) + + tsam = ltisys.tsam; + + if (tsam > 0) + disp (sprintf("Sampling time: %g s", tsam)); + endif + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/eig.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/eig.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/eig.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,31 @@ +## 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{p} =} eig (@var{sys}) +## Compute poles of LTI system. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function pol = eig (sys) + + pol = pole (sys); + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/feedback.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/feedback.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/feedback.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,159 @@ +## 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{sys} =} feedback (@var{sys1}) +## @deftypefnx {Function File} {@var{sys} =} feedback (@var{sys1}, @var{"+"}) +## @deftypefnx {Function File} {@var{sys} =} feedback (@var{sys1}, @var{sys2}) +## @deftypefnx {Function File} {@var{sys} =} feedback (@var{sys1}, @var{sys2}, @var{"+"}) +## @deftypefnx {Function File} {@var{sys} =} feedback (@var{sys1}, @var{sys2}, @var{feedin}, @var{feedout}) +## @deftypefnx {Function File} {@var{sys} =} feedback (@var{sys1}, @var{sys2}, @var{feedin}, @var{feedout}, @var{"+"}) +## Feedback connection of two LTI models. +## @example +## @group +## u + +--------+ y +## ------>(+)----->| sys1 |-------+-------> +## ^ - +--------+ | +## | | +## | +--------+ | +## +-------| sys2 |<------+ +## +--------+ +## @end group +## @end example +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function sys = feedback (sys1, sys2_or_fbsign, feedin_or_fbsign, feedout, fbsign) + + [p1, m1] = size (sys1); + + switch (nargin) + case 1 # sys = feedback (sys) + if (p1 != m1) + error ("feedback: argument must be a square system"); + endif + + fbsign = -1; + sys2 = eye (p1); + feedin = 1 : m1; + feedout = 1 : p1; + + case 2 + if (ischar (sys2_or_fbsign)) # sys = feedback (sys, "+") + if (p1 != m1) + error ("feedback: argument must be a square system"); + endif + + fbsign = checkfbsign (sys2_or_fbsign); + sys2 = eye (p1); + feedin = 1 : m1; + feedout = 1 : p1; + else # sys = feedback (sys1, sys2) + fbsign = -1; + sys2 = sys2_or_fbsign; + feedin = 1 : m1; + feedout = 1 : p1; + endif + + case 3 # sys = feedback (sys1, sys2, "+") + fbsign = checkfbsign (feedin_or_fbsign); + sys2 = sys2_or_fbsign; + feedin = 1 : m1; + feedout = 1 : p1; + + case 4 # sys = feedback (sys1, sys2, feedin, feedout) + fbsign = -1; + sys2 = sys2_or_fbsign; + feedin = feedin_or_fbsign; + + case 5 # sys = feedback (sys1, sys2, feedin, feedout, "+") + fbsign = checkfbsign (fbsign); + sys2 = sys2_or_fbsign; + feedin = feedin_or_fbsign; + + otherwise + print_usage (); + + endswitch + + [p2, m2] = size (sys2); + + l_feedin = length (feedin); + l_feedout = length (feedout); + + if (l_feedin != p2) + error ("feedback: feedin indices: %d, outputs sys2: %d", l_feedin, p2); + endif + + if (l_feedout != m2) + error ("feedback: feedout indices: %d, inputs sys2: %d", l_feedout, m2); + endif + + if (any (feedin > m1 | feedin < 1)) + error ("feedback: range of feedin indices exceeds dimensions of sys1"); + endif + + if (any (feedin > p1 | feedin < 1)) + error ("feedback: range of feedout indices exceeds dimensions of sys1"); + endif + + M11 = zeros (m1, p1); + M12 = zeros (m1, p2); + M21 = zeros (m2, p1); + M22 = zeros (m2, p2); + + for k = 1 : l_feedin + M12(feedin(k), k) = fbsign; + endfor + + for k = 1 : l_feedout + M21(k, feedout(k)) = 1; + endfor + + M = [M11, M12; + M21, M22]; + + in_idx = 1 : m1; + out_idx = 1 : p1; + + sys = __sysgroup__ (sys1, sys2); + sys = __sysconnect__ (sys, M); + sys = __sysprune__ (sys, out_idx, in_idx); + +endfunction + + +function fbsign = checkfbsign (fbsign) + + if (isnumeric (fbsign)) + fbsign = sign (fbsign); + elseif (ischar (fbsign)) + if (fbsign == "+") + fbsign = +1; + elseif (fbsign == "-") + fbsign = -1; + else + error ("feedback: invalid feedback sign"); + endif + else + error ("feedback: invalid feedback sign"); + endif + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/get.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/get.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/get.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,54 @@ +## 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{value} =} get (@var{sys}, @var{"property"}) +## Access property values of LTI objects. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function varargout = get (sys, varargin) + + if (nargin == 1) + [propv, valv] = __propnames__ (sys); + nrows = numel (propv); + str = strjust (strvcat (propv), "right"); + str = horzcat (repmat (" ", nrows, 1), str, repmat (": ", nrows, 1), strvcat (valv)); + disp (str); + else + for k = 1 : (nargin-1) + prop = lower (varargin{k}); + + switch (prop) + case {"inname", "inputname"} + val = sys.inname; + case {"outname", "outputname"} + val = sys.outname; + case {"tsam", "ts"} + val = sys.tsam; + otherwise + val = __get__ (sys, prop); + endswitch + + varargout{k} = val; + endfor + endif + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/horzcat.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/horzcat.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/horzcat.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,49 @@ +## 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 -*- +## Horizontal concatenation of LTI objects. If necessary, object conversion +## is done by sysgroup. Used by the parser for "[lti1, lti2]". + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function sys = horzcat (sys, varargin) + + for k = 1 : (nargin-1) + + sys1 = sys; + sys2 = varargin{k}; + + [p1, m1] = size (sys1); + [p2, m2] = size (sys2); + + if (p1 != p2) + error ("lti: horzcat: number of system outputs incompatible: [(%dx%d), (%dx%d)]", + p1, m1, p2, m2); + end + + sys = __sysgroup__ (sys1, sys2); + + out_scl = [eye(p1), eye(p2)]; + + sys = out_scl * sys; + + endfor + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/inv.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/inv.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/inv.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,37 @@ +## 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 -*- +## Inversion of LTI objects. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function retsys = inv (sys) + + [p, m] = size (sys); + + if (p != m) + error ("lti: inv: system must be square"); + endif + + retsys = __sysinv__ (sys); + + ## TODO: handle i/o names + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/isct.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/isct.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/isct.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,35 @@ +## 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{bool} =} isct (@var{sys}) +## Determine whether LTI model is a continuous-time system. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function bool = isct (ltisys) + + if (nargin != 1) + print_usage (); + endif + + bool = (ltisys.tsam == 0 || ltisys.tsam == -1); + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/isdt.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/isdt.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/isdt.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,35 @@ +## 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{bool} =} isdt (@var{sys}) +## Determine whether LTI model is a discrete-time system. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function bool = isdt (ltisys) + + if (nargin != 1) + print_usage (); + endif + + bool = (ltisys.tsam != 0); + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/issiso.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/issiso.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/issiso.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,37 @@ +## 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{bool} =} issiso (@var{sys}) +## Determine whether LTI model is single-input/single-output (SISO). +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function bool = issiso (sys) + + if (nargin != 1) + print_usage (); + endif + + [p, m] = size (sys); + + bool = (p*m == 1); + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/isstable.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/isstable.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/isstable.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,42 @@ +## 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{bool} =} isstable (@var{sys}) +## @deftypefn {Function File} {@var{bool} =} isstable (@var{sys}, @var{tol}) +## Determine whether LTI system is stable. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function bool = isstable (sys, tol = 0) + + if (nargin == 0 || nargin > 2) + print_usage (); + endif + + eigw = pole (sys); + + if (isct (sys)) + bool = all (real (eigw) < -tol*(1 + abs (eigw))); + else + bool = all (abs (eigw) < 1 - tol); + endif + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/lti.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/lti.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/lti.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,36 @@ +## 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 -*- +## Constructor for LTI objects. For internal use only. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function ltisys = lti (ny = 0, nu = 0, tsam = -1) + + inname = repmat ({""}, nu, 1); + outname = repmat ({""}, ny, 1); + + ltisys = struct ("tsam", tsam, + "inname", {inname}, + "outname", {outname}); + + ltisys = class (ltisys, "lti"); + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/mconnect.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/mconnect.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/mconnect.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,63 @@ +## 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{sys} =} mconnect (@var{sys}, @var{m}) +## @deftypefnx {Function File} {@var{sys} =} mconnect (@var{sys}, @var{m}, @var{inputs}, @var{outputs}) +## Arbitrary interconnections between the inputs and outputs of an LTI model. +## @example +## @group +## Solve the system equations of +## y(t) = G e(t) +## e(t) = u(t) + M y(t) +## in order to build +## y(t) = H u(t) +## The matrix M for a (p-by-m) system G has m rows and p columns. +## @end group +## @end example +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function sys = mconnect (sys, M, in_idx, out_idx) + + if (nargin != 2 || nargin != 4) + print_usage (); + endif + + [p, m] = size (sys); + [mrows, mcols] = size (M); + + if (p != mcols || m != mrows) + error ("mconnect: second argument must be a (%dx&d) matrix", m, p); + endif + + if (! isreal (M)) + error ("mconnect: second argument must be a matrix with real coefficients"); + endif + + M = round (M); + + sys = __sysconnect__ (sys, M); + + if (nargin == 4) + sys = __sysprune__ (sys, out_idx, in_idx); + endif + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/minus.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/minus.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/minus.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,43 @@ +## 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 -*- +## Binary subtraction of LTI objects. If necessary, object conversion +## is done by sysgroup. Used by the parser for "lti1 - lti2". + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function sys = minus (sys1, sys2) + + [p1, m1] = size (sys1); + [p2, m2] = size (sys2); + + if (p1 != p2 || m1 != m2) + error ("lti: minus: system dimensions incompatible: (%dx%d) - (%dx%d)", + p1, m1, p2, m2); + end + + sys = __sysgroup__ (sys1, sys2); + + in_scl = [eye(m1); eye(m2)]; + out_scl = [eye(p1), -eye(p2)]; + + sys = out_scl * sys * in_scl; + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/mldivide.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/mldivide.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/mldivide.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,41 @@ +## 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 -*- +## Matrix left division of LTI objects. If necessary, object conversion +## is done by sysgroup in mtimes. Used by the parser for "lti1 \ lti2". + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function sys = mldivide (sys1, sys2) + + sys1 = inv (sys1); # let octave decide which inv() it uses + + [p1, m1] = size (sys1); + [p2, m2] = size (sys2); + + if (m2 != p1) + error ("lti: mldivide: system dimensions incompatible: (%dx%d) \ (%dx%d)", + p1, m1, p2, m2); + end + + sys = sys1 * sys2; + +endfunction + Added: trunk/octave-forge/extra/control-oo/inst/@lti/mpower.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/mpower.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/mpower.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,55 @@ +## 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 -*- +## Matrix power of LTI objects. The exponent must be an integer. +## Used by the parser for "lti^int". + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function retsys = mpower (sys, e) + + if (! isnumeric (e) || ! isscalar (e) || ! isreal (e) || e != round (e)) + error ("lti: mpower: exponent must be an integer"); + endif + + [p, m] = size (sys); + + if (p != m) + error ("lti: mpower: system must be square"); + endif + + ex = round (abs (e)); # make sure ex is a positive integer + + switch (sign (e)) + case -1 # lti^-ex + sys = inv (sys); + retsys = sys; + case 0 # lti^0 + retsys = eye (p); + return; + case 1 # lti^ex + retsys = sys; + endswitch + + for k = 2 : ex + retsys = retsys * sys; + endfor + +endfunction Added: trunk/octave-forge/extra/control-oo/inst/@lti/mrdivide.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/mrdivide.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/mrdivide.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,41 @@ +## 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 -*- +## Matrix right division of LTI objects. If necessary, object conversion +## is done by sysgroup in mtimes. Used by the parser for "lti1 / lti2". + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function sys = mrdivide (sys1, sys2) + + sys2 = inv (sys2); # let octave decide which inv() it uses + + [p1, m1] = size (sys1); + [p2, m2] = size (sys2); + + if (m2 != p1) + error ("lti: mrdivide: system dimensions incompatible: (%dx%d) / (%dx%d)", + p1, m1, p2, m2); + end + + sys = sys1 * sys2; + +endfunction + Added: trunk/octave-forge/extra/control-oo/inst/@lti/mtimes.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/mtimes.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/mtimes.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,80 @@ +## 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 -*- +## Matrix multiplication of LTI objects. If necessary, object conversion +## is done by sysgroup. Used by the parser for "lti1 * lti2". + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function sys = mtimes (sys2, sys1) + + [p1, m1] = size (sys1); + [p2, m2] = size (sys2); + + if (m2 != p1) + error ("lti: mtimes: system dimensions incompatible: (%dx%d) * (%dx%d)", + p2, m2, p1, m1); + end + + M22 = zeros (m2, p2); + M21 = eye (m2, p1); + M12 = zeros (m1, p2); + M11 = zeros (m1, p1); + + M = [M22, M21; + M12, M11]; + + out_idx = 1 : p2; + in_idx = m2 + (1 : m1); + + sys = __sysgroup__ (sys2, sys1); + sys = __sysconnect__ (sys, M); + sys = __sysprune__ (sys, out_idx, in_idx); + +endfunction + + +%!shared sysmat, sysmat_exp +%! +%! sys1 = ss ([0, 1; -3, -2], [0; 1], [-5, 1], [2]); +%! sys2 = ss ([-10], [1], [-40], [5]); +%! +%! sysa = sys2 * sys1; +%! +%! [A, B, C, D] = ssdata (sysa); +%! sysmat = [A, B; C, D]; +%! +%! ## expected values computed by the "dark side" +%! A_exp = [ -10 -5 1 +%! 0 0 1 +%! 0 -3 -2 ]; +%! +%! B_exp = [ 2 +%! 0 +%! 1 ]; +%! +%! C_exp = [ -40 -25 5 ]; +%! +%! D_exp = [ 10 ]; +%! +%! sysmat_exp = [A_exp, B_exp; C_exp, D_exp]; +%! +%!assert (sysmat, sysmat_exp) + Added: trunk/octave-forge/extra/control-oo/inst/@lti/parallel.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/parallel.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/parallel.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,54 @@ +## 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{sys} =} parallel (@var{sys1}, @var{sys2}) +## Parallel connection of two LTI systems. +## @example +## @group +## .......................... +## : +--------+ : +## : +-->| sys1 |---+ : +## u : | +--------+ | + : y +## -------+ O---------> +## : | +--------+ | + : +## : +-->| sys2 |---+ : +## : +--------+ : +## :.........sys............: +## +## sys = parallel (sys1, sys2) +## @end group +## @end example +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function sys = parallel (sys1, sys2) + + if (nargin == 2) + sys = sys1 + sys2; + ## elseif (nargin == 6) + + ## TODO: implement "complicated" case sys = parallel (sys1, sys2, in1, in2, out1, out2) + + else + print_usage (); + endif + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/plus.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/plus.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/plus.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,44 @@ +## 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 -*- +## Binary addition of LTI objects. If necessary, object conversion +## is done by sysgroup. Used by the parser for "lti1 + lti2". +## Operation is also known as "parallel connection". + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function sys = plus (sys1, sys2) + + [p1, m1] = size (sys1); + [p2, m2] = size (sys2); + + if (p1 != p2 || m1 != m2) + error ("lti: plus: system dimensions incompatible: (%dx%d) + (%dx%d)", + p1, m1, p2, m2); + end + + sys = __sysgroup__ (sys1, sys2); + + in_scl = [eye(m1); eye(m2)]; + out_scl = [eye(p1), eye(p2)]; + + sys = out_scl * sys * in_scl; + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/pole.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/pole.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/pole.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,31 @@ +## 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{p} =} pole (@var{sys}) +## Compute poles of LTI system. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function pol = pole (sys) + + pol = __pole__ (sys); + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/series.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/series.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/series.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,106 @@ +## 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{sys} =} series (@var{sys1}, @var{sys2}) +## @deftypefnx {Function File} {@var{sys} =} series (@var{sys1}, @var{sys2}, @var{outputs1}, @var{inputs2}) +## Series connection of two LTI models. +## @example +## @group +## ..................................... +## u : +--------+ y1 u2 +--------+ : y +## ------>| sys1 |---------->| sys2 |-------> +## : +--------+ +--------+ : +## :................sys................. +## +## sys = series (sys1, sys2) +## +## +## ..................................... +## : v2 +--------+ : +## : ---------->| | : y +## : +--------+ y1 u2 | sys2 |-------> +## u : | |---------->| | : +## ------>| sys1 | z1 +--------+ : +## : | |----------> : +## : +--------+ : +## :................sys................. +## +## outputs1 = [1] +## inputs2 = [2] +## sys = series (sys1, sys2, outputs1, inputs2) +## @end group +## @end example +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function sys = series (sys1, sys2, out1, in2) + + if (nargin == 2) + sys = sys2 * sys1; + elseif (nargin == 4) + [p1, m1] = size (sys1); + [p2, m2] = size (sys2); + + if (! isvector (out1)) + error ("series: argument 3 (outputs1) invalid"); + endif + + if (! isvector (in2)) + error ("series: argument 4 (inputs2) invalid"); + endif + + l_out1 = length (out1); + l_in2 = length (in2); + + if (l_out1 > p1) + error ("series: outputs1 has too many indices for sys1"); + endif + + if (l_in2 > m2) + error ("series: inputs2 has too many indices for sys2"); + endif + + if (l_out1 != l_in2) + error ("series: number of outputs1 and inputs2 indices must be equal"); + endif + + if (any (out1 > m1 | out1 < 1)) + error ("series: range of outputs1 indices exceeds dimensions of sys1"); + endif + + if (any (in2 > p1 | in2 < 1)) + error ("series: range of inputs2 indices exceeds dimensions of sys2"); + endif + + out_scl = zeros (l_out1, p1); + in_scl = zeros (m2, l_in2); + + for k = 1 : l_out1 + out_scl(k, out1(k)) = 1; + in_scl(in2(k), k) = 1; + endfor + + sys = sys2 * in_scl * out_scl * sys1; + else + print_usage (); + endif + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/@lti/set.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/set.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@lti/set.m 2009-10-18 09:54:21 UTC (rev 6325) @@ -0,0 +1,90 @@ +## 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{sys} =} set (@var{sys}, @var{"property"}, @var{value}) +## Set or modify properties of LTI objects. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Cre... [truncated message content] |
From: <par...@us...> - 2009-10-18 10:23:43
|
Revision: 6327 http://octave.svn.sourceforge.net/octave/?rev=6327&view=rev Author: paramaniac Date: 2009-10-18 10:23:34 +0000 (Sun, 18 Oct 2009) Log Message: ----------- control-oo: add more files Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/@ss/__zero__.m Added Paths: ----------- trunk/octave-forge/extra/control-oo/inst/ocst/ trunk/octave-forge/extra/control-oo/inst/ocst/__swap__.m trunk/octave-forge/extra/control-oo/inst/ocst/__tzero__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgfmul__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgfslv__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zginit__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgpbal__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgreduce__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgrownorm__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgshsr__.m Modified: trunk/octave-forge/extra/control-oo/inst/@ss/__zero__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@ss/__zero__.m 2009-10-18 10:10:19 UTC (rev 6326) +++ trunk/octave-forge/extra/control-oo/inst/@ss/__zero__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -24,6 +24,6 @@ function [zer, gain] = __zero__ (sys) - [zer, gain] = tzero (sys.a, sys.b, sys.c, sys.d); + [zer, gain] = __tzero__ (sys.a, sys.b, sys.c, sys.d); endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-oo/inst/ocst/__swap__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__swap__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__swap__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -0,0 +1,41 @@ +## Copyright (C) 1996, 2000, 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} {} swap (@var{inputs}) +## @format +## [a1,b1] = swap(a,b) +## interchange a and b +## @end format +## @end deftypefn + +## Author: A. S. Hodel <a.s...@en...> +## Created: July 24, 1992 +## Conversion to Octave R. Bruce Tenison July 4, 1994 + +function [a1, b1] = __swap__ (a, b) + + if (nargin != 2) + print_usage (); + endif + + a1 = b; + b1 = a; + +endfunction + Added: trunk/octave-forge/extra/control-oo/inst/ocst/__tzero__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__tzero__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__tzero__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -0,0 +1,157 @@ +## Copyright (C) 1996, 2000, 2002, 2004, 2005, 2006, 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} {[@var{zer}, @var{gain}] =} tzero (@var{a}, @var{b}, @var{c}, @var{d}, @var{opt}) +## @deftypefnx {Function File} {[@var{zer}, @var{gain}] =} tzero (@var{sys}, @var{opt}) +## Compute transmission zeros of a continuous system: +## @iftex +## @tex +## $$ \dot x = Ax + Bu $$ +## $$ y = Cx + Du $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## . +## x = Ax + Bu +## y = Cx + Du +## @end example +## @end ifinfo +## or of a discrete one: +## @iftex +## @tex +## $$ x_{k+1} = Ax_k + Bu_k $$ +## $$ y_k = Cx_k + Du_k $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## x(k+1) = A x(k) + B u(k) +## y(k) = C x(k) + D u(k) +## @end example +## @end ifinfo +## +## @strong{Outputs} +## @table @var +## @item zer +## transmission zeros of the system +## @item gain +## leading coefficient (pole-zero form) of @acronym{SISO} transfer function +## returns gain=0 if system is multivariable +## @end table +## @strong{References} +## @enumerate +## @item Emami-Naeini and Van Dooren, Automatica, 1982. +## @item Hodel, @cite{Computation of Zeros with Balancing}, 1992 Lin. Alg. Appl. +## @end enumerate +## @end deftypefn + +## Author: R. Bruce Tenison <bte...@en...> +## Created: July 4, 1994 +## A. S. Hodel Aug 1995: allow for MIMO and system data structures +## Lukas Reichlin Oct 2009: adapted for LTI Syncope + +function [zer, gain] = __tzero__ (A, B, C, D) + + if (nargin != 4) + print_usage (); + endif + + Ao = A; # save for leading coefficient + Bo = B; + Co = C; + Do = D; + + [m, n, p] = __ssmatdim__ (A, B, C, D); + + if (m*p == 1) + siso = 1; + else + siso = 0; + endif + + ## see if it's a gain block + if (isempty (A)) + zer = []; + gain = D; + return; + endif + + ## First, balance the system via the zero computation generalized eigenvalue + ## problem balancing method (Hodel and Tiller, Linear Alg. Appl., 1992) + + ## balance coefficients + [A, B, C, D] = __zgpbal__ (A, B, C, D); + + if (isa ([A, B; C, D], "single")) + meps = 2*eps("single")*norm ([A, B; C, D], "fro"); + else + meps = 2*eps*norm ([A, B; C, D], "fro"); + endif + + ## ENVD algorithm + [A, B, C, D] = __zgreduce__ (A, B, C, D, meps); + + if (! isempty (A)) + ## repeat with dual system + [A, B, C, D] = __zgreduce__ (A', B', C', D', meps); + + ## transform back + A = A'; + B = B'; + C = C'; + D = D'; + endif + + zer = []; # assume none + + if (! isempty (C)) + [W, r, Pi] = qr ([C, D]'); + [nonz, ztmp] = __zgrownorm__ (r, meps); + + if (nonz) + ## We can now solve the generalized eigenvalue problem. + [pp, mm] = size (D); + nn = rows (A); + Afm = [A, B ; C, D] * W'; + Bfm = [eye(nn), zeros(nn,mm); zeros(pp,nn+mm)]*W'; + + jdx = (mm+1):(mm+nn); + Af = Afm(1:nn,jdx); + Bf = Bfm(1:nn,jdx); + zer = qz (Af, Bf); + endif + endif + + mz = length (zer); + + ## compute leading coefficient + if (nargout == 2 && siso) + n = rows (Ao); + if (mz == n) + gain = Do; + elseif (mz < n) + gain = Co*(Ao^(n-1-mz))*Bo; + endif + else + gain = []; + endif + +endfunction + Added: trunk/octave-forge/extra/control-oo/inst/ocst/__zgfmul__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgfmul__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgfmul__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -0,0 +1,103 @@ +## Copyright (C) 1996, 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} {@var{y} =} zgfmul (@var{a}, @var{b}, @var{c}, @var{d}, @var{x}) +## Compute product of @var{zgep} incidence matrix @math{F} with vector @var{x}. +## Used by @command{zgepbal} (in @command{zgscal}) as part of generalized conjugate gradient +## iteration. +## @end deftypefn + +## References: +## ZGEP: Hodel, "Computation of Zeros with Balancing," 1992, submitted to LAA +## Generalized CG: Golub and Van Loan, "Matrix Computations, 2nd ed" 1989 + +## Author: A. S. Hodel <a.s...@en...> +## Conversion to Octave July 3, 1994 + +function y = __zgfmul__ (a, b, c, d, x) + + if (nargin != 5) + print_usage (); + endif + + [n, m] = size (b); + [p, m1] = size (c); + nm = n+m; + y = zeros (nm+p, 1); + + ## construct F column by column + for jj = 1:n + Fj = zeros (nm+p, 1); + + ## rows 1:n: F1 + aridx = complement (jj, find (a(jj,:) != 0)); + acidx = complement (jj, find (a(:,jj) != 0)); + bidx = find (b(jj,:) != 0); + cidx = find (c(:,jj) != 0); + + Fj(aridx) = Fj(aridx) - 1; # off diagonal entries of F1 + Fj(acidx) = Fj(acidx) - 1; + ## diagonal entry of F1 + Fj(jj) = length (aridx) + length (acidx) + length (bidx) + length (cidx); + + ## B' incidence + if (! isempty (bidx)) + Fj(n+bidx) = 1; + endif + + ## -C incidence + if (! isempty (cidx)) + Fj(n+m+cidx) = -1; + endif + y = y + x(jj)*Fj; # multiply by corresponding entry of x + endfor + + for jj = 1:m + Fj = zeros (nm+p, 1); + bidx = find (b(:,jj) != 0); + ## B incidence + if (! isempty (bidx)) + Fj(bidx) = 1; + endif + didx = find (d(:,jj) != 0); + ## D incidence + if (! isempty (didx)) + Fj(n+m+didx) = 1; + endif + Fj(n+jj) = length(bidx) + length(didx); # F2 is diagonal + y = y + x(n+jj)*Fj; # multiply by corresponding entry of x + endfor + + for jj = 1:p + Fj = zeros (nm+p, 1); + cidx = find (c(jj,:) != 0); + ## -C' incidence + if (! isempty (cidx)) + Fj(cidx) = -1; + endif + didx = find(d(jj,:) != 0); + ## D' incidence + if (! isempty (didx)) + Fj(n+didx) = 1; + endif + Fj(n+m+jj) = length (cidx) + length (didx); # F2 is diagonal + y = y + x(n+m+jj)*Fj; # multiply by corresponding entry of x + endfor + +endfunction Added: trunk/octave-forge/extra/control-oo/inst/ocst/__zgfslv__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgfslv__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgfslv__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -0,0 +1,81 @@ +## Copyright (C) 1996, 1998, 2000, 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} {} zgfslv (@var{n}, @var{m}, @var{p}, @var{b}) +## Solve system of equations for dense zgep problem. +## @end deftypefn + +## Author: A. S. Hodel <a.s...@en...> +## Converted to Octave by R Bruce Tenison, July 3, 1994 + +function x = __zgfslv__ (n, m, p, b) + + if (nargin != 4) + print_usage (); + endif + + nmp = n+m+p; + gam1 = (2*n)+m+p; + gam2 = n+p; + gam3 = n+m; + + G1 = givens (sqrt (m), -sqrt (p))'; + G2 = givens (m+p, sqrt (n*(m+p)))'; + + x = b; + + ## 1) U1 e^n = sqrt(n)e_1^n + ## 2) U2 e^m = sqrt(m)e_1^m + ## 3) U3 e^p = sqrt(p)e_1^p + xdx1 = 1:n; + xdx2 = n+(1:m); + xdx3 = n+m+(1:p); + + x(xdx1,1) = __zgshsr__ (x(xdx1,1)); + x(xdx2,1) = __zgshsr__ (x(xdx2,1)); + x(xdx3,1) = __zgshsr__ (x(xdx3,1)); + + ## 4) Givens rotations to reduce stray non-zero elements + idx1 = [n+1, n+m+1]; + idx2 = [1, n+1]; + + x(idx1) = G1'*x(idx1); + x(idx2) = G2'*x(idx2); + + ## 6) Scale x, then back-transform to get x + en = ones (n, 1); + em = ones (m, 1); + ep = ones (p, 1); + lam = [gam1*en; gam2*em; gam3*ep]; + lam(1) = n+m+p; + lam(n+1) = 1; # dummy value to avoid divide by zero + lam(n+m+1) = n+m+p; + + x = x ./ lam; + x(n+1) = 0; # minimum norm solution + + ## back transform now. + x(idx2) = G2*x(idx2); + x(idx1) = G1*x(idx1); + x(xdx3,1) = __zgshsr__ (x(xdx3,1)); + x(xdx2,1) = __zgshsr__ (x(xdx2,1)); + x(xdx1,1) = __zgshsr__ (x(xdx1,1)); + +endfunction + Added: trunk/octave-forge/extra/control-oo/inst/ocst/__zginit__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zginit__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zginit__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -0,0 +1,97 @@ +## Copyright (C) 1996, 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} {@var{zz} =} zginit (@var{a}, @var{b}, @var{c}, @var{d}) +## Construct right hand side vector @var{zz} +## for the zero-computation generalized eigenvalue problem +## balancing procedure. Called by @command{zgepbal}. +## @end deftypefn + +## References: +## ZGEP: Hodel, "Computation of Zeros with Balancing," 1992, submitted to LAA +## Generalized CG: Golub and Van Loan, "Matrix Computations, 2nd ed" 1989 + +## Author: A. S. Hodel <a.s...@en...> +## Created: July 24, 1992 +## Conversion to Octave by R. Bruce Tenison, July 3, 1994 + +function zz = __zginit__ (a, b, c, d) + + if (nargin != 4) + print_usage (); + endif + + [nn, mm] = size (b); + [pp, mm] = size (d); + + nmp = nn+mm+pp; + + ## set up log vector zz + zz = zeros (nmp, 1); + + ## zz part 1: + for i = 1:nn + ## nonzero off diagonal entries of a + if (nn > 1) + nidx = complement (i, 1:nn); + a_row_i = a(i,nidx); + a_col_i = a(nidx,i); + arnz = a_row_i(find (a_row_i != 0)); + acnz = a_col_i(find (a_col_i != 0)); + else + arnz = acnz = []; + endif + + ## row of b + bidx = find (b(i,:) != 0); + b_row_i = b(i,bidx); + + ## column of c + cidx = find (c(:,i) != 0); + c_col_i = c(cidx,i); + + ## sum the entries + zz(i) = sum (log (abs (acnz))) - sum (log (abs (arnz))) ... + - sum (log (abs (b_row_i))) + sum (log (abs (c_col_i))); + endfor + + ## zz part 2: + bd = [b; d]; + for i = 1:mm + i1 = i+nn; + + ## column of [b;d] + bdidx = find (bd(:,i) != 0); + bd_col_i = bd(bdidx,i); + zz(i1) = sum (log (abs(bd_col_i))); + endfor + + ## zz part 3: + cd = [c, d]; + for i = 1:pp + i1 = i+nn+mm; + cdidx = find (cd(i,:) != 0); + cd_row_i = cd(i,cdidx); + zz(i1) = -sum (log (abs (cd_row_i))); + endfor + + ## now set zz as log base 2 + zz *= (1 / log (2)); + +endfunction Added: trunk/octave-forge/extra/control-oo/inst/ocst/__zgpbal__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgpbal__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgpbal__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -0,0 +1,104 @@ +## Copyright (C) 1996, 2000, 2002, 2004, 2005, 2006, 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/>. + +## Undocumented internal function. + +## -*- texinfo -*- +## @deftypefn {Function File} {} __zgpbal__ (@var{sys}) +## +## Used internally in @command{tzero}; minimal argument checking performed. +## +## Implementation of zero computation generalized eigenvalue problem +## balancing method (Hodel and Tiller, Allerton Conference, 1991) +## Based on Ward's balancing algorithm (@acronym{SIAM} J. Sci Stat. Comput., 1981). +## +## @command{__zgpbal__} computes a state/input/output weighting that attempts to +## reduced the range of the magnitudes of the nonzero elements of [@var{a}, @var{b}, +## @var{c}, @var{d}]. +## The weighting uses scalar multiplication by powers of 2, so no roundoff +## will occur. +## +## @command{__zgpbal__} should be followed by @command{zgpred}. +## @end deftypefn + +## References: +## ZGEP: Hodel, "Computation of Zeros with Balancing," 1992, submitted to LAA +## Generalized CG: Golub and Van Loan, "Matrix Computations, 2nd ed" 1989 + +## Author: A. S. Hodel <a.s...@en...> +## Created: July 24, 1992 +## Conversion to Octave by R. Bruce Tenison July 3, 1994 +## Adapted for LTI Syncope by Lukas Reichlin October 2009 + +function [a, b, c, d] = __zgpbal__ (a, b, c, d) + + if (nargin != 4) + print_usage (); + endif + + [mm, nn, pp] = __ssmatdim__ (a, b, c, d); + + np1 = nn+1; + nmp = nn+mm+pp; + + ## set up log vector zz, incidence matrix ff + zz = zginit (a, b, c, d); + + if (norm (zz)) + ## generalized conjugate gradient approach + xx = zgscal (a, b, c, d, zz, nn, mm, pp); + + for i = 1:nmp + xx(i) = floor (xx(i)+0.5); + xx(i) = 2.0^xx(i); + endfor + + ## now scale a + ## block 1: a = sigma a inv(sigma) + for i = 1:nn + a(i,1:nn) = a(i,1:nn)*xx(i); + a(1:nn,i) = a(1:nn,i)/xx(i); + endfor + ## block 2: b= sigma a phi + for j = 1:mm + j1 = j+nn; + b(1:nn,j) = b(1:nn,j)*xx(j1); + endfor + for i = 1:nn + b(i,1:mm) = b(i,1:mm)*xx(i); + endfor + for i = 1:pp + i1 = i+nn+mm; + ## block 3: c = psi C inv(sigma) + c(i,1:nn) = c(i,1:nn)*xx(i1); + endfor + for j = 1:nn + c(1:pp,j) = c(1:pp,j)/xx(j); + endfor + ## block 4: d = psi D phi + for j = 1:mm + j1 = j+nn; + d(1:pp,j) = d(1:pp,j)*xx(j1); + endfor + for i = 1:pp + i1 = i + nn + mm; + d(i,1:mm) = d(i,1:mm)*xx(i1); + endfor + endif + +endfunction Added: trunk/octave-forge/extra/control-oo/inst/ocst/__zgreduce__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgreduce__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgreduce__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -0,0 +1,104 @@ +## Copyright (C) 1996, 2000, 2005, 2007 +## Auburn University. All rights reserved. +## Copyright (C) 2009 Lukas F. 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 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} {} zgreduce (@var{sys}, @var{meps}) +## Implementation of procedure REDUCE in (Emami-Naeini and Van Dooren, +## Automatica, # 1982). +## @end deftypefn + +function [A, B, C, D] = __zgreduce__ (A, B, C, D, meps) + + if (nargin != 5) + print_usage (); + endif + + exit_1 = 0; # exit_1 = 1 or 2 on exit of loop + + [m, n, p] = __ssmatdim__ (A, B, C, D); + + if (n == 0) + exit_1 = 2; # there are no finite zeros + endif + + while (! exit_1) + [Q, R, Pi] = qr (D); # compress rows of D + D = Q'*D; + C = Q'*C; + + ## check row norms of D + [sig, tau] = __zgrownorm__ (D, meps); + + if (tau == 0) + exit_1 = 1; # exit_1 - reduction complete and correct + else + Cb = Db = []; + if (sig) + Cb = C(1:sig,:); + Db = D(1:sig,:); + endif + Ct = C(sig+(1:tau),:); + + ## compress columns of Ct + [pp, nn] = size (Ct); + rvec = nn:-1:1; + [V, Sj, Pi] = qr (Ct'); + V = V(:,rvec); + [rho, gnu] = __zgrownorm__ (Sj, meps); + + if (rho == 0) + exit_1 = 1; # exit_1 - reduction complete and correct + elseif (gnu == 0) + exit_1 = 2; # there are no zeros at all + else + mu = rho + sig; + + ## update system with Q + M = [A, B]; + [nn, mm] = size (B); + + pp = rows (D); + Vm = [V, zeros(nn,mm); zeros(mm,nn), eye(mm)]; + if (sig) + M = [M; Cb, Db]; + Vs = [V', zeros(nn,sig); zeros(sig,nn), eye(sig)]; + else + Vs = V'; + endif + + M = Vs*M*Vm; + + idx = 1:gnu; + jdx = nn + (1:mm); + sdx = gnu + (1:mu); + + A = M(idx,idx); + B = M(idx,jdx); + C = M(sdx,idx); + D = M(sdx,jdx); + + endif + endif + endwhile + + if (exit_1 == 2) + ## there are no zeros at all! + A = B = C = []; + endif + +endfunction Added: trunk/octave-forge/extra/control-oo/inst/ocst/__zgrownorm__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgrownorm__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgrownorm__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -0,0 +1,39 @@ +## Copyright (C) 1996, 2000, 2003, 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} {[@var{nonz}, @var{zer}] =} zgrownorm (@var{mat}, @var{meps}) +## Return @var{nonz} = number of rows of @var{mat} whose two norm +## exceeds @var{meps}, and @var{zer} = number of rows of mat whose two +## norm is less than @var{meps}. +## @end deftypefn + +function [sig, tau] = __zgrownorm__ (mat, meps) + + if (nargin != 2) + print_usage (); + endif + + rownorm = []; + for ii = 1:rows (mat) + rownorm(ii) = norm (mat(ii,:)); + endfor + sig = sum (rownorm > meps); + tau = sum (rownorm <= meps); + +endfunction Added: trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -0,0 +1,152 @@ +## Copyright (C) 1996, 1998, 2000, 2002, 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} {@var{x} =} zgscal (@var{f}, @var{z}, @var{n}, @var{m}, @var{p}) +## Generalized conjugate gradient iteration to +## solve zero-computation generalized eigenvalue problem balancing equation +## @math{fx=z}; called by @command{zgepbal}. +## @end deftypefn + +## References: +## ZGEP: Hodel, "Computation of Zeros with Balancing," 1992, submitted to LAA +## Generalized CG: Golub and Van Loan, "Matrix Computations, 2nd ed" 1989 + +## Author: A. S. Hodel <a.s...@en...> +## Created: July 24, 1992 +## Conversion to Octave R. Bruce Tenison July 3, 1994 + +function x = __zgscal__ (a, b, c, d, z, n, m, p) + + if (nargin != 8) + print_usage (); + endif + + ## initialize parameters: + ## Givens rotations, diagonalized 2x2 block of F, gcg vector initialization + + nmp = n+m+p; + + ## x_0 = x_{-1} = 0, r_0 = z + x = zeros (nmp, 1); + xk1 = x; + xk2 = x; + rk1 = z; + k = 0; + + ## construct balancing least squares problem + F = eye (nmp); + for kk = 1:nmp + F(1:nmp,kk) = __zgfmul__ (a, b, c, d, F(:,kk)); + endfor + + [U, H, k1] = krylov (F, z, nmp, 1e-12, 1); + if (! issquare (H)) + if (columns (H) != k1) + error ("zgscal(tzero): k1=%d, columns(H)=%d", k1, columns (H)); + elseif (rows (H) != k1+1) + error ("zgscal: k1=%d, rows(H) = %d", k1, rows (H)); + elseif (norm (H(k1+1,:)) > 1e-12*norm (H, "inf")) + zgscal_last_row_of_H = H(k1+1,:) + error ("zgscal: last row of H nonzero (norm(H)=%e)", norm (H, "inf")) + endif + H = H(1:k1,1:k1); + U = U(:,1:k1); + endif + + ## tridiagonal H can still be rank deficient, so do permuted qr + ## factorization + [qq, rr, pp] = qr (H); # H = qq*rr*pp' + nn = rank (rr); + qq = qq(:,1:nn); + rr = rr(1:nn,:); # rr may not be square, but "\" does least + xx = U*pp*(rr\qq'*(U'*z)); # squares solution, so this works + ## xx1 = pinv(F)*z; + ## zgscal_x_xx1_err = [xx,xx1,xx-xx1] + return; + + ## the rest of this is left from the original zgscal; + ## I've had some numerical problems with the GCG algorithm, + ## so for now I'm solving it with the krylov routine. + + ## initialize residual error norm + rnorm = norm (rk1, 1); + + xnorm = 0; + fnorm = 1e-12 * norm ([a, b; c, d], 1); + + gamk2 = 0; + omega1 = 0; + ztmz2 = 0; + + ## do until small changes to x + len_x = length(x); + while ((k < 2*len_x && xnorm > 0.5 && rnorm > fnorm) || k == 0) + k++; + + ## solve F_d z_{k-1} = r_{k-1} + zk1= __zgfslv__ (n, m, p, rk1); + + ## Generalized CG iteration + ## gamk1 = (zk1'*F_d*zk1)/(zk1'*F*zk1); + ztMz1 = zk1'*rk1; + gamk1 = ztMz1/(zk1'*__zgfmul__ (a, b, c, d, zk1)); + + if (rem (k, len_x) == 1) + omega = 1; + else + omega = 1/(1-gamk1*ztMz1/(gamk2*omega1*ztmz2)); + endif + + ## store x in xk2 to save space + xk2 = xk2 + omega*(gamk1*zk1 + xk1 - xk2); + + ## compute new residual error: rk = z - F xk, check end conditions + rk1 = z - __zgfmul__ (a, b, c, d, xk2); + rnorm = norm (rk1); + xnorm = max (abs (xk1 - xk2)); + + ## printf("zgscal: k=%d, gamk1=%e, gamk2=%e, \nztMz1=%e ztmz2=%e\n", ... + ## k,gamk1, gamk2, ztMz1, ztmz2); + ## xk2_1_zk1 = [xk2 xk1 zk1] + ## ABCD = [a,b;c,d] + ## prompt + + ## get ready for next iteration + gamk2 = gamk1; + omega1 = omega; + ztmz2 = ztMz1; + [xk1, xk2] = __swap__ (xk1, xk2); + endwhile + x = xk2; + + ## check convergence + if (xnorm> 0.5 && rnorm > fnorm) + warning ("zgscal(tzero): GCG iteration failed; solving with pinv"); + + ## perform brute force least squares; construct F + Am = eye (nmp); + for ii = 1:nmp + Am(:,ii) = __zgfmul__ (a, b, c, d, Am(:,ii)); + endfor + + ## now solve with qr factorization + x = pinv (Am) * z; + endif + +endfunction Added: trunk/octave-forge/extra/control-oo/inst/ocst/__zgshsr__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgshsr__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgshsr__.m 2009-10-18 10:23:34 UTC (rev 6327) @@ -0,0 +1,58 @@ +## Copyright (C) 1996, 2000, 2002, 2003, 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} {@var{x} =} zgshsr (@var{y}) +## Apply householder vector based on +## @iftex +## @tex +## $ e^m $ +## @end tex +## @end iftex +## @ifinfo +## @math{e^(m)} +## @end ifinfo +## to column vector @var{y}. +## Called by @command{zgfslv}. +## @end deftypefn + +## Author: A. S. Hodel <a.s...@en...> +## Created: July 24, 1992 +## Conversion to Octave by R. Bruce Tenison July 3, 1994 + +function x = __zgshsr__ (y) + + if (nargin != 1) + print_usage (); + endif + + if (! isvector (y)) + error ("y(%dx%d) must be a vector", rows (y), columns (y)); + endif + x = vec (y); + m = length (x); + if (m > 1) + beta = (1 + sqrt (m)) * x(1) + sum (x(2:m)); + beta /= (m + sqrt (m)); + x(1) -= (beta * (1 + sqrt (m))); + x(2:m) -= (beta * ones (m-1,1)); + else + x = -x; + endif + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-10-18 18:11:01
|
Revision: 6333 http://octave.svn.sourceforge.net/octave/?rev=6333&view=rev Author: paramaniac Date: 2009-10-18 18:10:54 +0000 (Sun, 18 Oct 2009) Log Message: ----------- control-oo: fix bug Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/@ss/__sys2tf__.m trunk/octave-forge/extra/control-oo/inst/@ss/ss.m trunk/octave-forge/extra/control-oo/inst/@tf/__sys2ss__.m trunk/octave-forge/extra/control-oo/inst/@tf/tf.m Modified: trunk/octave-forge/extra/control-oo/inst/@ss/__sys2tf__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@ss/__sys2tf__.m 2009-10-18 17:07:13 UTC (rev 6332) +++ trunk/octave-forge/extra/control-oo/inst/@ss/__sys2tf__.m 2009-10-18 18:10:54 UTC (rev 6333) @@ -22,8 +22,11 @@ ## Created: October 2009 ## Version: 0.1 -function retsys = __sys2tf__ (sys) +function [retsys, retlti] = __sys2tf__ (sys) error ("ss: ss2tf: not implemented yet"); + retsys = tf (num, den, get (sys, "tsam")); # tsam needed to set appropriate tfvar + retlti = sys.lti; # preserve lti properties + endfunction \ No newline at end of file Modified: trunk/octave-forge/extra/control-oo/inst/@ss/ss.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@ss/ss.m 2009-10-18 17:07:13 UTC (rev 6332) +++ trunk/octave-forge/extra/control-oo/inst/@ss/ss.m 2009-10-18 18:10:54 UTC (rev 6333) @@ -48,7 +48,8 @@ sys = a; return; elseif (isa (a, "lti")) # another lti object - sys = __sys2ss__ (a); + [sys, alti] = __sys2ss__ (a); + sys.lti = alti; # preserve lti properties return; elseif (isnumeric (a)) # static gain d = a; @@ -75,12 +76,12 @@ [b, c] = __gaincheck__ (b, c, d); argc = numel (varargin); - if (issample (varargin{1})) # sys = ss (a, b, c, d, tsam, "prop1, "val1", ...) + if (isscalar (varargin{1}) && (varargin{1} == abs (varargin{1}))) # sys = ss (a, b, c, d, tsam, "prop1, "val1", ...) tsam = varargin{1}; + argc -= 1; - if (argc > 1) + if (argc > 0) varargin = varargin(2:end); - argc -= 1; endif else # sys = ss (a, b, c, d, "prop1, "val1", ...) tsam = 0; Modified: trunk/octave-forge/extra/control-oo/inst/@tf/__sys2ss__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@tf/__sys2ss__.m 2009-10-18 17:07:13 UTC (rev 6332) +++ trunk/octave-forge/extra/control-oo/inst/@tf/__sys2ss__.m 2009-10-18 18:10:54 UTC (rev 6333) @@ -22,8 +22,11 @@ ## Created: October 2009 ## Version: 0.1 -function retsys = __sys2ss__ (sys) +function [retsys, retlti] = __sys2ss__ (sys) error ("tf: tf2ss: not implemented yet"); + retsys = ss (a, b, c, d); + retlti = sys.lti; # preserve lti properties + endfunction \ No newline at end of file Modified: trunk/octave-forge/extra/control-oo/inst/@tf/tf.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@tf/tf.m 2009-10-18 17:07:13 UTC (rev 6332) +++ trunk/octave-forge/extra/control-oo/inst/@tf/tf.m 2009-10-18 18:10:54 UTC (rev 6333) @@ -47,7 +47,8 @@ sys = num; return; elseif (isa (num, "lti")) # another lti object - sys = __sys2tf__ (num); + [sys, numlti] = __sys2tf__ (num); + sys.lti = numlti; # preserve lti properties return; elseif (isnumeric (num)) # static gain num = num2cell (num); @@ -83,13 +84,18 @@ den = __conv2tfpolycell__ (den); argc = numel (varargin); - if (issample (varargin{1})) # sys = tf (num, den, tsam, "prop1, "val1", ...) + if (isscalar (varargin{1}) && (varargin{1} == abs (varargin{1}))) # sys = tf (num, den, tsam, "prop1, "val1", ...) tsam = varargin{1}; - tfvar = "z"; - - if (argc > 1) + argc -= 1; + + if (varargin{1} != 0) + tfvar = "z"; + else + tfvar = "s"; + endif + + if (argc > 0) varargin = varargin(2:end); - argc -= 1; endif else # sys = tf (num, den, "prop1, "val1", ...) tsam = 0; @@ -98,6 +104,7 @@ endswitch + [p, m] = __tfnddim__ (num, den); tfdata = struct ("num", {num}, This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-10-21 06:54:27
|
Revision: 6351 http://octave.svn.sourceforge.net/octave/?rev=6351&view=rev Author: paramaniac Date: 2009-10-21 06:54:21 +0000 (Wed, 21 Oct 2009) Log Message: ----------- control-oo: update Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/control/place.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m Modified: trunk/octave-forge/extra/control-oo/inst/control/place.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/place.m 2009-10-21 06:53:09 UTC (rev 6350) +++ trunk/octave-forge/extra/control-oo/inst/control/place.m 2009-10-21 06:54:21 UTC (rev 6351) @@ -111,11 +111,11 @@ ## Check if the eigenvalues of (A-BK) are the same specified in P Pcalc = eig (A-B*K); - Pcalc = __sortcom__ (Pcalc); - P = __sortcom__ (P); + Pcalc = sort (Pcalc); + P = sort (P); if (max ((abs(Pcalc)-abs(P))./abs(P) ) > 0.1) - warning ("place: Pole placed at more than 10% relative error from specified"); + warning ("place: pole placed at more than 10% relative error from specified"); endif endfunction Modified: trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m 2009-10-21 06:53:09 UTC (rev 6350) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m 2009-10-21 06:54:21 UTC (rev 6351) @@ -150,3 +150,12 @@ endif endfunction + + +function [a1, b1] = __swap__ (a, b) + + a1 = b; + b1 = a; + +endfunction + This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-10-22 08:56:35
|
Revision: 6355 http://octave.svn.sourceforge.net/octave/?rev=6355&view=rev Author: paramaniac Date: 2009-10-22 08:56:22 +0000 (Thu, 22 Oct 2009) Log Message: ----------- control-oo: mark legacy files as adapted to reduce confusion Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/control/are.m trunk/octave-forge/extra/control-oo/inst/control/dare.m trunk/octave-forge/extra/control-oo/inst/control/dlyap.m trunk/octave-forge/extra/control-oo/inst/control/gram.m trunk/octave-forge/extra/control-oo/inst/control/isctrb.m trunk/octave-forge/extra/control-oo/inst/control/isobsv.m trunk/octave-forge/extra/control-oo/inst/control/issample.m trunk/octave-forge/extra/control-oo/inst/control/lqr.m trunk/octave-forge/extra/control-oo/inst/control/lyap.m trunk/octave-forge/extra/control-oo/inst/control/place.m trunk/octave-forge/extra/control-oo/inst/ocst/__tzero__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgfmul__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgfslv__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zginit__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgpbal__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgreduce__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgrownorm__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m trunk/octave-forge/extra/control-oo/inst/ocst/__zgshsr__.m Modified: trunk/octave-forge/extra/control-oo/inst/control/are.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/are.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/control/are.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -64,8 +64,11 @@ ## Author: A. S. Hodel <a.s...@en...> ## Created: August 1993 -## Adapted for LTI Syncope by Lukas Reichlin, October 2009 +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function x = are (a, b, c, opt) if (nargin == 3 || nargin == 4) Modified: trunk/octave-forge/extra/control-oo/inst/control/dare.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/dare.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/control/dare.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -75,6 +75,10 @@ ## Created: August 1993 ## Adapted-By: jwe +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function x = dare (a, b, q, r, opt) if (nargin == 4 || nargin == 5) Modified: trunk/octave-forge/extra/control-oo/inst/control/dlyap.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/dlyap.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/control/dlyap.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -87,6 +87,10 @@ ## Author: A. S. Hodel <a.s...@en...> ## Created: August 1993 +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function x = dlyap (a, b) if (nargin != 2) Modified: trunk/octave-forge/extra/control-oo/inst/control/gram.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/gram.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/control/gram.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -31,6 +31,10 @@ ## Author: A. S. Hodel <a.s...@en...> +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function W = gram (argin1, argin2) if (nargin != 2) Modified: trunk/octave-forge/extra/control-oo/inst/control/isctrb.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/isctrb.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/control/isctrb.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -59,8 +59,11 @@ ## Created: August 1993 ## Updated by A. S. Hodel (sc...@en...) Aubust, 1995 to use krylovb ## Updated by John Ingram (in...@en...) July, 1996 for packed systems -## Updated by Lukas Reichlin (luk...@gm...) October, 2009 for LTI objects +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function [retval, U] = isctrb (a, b, tol) deftol = 1; # assume default tolerance Modified: trunk/octave-forge/extra/control-oo/inst/control/isobsv.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/isobsv.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/control/isobsv.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -34,8 +34,11 @@ ## Author: A. S. Hodel <a.s...@en...> ## Created: August 1993 ## Updated by John Ingram (in...@en...) July 1996. -## Updated by Lukas Reichlin (luk...@gm...) October 2009. +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function [retval, U] = isobsv (a, c, tol) if (nargin < 1) Modified: trunk/octave-forge/extra/control-oo/inst/control/issample.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/issample.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/control/issample.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -23,6 +23,10 @@ ## Author: A. S. Hodel <a.s...@en...> ## Created: July 1995 +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: September 2009 +## Version: 0.1 + function bool = issample (tsam) if (nargin != 1) Modified: trunk/octave-forge/extra/control-oo/inst/control/lqr.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/lqr.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/control/lqr.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -113,8 +113,11 @@ ## Author: A. S. Hodel <a.s...@en...> ## Created: August 1993. -## Adapted for LTI Syncope by Lukas Reichlin, October 2009 +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function [k, p, e] = lqr (a, b, q, r, s) if (nargin != 4 && nargin != 5) Modified: trunk/octave-forge/extra/control-oo/inst/control/lyap.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/lyap.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/control/lyap.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -79,6 +79,10 @@ ## Created: August 1993 ## Adapted-By: jwe +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function x = lyap (a, b, c) if (nargin != 3 && nargin != 2) Modified: trunk/octave-forge/extra/control-oo/inst/control/place.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/place.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/control/place.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -42,8 +42,11 @@ ## ## code adaped by A.S.Hodel (a.s...@en...) for use in controls ## toolbox -## Adapted for LTI Syncope by Lukas Reichlin (luk...@gm...) +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function K = place (argin1, argin2, argin3) if (nargin == 3) Modified: trunk/octave-forge/extra/control-oo/inst/ocst/__tzero__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__tzero__.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__tzero__.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -65,8 +65,11 @@ ## Author: R. Bruce Tenison <bte...@en...> ## Created: July 4, 1994 ## A. S. Hodel Aug 1995: allow for MIMO and system data structures -## Lukas Reichlin Oct 2009: adapted for LTI Syncope +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function [zer, gain] = __tzero__ (A, B, C, D) if (nargin != 4) Modified: trunk/octave-forge/extra/control-oo/inst/ocst/__zgfmul__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgfmul__.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgfmul__.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -30,6 +30,10 @@ ## Author: A. S. Hodel <a.s...@en...> ## Conversion to Octave July 3, 1994 +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function y = __zgfmul__ (a, b, c, d, x) if (nargin != 5) Modified: trunk/octave-forge/extra/control-oo/inst/ocst/__zgfslv__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgfslv__.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgfslv__.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -24,6 +24,10 @@ ## Author: A. S. Hodel <a.s...@en...> ## Converted to Octave by R Bruce Tenison, July 3, 1994 +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function x = __zgfslv__ (n, m, p, b) if (nargin != 4) Modified: trunk/octave-forge/extra/control-oo/inst/ocst/__zginit__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zginit__.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zginit__.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -31,6 +31,10 @@ ## Created: July 24, 1992 ## Conversion to Octave by R. Bruce Tenison, July 3, 1994 +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function zz = __zginit__ (a, b, c, d) if (nargin != 4) Modified: trunk/octave-forge/extra/control-oo/inst/ocst/__zgpbal__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgpbal__.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgpbal__.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -43,8 +43,11 @@ ## Author: A. S. Hodel <a.s...@en...> ## Created: July 24, 1992 ## Conversion to Octave by R. Bruce Tenison July 3, 1994 -## Adapted for LTI Syncope by Lukas Reichlin October 2009 +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function [a, b, c, d] = __zgpbal__ (a, b, c, d) if (nargin != 4) Modified: trunk/octave-forge/extra/control-oo/inst/ocst/__zgreduce__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgreduce__.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgreduce__.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -20,8 +20,11 @@ ## @deftypefn {Function File} {} zgreduce (@var{sys}, @var{meps}) ## Implementation of procedure REDUCE in (Emami-Naeini and Van Dooren, ## Automatica, # 1982). -## @end deftypefn +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function [A, B, C, D] = __zgreduce__ (A, B, C, D, meps) if (nargin != 5) Modified: trunk/octave-forge/extra/control-oo/inst/ocst/__zgrownorm__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgrownorm__.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgrownorm__.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -23,6 +23,10 @@ ## norm is less than @var{meps}. ## @end deftypefn +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function [sig, tau] = __zgrownorm__ (mat, meps) if (nargin != 2) Modified: trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgscal__.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -31,6 +31,10 @@ ## Created: July 24, 1992 ## Conversion to Octave R. Bruce Tenison July 3, 1994 +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function x = __zgscal__ (a, b, c, d, z, n, m, p) if (nargin != 8) Modified: trunk/octave-forge/extra/control-oo/inst/ocst/__zgshsr__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__zgshsr__.m 2009-10-22 07:29:37 UTC (rev 6354) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__zgshsr__.m 2009-10-22 08:56:22 UTC (rev 6355) @@ -35,6 +35,10 @@ ## Created: July 24, 1992 ## Conversion to Octave by R. Bruce Tenison July 3, 1994 +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + function x = __zgshsr__ (y) if (nargin != 1) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-10-27 16:17:29
|
Revision: 6404 http://octave.svn.sourceforge.net/octave/?rev=6404&view=rev Author: paramaniac Date: 2009-10-27 16:17:16 +0000 (Tue, 27 Oct 2009) Log Message: ----------- control-oo: add some frequency response draft code which is needed by step Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/control/__timeresp__.m Added Paths: ----------- trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m trunk/octave-forge/extra/control-oo/inst/control/dcgain.m trunk/octave-forge/extra/control-oo/inst/control/freqresp.m Added: trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m 2009-10-27 16:17:16 UTC (rev 6404) @@ -0,0 +1,63 @@ +## 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 -*- + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function H = __freqresp__ (sys, w, resptype = 0) + + [p, m] = size (sys); + [A, B, C, D, Ts] = ssdata (sys); + I = eye (size (A)); + + % J = eye (m); + + if (Ts == 0 || Ts == -1) # continuous system + s = i * w; + else # discrete system + s = exp (i * w * Ts); + endif + + l_s = length (s); + H = zeros (p, m, l_s); + + switch (resptype) + case 0 # default system + for k = 1 : l_s + H(:, :, k) = C * inv (s(k)*I - A) * B + D; + endfor + + case 1 # inversed system + for k = 1 : l_s + H(:, :, k) = inv (C * inv (s(k)*I - A) * B + D); + endfor + + case 2 # inversed sensitivity + for k = 1 : l_s + H(:, :, k) = J + C * inv (s(k)*I - A) * B + D; + endfor + + case 3 # inversed complementary sensitivity + for k = 1 : l_s + H(:, :, k) = J + inv (C * inv (s(k)*I - A) * B + D); + endfor + endswitch + +endfunction \ No newline at end of file Modified: trunk/octave-forge/extra/control-oo/inst/control/__timeresp__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/__timeresp__.m 2009-10-27 09:01:18 UTC (rev 6403) +++ trunk/octave-forge/extra/control-oo/inst/control/__timeresp__.m 2009-10-27 16:17:16 UTC (rev 6404) @@ -63,6 +63,7 @@ switch (resptype) case "initial" str = "Response to Initial Conditions"; + yfinal = zeros (p, 1); ## preallocate memory y = zeros (l_t, p); @@ -84,6 +85,7 @@ case "step" str = "Step Response"; + yfinal = dcgain (sys); ## preallocate memory y = zeros (l_t, p, m); @@ -105,6 +107,7 @@ case "impulse" str = "Impulse Response"; + yfinal = zeros (p, m); ## preallocate memory y = zeros (l_t, p, m); @@ -168,7 +171,7 @@ for j = 1 : cols subplot (p, cols, (k-1)*cols+j) - stairs (t, y(:, k, j)) + stairs (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]) grid on if (k == 1) @@ -189,7 +192,7 @@ for j = 1 : cols subplot (p, cols, (k-1)*cols+j) - plot (t, y(:, k, j)) + plot (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]) grid on if (k == 1) Added: trunk/octave-forge/extra/control-oo/inst/control/dcgain.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/dcgain.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/control/dcgain.m 2009-10-27 16:17:16 UTC (rev 6404) @@ -0,0 +1,32 @@ +## 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 -*- + +## 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 Added: trunk/octave-forge/extra/control-oo/inst/control/freqresp.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/freqresp.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/control/freqresp.m 2009-10-27 16:17:16 UTC (rev 6404) @@ -0,0 +1,32 @@ +## 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 -*- + +## 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. |
From: <par...@us...> - 2009-10-30 16:31:41
|
Revision: 6412 http://octave.svn.sourceforge.net/octave/?rev=6412&view=rev Author: paramaniac Date: 2009-10-30 16:31:31 +0000 (Fri, 30 Oct 2009) Log Message: ----------- control-oo: update freqresp Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m trunk/octave-forge/extra/control-oo/inst/@tf/__freqresp__.m Modified: trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m 2009-10-30 14:52:35 UTC (rev 6411) +++ trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m 2009-10-30 16:31:31 UTC (rev 6412) @@ -26,8 +26,11 @@ [p, m] = size (sys); [A, B, C, D, Ts] = ssdata (sys); I = eye (size (A)); + J = eye (m); - % J = eye (m); + if (resptype != 0 && m != p) + error ("ss: freqresp: system must be square for response type %d", resptype); + endif if (Ts == 0 || Ts == -1) # continuous system s = i * w; @@ -58,6 +61,10 @@ for k = 1 : l_s H(:, :, k) = J + inv (C * inv (s(k)*I - A) * B + D); endfor + + otherwise + error ("ss: freqresp: invalid response type"); + endswitch endfunction \ No newline at end of file Modified: trunk/octave-forge/extra/control-oo/inst/@tf/__freqresp__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@tf/__freqresp__.m 2009-10-30 14:52:35 UTC (rev 6411) +++ trunk/octave-forge/extra/control-oo/inst/@tf/__freqresp__.m 2009-10-30 16:31:31 UTC (rev 6412) @@ -35,65 +35,44 @@ l_s = length (s); H = zeros (p, m, l_s); - switch (resptype) - case 0 # default system - for b = 1 : p - for a = 1 : m - num_pm = num{b, a}; - den_pm = den{b, a}; + for b = 1 : p + for a = 1 : m + num_pm = num{b, a}; + den_pm = den{b, a}; - for k = 1 : l_s - H(b, a, k) = polyval (num_pm, s(k)) / polyval (den_pm, s(k)); - endfor - endfor + for k = 1 : l_s + H(b, a, k) = polyval (num_pm, s(k)) / polyval (den_pm, s(k)); endfor -#{ - ## FIXME: response types 1-3 are wrong for MIMO systems! - case 1 # inversed system - for b = 1 : p - for a = 1 : m - num_pm = num{b, a}; - den_pm = den{b, a}; + endfor + endfor - for k = 1 : l_s - H(b, a, k) = polyval (den_pm, s(k)) / polyval (num_pm, s(k)); - endfor - endfor - endfor + if (resptype) + if (m != p) + error ("tf: freqresp: system must be square for response type %d", resptype); + endif - case 2 # inversed sensitivity - J = eye (p); + I = eye (p); - for b = 1 : p - for a = 1 : m - num_pm = num{b, a}; - den_pm = den{b, a}; - J_pm = J(b, a); + switch (resptype) + case 1 # inversed system + for k = 1 : l_s + H(:, :, k) = inv (H(:, :, k)); + endfor - for k = 1 : l_s - H(b, a, k) = J_pm + polyval (num_pm, s(k)) / polyval (den_pm, s(k)); - endfor + case 2 # inversed sensitivity + for k = 1 : l_s + H(:, :, k) = I + H(:, :, k); endfor - endfor - case 3 # inversed complementary sensitivity - J = eye (p); - - for b = 1 : p - for a = 1 : m - num_pm = num{b, a}; - den_pm = den{b, a}; - J_pm = J(b, a); - - for k = 1 : l_s - H(b, a, k) = J_pm + polyval (den_pm, s(k)) / polyval (num_pm, s(k)); - endfor + case 3 # inversed complementary sensitivity + for k = 1 : l_s + H(:, :, k) = I + inv(H(:, :, k)); endfor - endfor -#} - otherwise - error ("tf: freqresp: invalid response type"); - endswitch + otherwise + error ("tf: freqresp: invalid response type"); + endswitch + endif + endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-10-30 19:15:47
|
Revision: 6414 http://octave.svn.sourceforge.net/octave/?rev=6414&view=rev Author: paramaniac Date: 2009-10-30 19:15:28 +0000 (Fri, 30 Oct 2009) Log Message: ----------- control-oo: fix bug in tf Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/@tf/tf.m trunk/octave-forge/extra/control-oo/inst/control/issample.m Modified: trunk/octave-forge/extra/control-oo/inst/@tf/tf.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@tf/tf.m 2009-10-30 18:43:08 UTC (rev 6413) +++ trunk/octave-forge/extra/control-oo/inst/@tf/tf.m 2009-10-30 19:15:28 UTC (rev 6414) @@ -40,7 +40,7 @@ num = {}; den = {}; tsam = -1; - tfvar = "x"; % undefined + tfvar = "x"; # undefined case 1 if (isa (num, "tf")) # already in tf form @@ -56,7 +56,7 @@ [p, m] = size (num); den = tfpolyones (p, m); tsam = -1; - tfvar = "x"; % undefined + tfvar = "x"; # undefined elseif (ischar (num)) # s = tf ("s") tfvar = num; num = __conv2tfpolycell__ ([1, 0]); @@ -84,14 +84,16 @@ den = __conv2tfpolycell__ (den); argc = numel (varargin); - if (isscalar (varargin{1}) && (varargin{1} == abs (varargin{1}))) # sys = tf (num, den, tsam, "prop1, "val1", ...) + if (issample (varargin{1}, 1)) # sys = tf (num, den, tsam, "prop1, "val1", ...) tsam = varargin{1}; argc -= 1; - if (varargin{1} != 0) + if (varargin{1} == 0) + tfvar = "s"; + elseif (varargin{1} == -1) + tfvar = "x"; + else tfvar = "z"; - else - tfvar = "s"; endif if (argc > 0) Modified: trunk/octave-forge/extra/control-oo/inst/control/issample.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/issample.m 2009-10-30 18:43:08 UTC (rev 6413) +++ trunk/octave-forge/extra/control-oo/inst/control/issample.m 2009-10-30 19:15:28 UTC (rev 6414) @@ -27,12 +27,16 @@ ## Date: September 2009 ## Version: 0.1 -function bool = issample (tsam) +function bool = issample (tsam, flg = 0) - if (nargin != 1) + if (nargin < 1 || nargin > 2) print_usage (); endif - bool = (isscalar (tsam) && (tsam == abs (tsam)) && (tsam != 0)); + if (flg) + bool = (isscalar (tsam) && (tsam >= 0 || tsam == -1)); + else + bool = (isscalar (tsam) && (tsam == abs (tsam)) && (tsam != 0)); + endif endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-10-31 11:10:38
|
Revision: 6416 http://octave.svn.sourceforge.net/octave/?rev=6416&view=rev Author: paramaniac Date: 2009-10-31 11:10:28 +0000 (Sat, 31 Oct 2009) Log Message: ----------- control-oo: TF to SS conversion (SISO case) Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/@tf/__sys2ss__.m Added Paths: ----------- trunk/octave-forge/extra/control-oo/inst/ocst/__tf2ss__.m Modified: trunk/octave-forge/extra/control-oo/inst/@tf/__sys2ss__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@tf/__sys2ss__.m 2009-10-30 20:07:31 UTC (rev 6415) +++ trunk/octave-forge/extra/control-oo/inst/@tf/__sys2ss__.m 2009-10-31 11:10:28 UTC (rev 6416) @@ -24,8 +24,22 @@ function [retsys, retlti] = __sys2ss__ (sys) - error ("tf: tf2ss: not implemented yet"); + if (! issiso (sys)) + error ("tf: tf2ss: MIMO case not implemented yet"); + endif + [num, den] = tfdata (sys); + + num = num{1, 1}; + den = den{1, 1}; + + ## tfpoly ensures that there are no leading zeros + if (length (num) > length (den)) + error ("tf: tf2ss: system must be proper"); + endif + + [a, b, c, d] = __tf2ss__ (num, den); + retsys = ss (a, b, c, d); retlti = sys.lti; # preserve lti properties Added: trunk/octave-forge/extra/control-oo/inst/ocst/__tf2ss__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__tf2ss__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__tf2ss__.m 2009-10-31 11:10:28 UTC (rev 6416) @@ -0,0 +1,134 @@ +## Copyright (C) 1996, 1998, 2000, 2002, 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} {[@var{a}, @var{b}, @var{c}, @var{d}] =} tf2ss (@var{num}, @var{den}) +## Conversion from transfer function to state-space. +## The state space system: +## @iftex +## @tex +## $$ \dot x = Ax + Bu $$ +## $$ y = Cx + Du $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## . +## x = Ax + Bu +## y = Cx + Du +## @end example +## @end ifinfo +## is obtained from a transfer function: +## @iftex +## @tex +## $$ G(s) = { { \rm num }(s) \over { \rm den }(s) } $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## num(s) +## G(s)=------- +## den(s) +## @end example +## @end ifinfo +## +## The vector @var{den} must contain only one row, whereas the vector +## @var{num} may contain as many rows as there are outputs @var{y} of +## the system. The state space system matrices obtained from this function +## will be in controllable canonical form as described in @cite{Modern Control +## Theory}, (Brogan, 1991). +## @end deftypefn + +## Author: R. Bruce Tenison <bte...@en...> +## Created: June 22, 1994 +## mod A S Hodel July, Aug 1995 + +function [a, b, c, d] = __tf2ss__ (num, den) + + if (nargin != 2) + print_usage (); + elseif (isempty (num)) + error ("tf2ss: empty numerator"); + elseif (isempty (den)) + error ("tf2ss: empy denominator"); + elseif (! isvector (num)) + error ("num(%dx%d) must be a vector", rows (num), columns (num)); + elseif (! isvector (den)) + error ("den(%dx%d) must be a vector", rows (den), columns (den)); + endif + + ## strip leading zeros from num, den + nz = find (num != 0); + if (isempty (nz)) + num = 0; + else + num = num(nz(1):length(num)); + endif + nz = find (den != 0); + if (isempty (nz)) + error ("denominator is 0."); + else + den = den(nz(1):length(den)); + endif + + ## force num, den to be row vectors + num = vec (num)'; + den = vec (den)'; + nn = length (num); + nd = length (den); + if (nn > nd) + error ("deg(num)=%d > deg(den)= %d", nn, nd); + endif + + ## Check sizes + if (nd == 1) + a = b = c = []; + d = num(:,1) / den(1); + else + ## Pad num so that length(num) = length(den) + if (nd-nn > 0) + num = [zeros(1,nd-nn), num]; + endif + + ## Normalize the numerator and denominator vector w.r.t. the leading + ## coefficient + d1 = den(1); + num = num / d1; + den = den(2:nd)/d1; + sw = nd-1:-1:1; + + ## Form the A matrix + if (nd > 2) + a = [zeros(nd-2,1), eye(nd-2,nd-2); -den(sw)]; + else + a = -den(sw); + endif + + ## Form the B matrix + b = zeros (nd-1, 1); + b(nd-1,1) = 1; + + ## Form the C matrix + c = num(:,2:nd)-num(:,1)*den; + c = c(:,sw); + + ## Form the D matrix + d = num(:,1); + endif + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-11-05 15:58:14
|
Revision: 6438 http://octave.svn.sourceforge.net/octave/?rev=6438&view=rev Author: paramaniac Date: 2009-11-05 15:58:05 +0000 (Thu, 05 Nov 2009) Log Message: ----------- control-oo: use horzcat instead of sprinft for better performance Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/@ss/display.m trunk/octave-forge/extra/control-oo/inst/@tf/display.m trunk/octave-forge/extra/control-oo/inst/@tfpoly/tfpoly2str.m Modified: trunk/octave-forge/extra/control-oo/inst/@ss/display.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@ss/display.m 2009-11-05 14:29:25 UTC (rev 6437) +++ trunk/octave-forge/extra/control-oo/inst/@ss/display.m 2009-11-05 15:58:05 UTC (rev 6438) @@ -76,7 +76,7 @@ MAX_LEN = 12; # max length of row name and column name [mrows, mcols] = size (m); - disp (sprintf ("%s =", mname)); + disp ([mname, " ="]); row_name = strjust (strvcat (" ", rname), "left"); row_name = row_name(:, 1 : min (MAX_LEN, end)); Modified: trunk/octave-forge/extra/control-oo/inst/@tf/display.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@tf/display.m 2009-11-05 14:29:25 UTC (rev 6437) +++ trunk/octave-forge/extra/control-oo/inst/@tf/display.m 2009-11-05 15:58:05 UTC (rev 6438) @@ -44,7 +44,7 @@ disp (""); for m = 1 : nu - disp (sprintf ("Transfer function from input ""%s"" to output ...", inname{m})); + disp (["Transfer function from input \"", inname{m}, "\" to output ..."]); for p = 1 : ny dispfrac (sys.num{p, m}, sys.den{p, m}, sys.tfvar, outname{p}); endfor Modified: trunk/octave-forge/extra/control-oo/inst/@tfpoly/tfpoly2str.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@tfpoly/tfpoly2str.m 2009-11-05 14:29:25 UTC (rev 6437) +++ trunk/octave-forge/extra/control-oo/inst/@tfpoly/tfpoly2str.m 2009-11-05 15:58:05 UTC (rev 6438) @@ -42,16 +42,12 @@ endif if (lp == 1) - if (abs (a) != 1) - str = sprintf ("%s%s", cs, num2str (abs (a), 4)); - else - str = sprintf ("%s%s", cs, num2str (abs (a), 4)); - endif + str = [cs, num2str(abs (a), 4)]; else if (abs (a) != 1) - str = sprintf ("%s%s %s", cs, coeff(a), variab (tfvar, lp-1)); + str = [cs, coeff(a), " ", variab(tfvar, lp-1)]; else - str = sprintf ("%s%s%s", cs, coeff(a), variab (tfvar, lp-1)); + str = [cs, coeff(a), variab(tfvar, lp-1)]; endif endif @@ -68,9 +64,9 @@ endif if (abs (a) != 1) - str = sprintf ("%s%s%s %s", str, cs, coeff(a), variab (tfvar, lp-k)); + str = [str, cs, coeff(a), " ", variab(tfvar, lp-k)]; else - str = sprintf ("%s%s%s%s", str, cs, coeff(a), variab (tfvar, lp-k)); + str = [str, cs, coeff(a), variab(tfvar, lp-k)]; endif endif endfor @@ -85,7 +81,7 @@ cs = " + "; endif - str = sprintf ("%s%s%s", str, cs, num2str (abs (a), 4)); + str = [str, cs, num2str(abs (a), 4)]; endif endif endif @@ -111,7 +107,7 @@ if (n == 1) str = tfvar; else - str = sprintf ("%s^%d", tfvar, n); + str = [tfvar, "^", num2str(n)]; endif endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-11-11 08:43:55
|
Revision: 6459 http://octave.svn.sourceforge.net/octave/?rev=6459&view=rev Author: paramaniac Date: 2009-11-11 08:43:48 +0000 (Wed, 11 Nov 2009) Log Message: ----------- control-oo: add todo's Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/@lti/norm.m trunk/octave-forge/extra/control-oo/inst/@ss/__zero__.m trunk/octave-forge/extra/control-oo/inst/control/place.m Modified: trunk/octave-forge/extra/control-oo/inst/@lti/norm.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@lti/norm.m 2009-11-10 13:41:22 UTC (rev 6458) +++ trunk/octave-forge/extra/control-oo/inst/@lti/norm.m 2009-11-11 08:43:48 UTC (rev 6459) @@ -33,6 +33,8 @@ ## Date: November 2009 ## Version: 0.1 +## TODO: Use Fortran code from Slicot + function gain = norm (sys, ntype = "2", tol = 0.001) if (nargin > 3) # norm () is catched by built-in function Modified: trunk/octave-forge/extra/control-oo/inst/@ss/__zero__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@ss/__zero__.m 2009-11-10 13:41:22 UTC (rev 6458) +++ trunk/octave-forge/extra/control-oo/inst/@ss/__zero__.m 2009-11-11 08:43:48 UTC (rev 6459) @@ -22,6 +22,8 @@ ## Created: October 2009 ## Version: 0.1 +## TODO: Use Fortran code from Slicot + function [zer, gain] = __zero__ (sys) warning ("ss: zero: subroutine tzero is buggy, use results with caution"); Modified: trunk/octave-forge/extra/control-oo/inst/control/place.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/place.m 2009-11-10 13:41:22 UTC (rev 6458) +++ trunk/octave-forge/extra/control-oo/inst/control/place.m 2009-11-11 08:43:48 UTC (rev 6459) @@ -47,6 +47,8 @@ ## Date: October 2009 ## Version: 0.1 +## TODO: Support MIMO systems, use Fortran algorithm from Slicot + function K = place (argin1, argin2, argin3) if (nargin == 3) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-11-20 20:18:26
|
Revision: 6493 http://octave.svn.sourceforge.net/octave/?rev=6493&view=rev Author: paramaniac Date: 2009-11-20 20:18:16 +0000 (Fri, 20 Nov 2009) Log Message: ----------- control-oo: minor changes Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m trunk/octave-forge/extra/control-oo/inst/@ss/display.m trunk/octave-forge/extra/control-oo/inst/@tf/__freqresp__.m Modified: trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m 2009-11-20 15:11:45 UTC (rev 6492) +++ trunk/octave-forge/extra/control-oo/inst/@ss/__freqresp__.m 2009-11-20 20:18:16 UTC (rev 6493) @@ -16,6 +16,7 @@ ## along with this program. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- +## Frequency response of SS models. ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 @@ -32,10 +33,10 @@ error ("ss: freqresp: system must be square for response type %d", resptype); endif - if (Ts == 0 || Ts == -1) # continuous system + if (Ts > 0) # discrete system + s = exp (i * w * Ts); + else # continuous system s = i * w; - else # discrete system - s = exp (i * w * Ts); endif l_s = length (s); Modified: trunk/octave-forge/extra/control-oo/inst/@ss/display.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@ss/display.m 2009-11-20 15:11:45 UTC (rev 6492) +++ trunk/octave-forge/extra/control-oo/inst/@ss/display.m 2009-11-20 20:18:16 UTC (rev 6493) @@ -76,8 +76,6 @@ MAX_LEN = 12; # max length of row name and column name [mrows, mcols] = size (m); - disp ([mname, " ="]); - row_name = strjust (strvcat (" ", rname), "left"); row_name = row_name(:, 1 : min (MAX_LEN, end)); row_name = horzcat (repmat (" ", mrows+1, 3), row_name); @@ -98,6 +96,8 @@ endfor mat = horzcat (row_name, mat{:}); + + disp ([mname, " ="]); disp (mat); disp (""); Modified: trunk/octave-forge/extra/control-oo/inst/@tf/__freqresp__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@tf/__freqresp__.m 2009-11-20 15:11:45 UTC (rev 6492) +++ trunk/octave-forge/extra/control-oo/inst/@tf/__freqresp__.m 2009-11-20 20:18:16 UTC (rev 6493) @@ -16,6 +16,7 @@ ## along with this program. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- +## Frequency response of TF models. ## Author: Lukas Reichlin <luk...@gm...> ## Created: October 2009 @@ -26,10 +27,10 @@ [p, m] = size (sys); [num, den, Ts] = tfdata (sys); - if (Ts == 0 || Ts == -1) # continuous system + if (Ts > 0) # discrete system + s = exp (i * w * Ts); + else # continuous system s = i * w; - else # discrete system - s = exp (i * w * Ts); endif l_s = length (s); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-11-27 18:18:37
|
Revision: 6546 http://octave.svn.sourceforge.net/octave/?rev=6546&view=rev Author: paramaniac Date: 2009-11-27 18:18:29 +0000 (Fri, 27 Nov 2009) Log Message: ----------- control-oo: move test for ss zeros to ltimodels.m Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/control/ltimodels.m Removed Paths: ------------- trunk/octave-forge/extra/control-oo/inst/ocst/__test_tzero__.m Modified: trunk/octave-forge/extra/control-oo/inst/control/ltimodels.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/ltimodels.m 2009-11-27 18:13:05 UTC (rev 6545) +++ trunk/octave-forge/extra/control-oo/inst/control/ltimodels.m 2009-11-27 18:18:29 UTC (rev 6546) @@ -129,3 +129,41 @@ %! Hinf = norm (sys, inf); %!assert (H2, 1.2527, 1.5e-5); %!assert (Hinf, 2.7, 0.1); + +## transmission zeros of state-space models + +## Results from the "Dark Side" 7.5 and 7.8 +## +## -13.2759 +## 12.5774 +## -0.0155 + +## Results from Scilab 5.2.0b1 (trzeros) +## +## - 13.275931 +## 12.577369 +## - 0.0155265 + +%!shared z, z_exp +%! A = [ -0.7 -0.0458 -12.2 0 +%! 0 -0.014 -0.2904 -0.562 +%! 1 -0.0057 -1.4 0 +%! 1 0 0 0 ]; +%! +%! B = [ -19.1 -3.1 +%! -0.0119 -0.0096 +%! -0.14 -0.72 +%! 0 0 ]; +%! +%! C = [ 0 0 -1 1 +%! 0 0 0.733 0 ]; +%! +%! D = [ 0 0 +%! 0.0768 0.1134 ]; +%! +%! sys = ss (A, B, C, D); +%! z = sort (zero (sys)); +%! +%! z_exp = sort ([-13.2759; 12.5774; -0.0155]); +%! +%!assert (z, z_exp, 1e-4); Deleted: trunk/octave-forge/extra/control-oo/inst/ocst/__test_tzero__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__test_tzero__.m 2009-11-27 18:13:05 UTC (rev 6545) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__test_tzero__.m 2009-11-27 18:18:29 UTC (rev 6546) @@ -1,41 +0,0 @@ -## There's a bug in __tzero__.m as well as in the -## original function tzero.m from the OCST (control-1.0.17) - - -## Results from the "Dark Side" 7.5 and 7.8 -## -## -13.2759 -## 12.5774 -## -0.0155 - - -## Results from Scilab 5.2.0b1 (trzeros) -## -## - 13.275931 -## 12.577369 -## - 0.0155265 - - -%!shared z, z_exp -%! A = [ -0.7 -0.0458 -12.2 0 -%! 0 -0.014 -0.2904 -0.562 -%! 1 -0.0057 -1.4 0 -%! 1 0 0 0 ]; -%! -%! B = [ -19.1 -3.1 -%! -0.0119 -0.0096 -%! -0.14 -0.72 -%! 0 0 ]; -%! -%! C = [ 0 0 -1 1 -%! 0 0 0.733 0 ]; -%! -%! D = [ 0 0 -%! 0.0768 0.1134 ]; -%! -%! sys = ss (A, B, C, D); -%! z = sort (zero (sys)); -%! -%! z_exp = sort ([-13.2759; 12.5774; -0.0155]); -%! -%!assert (z, z_exp, 1e-4); \ 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-12-01 21:25:37
|
Revision: 6566 http://octave.svn.sourceforge.net/octave/?rev=6566&view=rev Author: paramaniac Date: 2009-12-01 21:25:26 +0000 (Tue, 01 Dec 2009) Log Message: ----------- control-oo: move tf2ss to control Added Paths: ----------- trunk/octave-forge/extra/control-oo/inst/control/__tf2ss__.m Removed Paths: ------------- trunk/octave-forge/extra/control-oo/inst/ocst/__tf2ss__.m Copied: trunk/octave-forge/extra/control-oo/inst/control/__tf2ss__.m (from rev 6565, trunk/octave-forge/extra/control-oo/inst/ocst/__tf2ss__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/__tf2ss__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/control/__tf2ss__.m 2009-12-01 21:25:26 UTC (rev 6566) @@ -0,0 +1,134 @@ +## Copyright (C) 1996, 1998, 2000, 2002, 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} {[@var{a}, @var{b}, @var{c}, @var{d}] =} tf2ss (@var{num}, @var{den}) +## Conversion from transfer function to state-space. +## The state space system: +## @iftex +## @tex +## $$ \dot x = Ax + Bu $$ +## $$ y = Cx + Du $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## . +## x = Ax + Bu +## y = Cx + Du +## @end example +## @end ifinfo +## is obtained from a transfer function: +## @iftex +## @tex +## $$ G(s) = { { \rm num }(s) \over { \rm den }(s) } $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## num(s) +## G(s)=------- +## den(s) +## @end example +## @end ifinfo +## +## The vector @var{den} must contain only one row, whereas the vector +## @var{num} may contain as many rows as there are outputs @var{y} of +## the system. The state space system matrices obtained from this function +## will be in controllable canonical form as described in @cite{Modern Control +## Theory}, (Brogan, 1991). +## @end deftypefn + +## Author: R. Bruce Tenison <bte...@en...> +## Created: June 22, 1994 +## mod A S Hodel July, Aug 1995 + +function [a, b, c, d] = __tf2ss__ (num, den) + + if (nargin != 2) + print_usage (); + elseif (isempty (num)) + error ("tf2ss: empty numerator"); + elseif (isempty (den)) + error ("tf2ss: empy denominator"); + elseif (! isvector (num)) + error ("num(%dx%d) must be a vector", rows (num), columns (num)); + elseif (! isvector (den)) + error ("den(%dx%d) must be a vector", rows (den), columns (den)); + endif + + ## strip leading zeros from num, den + nz = find (num != 0); + if (isempty (nz)) + num = 0; + else + num = num(nz(1):length(num)); + endif + nz = find (den != 0); + if (isempty (nz)) + error ("denominator is 0."); + else + den = den(nz(1):length(den)); + endif + + ## force num, den to be row vectors + num = vec (num)'; + den = vec (den)'; + nn = length (num); + nd = length (den); + if (nn > nd) + error ("deg(num)=%d > deg(den)= %d", nn, nd); + endif + + ## Check sizes + if (nd == 1) + a = b = c = []; + d = num(:,1) / den(1); + else + ## Pad num so that length(num) = length(den) + if (nd-nn > 0) + num = [zeros(1,nd-nn), num]; + endif + + ## Normalize the numerator and denominator vector w.r.t. the leading + ## coefficient + d1 = den(1); + num = num / d1; + den = den(2:nd)/d1; + sw = nd-1:-1:1; + + ## Form the A matrix + if (nd > 2) + a = [zeros(nd-2,1), eye(nd-2,nd-2); -den(sw)]; + else + a = -den(sw); + endif + + ## Form the B matrix + b = zeros (nd-1, 1); + b(nd-1,1) = 1; + + ## Form the C matrix + c = num(:,2:nd)-num(:,1)*den; + c = c(:,sw); + + ## Form the D matrix + d = num(:,1); + endif + +endfunction Deleted: trunk/octave-forge/extra/control-oo/inst/ocst/__tf2ss__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ocst/__tf2ss__.m 2009-12-01 17:30:41 UTC (rev 6565) +++ trunk/octave-forge/extra/control-oo/inst/ocst/__tf2ss__.m 2009-12-01 21:25:26 UTC (rev 6566) @@ -1,134 +0,0 @@ -## Copyright (C) 1996, 1998, 2000, 2002, 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} {[@var{a}, @var{b}, @var{c}, @var{d}] =} tf2ss (@var{num}, @var{den}) -## Conversion from transfer function to state-space. -## The state space system: -## @iftex -## @tex -## $$ \dot x = Ax + Bu $$ -## $$ y = Cx + Du $$ -## @end tex -## @end iftex -## @ifinfo -## @example -## . -## x = Ax + Bu -## y = Cx + Du -## @end example -## @end ifinfo -## is obtained from a transfer function: -## @iftex -## @tex -## $$ G(s) = { { \rm num }(s) \over { \rm den }(s) } $$ -## @end tex -## @end iftex -## @ifinfo -## @example -## num(s) -## G(s)=------- -## den(s) -## @end example -## @end ifinfo -## -## The vector @var{den} must contain only one row, whereas the vector -## @var{num} may contain as many rows as there are outputs @var{y} of -## the system. The state space system matrices obtained from this function -## will be in controllable canonical form as described in @cite{Modern Control -## Theory}, (Brogan, 1991). -## @end deftypefn - -## Author: R. Bruce Tenison <bte...@en...> -## Created: June 22, 1994 -## mod A S Hodel July, Aug 1995 - -function [a, b, c, d] = __tf2ss__ (num, den) - - if (nargin != 2) - print_usage (); - elseif (isempty (num)) - error ("tf2ss: empty numerator"); - elseif (isempty (den)) - error ("tf2ss: empy denominator"); - elseif (! isvector (num)) - error ("num(%dx%d) must be a vector", rows (num), columns (num)); - elseif (! isvector (den)) - error ("den(%dx%d) must be a vector", rows (den), columns (den)); - endif - - ## strip leading zeros from num, den - nz = find (num != 0); - if (isempty (nz)) - num = 0; - else - num = num(nz(1):length(num)); - endif - nz = find (den != 0); - if (isempty (nz)) - error ("denominator is 0."); - else - den = den(nz(1):length(den)); - endif - - ## force num, den to be row vectors - num = vec (num)'; - den = vec (den)'; - nn = length (num); - nd = length (den); - if (nn > nd) - error ("deg(num)=%d > deg(den)= %d", nn, nd); - endif - - ## Check sizes - if (nd == 1) - a = b = c = []; - d = num(:,1) / den(1); - else - ## Pad num so that length(num) = length(den) - if (nd-nn > 0) - num = [zeros(1,nd-nn), num]; - endif - - ## Normalize the numerator and denominator vector w.r.t. the leading - ## coefficient - d1 = den(1); - num = num / d1; - den = den(2:nd)/d1; - sw = nd-1:-1:1; - - ## Form the A matrix - if (nd > 2) - a = [zeros(nd-2,1), eye(nd-2,nd-2); -den(sw)]; - else - a = -den(sw); - endif - - ## Form the B matrix - b = zeros (nd-1, 1); - b(nd-1,1) = 1; - - ## Form the C matrix - c = num(:,2:nd)-num(:,1)*den; - c = c(:,sw); - - ## Form the D matrix - d = num(:,1); - endif - -endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-12-02 13:39:55
|
Revision: 6576 http://octave.svn.sourceforge.net/octave/?rev=6576&view=rev Author: paramaniac Date: 2009-12-02 13:39:43 +0000 (Wed, 02 Dec 2009) Log Message: ----------- control-oo: delete subdirectories inst/control and inst/examples. PKG_ADD caused problems with the Makefile when oct-files were present. Added Paths: ----------- trunk/octave-forge/extra/control-oo/inst/__axis2dlim__.m trunk/octave-forge/extra/control-oo/inst/__checkname__.m trunk/octave-forge/extra/control-oo/inst/__conv2tfpolycell__.m trunk/octave-forge/extra/control-oo/inst/__freqbounds__.m trunk/octave-forge/extra/control-oo/inst/__getfreqresp__.m trunk/octave-forge/extra/control-oo/inst/__markemptynames__.m trunk/octave-forge/extra/control-oo/inst/__ssmatdim__.m trunk/octave-forge/extra/control-oo/inst/__tf2ss__.m trunk/octave-forge/extra/control-oo/inst/__tfnddim__.m trunk/octave-forge/extra/control-oo/inst/__timeresp__.m trunk/octave-forge/extra/control-oo/inst/are.m trunk/octave-forge/extra/control-oo/inst/bode.m trunk/octave-forge/extra/control-oo/inst/bodemag.m trunk/octave-forge/extra/control-oo/inst/care.m trunk/octave-forge/extra/control-oo/inst/ctrb.m trunk/octave-forge/extra/control-oo/inst/dare.m trunk/octave-forge/extra/control-oo/inst/dcgain.m trunk/octave-forge/extra/control-oo/inst/dlqr.m trunk/octave-forge/extra/control-oo/inst/dlyap.m trunk/octave-forge/extra/control-oo/inst/estim.m trunk/octave-forge/extra/control-oo/inst/freqresp.m trunk/octave-forge/extra/control-oo/inst/gensig.m trunk/octave-forge/extra/control-oo/inst/gram.m trunk/octave-forge/extra/control-oo/inst/impulse.m trunk/octave-forge/extra/control-oo/inst/initial.m trunk/octave-forge/extra/control-oo/inst/isctrb.m trunk/octave-forge/extra/control-oo/inst/isdetectable.m trunk/octave-forge/extra/control-oo/inst/isobsv.m trunk/octave-forge/extra/control-oo/inst/issample.m trunk/octave-forge/extra/control-oo/inst/isstabilizable.m trunk/octave-forge/extra/control-oo/inst/kalman.m trunk/octave-forge/extra/control-oo/inst/lqr.m trunk/octave-forge/extra/control-oo/inst/lsim.m trunk/octave-forge/extra/control-oo/inst/ltimodels.m trunk/octave-forge/extra/control-oo/inst/lyap.m trunk/octave-forge/extra/control-oo/inst/margin.m trunk/octave-forge/extra/control-oo/inst/nichols.m trunk/octave-forge/extra/control-oo/inst/nyquist.m trunk/octave-forge/extra/control-oo/inst/obsv.m trunk/octave-forge/extra/control-oo/inst/optiPID.m trunk/octave-forge/extra/control-oo/inst/optiPIDfunction.m trunk/octave-forge/extra/control-oo/inst/place.m trunk/octave-forge/extra/control-oo/inst/pzmap.m trunk/octave-forge/extra/control-oo/inst/sigma.m trunk/octave-forge/extra/control-oo/inst/step.m trunk/octave-forge/extra/control-oo/inst/strseq.m trunk/octave-forge/extra/control-oo/inst/tfpolyones.m trunk/octave-forge/extra/control-oo/inst/tfpolyzeros.m Removed Paths: ------------- trunk/octave-forge/extra/control-oo/inst/control/ trunk/octave-forge/extra/control-oo/inst/examples/ Copied: trunk/octave-forge/extra/control-oo/inst/__axis2dlim__.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/__axis2dlim__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/__axis2dlim__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/__axis2dlim__.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -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 Copied: trunk/octave-forge/extra/control-oo/inst/__checkname__.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/__checkname__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/__checkname__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/__checkname__.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,38 @@ +## 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 -*- +## Check whether a cell contains the required number of strings. +## Used by set and __set__. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function name = __checkname__ (name, req_len) + + if (! iscell (name)) # catch the siso case, + name = {name}; # e.g. sys = set (sys, "inname", "u_1") + endif + + n = numel (name); + + if (n != req_len) + error ("lti: set: cell must contain %d strings", req_len); + endif + +endfunction \ No newline at end of file Copied: trunk/octave-forge/extra/control-oo/inst/__conv2tfpolycell__.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/__conv2tfpolycell__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/__conv2tfpolycell__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/__conv2tfpolycell__.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,42 @@ +## 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 -*- +## Convert a (cell of) row vector(s) to a cell of tfpoly objects. +## Used by tf and __set__. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function ndr = __conv2tfpolycell__ (nd) + + if (! iscell (nd)) + nd = {nd}; + endif + + [ndrows, ndcols] = size (nd); + + ndr = cell (ndrows, ndcols); + + for k = 1 : ndrows + for l = 1 : ndcols + ndr{k, l} = tfpoly (nd{k, l}); + endfor + endfor + +endfunction \ No newline at end of file Copied: trunk/octave-forge/extra/control-oo/inst/__freqbounds__.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/__freqbounds__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/__freqbounds__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/__freqbounds__.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,117 @@ +## Copyright (C) 1996, 2000, 2004, 2005, 2006, 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} {[@var{dec_min}, @var{dec_max}] =} __freqbounds__ (@var{sys}) +## Get default range of frequencies based on cutoff frequencies of system +## poles and zeros. +## Frequency range is the interval +## @iftex +## @tex +## $ [ 10^{w_{min}}, 10^{w_{max}} ] $ +## @end tex +## @end iftex +## @ifinfo +## [10^@var{wmin}, 10^@var{wmax}] +## @end ifinfo +## +## Used internally in @command{margin} (@command{bode}, @command{nyquist}) +## @end deftypefn + +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + +function [dec_min, dec_max] = __freqbounds__ (sys, wbounds = "std") + + zer = zero (sys); + pol = pole (sys); + tsam = get (sys, "tsam"); + digital = (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) + ## 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); + + ## avoid dividing empty matrices, it would work but looks nasty + if (! isempty (iiz)) + czer = log (zer(iiz))/tsam; + else + czer = []; + endif + + if (! isempty (iip)) + cpol = log (pol(iip))/tsam; + else + cpol = []; + endif + else + ## continuous + iip = find (abs(pol) > norm(pol)*eps); + iiz = find (abs(zer) > norm(zer)*eps); + + if (! isempty (zer)) + czer = zer(iiz); + else + czer = []; + endif + if (! isempty (pol)) + cpol = pol(iip); + else + cpol = []; + endif + endif + + if (isempty (iip) && isempty (iiz)) + ## no poles/zeros away from omega = 0; pick defaults + dec_min = 0; # -1 + dec_max = 2; # 3 + else + dec_min = floor (log10 (min (abs ([cpol, czer])))); + dec_max = ceil (log10 (max (abs ([cpol, czer])))); + endif + + ## expand to show the entirety of the "interesting" portion of the plot + switch (wbounds) + case "std" # standard + if (dec_min == dec_max) + dec_min -= 2; + dec_max += 2; + else + dec_min--; + dec_max++; + endif + case "ext" # extended (for nyquist) + dec_min -= 2; + dec_max += 2; + otherwise + error ("freqbounds: second argument invalid"); + endswitch + + ## run digital frequency all the way to pi + if (digital) + dec_max = log10 (pi/tsam); + endif + +endfunction Copied: trunk/octave-forge/extra/control-oo/inst/__getfreqresp__.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/__getfreqresp__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/__getfreqresp__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/__getfreqresp__.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,52 @@ +## 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 -*- +## Return frequency response H and frequency vector w. +## If w is empty, it will be calculated by __freqbounds__ + +## Author: Lukas Reichlin <luk...@gm...> +## Created: November 2009 +## Version: 0.1 + +function [H, w] = __getfreqresp__ (sys, w = [], mimoflag = 0, resptype = 0, wbounds = "std") + + ## check arguments + if(! isa (sys, "lti")) + error ("getfreqresp: first argument sys must be a LTI system"); + endif + + if (! isvector (w) && ! isempty (w)) + error ("getfreqresp: second argument w must be a vector of frequencies"); + endif + + if (! mimoflag && ! issiso (sys)) + error ("getfreqresp: require SISO system"); + endif + + ## find interesting frequency range w if not specified + if (isempty (w)) + ## begin plot at 10^dec_min, end plot at 10^dec_max [rad/s] + [dec_min, dec_max] = __freqbounds__ (sys, wbounds); + + w = logspace (dec_min, dec_max, 500); # [rad/s] + endif + + ## frequency response + H = __freqresp__ (sys, w, resptype); + +endfunction \ No newline at end of file Copied: trunk/octave-forge/extra/control-oo/inst/__markemptynames__.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/__markemptynames__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/__markemptynames__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/__markemptynames__.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,34 @@ +## 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 -*- +## Check whether a string of the cell "name" is empty and mark them +## with "?". Used by display routines. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function name = __markemptynames__ (name) + + for k = 1 : numel (name) + if (isempty (name{k})) + name{k} = "?"; + endif + endfor + +endfunction \ No newline at end of file Copied: trunk/octave-forge/extra/control-oo/inst/__ssmatdim__.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/__ssmatdim__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/__ssmatdim__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/__ssmatdim__.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,61 @@ +## 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 -*- +## Number of inputs, states and outputs of state space matrices. +## For internal use only. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2009 +## Version: 0.1 + +function [nu, nx, ny] = __ssmatdim__ (a, b, c, d) + + [arows, acols] = size (a); + [brows, bcols] = size (b); + [crows, ccols] = size (c); + [drows, dcols] = size (d); + + nu = bcols; # = dcols + nx = arows; # = acols + ny = crows; # = drows + + if (! issquare (a) && ! isempty (a)) + error ("ss: system matrix a(%dx%d) is not square", arows, acols); + endif + + if (brows != arows) + error ("ss: system matrices a(%dx%d) and b(%dx%d) are incompatible", + arows, acols, brows, bcols); + endif + + if (ccols != acols) + error ("ss: system matrices a(%dx%d) and c(%dx%d) are incompatible", + arows, acols, crows, ccols); + endif + + if (bcols != dcols) + error ("ss: system matrices b(%dx%d) and d(%dx%d) are incompatible", + brows, bcols, drows, dcols); + endif + + if (crows != drows) + error ("ss: system matrices c(%dx%d) and d(%dx%d) are incompatible", + crows, ccols, drows, dcols); + endif + +endfunction Copied: trunk/octave-forge/extra/control-oo/inst/__tf2ss__.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/__tf2ss__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/__tf2ss__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/__tf2ss__.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,134 @@ +## Copyright (C) 1996, 1998, 2000, 2002, 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} {[@var{a}, @var{b}, @var{c}, @var{d}] =} tf2ss (@var{num}, @var{den}) +## Conversion from transfer function to state-space. +## The state space system: +## @iftex +## @tex +## $$ \dot x = Ax + Bu $$ +## $$ y = Cx + Du $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## . +## x = Ax + Bu +## y = Cx + Du +## @end example +## @end ifinfo +## is obtained from a transfer function: +## @iftex +## @tex +## $$ G(s) = { { \rm num }(s) \over { \rm den }(s) } $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## num(s) +## G(s)=------- +## den(s) +## @end example +## @end ifinfo +## +## The vector @var{den} must contain only one row, whereas the vector +## @var{num} may contain as many rows as there are outputs @var{y} of +## the system. The state space system matrices obtained from this function +## will be in controllable canonical form as described in @cite{Modern Control +## Theory}, (Brogan, 1991). +## @end deftypefn + +## Author: R. Bruce Tenison <bte...@en...> +## Created: June 22, 1994 +## mod A S Hodel July, Aug 1995 + +function [a, b, c, d] = __tf2ss__ (num, den) + + if (nargin != 2) + print_usage (); + elseif (isempty (num)) + error ("tf2ss: empty numerator"); + elseif (isempty (den)) + error ("tf2ss: empy denominator"); + elseif (! isvector (num)) + error ("num(%dx%d) must be a vector", rows (num), columns (num)); + elseif (! isvector (den)) + error ("den(%dx%d) must be a vector", rows (den), columns (den)); + endif + + ## strip leading zeros from num, den + nz = find (num != 0); + if (isempty (nz)) + num = 0; + else + num = num(nz(1):length(num)); + endif + nz = find (den != 0); + if (isempty (nz)) + error ("denominator is 0."); + else + den = den(nz(1):length(den)); + endif + + ## force num, den to be row vectors + num = vec (num)'; + den = vec (den)'; + nn = length (num); + nd = length (den); + if (nn > nd) + error ("deg(num)=%d > deg(den)= %d", nn, nd); + endif + + ## Check sizes + if (nd == 1) + a = b = c = []; + d = num(:,1) / den(1); + else + ## Pad num so that length(num) = length(den) + if (nd-nn > 0) + num = [zeros(1,nd-nn), num]; + endif + + ## Normalize the numerator and denominator vector w.r.t. the leading + ## coefficient + d1 = den(1); + num = num / d1; + den = den(2:nd)/d1; + sw = nd-1:-1:1; + + ## Form the A matrix + if (nd > 2) + a = [zeros(nd-2,1), eye(nd-2,nd-2); -den(sw)]; + else + a = -den(sw); + endif + + ## Form the B matrix + b = zeros (nd-1, 1); + b(nd-1,1) = 1; + + ## Form the C matrix + c = num(:,2:nd)-num(:,1)*den; + c = c(:,sw); + + ## Form the D matrix + d = num(:,1); + endif + +endfunction Copied: trunk/octave-forge/extra/control-oo/inst/__tfnddim__.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/__tfnddim__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/__tfnddim__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/__tfnddim__.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,36 @@ +## 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 -*- +## Number of outputs and inputs of transfer function numerator and +## denominator. For internal use only. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function [nrows, ncols] = __tfnddim__ (num, den) + + [nrows, ncols] = size (num); + [drows, dcols] = size (den); + + if (nrows != drows || ncols != dcols) + error (sprintf ("tf: num(%dx%d) and den(%dx%d) must have equal dimensions", + nrows, ncols, drows, dcols)); + endif + +endfunction \ No newline at end of file Copied: trunk/octave-forge/extra/control-oo/inst/__timeresp__.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/__timeresp__.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/__timeresp__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/__timeresp__.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,309 @@ +## 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 -*- +## Common code for the time response functions step, impulse and initial. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function [y, t, x_arr] = __timeresp__ (sys, resptype, plotflag, tfinal, dt, x0) + + if (! isa (sys, "ss")) + sys = ss (sys); # sys must be proper + endif + + if (! isempty (tfinal) && ! isscalar (tfinal)) # time vector t passed + dt = tfinal(2) - tfinal(1); # assume that t is regularly spaced + tfinal = tfinal(end); + endif + + [A, B, C, D, tsam] = ssdata (sys); + + digital = ! isct (sys); # static gains are treated as analog systems + + if (digital) + if (! isempty (dt)) + warning ("timeresp: argument dt has no effect on sampling time of digital system"); + endif + + dt = tsam; + endif + + [tfinal, dt] = __simhorizon__ (A, digital, tfinal, dt); + + if (! digital) + sys = c2d (sys, dt, "zoh"); + endif + + [F, G] = ssdata (sys); # matrices C and D don't change + + n = rows (F); # number of states + m = columns (G); # number of inputs + p = rows (C); # number of outputs + + ## time vector + t = (0 : dt : tfinal)'; + l_t = length (t); + + switch (resptype) + case "initial" + str = "Response to Initial Conditions"; + yfinal = zeros (p, 1); + + ## preallocate memory + y = zeros (l_t, p); + x_arr = zeros (l_t, n); + + ## initial conditions + x = x0(:); # make sure that x is a column vector + + if (n != length (x0)) + error ("initial: x0 must be a vector with %d elements", n); + endif + + ## simulation + for k = 1 : l_t + y(k, :) = C * x; + x_arr(k, :) = x; + x = F * x; + endfor + + case "step" + str = "Step Response"; + yfinal = dcgain (sys); + + ## preallocate memory + y = zeros (l_t, p, m); + x_arr = zeros (l_t, n, m); + + for j = 1 : m # for every input channel + ## initial conditions + x = zeros (n, 1); + u = zeros (m, 1); + u(j) = 1; + + ## simulation + for k = 1 : l_t + y(k, :, j) = C * x + D * u; + x_arr(k, :, j) = x; + x = F * x + G * u; + endfor + endfor + + case "impulse" + str = "Impulse Response"; + yfinal = zeros (p, m); + + ## preallocate memory + y = zeros (l_t, p, m); + x_arr = zeros (l_t, n, m); + + for j = 1 : m # for every input channel + ## initial conditions + u = zeros (m, 1); + u(j) = 1; + + if (digital) + x = zeros (n, 1); # zero by definition + y(1, :, j) = D * u / dt; + x_arr(1, :, j) = x; + x = G * u / dt; + else + x = B * u; # B, not G! + y(1, :, j) = C * x; + x_arr(1, :, j) = x; + x = F * x; + endif + + ## simulation + for k = 2 : l_t + y (k, :, j) = C * x; + x_arr(k, :, j) = x; + x = F * x; + endfor + endfor + + if (digital) + y *= dt; + x_arr *= dt; + endif + + otherwise + error ("timeresp: invalid response type"); + + endswitch + + + if (plotflag) # display plot + + ## TODO: Set correct titles, especially for multi-input systems + ## Limitations probably due to the Gnuplot backend + + stable = isstable (sys); + outname = get (sys, "outname"); + + if (isempty (outname) || isequal ("", outname{:})) + outname = strseq ("y_", 1 : length (outname)); + else + outname = __markemptynames__ (outname); + endif + + if (strcmp (resptype, "initial")) + cols = 1; + else + cols = m; + endif + + if (digital) # discrete system + for k = 1 : p + for j = 1 : cols + + subplot (p, cols, (k-1)*cols+j); + + if (stable) + stairs (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]); + else + stairs (t, y(:, k, j)); + endif + + grid ("on"); + + if (k == 1 && j == 1) + title (str); + endif + + if (j == 1) + ylabel (sprintf ("Amplitude %s", outname{k})); + endif + + endfor + endfor + + xlabel ("Time [s]"); + + else # continuous system + for k = 1 : p + for j = 1 : cols + + subplot (p, cols, (k-1)*cols+j); + + if (stable) + plot (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]); + else + plot (t, y(:, k, j)); + endif + + grid ("on"); + + if (k == 1 && j == 1) + title (str); + endif + + if (j == 1) + ylabel (sprintf ("Amplitude %s", outname{k})); + endif + + endfor + endfor + + xlabel ("Time [s]"); + + endif + endif + +endfunction + + +function [tfinal, dt] = __simhorizon__ (A, digital, tfinal, Ts) + + ## code based on __stepimp__.m of Kai P. Mueller and A. Scottedward Hodel + + TOL = 1.0e-10; # values below TOL are assumed to be zero + N_MIN = 50; # min number of points + N_MAX = 2000; # max number of points + N_DEF = 1000; # default number of points + T_DEF = 10; # default simulation time + + n = rows (A); + eigw = eig (A); + + if (digital) + ## perform bilinear transformation on poles in z + for k = 1 : n + pol = eigw(k); + if (abs (pol + 1) < TOL) + eigw(k) = 0; + else + eigw(k) = 2 / Ts * (pol - 1) / (pol + 1); + endif + endfor + endif + + ## remove poles near zero from eigenvalue array eigw + nk = n; + for k = 1 : n + if (abs (real (eigw(k))) < TOL) + eigw(k) = 0; + nk -= 1; + endif + endfor + + if (nk == 0) + if (isempty (tfinal)) + tfinal = T_DEF; + endif + + if (! digital) + dt = tfinal / N_DEF; + endif + else + eigw = eigw(find (eigw)); + eigw_max = max (abs (eigw)); + + if (! digital) + dt = 0.2 * pi / eigw_max; + endif + + if (isempty (tfinal)) + eigw_min = min (abs (real (eigw))); + tfinal = 5.0 / eigw_min; + + ## round up + yy = 10^(ceil (log10 (tfinal)) - 1); + tfinal = yy * ceil (tfinal / yy); + endif + + if (! digital) + N = tfinal / dt; + + if (N < N_MIN) + dt = tfinal / N_MIN; + endif + + if (N > N_MAX) + dt = tfinal / N_MAX; + endif + endif + endif + + if (! isempty (Ts)) # catch case analog system with dt specified + dt = Ts; + endif + +endfunction Copied: trunk/octave-forge/extra/control-oo/inst/are.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/are.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/are.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/are.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,133 @@ +## Copyright (C) 1993, 1994, 1995, 2000, 2002, 2004, 2005, 2006, 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} {@var{x} =} are (@var{a}, @var{b}, @var{c}, @var{opt}) +## Solve the Algebraic Riccati Equation +## @iftex +## @tex +## $$ +## A^TX + XA - XBX + C = 0 +## $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## a' * x + x * a - x * b * x + c = 0 +## @end example +## @end ifinfo +## +## @strong{Inputs} +## @noindent +## for identically dimensioned square matrices +## @table @var +## @item a +## @var{n} by @var{n} matrix; +## @item b +## @var{n} by @var{n} matrix or @var{n} by @var{m} matrix; in the latter case +## @var{b} is replaced by @math{b:=b*b'}; +## @item c +## @var{n} by @var{n} matrix or @var{p} by @var{m} matrix; in the latter case +## @var{c} is replaced by @math{c:=c'*c}; +## @item opt +## (optional argument; default = @code{"B"}): +## String option passed to @code{balance} prior to ordered Schur decomposition. +## @end table +## +## @strong{Output} +## @table @var +## @item x +## solution of the @acronym{ARE}. +## @end table +## +## @strong{Method} +## Laub's Schur method (@acronym{IEEE} Transactions on +## Automatic Control, 1979) is applied to the appropriate Hamiltonian +## matrix. +## @seealso{balance, dare} +## @end deftypefn + +## Author: A. S. Hodel <a.s...@en...> +## Created: August 1993 + +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + +function x = are (a, b, c, opt = "B") + + if (nargin < 3 || nargin > 4) + print_usage (); + endif + + if (nargin == 4) + if (ischar (opt)) + opt = upper (opt(1)); + if (opt != "B" && opt != N && opt != "P" && opt != "S") + warning ("are: opt has invalid value ""%s""; setting to ""B""", opt); + opt = "B"; + endif + else + warning ("are: invalid argument opt, setting to ""B"""); + opt = "B"; + endif + endif + + n = issquare (a); + m = issquare (b); + p = issquare (c); + + if (! n) + error ("are: matrix a is not square"); + endif + + if (! isctrb (a, b)) + warning ("are: matrices a and b are not controllable"); + endif + + if (! m) + b = b * b'; # b must be symmetric + m = rows (b); + endif + + if (! isobsv (a, c)) + warning ("are: matrices a and c are not observable"); + endif + + ## to allow lqe design, don't force + ## b to be positive semi-definite + + if (! p) + c = c' * c; # c must be symmetric + p = rows (c); + endif + + if (n != m || n != p) + error ("are: matrices a, b and c not conformably dimensioned"); + endif + + [d, h] = balance ([a, -b; -c, -a'], opt); + [u, s] = schur (h, "A"); + + u = d * u; + n1 = n + 1; + n2 = 2 * n; + + x = u(n1:n2, 1:n) / u(1:n, 1:n); + +endfunction Copied: trunk/octave-forge/extra/control-oo/inst/bode.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/bode.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/bode.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/bode.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,74 @@ +## 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{mag}, @var{pha}, @var{w}] =} bode (@var{sys}) +## @deftypefnx {Function File} {[@var{mag}, @var{pha}, @var{w}] =} bode (@var{sys}, @var{w}) +## Bode diagram of LTI model's frequency response. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: November 2009 +## Version: 0.1 + +function [mag_r, pha_r, w_r] = bode (sys, w = []) + + if (nargin == 0 || nargin > 2) + print_usage (); + endif + + [H, w] = __getfreqresp__ (sys, w, false, 0, "std"); + + H = H(:); + mag = abs (H); + pha = unwrap (arg (H)) * 180 / pi; + + if (! nargout) + mag_db = 20 * log10 (mag); + + wv = [min(w), max(w)]; + ax_vec_mag = __axis2dlim__ ([w(:), mag_db(:)]); + ax_vec_mag(1:2) = wv; + ax_vec_pha = __axis2dlim__ ([w(:), pha(:)]); + ax_vec_pha(1:2) = wv; + + if (isct (sys)) + xl_str = "Frequency [rad/s]"; + else + xl_str = sprintf ("Frequency [rad/s] w_N = %g", pi / get (sys, "tsam")); + endif + + subplot (2, 1, 1) + semilogx (w, mag_db) + axis (ax_vec_mag) + grid ("on") + title ("Bode Diagram") + ylabel ("Magnitude [dB]") + + subplot (2, 1, 2) + semilogx (w, pha) + axis (ax_vec_pha) + grid ("on") + xlabel (xl_str) + ylabel ("Phase [deg]") + else + mag_r = mag; + pha_r = pha; + w_r = w; + endif + +endfunction \ No newline at end of file Copied: trunk/octave-forge/extra/control-oo/inst/bodemag.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/bodemag.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/bodemag.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/bodemag.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,63 @@ +## 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{mag}, @var{w}] =} bodemag (@var{sys}) +## @deftypefnx {Function File} {[@var{mag}, @var{w}] =} bodemag (@var{sys}, @var{w}) +## Bode magnitude diagram of LTI model's frequency response. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: November 2009 +## Version: 0.1 + +function [mag_r, w_r] = bodemag (sys, w = []) + + if (nargin == 0 || nargin > 2) + print_usage (); + endif + + [H, w] = __getfreqresp__ (sys, w, false, 0, "std"); + + H = H(:); + mag = abs (H); + + if (! nargout) + mag_db = 20 * log10 (mag); + + wv = [min(w), max(w)]; + ax_vec_mag = __axis2dlim__ ([w(:), mag_db(:)]); + ax_vec_mag(1:2) = wv; + + if (isct (sys)) + xl_str = "Frequency [rad/s]"; + else + xl_str = sprintf ("Frequency [rad/s] w_N = %g", pi / get (sys, "tsam")); + endif + + semilogx (w, mag_db) + axis (ax_vec_mag) + grid ("on") + title ("Bode Magnitude Diagram") + xlabel (xl_str) + ylabel ("Magnitude [dB]") + else + mag_r = mag; + w_r = w; + endif + +endfunction \ No newline at end of file Copied: trunk/octave-forge/extra/control-oo/inst/care.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/care.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/care.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/care.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,104 @@ +## 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 -*- + +## Author: Lukas Reichlin <luk...@gm...> +## Created: November 2009 +## Version: 0.1 + +function [x, l, g] = care (a, b, q, r, s = [], opt = "B") + + if (nargin < 4 || nargin > 6) + print_usage (); + endif + + if (nargin == 6) + if (ischar (opt)) + opt = upper (opt(1)); + if (opt != "B" && opt != N && opt != "P" && opt != "S") + warning ("dare: opt has invalid value ""%s""; setting to ""B""", opt); + opt = "B"; + endif + else + warning ("dare: invalid argument opt, setting to ""B"""); + opt = "B"; + endif + endif + + [n, m] = size (b); + p = issquare (q); + m1 = issquare (r); + + if (! m1) + error ("care: r is not square"); + elseif (m1 != m) + error ("care: b, r are not conformable"); + endif + + if (! p) + q = q' * q; + endif + + ## incorporate cross term into a and q + if (isempty (s)) + s = zeros (n, m); + ao = a; + qo = q; + else + [n2, m2] = size (s); + + if (n2 != n || m2 != m) + error ("care: s (%dx%d) must be identically dimensioned with b (%dx%d)", + n2, m2, n, m); + endif + + ao = a - (b/r)*s'; + qo = q - (s/r)*s'; + endif + + ## check stabilizability + if (! isstabilizable (ao, b, [], 0)) + error ("care: a and b not stabilizable"); + endif + + ## check detectability + dflag = isdetectable (ao, qo, [], 0); + + if (dflag == 0) + warning ("care: a and q not detectable"); + elseif (dflag == -1) + error ("care: a and q have poles on imaginary axis"); + endif + + ## to allow lqe design, don't force + ## qo to be positive semi-definite + + if (isdefinite (r) <= 0) + error ("care: r must be positive definite"); + endif + + ## unique stabilizing solution + x = are (ao, (b/r)*b', qo, opt); + + ## corresponding gain matrix + g = r \ (b'*x + s'); + + ## closed-loop poles + l = eig (a - b*g); + +endfunction \ No newline at end of file Copied: trunk/octave-forge/extra/control-oo/inst/ctrb.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/ctrb.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/ctrb.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/ctrb.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,56 @@ +## Copyright (C) 1997, 2000, 2002, 2004, 2005, 2006, 2007 Kai P. Mueller +## 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{co} =} ctrb (@var{sys}) +## @deftypefnx {Function File} {@var{co} =} ctrb (@var{a}, @var{b}) +## Controllability matrix. +## @end deftypefn + +## Author: Kai P. Mueller <mu...@if...> +## Created: November 4, 1997 +## based on is_controllable.m of Scottedward Hodel + +function co = ctrb (sys_or_a, b) + + if (nargin == 1) + [a, b] = ssdata (sys_or_a); + elseif (nargin == 2) + a = sys_or_a; + if (! isnumeric (a) || ! isnumeric (b) || + rows(a) != rows (b) || ! issquare (a)) + error ("ctrb: invalid arguments"); + endif + else + print_usage (); + endif + + [arows, acols] = size (a); + [brows, bcols] = size (b); + + co = zeros (arows, acols*bcols); + + for k = 1 : arows + co(:, ((k-1)*bcols + 1) : (k*bcols)) = b; + b = a * b; + endfor + +endfunction + + +%!assert (ctrb ([1, 0; 0, -0.5], [8; 8]), [8, 8; 8, -4]); Copied: trunk/octave-forge/extra/control-oo/inst/dare.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/dare.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/dare.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/dare.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,177 @@ +## Copyright (C) 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2006, 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} {@var{x} =} dare (@var{a}, @var{b}, @var{q}, @var{r}, @var{opt}) +## +## Return the solution, @var{x} of the discrete-time algebraic Riccati +## equation +## @iftex +## @tex +## $$ +## A^TXA - X + A^TXB (R + B^TXB)^{-1} B^TXA + Q = 0 +## $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## a' x a - x + a' x b (r + b' x b)^(-1) b' x a + q = 0 +## @end example +## @end ifinfo +## @noindent +## +## @strong{Inputs} +## @table @var +## @item a +## @var{n} by @var{n} matrix; +## +## @item b +## @var{n} by @var{m} matrix; +## +## @item q +## @var{n} by @var{n} matrix, symmetric positive semidefinite, or a @var{p} by @var{n} matrix, +## In the latter case @math{q:=q'*q} is used; +## +## @item r +## @var{m} by @var{m}, symmetric positive definite (invertible); +## +## @item opt +## (optional argument; default = @code{"B"}): +## String option passed to @code{balance} prior to ordered @var{QZ} decomposition. +## @end table +## +## @strong{Output} +## @table @var +## @item x +## solution of @acronym{DARE}. +## @end table +## +## @strong{Method} +## Generalized eigenvalue approach (Van Dooren; @acronym{SIAM} J. +## Sci. Stat. Comput., Vol 2) applied to the appropriate symplectic pencil. +## +## See also: Ran and Rodman, @cite{Stable Hermitian Solutions of Discrete +## Algebraic Riccati Equations}, Mathematics of Control, Signals and +## Systems, Vol 5, no 2 (1992), pp 165--194. +## @seealso{balance, are} +## @end deftypefn + +## Author: A. S. Hodel <a.s...@en...> +## Created: August 1993 +## Adapted-By: jwe + +## Adapted-By: Lukas Reichlin <luk...@gm...> +## Date: October 2009 +## Version: 0.1 + +function [x, l, g] = dare (a, b, q, r, s = [], opt = "B") + + if (nargin < 4 || nargin > 6) + print_usage (); + endif + + if (nargin == 6) + if (ischar (opt)) + opt = upper (opt(1)); + if (opt != "B" && opt != N && opt != "P" && opt != "S") + warning ("dare: opt has invalid value ""%s""; setting to ""B""", opt); + opt = "B"; + endif + else + warning ("dare: invalid argument opt, setting to ""B"""); + opt = "B"; + endif + endif + + [n, m] = size (b); + p = issquare (q); + m1 = issquare (r); + + if (! m1) + error ("dare: r is not square"); + elseif (m1 != m) + error ("dare: b, r are not conformable"); + endif + + if (! p) + q = q' * q; + endif + + ## incorporate cross term into a and q. + if (isempty (s)) + s = zeros (n, m); + ao = a; + qo = q; + else + [n2, m2] = size (s); + + if (n2 != n || m2 != m) + error ("dare: s (%dx%d) must be identically dimensioned with b (%dx%d)", + n2, m2, n, m); + endif + + ao = a - (b/r)*s'; + qo = q - (s/r)*s'; + endif + + ## check stabilizability + if (! isstabilizable (ao, b, [], 1)) + error ("dare: a and b not stabilizable"); + endif + + ## check detectability + dflag = isdetectable (ao, qo, [], 1); + + if (dflag == 0) + warning ("dare: (a,q) not detectable"); + elseif (dflag == -1) + error ("dare: (a,q) has non-minimal modes near unit circle"); + endif + + ## Checking positive definiteness + if (isdefinite (r) <= 0) + error ("dare: r must be positive definite"); + endif + + ## if (isdefinite (qo) < 0) + ## error ("dare: q not positive semidefinite"); + ## endif + + ## solve the riccati equation + s1 = [ ao, zeros(n); + -qo, eye(n) ]; + s2 = [ eye(n), (b/r)*b'; + zeros(n), ao' ]; + + [c, d, s1, s2] = balance (s1, s2, opt); + [aa, bb, u, lam] = qz (s1, s2, "S"); + + u = d*u; + n1 = n+1; + n2 = 2*n; + + ## unique stabilizing solution + x = u(n1:n2, 1:n) / u(1:n, 1:n); + + ## corresponding gain matrix + g = (r+b'*x*b) \ (b'*x*a + s'); + + ## closed-loop poles + l = eig (a - b*g); + +endfunction Copied: trunk/octave-forge/extra/control-oo/inst/dcgain.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/dcgain.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/dcgain.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/dcgain.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,35 @@ +## 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. +## @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 Copied: trunk/octave-forge/extra/control-oo/inst/dlqr.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/dlqr.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/dlqr.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/dlqr.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,47 @@ +## 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 -*- + +## Author: Lukas Reichlin <luk...@gm...> +## Created: November 2009 +## Version: 0.1 + +function [g, x, l] = dlqr (a, b, q, r = [], s = []) + + if (nargin < 3 || nargin > 5) + print_usage (); + endif + + if (isa (a, "lti")) + s = r; + r = q; + q = b; + [a, b, c, d, tsam] = ssdata (a); + elseif (nargin < 4) + print_usage (); + else + tsam = 1; # any value > 0 could be used here + endif + + if (tsam > 0) + [x, l, g] = dare (a, b, q, r, s); + else + [x, l, g] = care (a, b, q, r, s); + endif + +endfunction \ No newline at end of file Copied: trunk/octave-forge/extra/control-oo/inst/dlyap.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/dlyap.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/dlyap.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/dlyap.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,159 @@ +## Copyright (C) 1993, 1994, 1995, 2000, 2002, 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} {} dlyap (@var{a}, @var{b}) +## Solve the discrete-time Lyapunov equation +## +## @strong{Inputs} +## @table @var +## @item a +## @var{n} by @var{n} matrix; +## @item b +## Matrix: @var{n} by @var{n}, @var{n} by @var{m}, or @var{p} by @var{n}. +## @end table +## +## @strong{Output} +## @table @var +## @item x +## matrix satisfying appropriate discrete time Lyapunov equation. +## @end table +## +## Options: +## @itemize @bullet +## @item @var{b} is square: solve +## @iftex +## @tex +## $$ axa^T - x + b = 0 $$ +## @end tex +## @end iftex +## @ifinfo +## @code{a x a' - x + b = 0} +## @end ifinfo +## @item @var{b} is not square: @var{x} satisfies either +## @iftex +## @tex +## $$ axa^T - x + bb^T = 0 $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## a x a' - x + b b' = 0 +## @end example +## @end ifinfo +## @noindent +## or +## @iftex +## @tex +## $$ a^Txa - x + b^Tb = 0, $$ +## @end tex +## @end iftex +## @ifinfo +## @example +## a' x a - x + b' b = 0, +## @end example +## @end ifinfo +## @noindent +## whichever is appropriate. +## @end itemize +## +## @strong{Method} +## Uses Schur decomposition method as in Kitagawa, +## @cite{An Algorithm for Solving the Matrix Equation @math{X = F X F' + S}}, +## International Journal of Control, Volume 25, Number 5, pages 745--753 +## (1977). +## +## Column-by-column solution method as suggested in +## Hammarling, @cite{Numerical Solution of the Stable, Non-Negative +## Definite Lyapunov Equation}, @acronym{IMA} Journal of Numerical Analysis, Volume +## 2, pages 303--323 (1982). +## @end deftypefn + +## Author: A. S. Hodel <a.s...@en...> +## Created: August 1993 + +function x = dlyap (a, b) + + if (nargin != 2) + print_usage (); + endif + + if ((n = issquare (a)) == 0) + warning ("dlyap: a must be square"); + endif + + if ((m = issquare (b)) == 0) + [n1, m] = size (b); + if (n1 == n) + b = b*b'; + m = n1; + else + b = b'*b; + a = a'; + endif + endif + + if (n != m) + warning ("dlyap: a,b not conformably dimensioned"); + endif + + ## Solve the equation column by column. + + [u, s] = schur (a); + b = u'*b*u; + + j = n; + while (j > 0) + j1 = j; + + ## Check for Schur block. + + if (j == 1) + blksiz = 1; + elseif (s (j, j-1) != 0) + blksiz = 2; + j = j - 1; + else + blksiz = 1; + endif + + Ajj = kron (s(j:j1,j:j1), s) - eye (blksiz*n); + + rhs = reshape (b (:,j:j1), blksiz*n, 1); + + if (j1 < n) + rhs2 = s*(x(:,(j1+1):n) * s(j:j1,(j1+1):n)'); + rhs = rhs + reshape (rhs2, blksiz*n, 1); + endif + + v = - Ajj\rhs; + x(:,j) = v (1:n); + + if (blksiz == 2) + x (:, j1) = v ((n+1):blksiz*n); + endif + + j = j - 1; + + endwhile + + ## Back-transform to original coordinates. + + x = u*x*u'; + +endfunction Copied: trunk/octave-forge/extra/control-oo/inst/estim.m (from rev 6575, trunk/octave-forge/extra/control-oo/inst/control/estim.m) =================================================================== --- trunk/octave-forge/extra/control-oo/inst/estim.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/estim.m 2009-12-02 13:39:43 UTC (rev 6576) @@ -0,0 +1,57 @@ +## 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 -*- + +## Author: Lukas Reichlin <luk...@gm...> +## Created: November 2009 +## Version: 0.1 + +function est = estim (sys, l, sensors = [], known = []) + + if (nargin < 2 || nargin > 4) + print_usage (); + endif + + if (! isa (sys, "lti")) + error ("estim: first argument must be a LTI system"); + endif + + [a... [truncated message content] |
From: <par...@us...> - 2009-12-08 17:04:17
|
Revision: 6619 http://octave.svn.sourceforge.net/octave/?rev=6619&view=rev Author: paramaniac Date: 2009-12-08 17:04:04 +0000 (Tue, 08 Dec 2009) Log Message: ----------- control-oo: sync doc with code Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/augw.m trunk/octave-forge/extra/control-oo/inst/mixsyn.m Modified: trunk/octave-forge/extra/control-oo/inst/augw.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/augw.m 2009-12-08 10:54:03 UTC (rev 6618) +++ trunk/octave-forge/extra/control-oo/inst/augw.m 2009-12-08 17:04:04 UTC (rev 6619) @@ -17,7 +17,7 @@ ## -*- texinfo -*- ## @deftypefn{Function File} {@var{P} =} augw (@var{G}, @var{W1}, @var{W2}, @var{W3}) -## Extend plant for stacked S/T/KS problem. +## Extend plant for stacked S/KS/T problem. ## @example ## @group ## @@ -39,10 +39,10 @@ ## +----------------------------------------+ ## ## +--------+ -## | |-----> z1 (p1x1) -## r (px1) ----->| P(s) |-----> z2 (p2x1) -## | |-----> z3 (p3x1) -## u (mx1) ----->| |-----> e (px1) +## | |-----> z1 (p1x1) z1 = W1 e +## r (px1) ----->| P(s) |-----> z2 (p2x1) z2 = W2 u +## | |-----> z3 (p3x1) z3 = W3 y +## u (mx1) ----->| |-----> e (px1) e = r - y ## +--------+ ## ## +--------+ Modified: trunk/octave-forge/extra/control-oo/inst/mixsyn.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/mixsyn.m 2009-12-08 10:54:03 UTC (rev 6618) +++ trunk/octave-forge/extra/control-oo/inst/mixsyn.m 2009-12-08 17:04:04 UTC (rev 6619) @@ -17,9 +17,9 @@ ## -*- texinfo -*- ## @deftypefn{Function File} {[@var{K}, @var{N}, @var{gamma}] =} mixsyn (@var{G}, @var{W1}, @var{W2}, @var{W3}) -## Solve stacked S/T/KS H-inf problem, i.e. bound the largest singular values -## of S (for performance), T (for robustness and to avoid sensitivity to noise) -## and K S (to penalize large inputs). +## Solve stacked S/KS/T H-inf problem, i.e. bound the largest singular values +## of S (for performance), K S (to penalize large inputs) and +## T (for robustness and to avoid sensitivity to noise). ## @example ## @group ## @@ -39,10 +39,10 @@ ## +----------------------------------------+ ## ## +--------+ -## | |-----> z1 (p1x1) -## r (px1) ----->| P(s) |-----> z2 (p2x1) -## | |-----> z3 (p3x1) -## u (mx1) ----->| |-----> e (px1) +## | |-----> z1 (p1x1) z1 = W1 e +## r (px1) ----->| P(s) |-----> z2 (p2x1) z2 = W2 u +## | |-----> z3 (p3x1) z3 = W3 y +## u (mx1) ----->| |-----> e (px1) e = r - y ## +--------+ ## ## +--------+ @@ -56,8 +56,14 @@ ## +--------+ ## ## +--------+ -## r ------| N(s) |-----> z +## r ----->| N(s) |-----> z ## +--------+ +## +## Extended Plant: P = augw (G, W1, W2, W3) +## Entire System: N = lft (P, K) +## Open Loop: L = G * K +## Closed Loop: T = feedback (L) +## ## Reference: ## Skogestad, S. and Postlethwaite I. ## Multivariable Feedback Control: Analysis and Design This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-12-09 08:13:13
|
Revision: 6620 http://octave.svn.sourceforge.net/octave/?rev=6620&view=rev Author: paramaniac Date: 2009-12-09 08:13:06 +0000 (Wed, 09 Dec 2009) Log Message: ----------- control-oo: improve example optiPID Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/optiPID.m trunk/octave-forge/extra/control-oo/inst/optiPIDfunction.m Modified: trunk/octave-forge/extra/control-oo/inst/optiPID.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/optiPID.m 2009-12-08 17:04:04 UTC (rev 6619) +++ trunk/octave-forge/extra/control-oo/inst/optiPID.m 2009-12-09 08:13:06 UTC (rev 6620) @@ -1,118 +1,105 @@ -% Numerical Optimization of an A/H PID Controller -% Written by Lukas Reichlin July 2009 -% Required OCTAVE Packages: control, miscellaneous, optim -% Required MATLAB Toolboxes: Control, Optimization - -% Tabula Rasa -clear all; -close all; -clc; - -% Global Variables -global P w t dt mu_1 mu_2 mu_3 - -% Plant -numP = [1]; -denP = conv ([1, 1, 1], [1, 4, 6, 4, 1]); -P = tf (numP, denP); - -% Relative Weighting Factors: PLAY AROUND WITH THESE! -mu_1 = 1; % Minimize ITAE Criterion -mu_2 = 1; % Minimize Max Overshoot -mu_3 = 1; % Maximize Critical Distance - -% Simulation Settings: PLANT-DEPENDENT! -t_sim = 30; % Simulation Time [s] -dt = 0.05; % Sampling Time [s] -t = 0 : dt : t_sim; % Time Vector [s] -w = logspace (-1, 1, 100); % Frequency Range [rad/s] - -% A/H PID Controller: Ms = 2.0 (mu_min = 0.5) -[gamma, phi, w_gamma, w_phi] = margin (P); - -ku = gamma; -Tu = 2*pi / w_gamma; -kappa = inv (dcgain (P)); - -disp ('optiPID: Astrom/Hagglund PID controller parameters:'); - -kp_AH = ku * 0.72 * exp ( -1.60 * kappa + 1.20 * kappa^2 ) -Ti_AH = Tu * 0.59 * exp ( -1.30 * kappa + 0.38 * kappa^2 ) -Td_AH = Tu * 0.15 * exp ( -1.40 * kappa + 0.56 * kappa^2 ) -tau_AH = Td_AH / 10 - -numC_AH = kp_AH * [Ti_AH * Td_AH, Ti_AH, 1]; -denC_AH = conv ([Ti_AH, 0], [tau_AH^2, 2 * tau_AH, 1]); -C_AH = tf (numC_AH, denC_AH); - -% Initial Values -kp_0 = kp_AH; -Ti_0 = Ti_AH; -Td_0 = Td_AH; - -C_par_0 = [kp_0; Ti_0; Td_0]; - -% Optimization -if (exist ('fminsearch')) - warning ('optiPID: optimization starts, please be patient ...'); -else - error ('optiPID: please install optim package to proceed'); -end - -C_par_opt = fminsearch (@optiPIDfunction, C_par_0); - -% Optimized Controller -disp ('optiPID: optimized PID controller parameters:'); - -kp_opt = C_par_opt(1) -Ti_opt = C_par_opt(2) -Td_opt = C_par_opt(3) -tau_opt = Td_opt / 10 - -numC_opt = kp_opt * [Ti_opt * Td_opt, Ti_opt, 1]; -denC_opt = conv ([Ti_opt, 0], [tau_opt^2, 2 * tau_opt, 1]); -C_opt = tf (numC_opt, denC_opt); - -% Open Loop -L_AH = P * C_AH; -L_opt = P * C_opt; - -% Closed Loop -T_AH = feedback (L_AH, 1); -T_opt = feedback (L_opt, 1); - -% A Posteriori Stability Check -disp ('optiPID: for stability, all eigenvalues should have negative real parts:'); - -eigw_AH = eig (T_AH) -eigw_opt = eig (T_opt) - -% Stability Margins -disp ('optiPID: gain margin gamma [-] and phase margin phi [deg]:'); - -[gamma_AH, phi_AH] = margin (L_AH) -[gamma_opt, phi_opt] = margin (L_opt) - -% Plot Step Response -[y_AH, t_AH] = step (T_AH, t); -[y_opt, t_opt] = step (T_opt, t); - -figure (1) -plot (t_AH, y_AH, 'b', t_opt, y_opt, 'r') -grid ('on') -title ('Step Response') -xlabel ('Time [s]') -ylabel ('Output [-]') -legend ('A/H', 'Optimized', 'Location', 'SouthEast') - -% Plot Nyquist Diagram -[re_AH, im_AH] = nyquist (L_AH, w); -[re_opt, im_opt] = nyquist (L_opt, w); - -figure (2) -plot (re_AH(:), im_AH(:), 'b', re_opt(:), im_opt(:), 'r') -grid ('on') -title ('Nyquist Diagram') -xlabel ('Re\{L(jw)\}') -ylabel ('Im\{L(jw)\}') -legend ('A/H', 'Optimized', 'Location', 'SouthEast') +% Numerical Optimization of an A/H PID Controller +% Written by Lukas Reichlin July 2009 +% Required OCTAVE Packages: control, miscellaneous, optim +% Required MATLAB Toolboxes: Control, Optimization + +% Tabula Rasa +clear all; +close all; +clc; + +% Global Variables +global P t dt mu_1 mu_2 mu_3 + +% Plant +numP = [1]; +denP = conv ([1, 1, 1], [1, 4, 6, 4, 1]); +P = tf (numP, denP); + +% Relative Weighting Factors: PLAY AROUND WITH THESE! +mu_1 = 1; % Minimize ITAE Criterion +mu_2 = 1; % Minimize Max Overshoot +mu_3 = 1; % Maximize Critical Distance + +% Simulation Settings: PLANT-DEPENDENT! +t_sim = 30; % Simulation Time [s] +dt = 0.05; % Sampling Time [s] +t = 0 : dt : t_sim; % Time Vector [s] + +% A/H PID Controller: Ms = 2.0 (mu_min = 0.5) +[gamma, phi, w_gamma, w_phi] = margin (P); + +ku = gamma; +Tu = 2*pi / w_gamma; +kappa = inv (dcgain (P)); + +disp ('optiPID: Astrom/Hagglund PID controller parameters:'); + +kp_AH = ku * 0.72 * exp ( -1.60 * kappa + 1.20 * kappa^2 ) +Ti_AH = Tu * 0.59 * exp ( -1.30 * kappa + 0.38 * kappa^2 ) +Td_AH = Tu * 0.15 * exp ( -1.40 * kappa + 0.56 * kappa^2 ) +tau_AH = Td_AH / 10 + +numC_AH = kp_AH * [Ti_AH * Td_AH, Ti_AH, 1]; +denC_AH = conv ([Ti_AH, 0], [tau_AH^2, 2 * tau_AH, 1]); +C_AH = tf (numC_AH, denC_AH); + +% Initial Values +kp_0 = kp_AH; +Ti_0 = Ti_AH; +Td_0 = Td_AH; + +C_par_0 = [kp_0; Ti_0; Td_0]; + +% Optimization +if (exist ('fminsearch')) + warning ('optiPID: optimization starts, please be patient ...'); +else + error ('optiPID: please install optim package to proceed'); +end + +C_par_opt = fminsearch (@optiPIDfunction, C_par_0); + +% Optimized Controller +disp ('optiPID: optimized PID controller parameters:'); + +kp_opt = C_par_opt(1) +Ti_opt = C_par_opt(2) +Td_opt = C_par_opt(3) +tau_opt = Td_opt / 10 + +numC_opt = kp_opt * [Ti_opt * Td_opt, Ti_opt, 1]; +denC_opt = conv ([Ti_opt, 0], [tau_opt^2, 2 * tau_opt, 1]); +C_opt = tf (numC_opt, denC_opt); + +% Open Loop +L_AH = P * C_AH; +L_opt = P * C_opt; + +% Closed Loop +T_AH = feedback (L_AH, 1); +T_opt = feedback (L_opt, 1); + +% A Posteriori Stability Check +disp ('optiPID: for stability, all eigenvalues should have negative real parts:'); + +eigw_AH = eig (T_AH) +eigw_opt = eig (T_opt) + +% Stability Margins +disp ('optiPID: gain margin gamma [-] and phase margin phi [deg]:'); + +[gamma_AH, phi_AH] = margin (L_AH) +[gamma_opt, phi_opt] = margin (L_opt) + +% Plot Step Response +[y_AH, t_AH] = step (T_AH, t); +[y_opt, t_opt] = step (T_opt, t); + +figure (1) +plot (t_AH, y_AH, 'b', t_opt, y_opt, 'r') +grid ('on') +title ('Step Response') +xlabel ('Time [s]') +ylabel ('Output [-]') +legend ('A/H', 'Optimized', 'Location', 'SouthEast') Modified: trunk/octave-forge/extra/control-oo/inst/optiPIDfunction.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/optiPIDfunction.m 2009-12-08 17:04:04 UTC (rev 6619) +++ trunk/octave-forge/extra/control-oo/inst/optiPIDfunction.m 2009-12-09 08:13:06 UTC (rev 6620) @@ -1,54 +1,54 @@ -function J = optiPIDfunction (C_par) - -% Objective Function -% written by Lukas Reichlin -% See "Analysis and Synthesis of SISO Control Systems" -% by Lino Guzzella for further details - -global P w t dt mu_1 mu_2 mu_3 - -% Function Argument -> Controller Parameters -kp = C_par(1); -Ti = C_par(2); -Td = C_par(3); -tau = Td / 10; - -% PID Controller with Roll-Off -numC = kp * [Ti * Td, Ti, 1]; -denC = conv ([Ti, 0], [tau^2, 2 * tau, 1]); -C = tf (numC, denC); - -% Open Loop -L = P * C; - -% Sum Block: e = r - y -SUM = ss ([1, -1]); % Matlab converts to SS (and back) for MIMO TF connections - -% Group Sum Block and Open Loop -SUML = append (SUM, L); - -% Build System Interconnections -CM = [3, 1; % Controller Input with Sum Block Output - 2, 2]; % Sum Block Negative Input with Plant Output - -inputs = [1]; % Input 1: reference r(t) -outputs = [1, 2]; % Output 1: error e(t), Output 2: output y(t) - -SUML = connect (SUML, CM, inputs, outputs); - -% Simulation -[y, t_y] = step (SUML, t); - -% ITAE Criterion - use for-loop instead of lsim for performance reasons -itae = 0; - -for k = 1 : length (y) - itae = itae + t_y(k) * abs (y(k, 1)) * dt; -end - -% Return Difference -Q = 1 + L; -mag = bode (Q, w); - -% Objective Function -J = mu_1 * itae + mu_2 * (max (y(:, 2)) - 1) + mu_3 * (1 - min (mag)); +function J = optiPIDfunction (C_par) + +% Objective Function +% written by Lukas Reichlin +% See "Analysis and Synthesis of SISO Control Systems" +% by Lino Guzzella for further details + +global P t dt mu_1 mu_2 mu_3 + +% Function Argument -> Controller Parameters +kp = C_par(1); +Ti = C_par(2); +Td = C_par(3); +tau = Td / 10; + +% PID Controller with Roll-Off +numC = kp * [Ti * Td, Ti, 1]; +denC = conv ([Ti, 0], [tau^2, 2 * tau, 1]); +C = tf (numC, denC); + +% Open Loop +L = P * C; + +% Sum Block: e = r - y +SUM = ss ([1, -1]); % Matlab converts to SS (and back) for MIMO TF connections + +% Group Sum Block and Open Loop +SUML = append (SUM, L); + +% Build System Interconnections +CM = [3, 1; % Controller Input with Sum Block Output + 2, 2]; % Sum Block Negative Input with Plant Output + +inputs = [1]; % Input 1: reference r(t) +outputs = [1, 2]; % Output 1: error e(t), Output 2: output y(t) + +SUML = connect (SUML, CM, inputs, outputs); + +% Simulation +[y, t_y] = step (SUML, t); + +% ITAE Criterion - use for-loop instead of lsim for performance reasons +itae = 0; + +for k = 1 : length (y) + itae = itae + t_y(k) * abs (y(k, 1)) * dt; +end + +% Critical Distance mu = 1/Ms +S = inv (1 + L); +mu = 1 / norm (S, inf); + +% Objective Function +J = mu_1 * itae + mu_2 * (max (y(:, 2)) - 1) + mu_3 * (1 - mu); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-12-11 10:48:29
|
Revision: 6632 http://octave.svn.sourceforge.net/octave/?rev=6632&view=rev Author: paramaniac Date: 2009-12-11 10:48:15 +0000 (Fri, 11 Dec 2009) Log Message: ----------- control-oo: display system names of tf and ss objects Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/@ss/display.m trunk/octave-forge/extra/control-oo/inst/@tf/display.m Modified: trunk/octave-forge/extra/control-oo/inst/@ss/display.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@ss/display.m 2009-12-10 14:15:33 UTC (rev 6631) +++ trunk/octave-forge/extra/control-oo/inst/@ss/display.m 2009-12-11 10:48:15 UTC (rev 6632) @@ -47,16 +47,15 @@ stname = __markemptynames__ (stname); endif - ## fprintf ("SS object %s:\n", inputname(1)); disp (""); if (! isempty (sys.a)) - dispmat (sys.a, "a", stname, stname); - dispmat (sys.b, "b", stname, inname); - dispmat (sys.c, "c", outname, stname); + dispmat (sys.a, [inputname(1), ".a"], stname, stname); + dispmat (sys.b, [inputname(1), ".b"], stname, inname); + dispmat (sys.c, [inputname(1), ".c"], outname, stname); endif - dispmat (sys.d, "d", outname, inname); + dispmat (sys.d, [inputname(1), ".d"], outname, inname); display (sys.lti); # display sampling time Modified: trunk/octave-forge/extra/control-oo/inst/@tf/display.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/@tf/display.m 2009-12-10 14:15:33 UTC (rev 6631) +++ trunk/octave-forge/extra/control-oo/inst/@tf/display.m 2009-12-11 10:48:15 UTC (rev 6632) @@ -40,11 +40,10 @@ outname = __markemptynames__ (outname); endif - ## fprintf ("TF object %s:\n", inputname(1)); disp (""); for m = 1 : nu - disp (["Transfer function from input \"", inname{m}, "\" to output ..."]); + disp (["Transfer function \"", inputname(1), "\" from input \"", inname{m}, "\" to output ..."]); for p = 1 : ny dispfrac (sys.num{p, m}, sys.den{p, m}, sys.tfvar, outname{p}); endfor This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-12-11 17:30:03
|
Revision: 6634 http://octave.svn.sourceforge.net/octave/?rev=6634&view=rev Author: paramaniac Date: 2009-12-11 17:29:52 +0000 (Fri, 11 Dec 2009) Log Message: ----------- control-oo: add/update texinfo strings Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/care.m trunk/octave-forge/extra/control-oo/inst/dare.m trunk/octave-forge/extra/control-oo/inst/dlqr.m trunk/octave-forge/extra/control-oo/inst/lqr.m Modified: trunk/octave-forge/extra/control-oo/inst/care.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/care.m 2009-12-11 14:24:53 UTC (rev 6633) +++ trunk/octave-forge/extra/control-oo/inst/care.m 2009-12-11 17:29:52 UTC (rev 6634) @@ -16,6 +16,12 @@ ## along with this program. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- +## @deftypefn {Function File} {[@var{x}, @var{l}, @var{g}] =} care (@var{a}, @var{b}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{x}, @var{l}, @var{g}] =} care (@var{a}, @var{b}, @var{q}, @var{r}, @var{s}) +## Return unique stabilizing solution x of the continuous-time +## riccati equation as well as the closed-loop poles l and the +## corresponding gain matrix g. +## @end deftypefn ## Author: Lukas Reichlin <luk...@gm...> ## Created: November 2009 Modified: trunk/octave-forge/extra/control-oo/inst/dare.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/dare.m 2009-12-11 14:24:53 UTC (rev 6633) +++ trunk/octave-forge/extra/control-oo/inst/dare.m 2009-12-11 17:29:52 UTC (rev 6634) @@ -17,8 +17,9 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{x} =} dare (@var{a}, @var{b}, @var{q}, @var{r}, @var{opt}) -## +## @deftypefn {Function File} {[@var{x}, @var{l}, @var{g}] =} dare (@var{a}, @var{b}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{x}, @var{l}, @var{g}] =} dare (@var{a}, @var{b}, @var{q}, @var{r}, @var{s}) +## @deftypefnx {Function File} {[@var{x}, @var{l}, @var{g}] =} dare (@var{a}, @var{b}, @var{q}, @var{r}, @var{s}, @var{opt}) ## Return the solution, @var{x} of the discrete-time algebraic Riccati ## equation ## @iftex Modified: trunk/octave-forge/extra/control-oo/inst/dlqr.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/dlqr.m 2009-12-11 14:24:53 UTC (rev 6633) +++ trunk/octave-forge/extra/control-oo/inst/dlqr.m 2009-12-11 17:29:52 UTC (rev 6634) @@ -16,6 +16,13 @@ ## along with this program. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- +## @deftypefn {Function File} {[@var{g}, @var{x}, @var{l}] =} dlqr (@var{sys}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} dlqr (@var{sys}, @var{q}, @var{r}, @var{s}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} dlqr (@var{a}, @var{b}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} dlqr (@var{a}, @var{b}, @var{q}, @var{r}, @var{s}) +## Return linear-quadratic state-feedback gain matrix g for a LTI system as well as +## the solution x of the associated riccati equation and the closed-loop poles l. +## @end deftypefn ## Author: Lukas Reichlin <luk...@gm...> ## Created: November 2009 Modified: trunk/octave-forge/extra/control-oo/inst/lqr.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/lqr.m 2009-12-11 14:24:53 UTC (rev 6633) +++ trunk/octave-forge/extra/control-oo/inst/lqr.m 2009-12-11 17:29:52 UTC (rev 6634) @@ -16,6 +16,13 @@ ## along with this program. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- +## @deftypefn {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{sys}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{sys}, @var{q}, @var{r}, @var{s}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}, @var{s}) +## Return linear-quadratic state-feedback gain matrix g for a LTI system as well as +## the solution x of the associated riccati equation and the closed-loop poles l. +## @end deftypefn ## Author: Lukas Reichlin <luk...@gm...> ## Created: November 2009 This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-12-12 14:05:02
|
Revision: 6639 http://octave.svn.sourceforge.net/octave/?rev=6639&view=rev Author: paramaniac Date: 2009-12-12 14:04:48 +0000 (Sat, 12 Dec 2009) Log Message: ----------- control-oo: add missing texinfo strings Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/estim.m trunk/octave-forge/extra/control-oo/inst/kalman.m Modified: trunk/octave-forge/extra/control-oo/inst/estim.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/estim.m 2009-12-11 22:16:38 UTC (rev 6638) +++ trunk/octave-forge/extra/control-oo/inst/estim.m 2009-12-12 14:04:48 UTC (rev 6639) @@ -16,6 +16,10 @@ ## along with this program. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- +## @deftypefn {Function File} {@var{est} =} estim (@var{sys}, @var{l}) +## @deftypefnx {Function File} {@var{est} =} estim (@var{sys}, @var{l}, @var{sensors}, @var{known}) +## Return state estimator for a given estimator gain +## @end deftypefn ## Author: Lukas Reichlin <luk...@gm...> ## Created: November 2009 Modified: trunk/octave-forge/extra/control-oo/inst/kalman.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/kalman.m 2009-12-11 22:16:38 UTC (rev 6638) +++ trunk/octave-forge/extra/control-oo/inst/kalman.m 2009-12-12 14:04:48 UTC (rev 6639) @@ -16,6 +16,10 @@ ## along with this program. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- +## @deftypefn {Function File} {[@var{est}, @var{g}, @var{x}] =} kalman (@var{sys}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{est}, @var{g}, @var{x}] =} kalman (@var{sys}, @var{q}, @var{r}, @var{s}) +## Design kalman estimator for LTI systems. +## @end deftypefn ## Author: Lukas Reichlin <luk...@gm...> ## Created: November 2009 This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-12-14 17:21:29
|
Revision: 6644 http://octave.svn.sourceforge.net/octave/?rev=6644&view=rev Author: paramaniac Date: 2009-12-14 17:21:21 +0000 (Mon, 14 Dec 2009) Log Message: ----------- control-oo: fix bug Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/h2syn.m trunk/octave-forge/extra/control-oo/inst/hinfsyn.m Modified: trunk/octave-forge/extra/control-oo/inst/h2syn.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/h2syn.m 2009-12-13 20:51:27 UTC (rev 6643) +++ trunk/octave-forge/extra/control-oo/inst/h2syn.m 2009-12-14 17:21:21 UTC (rev 6644) @@ -82,7 +82,7 @@ N = lft (P, K); varargout{1} = N; if (nargout > 2) - varargout{2} = norm (N); + varargout{2} = norm (N, 2); endif endif Modified: trunk/octave-forge/extra/control-oo/inst/hinfsyn.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/hinfsyn.m 2009-12-13 20:51:27 UTC (rev 6643) +++ trunk/octave-forge/extra/control-oo/inst/hinfsyn.m 2009-12-14 17:21:21 UTC (rev 6644) @@ -81,7 +81,7 @@ N = lft (P, K); varargout{1} = N; if (nargout > 2) - varargout{2} = norm (N); + varargout{2} = norm (N, inf); endif endif This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2009-12-18 21:11:05
|
Revision: 6654 http://octave.svn.sourceforge.net/octave/?rev=6654&view=rev Author: paramaniac Date: 2009-12-18 21:10:55 +0000 (Fri, 18 Dec 2009) Log Message: ----------- control-oo: simplify optiPID Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/optiPID.m trunk/octave-forge/extra/control-oo/inst/optiPIDfunction.m Modified: trunk/octave-forge/extra/control-oo/inst/optiPID.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/optiPID.m 2009-12-17 12:14:46 UTC (rev 6653) +++ trunk/octave-forge/extra/control-oo/inst/optiPID.m 2009-12-18 21:10:55 UTC (rev 6654) @@ -19,14 +19,14 @@ % Relative Weighting Factors: PLAY AROUND WITH THESE! mu_1 = 1; % Minimize ITAE Criterion mu_2 = 1; % Minimize Max Overshoot -mu_3 = 1; % Maximize Critical Distance +mu_3 = 1; % Minimize Sensitivity Ms % Simulation Settings: PLANT-DEPENDENT! -t_sim = 30; % Simulation Time [s] -dt = 0.05; % Sampling Time [s] -t = 0 : dt : t_sim; % Time Vector [s] +t_sim = 30; % Simulation Time [s] +dt = 0.05; % Sampling Time [s] +t = 0 : dt : t_sim; % Time Vector [s] -% A/H PID Controller: Ms = 2.0 (mu_min = 0.5) +% A/H PID Controller: Ms = 2.0 [gamma, phi, w_gamma, w_phi] = margin (P); ku = gamma; @@ -34,7 +34,6 @@ kappa = inv (dcgain (P)); disp ('optiPID: Astrom/Hagglund PID controller parameters:'); - kp_AH = ku * 0.72 * exp ( -1.60 * kappa + 1.20 * kappa^2 ) Ti_AH = Tu * 0.59 * exp ( -1.30 * kappa + 0.38 * kappa^2 ) Td_AH = Tu * 0.15 * exp ( -1.40 * kappa + 0.56 * kappa^2 ) @@ -45,12 +44,8 @@ C_AH = tf (numC_AH, denC_AH); % Initial Values -kp_0 = kp_AH; -Ti_0 = Ti_AH; -Td_0 = Td_AH; +C_par_0 = [kp_AH; Ti_AH; Td_AH]; -C_par_0 = [kp_0; Ti_0; Td_0]; - % Optimization if (exist ('fminsearch')) warning ('optiPID: optimization starts, please be patient ...'); @@ -62,7 +57,6 @@ % Optimized Controller disp ('optiPID: optimized PID controller parameters:'); - kp_opt = C_par_opt(1) Ti_opt = C_par_opt(2) Td_opt = C_par_opt(3) @@ -81,14 +75,12 @@ T_opt = feedback (L_opt, 1); % A Posteriori Stability Check -disp ('optiPID: for stability, all poles should have negative real parts:'); +disp ('optiPID: closed-loop stability check:'); +st_AH = isstable (T_AH) +st_opt = isstable (T_opt) -p_AH = pole (T_AH) -p_opt = pole (T_opt) - % Stability Margins disp ('optiPID: gain margin gamma [-] and phase margin phi [deg]:'); - [gamma_AH, phi_AH] = margin (L_AH) [gamma_opt, phi_opt] = margin (L_opt) Modified: trunk/octave-forge/extra/control-oo/inst/optiPIDfunction.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/optiPIDfunction.m 2009-12-17 12:14:46 UTC (rev 6653) +++ trunk/octave-forge/extra/control-oo/inst/optiPIDfunction.m 2009-12-18 21:10:55 UTC (rev 6654) @@ -46,9 +46,9 @@ itae = itae + t_y(k) * abs (y(k, 1)) * dt; end -% Critical Distance mu = 1/Ms +% Sensitivity S = inv (1 + L); -mu = 1 / norm (S, inf); +Ms = norm (S, inf); % Objective Function -J = mu_1 * itae + mu_2 * (max (y(:, 2)) - 1) + mu_3 * (1 - mu); +J = mu_1 * itae + mu_2 * (max (y(:, 2)) - 1) + mu_3 * Ms; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-01-10 18:15:05
|
Revision: 6732 http://octave.svn.sourceforge.net/octave/?rev=6732&view=rev Author: paramaniac Date: 2010-01-10 18:14:58 +0000 (Sun, 10 Jan 2010) Log Message: ----------- control-oo: fix silly mistake Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/dlyap.m trunk/octave-forge/extra/control-oo/inst/lyap.m Modified: trunk/octave-forge/extra/control-oo/inst/dlyap.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/dlyap.m 2010-01-10 17:51:13 UTC (rev 6731) +++ trunk/octave-forge/extra/control-oo/inst/dlyap.m 2010-01-10 18:14:58 UTC (rev 6732) @@ -55,7 +55,7 @@ [x, scale] = slsb03md (a, -b, true); # AXA' - X = -B - x *= scale; + x /= scale; # 0 < scale <= 1 elseif (nargin == 3) # Sylvester equation Modified: trunk/octave-forge/extra/control-oo/inst/lyap.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/lyap.m 2010-01-10 17:51:13 UTC (rev 6731) +++ trunk/octave-forge/extra/control-oo/inst/lyap.m 2010-01-10 18:14:58 UTC (rev 6732) @@ -55,7 +55,7 @@ [x, scale] = slsb03md (a, -b, false); # AX + XA' = -B - x *= scale; + x /= scale; # 0 < scale <= 1 elseif (nargin == 3) # Sylvester equation This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-01-11 11:28:33
|
Revision: 6743 http://octave.svn.sourceforge.net/octave/?rev=6743&view=rev Author: paramaniac Date: 2010-01-11 11:28:08 +0000 (Mon, 11 Jan 2010) Log Message: ----------- control-oo: prepare update Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/dlyap.m trunk/octave-forge/extra/control-oo/inst/lyap.m Modified: trunk/octave-forge/extra/control-oo/inst/dlyap.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/dlyap.m 2010-01-11 10:17:20 UTC (rev 6742) +++ trunk/octave-forge/extra/control-oo/inst/dlyap.m 2010-01-11 11:28:08 UTC (rev 6743) @@ -18,14 +18,17 @@ ## -*- texinfo -*- ## @deftypefn{Function File} {@var{x} =} dlyap (@var{a}, @var{b}) ## @deftypefnx{Function File} {@var{x} =} dlyap (@var{a}, @var{b}, @var{c}) +## @deftypefnx{Function File} {@var{x} =} dlyap (@var{a}, @var{b}, @var{[]}, @var{e}) ## Solve discrete-time Lyapunov or Sylvester equations. -## Uses SLICOT SB03MD and SB04QD by courtesy of NICONET e.V. +## Uses SLICOT SB03MD, SB04QD and SG03AD by courtesy of NICONET e.V. ## <http://www.slicot.org> ## @example ## @group -## AXA' - X + B = 0 (Lyapunov Equation) +## AXA' - X + B = 0 (Lyapunov Equation) ## -## AXB' - X + C = 0 (Sylvester Equation) +## AXB' - X + C = 0 (Sylvester Equation) +## +## AXA' - EXE' + B = 0 (Generalized Lyapunov Equation) ## @end group ## @end example ## @end deftypefn @@ -34,53 +37,58 @@ ## Created: January 2010 ## Version: 0.1 -function x = dlyap (a, b, c) +function x = dlyap (a, b, c, e) - if (nargin == 2) # Lyapunov equation + switch (nargin) + case 2 # Lyapunov equation - na = issquare (a); - nb = issquare (b); + na = issquare (a); + nb = issquare (b); - if (! na) - error ("lyap: a must be square"); - endif + if (! na) + error ("lyap: a must be square"); + endif - if (! nb) - error ("lyap: b must be square") - endif + if (! nb) + error ("lyap: b must be square") + endif - if (na != nb) - error ("lyap: a and b must be of identical size"); - endif + if (na != nb) + error ("lyap: a and b must be of identical size"); + endif + + [x, scale] = slsb03md (a, -b, true); # AXA' - X = -B - [x, scale] = slsb03md (a, -b, true); # AXA' - X = -B - - x /= scale; # 0 < scale <= 1 - - elseif (nargin == 3) # Sylvester equation - - n = issquare (a); - m = issquare (b); - [crows, ccols] = size (c); - - if (! n) - error ("dlyap: a must be square"); - endif - - if (! m) - error ("dlyap: b must be square"); - endif - - if (crows != n || ccols != m) - error ("dlyap: c must be a (%dx%d) matrix", n, m); - endif + x /= scale; # 0 < scale <= 1 - x = slsb04qd (-a, b, c); # AXB' - X = -C + case 3 # Sylvester equation + + n = issquare (a); + m = issquare (b); + [crows, ccols] = size (c); - else - print_usage (); - endif + if (! n) + error ("dlyap: a must be square"); + endif + if (! m) + error ("dlyap: b must be square"); + endif + + if (crows != n || ccols != m) + error ("dlyap: c must be a (%dx%d) matrix", n, m); + endif + + x = slsb04qd (-a, b, c); # AXB' - X = -C + + case 4 + error ("dlyap: case not implemented yet"); + + otherwise + print_usage (); + + endswitch + endfunction Modified: trunk/octave-forge/extra/control-oo/inst/lyap.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/lyap.m 2010-01-11 10:17:20 UTC (rev 6742) +++ trunk/octave-forge/extra/control-oo/inst/lyap.m 2010-01-11 11:28:08 UTC (rev 6743) @@ -18,14 +18,17 @@ ## -*- texinfo -*- ## @deftypefn{Function File} {@var{x} =} lyap (@var{a}, @var{b}) ## @deftypefnx{Function File} {@var{x} =} lyap (@var{a}, @var{b}, @var{c}) +## @deftypefnx{Function File} {@var{x} =} lyap (@var{a}, @var{b}, @var{[]}, @var{e}) ## Solve continuous-time Lyapunov or Sylvester equations. -## Uses SLICOT SB03MD and SB04MD by courtesy of NICONET e.V. +## Uses SLICOT SB03MD, SB04MD and SG03AD by courtesy of NICONET e.V. ## <http://www.slicot.org> ## @example ## @group -## AX + XA' + B = 0 (Lyapunov Equation) +## AX + XA' + B = 0 (Lyapunov Equation) ## -## AX + XB + C = 0 (Sylvester Equation) +## AX + XB + C = 0 (Sylvester Equation) +## +## AXE' + EXA' + B = 0 (Generalized Lyapunov Equation) ## @end group ## @end example ## @end deftypefn @@ -34,53 +37,58 @@ ## Created: January 2010 ## Version: 0.1 -function x = lyap (a, b, c) +function x = lyap (a, b, c, e) - if (nargin == 2) # Lyapunov equation + switch (nargin) + case 2 # Lyapunov equation - na = issquare (a); - nb = issquare (b); + na = issquare (a); + nb = issquare (b); - if (! na) - error ("lyap: a must be square"); - endif + if (! na) + error ("lyap: a must be square"); + endif - if (! nb) - error ("lyap: b must be square") - endif + if (! nb) + error ("lyap: b must be square") + endif - if (na != nb) - error ("lyap: a and b must be of identical size"); - endif - - [x, scale] = slsb03md (a, -b, false); # AX + XA' = -B + if (na != nb) + error ("lyap: a and b must be of identical size"); + endif + + [x, scale] = slsb03md (a, -b, false); # AX + XA' = -B + + x /= scale; # 0 < scale <= 1 - x /= scale; # 0 < scale <= 1 - - elseif (nargin == 3) # Sylvester equation + case 3 # Sylvester equation - n = issquare (a); - m = issquare (b); - [crows, ccols] = size (c); + n = issquare (a); + m = issquare (b); + [crows, ccols] = size (c); - if (! n) - error ("lyap: a must be square"); - endif - - if (! m) - error ("lyap: b must be square"); - endif - - if (crows != n || ccols != m) - error ("lyap: c must be a (%dx%d) matrix", n, m); - endif - - x = slsb04md (a, b, -c); # AX + XB = -C - - else - print_usage (); - endif + if (! n) + error ("lyap: a must be square"); + endif + if (! m) + error ("lyap: b must be square"); + endif + + if (crows != n || ccols != m) + error ("lyap: c must be a (%dx%d) matrix", n, m); + endif + + x = slsb04md (a, b, -c); # AX + XB = -C + + case 4 + error ("lyap: case not implemented yet"); + + otherwise + print_usage (); + + endswitch + endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2010-01-12 10:31:57
|
Revision: 6751 http://octave.svn.sourceforge.net/octave/?rev=6751&view=rev Author: paramaniac Date: 2010-01-12 10:31:50 +0000 (Tue, 12 Jan 2010) Log Message: ----------- control-oo: add TODOs Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/dlyap.m trunk/octave-forge/extra/control-oo/inst/hinfsyn.m Modified: trunk/octave-forge/extra/control-oo/inst/dlyap.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/dlyap.m 2010-01-11 22:13:44 UTC (rev 6750) +++ trunk/octave-forge/extra/control-oo/inst/dlyap.m 2010-01-12 10:31:50 UTC (rev 6751) @@ -161,4 +161,7 @@ %! 4.0000 7.0000 1.0000 %! 5.0000 3.0000 2.0000]; %! -%!assert (X, X_exp, 1e-4); \ No newline at end of file +%!assert (X, X_exp, 1e-4); + +## Generalized Lyapunov +TODO: add a test \ No newline at end of file Modified: trunk/octave-forge/extra/control-oo/inst/hinfsyn.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/hinfsyn.m 2010-01-11 22:13:44 UTC (rev 6750) +++ trunk/octave-forge/extra/control-oo/inst/hinfsyn.m 2010-01-12 10:31:50 UTC (rev 6751) @@ -27,6 +27,9 @@ ## Created: December 2009 ## Version: 0.1 +## TODO: find optimal instead of suboptimal controller +## TODO: improve compatibility for nargin >= 4 + function [K, varargout] = hinfsyn (P, nmeas, ncon, gmax = 1e15) ## check input arguments This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |