Demo of the "ident" package

1. Installation
2. Usage
3. Equivalence with PEM in the OE case
4. Permutation of variables
5. Unstable system
6. Trajectory of minimal length
7. Multiple short trajectories
8. Missing data
9. Step response identification
10. References

1 Installation

The software is hosted in GitHub. You can download it as a single zip file, or install it in a sub-directory \(\mathtt{slra}\) to your current directory by pasting the following commands into MATLAB:

%% Installation
unzip('https://github.com/slra/slra/zipball/master')
addpath(fullfile(cd, 'slra-slra-3e77deb'))

2 Usage

The package has two main functions:

[sys, info, wh] = ident(w, m, ell, opt); % identification 
[M, wh] = misfit(w, sys, opt);           % validation

data \(\mathtt{w}\) — set of time series \(\{\, w^1, \ldots, w^N \,\}\), \(w^k(t) \in \mathbb{R}^q\) specified by
- \(T\times q\) matrix — \(q\)-variate time series with \(T\) samples
- \(T\times q \times N\) array — \(N\), \(q\)-variate time series with \(T\) samples
- cell array of \(T_k\times q\) matrices — vector time series with different number of samples
- missing data — \(\mathtt{NaN}\) elements of \(\mathtt{w}\)
model \(\mathtt{sys}\) — LTI system in \(\mathtt{ss}\) object format
model class — \(\mathcal{L}_{m,\ell}\) is LTI systems with at most \(\mathtt{m}\) inputs and lag at most \(\mathtt{ell}\)
identification criterion — misfit \(\mathtt{M}\), defined as \begin{equation*} M(w,\mathcal{B}) := \min_{\hat w^1,\ldots,\hat w^N\in\mathcal{B}} \sqrt{\textstyle\sum_{k=1}^N \|w^k - \hat w^k\|^2_2 } \end{equation*}
identification problem, solved by \(\mathtt{ident}\) — \(\min_{\mathcal{B}\in\mathcal{L}_{m,\ell}} M(w,\mathcal{B})\)
extra options \(\mathtt{opt}\)
- \(\mathtt{exct}\) — exact/fixed variables (e.g., output error setup)
- \(\mathtt{wini}\) — fixed initial conditions (e.g., step response data)
- \(\mathtt{sys0}\), \(\mathtt{maxiter}\) — initial model and maximum # of iterations for the optimization method
online help — \(\mathtt{help ident}\)

3 Equivalence with PEM in the OE case

%% simulation parameters
clear all, randn('seed', 0), rand('seed', 0), warning off
m = 1; p = 1; ell = 2; T = 300; s = 0.1;     
opt_oe.exct = 1:m; % fixed inputs = output error identification
opt_eiv.exct = []; % errors-in-variables setup 

%% generate data
n = ell * p; q = m + p; sys0 = drss(n, p, m); % random "true" system
u0 = rand(T, m); xini0 = rand(n, 1);          % random "true" trajectory
y0 = lsim(sys0, u0, [], xini0); u = u0;       
yt = randn(T, p); y = y0 + s * norm(y0) * yt / norm(yt); % output error
w0 = [u0 y0]; w = [u y];

%% compare ident and pem
tic, sysh_ident = ident(w, m, ell, opt_oe); t_ident = toc;
tic, sysh_pem = pem(iddata(y, u), n); t_pem = toc;
[Yh, M] = compare(iddata(y, u), idss(sysh_ident), sysh_pem); 
ans = [[mean(M{1}); mean(M{2})] [t_ident; t_pem]]

86.831	0.0291
86.806	1.4404

4 Permutation of variables

%% the model for [y u] is the inverse of the model for [u y]
sysh1 = ident(w, m, ell, opt_eiv);
sysh2 = ident(fliplr(w), m, ell, opt_eiv);
ans = [eig(sysh1) eig(inv(sysh2))]

0.01129	0.011284
0.69769	0.69767

5 Unstable system

noisy data

%% the model for the reversed time data has reciprocal poles
sysh3 = ident(flipud(w), m, ell, opt_eiv); 
ans = [sort(1 ./ (eig(sysh1))) eig(sysh3)]

1.433	1.433
88.573	88.566

exact data

%% pem fails to identify unstable model
sysh_ident = ident(flipud(w0), m, ell, opt_oe); % note that now we use the true data w0
opt_pem = ssestOptions; opt_pem.Advanced.StabilityThreshold.z = inf; % disable the stability constraint
try 
  sysh_pem = pem(iddata(flipud(y0), flipud(u0)), n, opt_pem, 'dist', 'none');
catch
  sysh_pem = zeros(n); % zero indicates failure 
end
ans = [sort(1 ./ eig(sys0)) eig(sysh_ident) sort(eig(sysh_pem))]

1.510	1.510	0
25.751	25.751	0

6 Trajectory of minimal length

%% pem needs more data to identify a model
Tmin = ell + q * (ell + 1) - 1; TT = (Tmin - 1):15; np = length(TT);
work_ident = zeros(1, np); correct_ident = zeros(1, np); 
work_pem   = zeros(1, np); correct_pem   = zeros(1, np);
for i = 1:np
  try 
    sysh_ident = ident(w0(1:TT(i), :), m, ell, opt_oe); work_ident(i) = 1; 
    correct_ident(i) = norm(sys0 - sysh_ident) < 1e-5; 
  end
  try 
    sysh_pem = pem(iddata(y0(1:TT(i), :), u0(1:TT(i), :)), n, 'dist', 'none'); work_pem(i) = 1; 
    correct_pem(i) = norm(sys0 - sysh_pem) < 1e-5; 
  end
end 
ans = [TT; work_ident; correct_ident; work_pem; correct_pem]

6	7	8	9	10	11	12	13	14	15
0	0	1	1	1	1	1	1	1	1
0	0	1	1	1	1	1	1	1	1
0	0	0	0	0	0	0	1	1	1
0	0	0	0	0	0	0	1	1	1

7 Multiple short trajectories

%% N short trajectories equivalent to one long trajectory (N * Tshort = T)
Tshort = 13; N = round(T / Tshort); 
u_mult = {}; y_mult = {}; w_mult = []; 
for k = 1:N,
  u0 = rand(Tshort, m); xini0 = rand(n, 1); 
  y0 = lsim(sys0, u0, [], xini0); 
  yt = randn(Tshort, p); y_mult{k} = y0 + s * norm(y0) * yt / norm(yt); 
  u_mult{k} = u0; w_mult(:, :, k) = [u0 y_mult{k}]; 
end
tic, sysh = ident(w_mult, m, ell, opt_oe); t_ident = toc;
tic, sysh_pem = pem(iddata(y_mult, u_mult), n, 'dist', 'none'); t_pem = toc;
[Yh, M] = compare(iddata(y_mult, u_mult), idss(sysh), sysh_pem); 
ans = [[mean(M{1}); mean(M{2})] [t_ident; t_pem]]

88.303	0.0083
88.462	1.1033

8 Missing data

%% randomly distributed missing data in time and variables
Tp = 100; Im = randperm(q * Tp); Tm = round(0.2 * q * Tp); 
wm = w(1:Tp, :); wm(Im(1:Tm)) = NaN; um = wm(:, 1:m); ym = wm(:, m + 1:end);

%% use misdata for the comparison
tic, sysh_n4sid = ss(n4sid(misdata(iddata(ym, um), 10), n, 'dist', 'none')); t_n4sid = toc;
tic, opt_oe.sys0 = sysh_n4sid; sysh_ident = ident(wm, m, ell, opt_oe); t_ident = toc; opt_oe.sys0 = []; 
tic, sysh_pem = pem(misdata(iddata(ym, um), 10), sysh_n4sid, 'dist', 'none'); t_pem = toc;
[Yh, M] = compare(iddata(y0, u0), sysh_n4sid, idss(sysh_ident), sysh_pem); 
ans = [[mean(M{1}); mean(M{2}); mean(M{3})] [t_n4sid; t_ident; t_pem]]

58.51	6.8132
83.93	0.0726
60.36	7.3597

9 Step response identification

%% simulate step response data
ys0 = step(sys0); Ts = length(ys0); us = ones(Ts, m);
yt = randn(Ts, p); ys = ys0 + s * norm(ys0) * yt / norm(yt); 
opt_oe.w0 = 0; opt_oe.sys0 = sys0; [sysh, info, wh] = ident([us ys], m, ell, opt_oe);
figure('visible', 'off');
plot(1:Ts, ys0, '-r', 1:Ts, ys, ':k', 1:Ts, wh(:, end), '--b')
legend('true', 'noisy', 'approx.   .')
plot2svg demo_f.svg, set(gca, 'fontsize', 25), print -dpdf demo_f.pdf;

10 References

I. Markovsky and K. Usevich. Software for weighted structured low-rank approximation. J. Comput. Appl. Math., 256:278-292, 2014. [ bib | DOI | pdf | .html | Abstract ]
I. Markovsky. A software package for system identification in the behavioral setting. Control Engineering Practice, 21:1422-1436, 2013. [ bib | DOI | pdf | .html | Abstract ]