Code and simulation examples for "An application of system identification in metrology"

Instructions

  1. Download and unpack this archive.
  2. Install the IDENT package (needed by ident_io, ident_fd, and ident_fp).
  3. Execute the code listed in section "Examples".

Examples

Effect of the noise variance (section 4.1, simulation study)

clear all, a = 0.35; xini = 1; ub = -1; T = 30; N = 100;
sensor = @(a) c2d(ss(-a, a, 1, 0), 1); sys = sensor(a);
s = 0.02; name = ['cep-s' num2str(s)]; test_cep 
s = 0.03; name = ['cep-s' num2str(s)]; test_cep 
s = 0.04; name = ['cep-s' num2str(s)]; test_cep

Estimation with real data (section 4.1, real-life measurements)

clear all, load y-tea, name = 'cep-tea';
n = 1; g = 1; opt = []; 

[ub, sys, info, yh] = stepid_io(y, g, n, opt); sn = std(y - yh); 
gp = dcgain(sys); sys.b = sys.b / gp; sys.d = sys.d / gp; 

figure(3), clf, hold on % plot the estimation errors
T = length(y); t = 1:T; plot(t, y, ':k'), plot(t, yh, '-.r')
legend('data y   .', 'best fit', 'location', 'SouthEast'), 
ax = axis; axis([t(1) t(end) 35 ax(4)]), 
%print_fig('cep-tea-fit', 30) 

<<methods>> est_error = @(uh) sum(abs(ub - uh), 2); i = 1; 
<<run-methods>>
<<plot-results>>

Weight measurement (section 4.2)

clear all, 
m = 1; M = 1; k = 1; d = 1.1; g = 9.81; 
A = [0, 1; -k / (m + M), -d / (m + M)]; 
B = [0; -g / (m + M)]; C = [1, 0]; D = 0; 
sys = c2d(ss(A, B, C, D), 1); 
ub = M; xini = [0; 0]; T = 25; N = 100;
s = 0.01; name = ['cep-mass-s' num2str(s)]; test_cep 
s = 0.02; name = ['cep-mass-s' num2str(s)]; test_cep 
s = 0.03; name = ['cep-mass-s' num2str(s)]; test_cep

Thermometer and pressure sensors (not used in the paper)   P

clear all, A = [0.5 0.8]; tini = 1; ub = -1; T = 25; N = 50; S = [0.02 0.05];
sensor = @(a) c2d(ss(-diag(a), a(:), eye(length(a)), zeros(length(a), 1)), 1);
a = A(1); s = S(1); name = 'cep-thermometer'; test_cep 
a = A(2); s = S(2); name = 'cep-pressure'; test_cep 

test_methods = {'nv' 'dd' 'kf' 'fp' 'fd'};
a = A(1:2); s = S(1:2); name = 'cep-thermometer-pressure'; test_cep

Decay of the error when using multiple identical sensors (not used in the paper)

clear all, a = 0.2; tini = 1; ub = -1; T = 50; N = 10; s = 0.05;
sensor = @(a) c2d(ss(-a, a, 1, 0), 1);
p_ = 2; name = ['cep-p' num2str(s)]; test_cep 
p_ = 3; name = ['cep-p' num2str(s)]; test_cep 
p_ = 4; name = ['cep-p' num2str(s)]; test_cep

Decay of the error when using different sensors (not used in the paper)

clear all, A = [0.35 0.4 0.45 .5]; tini = 1; ub = -1; T = 50; N = 50; s = 0.05;
sensor = @(a) c2d(ss(-diag(a), a(:), eye(length(a)), zeros(length(a), 1)), 1);
a = A(1:2); name = ['cep-n' num2str(s)]; test_cep 
a = A(1:3); name = ['cep-n' num2str(s)]; test_cep 
a = A(1:4); name = ['cep-n' num2str(s)]; test_cep

Simulation function   P

test_cep

<<methods>>
if exist('p_'), 
  sys = ones(p_, 1) * sys; xini = kron(ones(p_, 1), xini);
end
g = dcgain(sys); p = size(g, 1); n = size(sys, 'order'); 

y0 = lsim(sys, ones(T, 1) * ub, [], xini);
est_error = @(uh) sum(abs(ub - uh), 2); opt = [];

for i = 1:N, i = i % Monte-Carlo simulation
  yn = randn(size(y0)); sn = s / norm(yn) * norm(y0); 
  if length(sn) ~= p, sn = sn * ones(1, p), end
  y = y0 + yn * diag(sn);
  <<run-methods>>
end
<<plot-results>>

<<methods>>

all_methods = {'nv' 'dd' 'kf' 'io' 'fp' 'fd'};
if ~exist('test_methods'), test_methods = all_methods; end, 
call_str = containers.Map(all_methods, ... % calling command as a string
  {'uh = stepid_nv(y, g);'
  'uh = stepid_dd(y, g, n);'
  'uh = stepid_kf(y, sys, diag(sn .^ 2));'
  '[uh, iter_io(i, :)] = recursify(@stepid_io, y, g, n, opt, 8);' 
  '[uh, iter_fp(i, :)] = recursify(@stepid_fp, y, g, n, opt, 8);'
  '[uh, iter_fd(i, :)] = recursify(@stepid_fd, y, g, n, opt, 8);'});
line_style = containers.Map(all_methods, ... % plotting line styles
  {''':k''', '''--b''', '''-.r''', '''-ob''', '''-sg''', '''-^g'''});

<<run-methods>>

for m = test_methods 
  %try, 
  eval([call_str(m{1}) ' e_' m{1} '(:, i) = est_error(uh);']), 
  %catch, eval(['e_' m '(:, i) = NaN; iter_' m '(i, :) = NaN;']), end
end

<<plot-results>>

figure(1), clf, hold on % plot the estimation errors
Tmin = (n + length(ub)) * 2 + 1; t = Tmin:T;
for m = test_methods, eval(['plot(t, mean(e_' m{1} '(t, :), 2),' line_style(m{1}) ')']), end
legend(test_methods), ax = axis; 
axis([Tmin T 0 mean(e_nv(Tmin, :))]), %print_fig(name, 30), pause(1) 
if 0
  Tmid = round(T / 2);
  axis([Tmid T 0 mean(e_dd(Tmid, :))]), legend off, print_fig([name '-close'], 30) 
end

figure(2), clf, hold on % plot the number of iterations
opt_test_methods = intersect(test_methods, {'io' 'fp' 'fd'});
for m = opt_test_methods, eval(['plot(t, mean(iter_' m{1} '(:, t), 1), ' line_style(m{1}) ')']), end
ax = axis; axis([Tmin T 0 ax(4)]), legend(opt_test_methods, 'location', 'SouthWest')
%print_fig([name '-iter'], 30), pause(1)

Estimation functions

stepid_nv

function uh = stepid_nv(y, g, n), uh = y * pinv(g');

stepid_dd

function uh = stepid_dd(y, g, n, ff)
if nargin == 3, ff = 1; end
T  = size(y, 1); [p, m] = size(g); 
dy = diff(y); yt = y((n + 1):T, :)'; Tt = T - n; 
if n == 0, dy = [0; dy]; end 
x  = rls([kron(ones(Tt, 1), g), blkhank(dy, Tt)], yt(:), ff);
uh = x(1:m, p:p:end)'; uh = [NaN * ones(n + 1, m); uh];

rls

function x = rls(a, b, f)
[m, n] = size(a); finv = 1 / f;
ai = a(1:n, 1:n); x = zeros(n, m); 
x(:, n) = ai \ b(1:n); p = inv(ai' * ai); 
for i = (n + 1):m
  ai = a(i, :);
  k  = finv * p * ai' / (1 + finv * ai * p * ai');
  x(:, i) = x(:, i - 1) + k * (b(i) - ai * x(:, i - 1));
  p  = 1 / f * (p - k * ai * p);
end

stepid_kf

function uh = stepid_kf(y, sys, v, x0, p0)
[a, b, c, d] = ssdata(sys); [p, m] = size(d); n = size(a, 1);
if ~exist('x0'), x0 = zeros(n + m, 1); end
if ~exist('p0'), p0 = 1e8 * eye(n + m); end
a_aut = [a b; zeros(m, n) eye(m)]; c_aut = [c d];
xeh = tvkf_oe(y, a_aut, c_aut, v, x0, p0); 
uh = xeh((n + 1):end, :)';

tvkf_oe

function x = tvkf_oe(y, a, c, v, x0, p)
T = size(y, 1); y = y'; 
x = zeros(size(a, 1), T); x(:, 1) = x0; 
for t = 1:(T-1)
  k  = (a * p * c') * pinv(v + c * p * c');
  x(:, t + 1) = a * x(:, t) + k * (y(:, t) - c * x(:, t));
  p  = a * p * a' - k * (v + c * p * c') * k';
end

ident_io   P

function [uh, sysh, info, yh] = stepid_io(y, g, n, opt)
[p, m] = size(g); T = size(y, 1); opt.exct = 1:m; ell = ceil(n / p); 
[sysh, info, wh] = ident([NaN * ones(ell, m + p); ones(T, m) y], m, ell, opt);
uh = sum(g \ dcgain(sysh), 2); yh = wh((ell + 1):end, (m + 1):end);

ident_fd   P

function [uh, sysh, info] = stepid_fd(y, g, n, opt)
T = size(y, 1); [p, m] = size(g); opt.n = n;
if isfield(opt, 'sys0') && ~isempty(opt.sys0.b)
  [a, b, c, d] = ssdata(opt.sys0);
  opt.sys0 = ss(a, [], c, [], 1);
end
[sysh, info] = ident(diff(y), 0, ceil(n / p), opt);
O = sysh.c; %ell = ceil((n + m) / p);
for t = 2:T, O = [O; O(end - p + 1:end, :) * sysh.a]; end 
yt = y(1:T, :)'; xe = [kron(ones(T, 1), g) O] \ yt(:); uh = xe(1:m);

ident_fp   P

function [uh, sysh, info] = stepid_fp(y, g, n, opt)
T = size(y, 1); [p, m] = size(g); 
opt.psi = toeplitz([1; -1; zeros(n, 1)], [1 zeros(1, n)])'; 
if isfield(opt, 'sys0') && ~isempty(opt.sys0.b)
  [a, b, c, d] = ssdata(opt.sys0); n = size(a, 1); 
  if ~isfield(opt, 'ub0'), opt.ub0 = ones(size(b, 2), 1); end
  opt.sys0 = ss([a b * opt.ub0; zeros(1, n) 1], [], [c d * opt.ub0], [], 1);
end
opt.n = n + 1;
[sysh, info, yh] = ident(y, 0, ceil((n + 1) / p), opt); 
sysh = canon(sysh, 'modal'); 
ind = find(abs(diag(sysh.a) - 1) < 1e-12); 
O = sysh.c; ell = ceil((n + m) / p);
for t = 2:ell, O = [O; O(end - p + 1:end, :) * sysh.a]; end 
yht = yh(1:ell, :)'; xini = O \ yht(:); uh = pinv(g) * sysh.c(:, ind) * xini(ind);

recursify   P

function [uh, iter] = recursify(method, y, g, n, opt, Tmin)
opt.maxiter = 100; T = size(y, 1);
if ~exist('Tmin'), Tmin = 2 * (n + 1) + 1; end
uh = zeros(T, 1); uh(1:Tmin - 1) = NaN; iter(1:Tmin - 1) = NaN;
for t = Tmin:T
  [uh(t), sysh, info] = method(y(1:t, :), g, n, opt);
  opt.sys0 = sysh; 
  iter(t) = info.iter;
end

Created: 2015-06-12 Fri 12:24

Emacs 24.3.1 (Org mode 8.2.10)

Validate