Mercurial > forge
changeset 10888:256b88bc1861 octave-forge
control: move multiplot time response into place
| author | paramaniac |
|---|---|
| date | Sat, 22 Sep 2012 10:23:23 +0000 |
| parents | 5af195699f20 |
| children | 63e17eb83383 |
| files | main/control/inst/__time_response__.m main/control/inst/step.m |
| diffstat | 2 files changed, 334 insertions(+), 212 deletions(-) [+] |
line wrap: on
line diff
--- a/main/control/inst/__time_response__.m Sat Sep 22 10:14:12 2012 +0000 +++ b/main/control/inst/__time_response__.m Sat Sep 22 10:23:23 2012 +0000 @@ -1,4 +1,4 @@ -## Copyright (C) 2009, 2010 Lukas F. Reichlin +## Copyright (C) 2009, 2010, 2012 Lukas F. Reichlin ## ## This file is part of LTI Syncope. ## @@ -20,241 +20,350 @@ ## Author: Lukas Reichlin <lukas.reichlin@gmail.com> ## Created: October 2009 -## Version: 0.2 +## Version: 0.3 -function [y, t, x_arr] = __time_response__ (sys, resptype, plotflag, tfinal, dt, x0, sysname) +function [y, t, x] = __time_response__ (response, args, sysname, plotflag) - if (! isa (sys, "ss")) - sys = ss (sys); # sys must be proper - endif + sys_idx = find (cellfun (@isa, args, {"lti"})); # look for LTI models, 'find' needed for plot styles + sys_cell = cellfun (@ss, args(sys_idx), "uniformoutput", false); # convert to state-space - 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); + if (! size_equal (sys_cell{:})) + error ("%s: models must have equal sizes", response); endif - [A, B, C, D, tsam] = ssdata (sys); + vec_idx = find (cellfun (@is_real_matrix, args)); # indices of vector arguments + n_vec = length (vec_idx); # number of vector arguments + n_sys = length (sys_cell); # number of LTI systems - discrete = ! isct (sys); # static gains are treated as analog systems - tsam = abs (tsam); # use 1 second if tsam is unspecified (-1) + tfinal = []; + dt = []; + x0 = []; - if (discrete) - if (! isempty (dt)) - warning ("time_response: argument dt has no effect on sampling time of discrete system"); + ## extract tfinal/t, dt, x0 from args + if (strcmpi (response, "initial")) + if (n_vec < 1) + error ("initial: require initial state vector 'x0'"); + else # initial state vector x0 specified + arg = args{vec_idx(1)}; + if (is_real_vector (arg)) + x0 = arg; + else + error ("initial: initial state vector 'x0' must be a vector of real values"); + endif + if (n_vec > 1) # tfinal or time vector t specified + arg = args{vec_idx(2)}; + if (issample (arg)) + tfinal = arg; + elseif (isempty (arg)) + ## tfinal = []; # nothing to do here + elseif (is_real_vector (arg)) + dt = abs (arg(2) - arg(1)); # assume that t is regularly spaced + tfinal = arg(end); + else + warning ("initial: argument number %d ignored", vec_idx(2)); + endif + if (n_vec > 2) # sampling time dt specified + arg = args{vec_idx(3)}; + if (issample (arg)) + dt = arg; + else + warning ("initial: argument number %d ignored", vec_idx(3)); + endif + if (n_vec > 3) + warning ("initial: ignored"); + endif + endif + endif + endif + else # step or impulse response + if (n_vec > 0) # tfinal or time vector t specified + arg = args{vec_idx(1)}; + if (issample (arg)) + tfinal = arg; + elseif (isempty (arg)) + ## tfinal = []; # nothing to do here + elseif (is_real_vector (arg)) + dt = abs (arg(2) - arg(1)); # assume that t is regularly spaced + tfinal = arg(end); + else + warning ("%s: argument number %d ignored", response, vec_idx(1)); + endif + if (n_vec > 1) # sampling time dt specified + arg = args{vec_idx(2)}; + if (issample (arg)) + dt = arg; + else + warning ("%s: argument number %d ignored", response, vec_idx(2)); + endif + if (n_vec > 2) + warning ("%s: ignored", response); + endif + endif endif - - dt = tsam; endif + ## TODO: share common code between initial and step/impulse - [tfinal, dt] = __sim_horizon__ (A, discrete, tfinal, dt); + [tfinal, dt] = cellfun (@__sim_horizon__, sys_cell, {tfinal}, {dt}, "uniformoutput", false); + tfinal = max ([tfinal{:}]); - 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 + ct_idx = cellfun (@isct, sys_cell); + sys_dt_cell = sys_cell; + tmp = cellfun (@c2d, sys_cell(ct_idx), dt(ct_idx), {"zoh"}, "uniformoutput", false); + sys_dt_cell(ct_idx) = tmp; ## time vector - t = reshape (0 : dt : tfinal, [], 1); - l_t = length (t); - - switch (resptype) - case "initial" - str = ["Response of ", sysname, " to Initial Conditions"]; - yfinal = zeros (p, 1); - - ## 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 + t = @cellfun (@(dt) reshape (0 : dt : tfinal, [], 1), dt, "uniformoutput", false); - ## simulation - for k = 1 : l_t - y(k, :) = C * x; - x_arr(k, :) = x; - x = F * x; - endfor - - case "step" - str = ["Step Response of ", sysname]; - 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; + ## function [y, x_arr] = __initial_response__ (sys, sys_dt, t, x0) + ## function [y, x_arr] = __step_response__ (sys_dt, t) + ## function [y, x_arr] = __impulse_response__ (sys, sys_dt, t) - ## 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 - + switch (response) + case "initial" + [y, x] = cellfun (@__initial_response__, sys_dt_cell, t, {x0}, "uniformoutput", false); + case "step" + [y, x] = cellfun (@__step_response__, sys_dt_cell, t, "uniformoutput", false); case "impulse" - str = ["Impulse Response of ", sysname]; - 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 (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 - + [y, x] = cellfun (@__impulse_response__, sys_cell, sys_dt_cell, t, "uniformoutput", false); otherwise error ("time_response: invalid response type"); - endswitch + + if (plotflag) # display plot + [p, m] = size (sys_cell{1}); + switch (response) + case "initial" + str = "Response to Initial Conditions"; + cols = 1; + ## yfinal = zeros (p, 1); + case "step" + str = "Step Response"; + cols = m; + ## yfinal = dcgain (sys_cell{1}); + case "impulse" + str = "Impulse Response"; + cols = m; + ## yfinal = zeros (p, m); + otherwise + error ("time_response: invalid response type"); + endswitch + + style_idx = find (cellfun (@ischar, args)); + outname = get (sys_cell{end}, "outname"); + outname = __labels__ (outname, "y"); + colororder = get (gca, "colororder"); + rc = rows (colororder); - 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; + for k = 1 : n_sys # for every system + if (k == n_sys) + lim = numel (args); + else + lim = sys_idx(k+1); + endif + style = args(style_idx(style_idx > sys_idx(k) & style_idx <= lim)); + if (isempty (style)) + color = colororder(1+rem (k-1, rc), :); + style = {"color", color}; + endif + discrete = ! ct_idx(k); + if (discrete) # discrete-time system + for i = 1 : p # for every output + for j = 1 : cols # for every input (except for initial where cols=1) + subplot (p, cols, (i-1)*cols+j); + stairs (t{k}, y{k}(:, i, j), style{:}); + hold on; + grid on; + if (k == n_sys) + axis tight; + ylim (__axis_margin__ (ylim)) + if (j == 1) + ylabel (outname{i}); + if (i == 1) + title (str); + endif + endif + endif + endfor + endfor + else # continuous-time system + for i = 1 : p # for every output + for j = 1 : cols # for every input (except for initial where cols=1) + subplot (p, cols, (i-1)*cols+j); + ##if (n_sys == 1 && isstable (sys_cell{1})) + ## plot (t{k}, y{k}(:, i, j), style{:}, [t{k}(1), t{k}(end)], repmat (yfinal(i,j), 1, 2)); + ## ## TODO: plot final value first such that its line doesn't overprint the response + ##else + plot (t{k}, y{k}(:, i, j), style{:}); + ##endif + hold on; + grid on; + if (k == n_sys) + axis tight + ylim (__axis_margin__ (ylim)) + if (j == 1) + ylabel (outname{i}); + if (i == 1) + title (str); + endif + endif + endif + endfor + endfor + endif + endfor + xlabel ("Time [s]"); + if (p == 1 && m == 1) + legend (sysname) 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 + hold off; endif endfunction -function [tfinal, dt] = __sim_horizon__ (A, discrete, tfinal, Ts) +function [y, x_arr] = __initial_response__ (sys_dt, t, x0) + + [F, G, C, D] = ssdata (sys_dt); # system must be proper + + 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); # system must be proper + + 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) + + [~, B] = ssdata (sys); + [F, G, C, D, dt] = ssdata (sys_dt); # system must be proper + dt = abs (dt); # use 1 second if tsam is unspecified (-1) + 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 + + +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 + 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); + ev = pole (sys); + n = length (ev); # number of states/poles + continuous = isct (sys); + discrete = ! continuous; if (discrete) + dt = Ts = abs (get (sys, "tsam")); ## perform bilinear transformation on poles in z for k = 1 : n - pol = eigw(k); + pol = ev(k); if (abs (pol + 1) < TOL) - eigw(k) = 0; + ev(k) = 0; else - eigw(k) = 2 / Ts * (pol - 1) / (pol + 1); + ev(k) = 2 / Ts * (pol - 1) / (pol + 1); endif endfor endif - ## remove poles near zero from eigenvalue array eigw + ## remove poles near zero from eigenvalue array ev nk = n; for k = 1 : n - if (abs (real (eigw(k))) < TOL) - eigw(k) = 0; + if (abs (real (ev(k))) < TOL) + ev(k) = 0; nk -= 1; endif endfor @@ -264,27 +373,27 @@ tfinal = T_DEF; endif - if (! discrete) + if (continuous) dt = tfinal / N_DEF; endif else - eigw = eigw(find (eigw)); - eigw_max = max (abs (eigw)); + ev = ev(find (ev)); + ev_max = max (abs (ev)); - if (! discrete) - dt = 0.2 * pi / eigw_max; + if (continuous) + dt = 0.2 * pi / ev_max; endif if (isempty (tfinal)) - eigw_min = min (abs (real (eigw))); - tfinal = 5.0 / eigw_min; + 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) + if (continuous) N = tfinal / dt; if (N < N_MIN) @@ -297,7 +406,7 @@ endif endif - if (! isempty (Ts)) # catch case cont. system with dt specified + if (continuous && ! isempty (Ts)) # catch case cont. system with dt specified dt = Ts; endif
--- a/main/control/inst/step.m Sat Sep 22 10:14:12 2012 +0000 +++ b/main/control/inst/step.m Sat Sep 22 10:23:23 2012 +0000 @@ -1,4 +1,4 @@ -## Copyright (C) 2009 Lukas F. Reichlin +## Copyright (C) 2009, 2012 Lukas F. Reichlin ## ## This file is part of LTI Syncope. ## @@ -54,22 +54,35 @@ ## Author: Lukas Reichlin <lukas.reichlin@gmail.com> ## Created: October 2009 -## Version: 0.1 - -function [y_r, t_r, x_r] = step (sys, tfinal = [], dt = []) +## Version: 0.2 - ## TODO: multiplot feature: step (sys1, "b", sys2, "r", ...) +function [y_r, t_r, x_r] = step2 (varargin) - if (nargin == 0 || nargin > 3) + if (nargin == 0) print_usage (); endif - [y, t, x] = __time_response__ (sys, "step", ! nargout, tfinal, dt, [], inputname (1)); + if (nargout) + sysname = {}; + else + sys_idx = find (cellfun (@isa, varargin, {"lti"})); + len = length (sys_idx); + sysname = cell (len, 1); + for k = 1 : len + try + sysname{k} = inputname(sys_idx(k)); + catch + sysname{k} = ""; + end_try_catch + endfor + endif + + [y, t, x] = __time_response_2__ ("step", varargin, sysname, ! nargout); if (nargout) - y_r = y; - t_r = t; - x_r = x; + y_r = y{1}; + t_r = t{1}; + x_r = x{1}; endif endfunction