Example script to obtain the parameters of a MIMO ARX model (VARX) using least squares.
N = 1000; % number of samples
p = 5; % past window of VARX model
directfeedthrough = 1; % VARX model includes direct feedthrough
r = 2; % number of inputs
l = 3; % number of outputs
u = randn(r,N); % input data
y = randn(l,N); % output data
% concatenate the past data vectors
m = r+l;
z = [u; y];
Z = zeros(p*m,N-p);
for i = 1:p
Z((i-1)*m+1:i*m,:) = z(:,i:N+i-p-1);
end
% solve VARX problem
Y = y(:,p+1:N);
U = u(:,p+1:N);
if directfeedthrough
Z = [Z; U];
end
VARX = Y/Z; % least squares estimate
% convert solution to idpoly object
A = cell(l,l);
B = cell(l,r);
if directfeedthrough
VARX0 = [VARX eye(l)];
else
VARX0 = [VARX zeros(l,r) eye(l)];
end
for i = 1:l
for j = 1:l
A{i,j} = fliplr(VARX0(i,r+j:m:m*p+r+j));
end
for k = 1:r
B{i,k} = fliplr(VARX0(i,k:m:m*p+k));
end
end
Ts = 1; % sample time
E = Y - VARX*Z;
NoiseVariance = cov(E');
sys = idpoly(A,B,[],[],[],NoiseVariance,Ts); % MIMO ARX model object
For larger models, the least squares problem can become ill conditioned and would require regularization to get good estimate. If regularization is needed, I recommend to look at source code of the following matlab function.
