From: <par...@us...> - 2009-10-26 17:18:59
|
Revision: 6392 http://octave.svn.sourceforge.net/octave/?rev=6392&view=rev Author: paramaniac Date: 2009-10-26 17:18:48 +0000 (Mon, 26 Oct 2009) Log Message: ----------- control-oo: add draft code for time responses (step, impulse and initial) Modified Paths: -------------- trunk/octave-forge/extra/control-oo/INDEX Added Paths: ----------- trunk/octave-forge/extra/control-oo/inst/control/__timeresp__.m trunk/octave-forge/extra/control-oo/inst/control/impulse.m trunk/octave-forge/extra/control-oo/inst/control/initial.m trunk/octave-forge/extra/control-oo/inst/control/step.m Modified: trunk/octave-forge/extra/control-oo/INDEX =================================================================== --- trunk/octave-forge/extra/control-oo/INDEX 2009-10-26 11:23:39 UTC (rev 6391) +++ trunk/octave-forge/extra/control-oo/INDEX 2009-10-26 17:18:48 UTC (rev 6392) @@ -25,6 +25,9 @@ zero Time Domain Analysis gensig + impulse + initial + step Frequency Domain Analysis Model Simplification minreal Added: trunk/octave-forge/extra/control-oo/inst/control/__timeresp__.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/__timeresp__.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/control/__timeresp__.m 2009-10-26 17:18:48 UTC (rev 6392) @@ -0,0 +1,242 @@ +## 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 [y, t, x_arr] = __timeresp__ (sys, resptype, plotflag, tfinal, dt, x0) + + if (! isa (sys, "ss")) + sys = ss (sys); # sys must be proper + 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) # don't transform static gains + sys = c2d (sys, dt); + 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); + + ## preallocate memory + y = zeros (l_t, p); + x_arr = zeros (l_t, n); + + switch (resptype) + case "initial" + str = "Response to Initial Conditions"; + + ## initial conditions + x = x0(:); # make sure that x is a row 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"; + + ## initial conditions + x = zeros (n, 1); + + ## simulation + for k = 1 : l_t + y(k, :) = C * x + D; + x_arr(k, :) = x; + x = F * x + G; + endfor + + case "impulse" + str = "Impulse Response"; + + ## initial conditions + if (digital) + x = G / dt; + y(1, :) = D / dt; + x_arr(1, :) = x; + else + if (D'*D > 0) + warning ("impulse: system is not strictly proper"); + endif + + x = B; # B, not G! + y(1, :) = C * x; + x_arr(1, :) = x; + x = F * x; + endif + + ## simulation + for k = 2 : l_t + y (k, :) = C * x; + x_arr(k, :) = x; + x = F * x; + endfor + + if (digital) + y *= dt; + endif + + otherwise + error ("timeresp: invalid response type"); + + endswitch + + + if (plotflag) + outname = get (sys, "outname"); + + if (isempty (outname) || isequal ("", outname{:})) + outname = strseq ("y_", 1 : length (outname)); + else + outname = __markemptynames__ (outname); + endif + + if (digital) # discrete system + for k = 1 : p + subplot (p, 1, k) + stairs (t, y(:, k)) + grid on + if (k == 1) + title (str) + endif + ylabel (sprintf ("Amplitude %s", outname{k})) + endfor + xlabel ("Time [s]") + else # continuous system + for k = 1 : p + subplot (p, 1, k) + plot (t, y(:, k)) + grid on + if (k == 1) + title (str) + endif + ylabel (sprintf ("Amplitude %s", outname{k})) + 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 Added: trunk/octave-forge/extra/control-oo/inst/control/impulse.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/impulse.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/control/impulse.m 2009-10-26 17:18:48 UTC (rev 6392) @@ -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 -*- + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function [y_r, t_r, x_r] = impulse (sys, tfinal = [], dt = []) + + [y, t, x] = __timeresp__ (sys, "impulse", ! nargout, tfinal, dt); + + if (nargout) + y_r = y; + t_r = t; + x_r = x; + endif + +endfunction Added: trunk/octave-forge/extra/control-oo/inst/control/initial.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/initial.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/control/initial.m 2009-10-26 17:18:48 UTC (rev 6392) @@ -0,0 +1,98 @@ +## 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 [y_r, t_r, x_r] = initial (sys, x0, tfinal = [], dt = []) + + [y, t, x] = __timeresp__ (sys, "initial", ! nargout, tfinal, dt, x0); + + if (nargout) + y_r = y; + t_r = t; + x_r = x; + endif + +endfunction + + +%!shared initial_c, initial_c_exp, initial_d, initial_d_exp +%! +%! A = [ -2.8 2.0 -1.8 +%! -2.4 -2.0 0.8 +%! 1.1 1.7 -1.0 ]; +%! +%! B = [ -0.8 0.5 0 +%! 0 0.7 2.3 +%! -0.3 -0.1 0.5 ]; +%! +%! C = [ -0.1 0 -0.3 +%! 0.9 0.5 1.2 +%! 0.1 -0.1 1.9 ]; +%! +%! D = [ -0.5 0 0 +%! 0.1 0 0.3 +%! -0.8 0 0 ]; +%! +%! x_0 = [1, 2, 3]; +%! +%! sysc = ss (A, B, C, D); +%! +%! [yc, tc, xc] = initial (sysc, x_0, 0.2, 0.1); +%! initial_c = round (1e4 * [yc, tc, xc]) / 1e4; +%! +%! sysd = c2d (sysc, 2); +%! +%! [yd, td, xd] = initial (sysd, x_0, 4); +%! initial_d = round (1e4 * [yd, td, xd]) / 1e4; +%! +%! ## expected values computed by the "dark side" +%! +%! yc_exp = [ -1.0000 5.5000 5.6000 +%! -0.9872 5.0898 5.7671 +%! -0.9536 4.6931 5.7598 ]; +%! +%! tc_exp = [ 0.0000 +%! 0.1000 +%! 0.2000 ]; +%! +%! xc_exp = [ 1.0000 2.0000 3.0000 +%! 0.5937 1.6879 3.0929 +%! 0.2390 1.5187 3.0988 ]; +%! +%! initial_c_exp = [yc_exp, tc_exp, xc_exp]; +%! +%! yd_exp = [ -1.0000 5.5000 5.6000 +%! -0.6550 3.1673 4.2228 +%! -0.5421 2.6186 3.4968 ]; +%! +%! td_exp = [ 0 +%! 2 +%! 4 ]; +%! +%! xd_exp = [ 1.0000 2.0000 3.0000 +%! -0.4247 1.5194 2.3249 +%! -0.3538 1.2540 1.9250 ]; +%! +%! initial_d_exp = [yd_exp, td_exp, xd_exp]; +%! +%!assert (initial_c, initial_c_exp) +%!assert (initial_d, initial_d_exp) Added: trunk/octave-forge/extra/control-oo/inst/control/step.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/control/step.m (rev 0) +++ trunk/octave-forge/extra/control-oo/inst/control/step.m 2009-10-26 17:18:48 UTC (rev 6392) @@ -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 -*- + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +function [y_r, t_r, x_r] = step (sys, tfinal = [], dt = []) + + [y, t, x] = __timeresp__ (sys, "step", ! nargout, tfinal, dt); + + if (nargout) + y_r = y; + t_r = t; + x_r = x; + endif + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |