From: <par...@us...> - 2012-09-15 14:36:28
|
Revision: 11027 http://octave.svn.sourceforge.net/octave/?rev=11027&view=rev Author: paramaniac Date: 2012-09-15 14:36:20 +0000 (Sat, 15 Sep 2012) Log Message: ----------- control: add draft code for multiplot time responses Added Paths: ----------- trunk/octave-forge/main/control/devel/__time_response_2__.m trunk/octave-forge/main/control/devel/step2.m Added: trunk/octave-forge/main/control/devel/__time_response_2__.m =================================================================== --- trunk/octave-forge/main/control/devel/__time_response_2__.m (rev 0) +++ trunk/octave-forge/main/control/devel/__time_response_2__.m 2012-09-15 14:36:20 UTC (rev 11027) @@ -0,0 +1,363 @@ +## Copyright (C) 2009, 2010 Lukas F. Reichlin +## +## This file is part of LTI Syncope. +## +## LTI Syncope is free software: you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## LTI Syncope is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with LTI Syncope. 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.2 + +% function [y, t, x_arr] = __time_response_2__ (sys, resptype, plotflag, tfinal, dt, x0, sysname) +function [y, t, x_arr] = __time_response_2__ (resptype, args) + + sys_idx = cellfun (@isa, args, {"lti"}); # look for LTI models + sys_cell = cellfun (@ss, args(sys_idx)); # system must be proper + + if (! size_equal (sys_cell{:})) + error ("%s: models must have equal sizes", resptype); + endif + + tmp = cellfun (@is_real_matrix, args); + vec_idx = find (tmp); + n_vec = length (vec_idx) + + if (n_vec >= 1) + arg = args{vec_idx(1)}; + + endif + + ## extract tfinal/t, dt, x0 + + [tfinal, dt] = cellfun (@__sim_horizon__, sys_cell, {tfinal}, dt); + + tfinal = max (tfinal); + + % __sim_horizon__ (sys, tfinal, dt); + + + hier sim_horizon + + ct_idx = cellfun (@isct, sys_cell) + sys_dt_cell = sys_cell; + tmp = (@c2d, sys_cell(ct_idx), dt, {"zoh"}, "uniformoutput", false) + sys_dt_cell(ct_idx) = tmp; + + +endfunction +%{ + if (! isa (sys, "ss")) + sys = ss (sys); # sys must be proper + endif + + if (is_real_vector (tfinal) && length (tfinal) > 1) # 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); + + discrete = ! isct (sys); # static gains are treated as analog systems + tsam = abs (tsam); # use 1 second if tsam is unspecified (-1) + + if (discrete) + if (! isempty (dt)) + warning ("time_response: argument dt has no effect on sampling time of discrete system"); + endif + + dt = tsam; + endif + + [tfinal, dt] = __sim_horizon__ (A, discrete, tfinal, dt); + + if (! discrete) + 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 = reshape (0 : dt : tfinal, [], 1); + l_t = length (t); +%} + +function [y, x_arr] = __initial_response__ (sys, sys_dt, t, x0) + + [A, B, C, D] = ssdata (sys); + [F, G] = ssdata (sys_dt); + + n = rows (F); # number of states + m = columns (G); # number of inputs + p = rows (C); # number of outputs + l_t = length (t); + + ## preallocate memory + y = zeros (l_t, p); + x_arr = zeros (l_t, n); + + ## initial conditions + x = reshape (x0, [], 1); # make sure that x is a column vector + + if (n != length (x0) || ! is_real_vector (x0)) + error ("initial: x0 must be a real vector with %d elements", n); + endif + + ## simulation + for k = 1 : l_t + y(k, :) = C * x; + x_arr(k, :) = x; + x = F * x; + endfor + +endfunction + + +function [y, x_arr] = __step_response__ (sys_dt, t) + + [F, G, C, D] = ssdata (sys_dt); + + n = rows (F); # number of states + m = columns (G); # number of inputs + p = rows (C); # number of outputs + l_t = length (t); + + ## 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 + +endfunction + + +function [y, x_arr] = __impulse_response__ (sys, sys_dt, t) + + [A, B, C, D, dt] = ssdata (sys); + [F, G] = ssdata (sys_dt); + discrete = ! isct (sys); + + n = rows (F); # number of states + m = columns (G); # number of inputs + p = rows (C); # number of outputs + l_t = length (t); + + ## 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 (discrete) + 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 (discrete) + y *= dt; + x_arr *= dt; + endif + +endfunction + +%{ + if (plotflag) # display plot + + ## TODO: Set correct titles, especially for multi-input systems + + stable = isstable (sys); + outname = get (sys, "outname"); + outname = __labels__ (outname, "y_"); + + if (strcmp (resptype, "initial")) + cols = 1; + else + cols = m; + endif + + if (discrete) # 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] = __sim_horizon__ (A, discrete, tfinal, Ts) +function [tfinal, dt] = __sim_horizon__ (sys, 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 + + ev = pole (sys); + n = length (ev); + + if (discrete) + ## perform bilinear transformation on poles in z + for k = 1 : n + pol = ev(k); + if (abs (pol + 1) < TOL) + ev(k) = 0; + else + ev(k) = 2 / Ts * (pol - 1) / (pol + 1); + endif + endfor + endif + + ## remove poles near zero from eigenvalue array ev + nk = n; + for k = 1 : n + if (abs (real (ev(k))) < TOL) + ev(k) = 0; + nk -= 1; + endif + endfor + + if (nk == 0) + if (isempty (tfinal)) + tfinal = T_DEF; + endif + + if (! discrete) + dt = tfinal / N_DEF; + endif + else + ev = ev(find (ev)); + ev_max = max (abs (ev)); + + if (! discrete) + dt = 0.2 * pi / ev_max; + endif + + if (isempty (tfinal)) + ev_min = min (abs (real (ev))); + tfinal = 5.0 / ev_min; + + ## round up + yy = 10^(ceil (log10 (tfinal)) - 1); + tfinal = yy * ceil (tfinal / yy); + endif + + if (! discrete) + 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 cont. system with dt specified + dt = Ts; + endif + +endfunction Added: trunk/octave-forge/main/control/devel/step2.m =================================================================== --- trunk/octave-forge/main/control/devel/step2.m (rev 0) +++ trunk/octave-forge/main/control/devel/step2.m 2012-09-15 14:36:20 UTC (rev 11027) @@ -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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} {[@var{y}, @var{t}, @var{x}] =} step (@var{sys}) +## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} step (@var{sys}, @var{t}) +## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} step (@var{sys}, @var{tfinal}) +## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} step (@var{sys}, @var{tfinal}, @var{dt}) +## Step response of LTI system. +## If no output arguments are given, the response is printed on the screen. +## +## @strong{Inputs} +## @table @var +## @item sys +## LTI model. +## @item t +## Time vector. Should be evenly spaced. If not specified, it is calculated by +## the poles of the system to reflect adequately the response transients. +## @item tfinal +## Optional simulation horizon. If not specified, it is calculated by +## the poles of the system to reflect adequately the response transients. +## @item dt +## Optional sampling time. Be sure to choose it small enough to capture transient +## phenomena. If not specified, it is calculated by the poles of the system. +## @end table +## +## @strong{Outputs} +## @table @var +## @item y +## Output response array. Has as many rows as time samples (length of t) +## and as many columns as outputs. +## @item t +## Time row vector. +## @item x +## State trajectories array. Has @code{length (t)} rows and as many columns as states. +## @end table +## +## @seealso{impulse, initial, lsim} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +% function [y_r, t_r, x_r] = step2 (sys, tfinal = [], dt = []) +function [y_r, t_r, x_r] = step2 (varargin) + + if (nargin == 0) + print_usage (); + endif + + %tmp = cellfun (@isa, varargin, {"lti"}); + %sys_idx = find (tmp); + + + %sys_names = arrayfun (@inputname, sys_idx); + + [y, t, x] = __time_response__ ("step", varargin, ! nargout); + + 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. |