Code covered by the BSD License  

Highlights from
MIMOtool

image thumbnail

MIMOtool

by

 

12 Nov 2001 (Updated )

Multi Input Multi Output Systems Toolbox

[th,datan]=rarx4(y,u,mxs,adm,adg,data0)
function [th,datan]=rarx4(y,u,mxs,adm,adg,data0)

% [th,datan]=rarx4(y,u,mxs,adm,adg,data0) MIMO model obtained joining 
% all the MISO recursive arx models between inputs and each output.
% Each rarx model is chosen on the first time, with the structure that achieve 
% the smallest AIC, according to a simple search algorithm, in which the first 
% half of data is used for estimation, the second for cross validation.
% Each column of y,(u) contains the story of the related output,(input).
% The stable part of the whole system is balanced and truncated in order
% to reduce the states to mxs if possible.
% mxs=inf simply skips balance and truncation, it is advisable to skip it
% during recursive use, since it requires many computations.
% data0 should be zero in the first execution of this function, then, 
% once the structure is selected, use data0 as the output datan of the last
% execution, it contains parameters, covariance matrix, regression vectors, for
% each of the transfer functions.
% adm: Adaptation mechanism. adg: Adaptation gain
% adm='ff', adg=lam:  Forgetting factor algorithm, with forg factor lam
% adm='kf', adg=R1: The Kalman filter algorithm with R1 as covariance 
% matrix of the parameter changes per time step
% adm='ng', adg=gam: A normalized gradient algorithm, with gain gam
% adm='ug', adg=gam: An Unnormalized gradient algorithm with gain gam
%
% Example:
% [a,b,c,d]=unpck(bilins2z(sysrand(4,2,3,1),1));
% th0=ms2th(modstruc(a,b,c,d,zeros(size(c'))),'d');
% e=rand(200,2);u=rand(200,3);y=idsim([u e],th0);
% [th4,datan]=rarx4(y(1:100,:),u(1:100,:),inf,'ff',0.8,0);
% for i=1:100, [th4,datan]=rarx4(y(i,:),u(i,:),inf,'ff',0.8,datan); end
% vid=ver('ident'); if vid.Version(1)=='4', [a4,b4,c4,d4,k4,x04]=th2ss(th4); 
% else [a4,b4,c4,d4,k4,x04]=ssdata(th2ido(th4)); end
% subplot(2,1,1),dsigma(a,b,c,d,1);
% subplot(2,1,2),dsigma(a4,b4,c4,d4,1);

% This is the jtools version (with mxs).
% It requires the files sys2sys,sysbal3 and infnan2x. 

% Giampiero Campa 12-12-95



% initialization

no=size(y,2);
ni=size(u,2);

dt1=[];
dt2=[];

fex=size(data0,1);
f=size(y,1);
l=floor(f/2);
cnt=0;

% outer cycle (constructs matrices row by row)
A=[];B=[];C=[];D=[];K=[];X0=[];
for i=1:no,

cnt=cnt+1;

if fex<2, % first time: computes optimal structure

n0=struc(1:3,1:3,1:3);
n0=[n0(:,1) n0(:,2)*ones(1,ni) n0(:,3)*ones(1,ni)];
n1=sls4(arxstruc([y(1:l,i) u(1:l,:)],[y(l+1:f,i) u(l+1:f,:)],n0),'AIC');
vtx=(n1==3);
c=0;

while any(vtx)&(c<9),
n0=struc(n1(1)-vtx(1):n1(1)+vtx(1),n1(2)-vtx(2):n1(2)+vtx(2),n1(2+ni)-vtx(2+ni):n1(2+ni)+vtx(2+ni));
n0=[n0(:,1) n0(:,2)*ones(1,ni) n0(:,3)*ones(1,ni)];
n2=sls4(arxstruc([y(1:l,i) u(1:l,:)],[y(l+1:f,i) u(l+1:f,:)],n0),'AIC');
vtx=(n2==(n1+1));
n1=n2;c=c+1;
end

% execution of algorithm for first time
[tmn,yh,Pn,Phn]=rarx([y(:,i) u],n1,adm,adg);

else
% NOT FIRST TIME, it is all in data0

n1=data0(1,(8+2*ni)*(cnt-1)+[1:1+2*ni]);

la1=data0(1,(8+2*ni)*(cnt-1)+2+2*ni);
la2=data0(1,(8+2*ni)*(cnt-1)+3+2*ni);
lb1=data0(1,(8+2*ni)*(cnt-1)+4+2*ni);
lb2=data0(1,(8+2*ni)*(cnt-1)+5+2*ni);
lc1=data0(1,(8+2*ni)*(cnt-1)+6+2*ni);
lc2=data0(1,(8+2*ni)*(cnt-1)+7+2*ni);
ldt=data0(1,(8+2*ni)*(cnt-1)+8+2*ni);

% construction of matrices from second row of data0
tm0=reshape(data0(2,ldt+[1:la1*la2]),la1,la2);
P0=reshape(data0(2,ldt+la1*la2+[1:lb1*lb2]),lb1,lb2);
Ph0=reshape(data0(2,ldt+la1*la2+lb1*lb2+[1:lc1*lc2]),lc1,lc2);

% execution of recursive algorithm
[tmn,yh,Pn,Phn]=rarx([y(:,i) u],n1,adm,adg,tm0(la1,:),P0,Ph0);

% end if first time
end

% takes only last row of tmn to avoid data explosion 
% anyway could be skipped without many consequences
tmn=tmn(size(tmn,1),:);

% new data construction
la1=size(tmn,1);
la2=size(tmn,2);
lb1=size(Pn,1);
lb2=size(Pn,2);
lc1=size(Phn,1);
lc2=size(Phn,2);
ldt=size(dt2,2);

dt1=[dt1 n1 la1 la2 lb1 lb2 lc1 lc2 ldt];
dt2=[dt2 reshape(tmn,1,la1*la2) reshape(Pn,1,lb1*lb2) reshape(Phn,1,lc1*lc2)];

% puts theta in state space form
BB=[];
for j=1:ni
BB=[BB zeros(1,n1(1+ni+j)) tmn(la1,n1(1)+[1:n1(1+j)])];
end
AA=[1 tmn(la1,1:n1(1))];
th=arx2th(AA,BB,1,ni);
[a,b,c,d,k,x0]=th2ss(th);

% state space matrices construction
A=[A zeros(size(A,1),size(a,2)); zeros(size(a,1),size(A,2)) a];
B=[B;b];
C=[C zeros(size(C,1),size(c,2)); zeros(size(c,1),size(C,2)) c];
D=[D;d];
K=[K zeros(size(K,1),size(k,2)); zeros(size(k,1),size(K,2)) k];
X0=[X0; x0];

% outer cycle
end

% datan from dt1 and dt2
l1=size(dt1,2);
l2=size(dt2,2);
if l2>l1,
datan=[dt1 zeros(1,l2-l1);dt2];
else
datan=[dt1; dt2 zeros(1,l1-l2)];
end

% balance and truncation if required
if ~isinf(mxs),

  % system in continuous time domain
  s0=bilinz2s(pck(A,[B K X0(1:size(A,1))],C,[D eye(size(K,2)) zeros(size(D,1),1)]),1);

  % balance and truncation
  s2=sys2sys(s0,'r',mxs);
  [ty2,no2,ni2,ns2]=minfo(s2);

  % final discrete time system
  s3=bilins2z(s2,1);

  A=s3(1:ns2,1:ns2);
  B=s3(1:ns2,ns2+[1:ni]);
  C=s3(ns2+[1:no],1:ns2);
  D=s3(ns2+[1:no],ns2+[1:ni]);
  K=s3(1:ns2,ns2+ni+[1:no]);
  X0=s3(1:ns2,ns2+ni+no+1);

% end if
end 

% final theta format
th=ss4th(A,B,C,D,K,X0);

% Comparation with spectral analysis if is the first time or a reduction is required
%if (~isinf(mxs) | fex<2),
%bodeplot([th2ff(th) spa([y u],30*ones(1,no))]);
%end

% subfunctions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function theta=ss4th(a,b,c,d,k,x0,cd,parval,lambda,T)

% THETA = ss4th(A,B,C,D,K,X0,CD,PARVAL,LAMBDA,T)
% Pack standard ss parameterizations into THETA format

if nargin<10, T=[];end
if nargin<9, lambda=[];end
if nargin<8, parval=[];end 
if nargin<7, cd=[];end

if isempty(T),T=1;end

if nargin<5,error('All of the matrices A, B, C, D, K, must be specified!'),end
[nx,nd]=size(a);
if nx ~= nd, error('The A-matrix must be square!'),end
if nargin<6, x0=zeros(nx,1);end,if isempty(x0),x0=zeros(nx,1);end
[nd,nu]=size(b);
if nx ~= nd & nu~=0, error('A and B must have the same number of rows!'),
end
[ny,nd]=size(c);
if nx ~= nd, error('A and C must have the same number of columns!'),
end
[nd1,nd2]=size(d);
if nd1 ~= ny & nu~=0, error('D and C must have the same number of rows!'),end
if nd2 ~= nu, error('D and B must have the same number of columns!'),
end
[nd1,nd2]=size(k);
if nd1 ~= nx, error('K and A must have the same numb er of rows!'),
end
if nd2 ~= ny, error('The number of columns in K must equal the number of coulumns in C!'),end
[nd1,nd2]=size(x0);
if nd1~=nx, error('X0 and A must have the same number of rows!'),end
if nd2 ~= 1, error('X0 must be a coulmn vector!'),end
ms(1:nx,1:nx)=a;
if nu>0, ms(1:nx,nx+1:nx+nu)=b;end
ms(1:nx,nx+nu+1:nx+nu+ny)=c';
if nu>0, ms(1:ny,nx+nu+ny+1:nx+2*nu+ny)=d;end
nn=nx+2*(nu+ny);
ms(1:nx,nx+2*nu+ny+1:nn)=k;
ms(1:nx,nn+1:nn+1)=x0;
ms(1,nn+2)=ny;ms(2,nn+2)=nx;

d1=sum(sum(isnan(ms))');

if isempty(parval),parval=zeros(1,d1);end
if isempty(cd), cd='d';end, cd=cd(1);
if cd~='c' & cd~='d', 
  error('CD must be one of ''c''(ontinuous) or ''d''(iscrete)')
end
[nx,nn]=size(ms);nyy=ms(1,nn);
if nyy<0,arx=1;else arx=0;end
d=length(parval);
if d~=d1 & ~arx, 
  error('Incorrect number of parameter values  have been specified; must be equal to the number of ''NaN'':s in the model structure!')
end

if T<0,
  error('The sampling interval T must be a positive number. Use ''c'' in the second argument to indicate continuous time model!')
end

[rarg,carg]=size(ms);
if cd=='c',theta(1,2)=-abs(T);else theta(1,2)=abs(T);end
if cd=='d',Tmod=-1;else Tmod=abs(T);end
theta(1,6:7)=[rarg,carg];

[a,b,c,Dd,k,x0]=ssmodx4(parval,Tmod,ms);

if ms(1,carg)>0,
   if max(abs(eig(a-k*c)))>1,disp('WARNING: Predictor unstable!'),end
end
[ny,nx]=size(c);[nx,nu]=size(b);
if isempty(lambda),lambda=eye(ny);end
theta(1,1)=det(lambda);
ti=fix(clock);
ti(1)=ti(1)/100;
theta(1,3:5)=[nu ny d];theta(2,1)=theta(1,1);
theta(2,2:6)=ti(1:5);
theta(2,7)=21;
if cd=='c', theta(2,8)=1;else theta(2,8)=2;end
if d>0,
	theta(3,1:d)=parval;
end
theta(4+d:3+d+rarg,1:carg)=ms;
theta(4+d+rarg:3+ny+d+rarg,1:ny)=lambda;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [nn,Vmod]=sls4(V,c)

%
% nn = sls4(V,c)    or   [nn,Vm] = sls(V,c)
% copy of selstruc.m with graphic mode disabled,
% meant to be used as a subroutine of arx4 and rarx4
%

nn=[];
if nargin<2,c='a';,end
if c<0,c='a';end
if c<0,c='a';end
[nl1,nm1]=size(V);
nu=floor((nl1-2)/2);
Nc=V(1,nm1);
if c(1)=='a' | c(1)=='A',alpha=2;end
if c(1)=='m' | c(1)=='M', alpha=log(Nc);end

for kj=1:nm1-1
Vmod(1,kj)=V(1,kj)*(1+(alpha/Nc)*sum(V(2:nu+2,kj)));
end
Vmod(2:nl1,1:nm1-1)=V(2:nl1,1:nm1-1);

Vmod(1,1:nm1-1)=log(Vmod(1,1:nm1-1));
if length(nn)==0,[vm,sel]=min(Vmod(1,1:nm1-1));
nn=V(2:2+2*nu,sel)';end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [A,B,C,D,K,X0]=ssmodx4(th,T,mod)

%	[A,B,C,D,K,X0] = ssmodx4(PARVAL,T,MS)
%
% copy of ssmodx9.m with a few changes,
% meant to be used as a subroutine of arx4 and rarx4
%
	
[dum,nn]=size(mod);nyy=mod(1,nn);nx=mod(2,nn);
ny=abs(nyy);nu=(nn-nx-2*ny-2)/2;
s=1;
if nyy<0,arx=1;else arx=0;end
na=0;
if arx,
as1=mod(1:nx,1:nx);
sumas=sum4vms(as1');
nr=find(sumas==0&~isnan(sumas));
if isempty(nr),na=nx/ny;else na=(nr(1)-1)/ny;end
end
A=mod(1:nx,1:nx);
for kr=1:nx
	for kc=1:nx
	if isnan(A(kr,kc)), A(kr,kc)=th(s);s=s+1;end
	end
end
B=mod(1:nx,nx+1:nx+nu);
for kr=1:nx
	for kc=1:nu
	if isnan(B(kr,kc)),B(kr,kc)=th(s);s=s+1;end
	end
end
if na==0
C=mod(1:nx,nx+nu+1:nx+nu+ny)';
for kr=1:ny
	for kc=1:nx
	if isnan(C(kr,kc)),C(kr,kc)=th(s);s=s+1;end
	end
end
if nu==0, D=[];else D=mod(1:ny,nx+nu+ny+1:nx+2*nu+ny);end
for kr=1:ny
	for kc=1:nu
	if isnan(D(kr,kc)),D(kr,kc)=th(s);s=s+1;end
	end
end
else
C=A(1:ny,1:nx);
if nu>0,D=B(1:ny,:);else D=[];end
end
K=mod(1:nx,nx+2*nu+ny+1:nx+2*nu+2*ny);
for kr=1:nx
	for kc=1:ny
	if isnan(K(kr,kc)),K(kr,kc)=th(s);s=s+1;end
	end
end

X0=mod(1:nx,nx+2*nu+2*ny+1);
for kr=1:nx
	if isnan(X0(kr)),X0(kr)=th(s);s=s+1;end
end
if T>0 % We shall in this case sample the model with sampling interval T
s = expm([[A [B K]]*T; zeros(nu+ny,nx+nu+ny)]);
A = s(1:nx,1:nx);
B = s(1:nx,nx+1:nx+nu);
K = s(1:nx,nx+nu+1:nx+nu+ny);
end

Contact us