Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[11074] trunk/octave-forge/main/control/inst
|
From: <par...@us...> - 2012-09-22 10:23:30
|
Revision: 11074
http://octave.svn.sourceforge.net/octave/?rev=11074&view=rev
Author: paramaniac
Date: 2012-09-22 10:23:23 +0000 (Sat, 22 Sep 2012)
Log Message:
-----------
control: move multiplot time response into place
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/__time_response__.m
trunk/octave-forge/main/control/inst/step.m
Modified: trunk/octave-forge/main/control/inst/__time_response__.m
===================================================================
--- trunk/octave-forge/main/control/inst/__time_response__.m 2012-09-22 10:14:12 UTC (rev 11073)
+++ trunk/octave-forge/main/control/inst/__time_response__.m 2012-09-22 10:23:23 UTC (rev 11074)
@@ -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 <luk...@gm...>
## 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
+ 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;
- [F, G] = ssdata (sys); # matrices C and D don't change
+ ## time vector
+ t = @cellfun (@(dt) reshape (0 : dt : tfinal, [], 1), dt, "uniformoutput", false);
- n = rows (F); # number of states
- m = columns (G); # number of inputs
- p = rows (C); # number of outputs
+ ## 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)
- ## time vector
- t = reshape (0 : dt : tfinal, [], 1);
- l_t = length (t);
-
- switch (resptype)
+ switch (response)
case "initial"
- str = ["Response of ", sysname, " to Initial Conditions"];
- yfinal = zeros (p, 1);
+ [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"
+ [y, x] = cellfun (@__impulse_response__, sys_cell, sys_dt_cell, t, "uniformoutput", false);
+ otherwise
+ error ("time_response: invalid response type");
+ endswitch
- ## 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);
+ 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);
+
+ 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
+ hold off;
+ endif
- ## simulation
- for k = 1 : l_t
- y(k, :) = C * x;
- x_arr(k, :) = x;
- x = F * x;
- endfor
+endfunction
- 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);
+function [y, x_arr] = __initial_response__ (sys_dt, t, x0)
- for j = 1 : m # for every input channel
- ## initial conditions
- x = zeros (n, 1);
- u = zeros (m, 1);
- u(j) = 1;
+ [F, G, C, D] = ssdata (sys_dt); # system must be proper
- ## 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
+ n = rows (F); # number of states
+ m = columns (G); # number of inputs
+ p = rows (C); # number of outputs
+ l_t = length (t);
- case "impulse"
- str = ["Impulse Response of ", sysname];
- yfinal = zeros (p, m);
+ ## preallocate memory
+ y = zeros (l_t, p);
+ x_arr = zeros (l_t, n);
- ## preallocate memory
- y = zeros (l_t, p, m);
- x_arr = zeros (l_t, n, m);
+ ## initial conditions
+ x = reshape (x0, [], 1); # make sure that x is a column vector
- for j = 1 : m # for every input channel
- ## initial conditions
- u = zeros (m, 1);
- u(j) = 1;
+ if (n != length (x0) || ! is_real_vector (x0))
+ error ("initial: x0 must be a real vector with %d elements", n);
+ endif
- 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 = 1 : l_t
+ y(k, :) = C * x;
+ x_arr(k, :) = x;
+ x = F * x;
+ endfor
- ## 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
-
- otherwise
- error ("time_response: invalid response type");
-
- endswitch
-
+endfunction
- if (plotflag) # display plot
- ## TODO: Set correct titles, especially for multi-input systems
+function [y, x_arr] = __step_response__ (sys_dt, t)
- stable = isstable (sys);
- outname = get (sys, "outname");
- outname = __labels__ (outname, "y_");
+ [F, G, C, D] = ssdata (sys_dt); # system must be proper
- if (strcmp (resptype, "initial"))
- cols = 1;
- else
- cols = m;
- endif
+ n = rows (F); # number of states
+ m = columns (G); # number of inputs
+ p = rows (C); # number of outputs
+ l_t = length (t);
- if (discrete) # discrete system
- for k = 1 : p
- for j = 1 : cols
+ ## preallocate memory
+ y = zeros (l_t, p, m);
+ x_arr = zeros (l_t, n, m);
- subplot (p, cols, (k-1)*cols+j);
+ for j = 1 : m # for every input channel
+ ## initial conditions
+ x = zeros (n, 1);
+ u = zeros (m, 1);
+ u(j) = 1;
- if (stable)
- stairs (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]);
- else
- stairs (t, y(:, k, j));
- endif
+ ## 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
- grid ("on");
+endfunction
- if (k == 1 && j == 1)
- title (str);
- endif
- if (j == 1)
- ylabel (sprintf ("Amplitude %s", outname{k}));
- endif
+function [y, x_arr] = __impulse_response__ (sys, sys_dt, t)
- endfor
- endfor
+ [~, 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);
- xlabel ("Time [s]");
+ n = rows (F); # number of states
+ m = columns (G); # number of inputs
+ p = rows (C); # number of outputs
+ l_t = length (t);
- else # continuous system
- for k = 1 : p
- for j = 1 : cols
+ ## preallocate memory
+ y = zeros (l_t, p, m);
+ x_arr = zeros (l_t, n, m);
- subplot (p, cols, (k-1)*cols+j);
+ for j = 1 : m # for every input channel
+ ## initial conditions
+ u = zeros (m, 1);
+ u(j) = 1;
- if (stable)
- plot (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]);
- else
- plot (t, y(:, k, j));
- endif
+ 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
- grid ("on");
+ ## simulation
+ for k = 2 : l_t
+ y (k, :, j) = C * x;
+ x_arr(k, :, j) = x;
+ x = F * x;
+ endfor
+ endfor
- if (k == 1 && j == 1)
- title (str);
- endif
-
- if (j == 1)
- ylabel (sprintf ("Amplitude %s", outname{k}));
- endif
-
- endfor
- endfor
-
- xlabel ("Time [s]");
-
- endif
+ if (discrete)
+ y *= dt;
+ x_arr *= dt;
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
+ 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
Modified: trunk/octave-forge/main/control/inst/step.m
===================================================================
--- trunk/octave-forge/main/control/inst/step.m 2012-09-22 10:14:12 UTC (rev 11073)
+++ trunk/octave-forge/main/control/inst/step.m 2012-09-22 10:23:23 UTC (rev 11074)
@@ -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 <luk...@gm...>
## Created: October 2009
-## Version: 0.1
+## Version: 0.2
-function [y_r, t_r, x_r] = step (sys, tfinal = [], dt = [])
+function [y_r, t_r, x_r] = step2 (varargin)
- ## TODO: multiplot feature: step (sys1, "b", sys2, "r", ...)
-
- 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
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
