## Copyright (C) 2009-2016 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 .
## -*- texinfo -*-
## Common code for the time response functions step, impulse and initial.
## Author: Lukas Reichlin
## Created: October 2009
## Version: 0.5
function [y, t, x] = __time_response__ (response, args, nout)
idx = cellfun (@islogical, args);
tmp = cellfun (@double, args(idx), "uniformoutput", false);
args(idx) = tmp;
sys_idx = cellfun (@isa, args, {"lti"}); # LTI models
mat_idx = cellfun (@is_real_matrix, args); # matrices
sty_idx = cellfun (@ischar, args); # strings (style arguments)
inv_idx = ! (sys_idx | mat_idx | sty_idx); # invalid arguments
if (any (inv_idx))
warning ("%s: arguments number %s are invalid and are being ignored", ...
response, mat2str (find (inv_idx)(:).'));
endif
if (nnz (sys_idx) == 0)
error ("%s: require at least one LTI model", response);
endif
if (nout > 0 && (nnz (sys_idx) > 1 || any (sty_idx)))
evalin ("caller", "print_usage ()");
endif
if (! size_equal (args{sys_idx}))
error ("%s: all LTI models must have equal size", response);
endif
if (any (find (sty_idx) < find (sys_idx)(1)))
warning ("%s: strings in front of first LTI model are being ignored", response);
endif
tfinal = []; dt = []; x0 = []; # default arguments
switch (response)
case "initial"
switch (nnz (mat_idx))
case 0
error ("initial: require initial state vector 'x0'");
case 1
x0 = args{mat_idx};
case 2
[x0, tfinal] = args{mat_idx};
case 3
[x0, tfinal, dt] = args{mat_idx};
otherwise
evalin ("caller", "print_usage ()");
endswitch
if (! is_real_vector (x0))
error ("initial: initial state vector 'x0' must be a real-valued vector");
endif
case {"step", "impulse", "ramp"}
switch (nnz (mat_idx))
case 0
## nothing to here, just prevent case 'otherwise'
case 1
tfinal = args{mat_idx};
case 2
[tfinal, dt] = args{mat_idx};
otherwise
evalin ("caller", "print_usage ()");
endswitch
otherwise
error ("time_response: invalid response type '%s'", response);
endswitch
switch (response)
case "step"
response1 = "zoh";
case "impulse"
response1 = "impulse";
otherwise
response1 = "zoh";
endswitch
if (issample (tfinal) || isempty (tfinal))
## nothing to do here
elseif (is_real_vector (tfinal))
dt = abs (tfinal(end) - tfinal(1)) / (length (tfinal) - 1);
tfinal = abs (tfinal(end));
else
evalin ("caller", "print_usage ()");
endif
if (isempty (dt))
## nothing to do here
elseif (issample (dt))
## nothing to do here
else
evalin ("caller", "print_usage ()");
endif
[tfinal, dt] = cellfun (@__sim_horizon__, args(sys_idx), {tfinal}, {dt}, "uniformoutput", false);
tfinal = max ([tfinal{:}]);
ct_idx = cellfun (@isct, args(sys_idx));
sys_dt = args(sys_idx);
tmp = cellfun (@c2d, args(sys_idx)(ct_idx), dt(ct_idx), {response1}, "uniformoutput", false);
sys_dt(ct_idx) = tmp;
## time vector
t = cellfun (@(dt) vec (linspace (0, tfinal, tfinal/dt+1)), dt, "uniformoutput", false);
## alternative code
## t = cellfun (@(dt) vec (0 : dt : tfinal), dt, "uniformoutput", false);
## function [y, x_arr] = __initial_response__ (sys_dt, t, x0)
## function [y, x_arr] = __step_response__ (sys_dt, t)
## function [y, x_arr] = __impulse_response__ (sys, sys_dt, t)
## function [y, x_arr] = __ramp_response__ (sys_dt, t)
switch (response)
case "initial"
[y, x] = cellfun (@__initial_response__, sys_dt, t, {x0}, "uniformoutput", false);
case "step"
[y, x] = cellfun (@__step_response__, sys_dt, t, "uniformoutput", false);
case "impulse"
[y, x] = cellfun (@__impulse_response__, args(sys_idx), sys_dt, t, "uniformoutput", false);
case "ramp"
[y, x] = cellfun (@__ramp_response__, sys_dt, t, "uniformoutput", false);
otherwise
error ("time_response: invalid response type");
endswitch
if (nout == 0) # display plot
## extract plotting styles
tmp = cumsum (sys_idx);
tmp(sys_idx | ! sty_idx) = 0;
n_sys = nnz (sys_idx);
sty = arrayfun (@(x) args(tmp == x), 1:n_sys, "uniformoutput", false);
## default plotting styles if empty
colororder = get (gca, "colororder");
rc = rows (colororder);
def = arrayfun (@(k) {"color", colororder(1+rem (k-1, rc), :)}, 1:n_sys, "uniformoutput", false);
idx = cellfun (@isempty, sty);
sty(idx) = def(idx);
## get system names for legend
leg = cell (1, n_sys);
idx = find (sys_idx);
for k = 1 : n_sys
leg{k} = evalin ("caller", sprintf ("inputname(%d)", idx(k)), "''");
endfor
outname = get (args(sys_idx){end}, "outname");
outname = __labels__ (outname, "y");
[p, m] = size (args(sys_idx){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);
case "ramp"
str = "Ramp Response";
cols = m;
otherwise
error ("time_response: invalid response type");
endswitch
for k = 1 : n_sys # for every system
if (ct_idx(k)) # continuous-time system
for i = 1 : p # for every output
for j = 1 : cols # for every input (except for initial where cols=1)
if (p != 1 || cols != 1)
subplot (p, cols, (i-1)*cols+j);
endif
plot (t{k}, y{k}(:, i, j), sty{k}{:});
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 # discrete-time system
for i = 1 : p # for every output
for j = 1 : cols # for every input (except for initial where cols=1)
if (p != 1 || cols != 1)
subplot (p, cols, (i-1)*cols+j);
endif
stairs (t{k}, y{k}(:, i, j), sty{k}{:});
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 (leg)
endif
hold off;
endif
endfunction
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_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
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 = G * u; #NO NO 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 [y, x_arr] = __ramp_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, l_t);
u(j, :) = t;
## simulation
for k = 1 : l_t
y(k, :, j) = C * x + D * u(:, k);
x_arr(k, :, j) = x;
x = F * x + G * u(:, k);
endfor
endfor
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
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 = 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
#prc_1
#Problem : pure imaginary poles are removed for calculating tfinal and dt. This could
# lead to choose tfinal = T_DEF and dt = tfinal / N_DEF and could cause
# aliasing (oscillating system)
#Test : sys=tf([0 0 1], [1e-10 0 1]); step(sys); t=0:5e-6:3e-3; hold on; step(sys, t,'r'); axis([0 0.0035 0 2])
#Solution: remove poles according to module and not real part
# if (abs (real (ev(k))) < TOL)
if (abs (ev(k)) < TOL)
ev(k) = 0;
nk -= 1;
endif
endfor
if (nk == 0)
if (isempty (tfinal))
tfinal = T_DEF;
endif
if (continuous)
dt = tfinal / N_DEF;
endif
else
ev = ev(find (ev));
ev_max = max (abs (ev));
if (continuous)
dt = 0.2 * pi / ev_max;
endif
if (isempty (tfinal))
ev_min = min (abs (real (ev)));
tfinal = 5.0 / ev_min;
## round up
#prc_2
#Problem: in the case of imaginary poles only (with modif prc_1), tfinal is Inf but was
# converted here to NaN => problem below with N and range tests on N.
#Solution: check if tfinal is Infinite before rounding up
if ~isinf(tfinal)
yy = 10^(ceil (log10 (tfinal)) - 1);
tfinal = yy * ceil (tfinal / yy);
endif
endif
if (continuous)
N = tfinal / dt;
if (N < N_MIN)
dt = tfinal / N_MIN;
endif
if (N > N_MAX)
#prc_3
#Problem : This was increasing dt: aliasing could occurs => could lead to wrong result
#Solution: modify tfinal rather than ts
#Test : sys=tf([0 0 1], [1e-10 1e-8 1]); step(sys); t=0:5e-6:3e-3; hold on; step(sys, t,'r'); axis([0.0015 0.0032 0 2])
#
# dt = tfinal / N_MAX;
tfinal=N_MAX*dt;
endif
endif
endif
if (continuous && ! isempty (Ts)) # catch case cont. system with dt specified
dt = Ts;
endif
#prc_4
#Problem : Generation of time vector at line 129 produces vector with modified dt if
# tfinal is not an entire multiple of dt. The consequence is that the resulting
# response is not correct (slow drift in time)
#Solution: Make tfinal an entire multiple of dt.
#Test: sys=tf([0 0 1], [1e-10 3e-7 1]); step(sys); t=0:5e-6:3e-3; hold on; step(sys, t,'r'); axis([0.0029 0.00302 0.985 1.015])
tfinal=round(tfinal/dt)*dt;
endfunction