Demo of the ident package
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.39396 |
| 86.806 |
6.2899 |
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 ]
|
