Global Index (short | long) | Local contents | Local Index (short | long)
[TH,bestchoice,nchoice,failflag] = ...
N4SID Estimates a state-space model using a sub-space method.
TH=N4SID(Z) or [TH,AO]=N4SID(Z,ORDER,NY,AUXORD,DKX,MAXSIZE,TSAMP)
TH: Returned as the estimated state-space model in the THETA format.
No model covariances are given.
Z : The output input data [y u], with y and u as column vectors
For multi-variable systems, Z=[y1 y2 ... yp u1 u2 ... un]
ORDER: The order of the model (Dimension of state vector). If entered
as a vector (e.g. 3:10) information about all these orders will be
given in a plot, Default; ORDER=1:10;
If ORDER is entered as 'best', the default order among 1:10 is
chosen.
NY: The number of outputs in the data matrix. Default NY =1.
AUXORD: An auxiliary order, that is used for the selection of state
variables. Default 1.2*ORDER+3. If AUXORD is entered as a row vector
the best value (min pred error) in this vector will be selected.
DKX: This is a vector defining the structure: DKX =[D,K,X]
D=1 indicates that a direct term from input to output will be
estimated, while D=0 means that a delay from input to output
is postulated.
K=1 indicates that the K-matrix is estimated, while K=0 means that
K will be fixed to zero.
X=1 indicates that the initial state is estimated, X=0 that the
initial state is set to zero.
To define an arbitrary input delay structure NK, where NK(ku) is
the delay from input number ku to any of the outputs, let
DKX=[D,K,X,NK]. NK is thus a row vector of length=no of input
channels. When NK is specified, it overrides the value of D.
Default: DKX = [0, 1, 1]
TRACE: Letting the last given argument be 'trace' gives info to screen
about fit and choice of AUXORD
MAXSIZE: See also AUXVAR.
AO: The chosen value of AUXORD.
The algorithm implements Van Overschee's and De Moor's method for
identification of general multivariable linear systems in state space.
See also CANSTART, PEM.
| This function calls | This function is called by |
|---|---|
function [TH,bestchoice,nchoice,failflag] = ...
n4sid(z,order,l,auxord,dkx,maxsize,Tsamp,refine,arg,trace)
% M. Viberg, 8-13-1992, T. McKelvey, L. Ljung 9-26-1993.
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 3.5 $ $Date: 1997/12/02 03:40:05 $
global XIDplotw
TH=[];bestchoice=[];nchoice=[];failflag=[];
if nargin<1
disp('Usage: TH = N4SID(Z);')
disp(' TH = N4SID(Z,ORDER,NY,AUXORD,DKX,MAXSIZE,TSAMP,''trace'');')
return
end
[Ncap,NN] = size(z);
if NN>Ncap,error('Data should be organized in column vectors.'),end
failflag=0;
trace=0;
if nargin==10,if strcmp(lower(trace),'trace'),trace=1;end,end
if nargin==9,if strcmp(lower(arg),'trace'),trace=1;arg=[];end,end
if nargin==8,if strcmp(lower(refine),'trace'),trace=1;refine=[];end,end
if nargin==7,if strcmp(lower(Tsamp),'trace'),trace=1;Tsamp=[];end,end
if nargin==6,if strcmp(lower(maxsize),'trace'),trace=1;maxsize=[];end,end
if nargin==5,if strcmp(lower(dkx),'trace'),trace=1;dkx=[];end,end
if nargin==4,if strcmp(lower(auxord),'trace'),trace=1;auxord=[];end,end
if nargin==3,if strcmp(lower(l),'trace'),trace=1;l=[];end,end
if nargin==2,if strcmp(lower(order),'trace'),trace=1;order=[];end,end
if nargin < 9, arg=[];end,if isempty(arg),arg='nogui';end
if nargin < 8, refine=[];end,if isempty(refine),refine=1;end
if nargin < 7, Tsamp=[];end,if isempty(Tsamp),Tsamp=1;end
if nargin < 6, maxsize=[];end
if nargin < 5, dkx=[];end,if isempty(dkx),dkx=[0,1,1];end
if nargin < 4, auxord=[];end
if nargin < 3, l=[];end, if isempty(l),l=1;end
if nargin < 2, order=[];end,if isempty(order),order=1:10;end
m=NN-l;
if isstr(order)
if Ncap<19*l+30*(m+1)+1
error('There are too few data points to support this order choice.')
end
def_order=1;order=1:10;
else
def_order=0;
end
if dkx(2),oe='no';else oe='oe';end
if dkx(1),dmat='d';else dmat='no';end
if length(dkx)>3
nk=dkx(4:length(dkx));
if length(nk)<m
error(['If you specify a delay structure in DKX, you must give nk',...
'for each input, i.e. DKX must have length 3+no of inputs.']);
end
if any(nk~=fix(nk))|any(nk<0)
error(['The delays must be non-negative integers.'])
end
if any(nk==0),dmat='d';else dmat='no';end
elseif dkx(1), nk=zeros(1,m);
else nk=ones(1,m);
end
delayu=find(nk>1);
if ~isempty(delayu)
refine=1; % In this case we override the refine value
zt=z(max(nk):Ncap,:);
for ku=delayu(:)'
zt(:,l+ku)=z(max(nk)-nk(ku)+1:Ncap+1-nk(ku),l+ku);
end
z=zt;[Ncap,dum]=size(z);
end
if length(auxord)>1,tryaux=1;else tryaux=0;end
if isempty(trace)
if tryaux
trace=1;
else
trace=0;
end
end
if m==0,timeseries=1;else timeseries=0;end
if timeseries&dkx(2)==0
error('For a time-series model, the K-matrix has to be estimated.')
end
auxord=auxord(:)';
if l>NN
error('You have specified more outputs than the number of data columns.')
end
if length(order)>1,
n=max(order);
ordchoice=1;
if tryaux
error('The chosen ORDER and the chosen AUXORD cannot both be vectors.')
end
else
n=order;
ordchoice=0;
laux=length(auxord);
if laux>1
if strcmp(lower(oe),'oe'),ncheck=n+1;else ncheck=n+2;end
auxord=auxord(find(auxord>=ncheck&auxord<=(Ncap+1)/6));
if laux>length(auxord)&strcmp(arg,'nogui'), disp(['AUXORD changed to ',...
int2str(min(auxord)),':',int2str(max(auxord))]),end
end
end
if isempty(auxord)
auxord=fix(1.2*n)+3;
end
if isempty(maxsize),maxsize=idmsize(Ncap,2*auxord(1)*(l+m));end
y = z(:,1:l);
u = z(:,l+1:l+m);
bestloss=[];
if length(auxord)>1&strcmp(arg,'gui')
handax=get(XIDplotw(11,1),'userdata');handax=handax(3);
axes(handax),cla,xlabel('Auxiliary order'),ylabel('Fit to working data')
set(gca,'xlim',[min(auxord)-1,max(auxord)+1])
set(XIDplotw(11,1),'vis','on')
end
for icount=auxord %
if strcmp(arg,'guichoice')
eval('hh=findobj(XIDplotw(10,1),''label'',menulabel(''&Help''));')
R=get(hh,'userdata');
[nR,cR]=size(R);
l=R(1,1);i=R(1,2);m=R(1,3);Ncap=R(1,4);
j = Ncap+1-2*i;
if m==0,timeseries=1;else timeseries=0;end
R=R(2:nR,:);
n=order;
else
i=icount;
if strcmp(lower(oe),'oe'),ncheck=n+1;else ncheck=n+2;end
j = Ncap+1-2*i;
if j<2*(m+1)*i
itemp = fix((Ncap-l)/(1.5*(l+m)+m));
if itemp<ncheck
error('Too few data points for this choice of orders.')
end
if itemp<i
i = itemp;
disp(['AUXORD has been changed to ' num2str(i)])
j = Ncap+1-2*i;
end
end
if i<ncheck
i=ncheck;
disp(['AUXORD has been changed to ',int2str(i)])
end
nrr=2*i*(l+m);
if nrr*1.2*nrr>maxsize,
maxsize=nrr*1.2*nrr;
disp(['MAXSIZE has been changed to ',int2str(maxsize)]);
end
%%% Hankel matrices:
nele=floor(maxsize/(2*i*(l+m)));
if nele==0,error(['MAXSIZE must be larger than ',int2str(maxsize) '.']),end
nloop=floor((Ncap-2*i)/nele-eps)+1;
R=[];H1=zeros(nrr,nrr+nele+2*i-1);
for kk=1:nloop
jj=[1+(kk-1)*nele:min(Ncap,kk*nele+2*i-1)];
Y = ssssaux('idblockh',y(jj,:)/sqrt(j),2*i);
if ~timeseries
U = ssssaux('idblockh',u(jj,:)/sqrt(j),2*i);
else
U=[];
end
H = [U;Y];
if ~isempty(R)
H1(:,1:min(ncR,nrr)+length(jj)-2*i+1) =...
[H,R(1:min(nrR,2*i*(l+m)),1:min(ncR,2*i*(l+m)))];
else
H1(:,1:length(jj)-2*i+1)=H;ncR=0;
end
R=triu(qr(H1(:,1:min(ncR,nrr)+length(jj)-2*i+1)'))';[nrR,ncR]=size(R);
end %for kk
if ncR<i*(2*m+l)+l
failflag=2;
disp('Error: Too few data points for this choice of orders.')
return
end
end % if strcmp(arg,'guichoice')
%%% Z_i:
if ordchoice&~def_order
if strcmp(arg,'gui')
handax=get(XIDplotw(10,1),'userdata');handax=handax(3);
set(handax,'vis','off');axes(handax),cla
[nR,nC]=size(R);
eval('hh=findobj(XIDplotw(10,1),''label'',menulabel(''&Help''));')
set(hh,'userdata',[[l,i,m,Ncap,zeros(1,nC-4)];R]);
else
hfig=figure;
end
end
if timeseries,
%%% Split R
R41t3 = R(l*(i+1)+1:l*2*i,1:l*(i+1));
R44 = R(l*(i+1)+1:l*2*i,l*(i+1)+1:l*2*i);
R3t41t2 = R(l*i+1:l*2*i,1:l*i);
R3t43t4 = R(l*i+1:l*2*i,l*i+1:l*2*i);
R2t41 = R(l*(i-1)+1:l*2*i,1:l*(i-1));
R2t42t4 = R(l*(i-1)+1:l*2*i,l*(i-1)+1:l*2*i);
%%% Calculate QSVD's
[Um1,Sm1,Xm1,Vm1,Tm1] = ssssaux('gsvd',R41t3',R44');
[U,S,X,V,T] = ssssaux('gsvd',R3t41t2',R3t43t4');
[Up1,Sp1,Xp1,Vp1,Tp1] = ssssaux('gsvd',R2t41',R2t42t4');
if ordchoice
dS=diag(S);
testo=log(diag(S));
ndef=min(n,max(find(testo>(max(testo)+min(testo))/2)));
if ~def_order
[dum,dum,dum,xx,yy]=makebars(1:length(dS),log(dS)');
mdS=floor(min(log(dS)));
zer=find(yy==0);yy(zer)=mdS*ones(size(zer));
if ~strcmp(arg,'gui')
handax=gca;
else
hfig=XIDplotw(10,1);
end
axes(handax)
line(xx,yy,'color','y');%%
ylabel('Log of Singular values');xlabel('Model order')
title('Model singular values vs order')
text(0.97,0.95,'Red: Default Choice','units','norm','fontsize',10,...
'HorizontalAlignment','right');
set(hfig,'vis','on')
patch(xx(1:ndef*5-3),yy(1:ndef*5-3),'y');%%
patch(xx(5*ndef-3:5*ndef),yy(5*ndef-3:5*ndef),'r');%%%
patch(xx(5*ndef:n*5+1),yy(5*ndef:n*5+1),'y');%%
set(handax,'vis','on')
if strcmp(arg,'gui')
set(XIDplotw(10,3),'userdata',[[ndef,1:length(dS)];[dS(ndef),dS']]);
return
end
title('Select model order in Command Window.')
n=input('Select model order:(''Return'' gives default) ');
if isempty(n)
n=ndef;
disp(['Order chosen to ',int2str(n)]);
end
else
n=ndef;
end % if def_order
nchoice=n;
end
%%% Find matrices
A = inv(sqrt(S(1:n,1:n)))*pinv(X(1:l*(i-1),1:n))*Xm1(:,1:n)*...
Sm1(1:n,1:n)*Um1(1:l*i,1:n)'*U(:,1:n)*inv(sqrt(S(1:n,1:n)));
CC = X(:,1:n)*sqrt(S(1:n,1:n));
C = CC(1:l,:);
GG = R(1:l*i,1:l*i)*U(:,1:n)*sqrt(S(1:n,1:n));
G = GG(l*(i-1)+1:l*i,:)';
Lam0 = R(l*i+1:l*(i+1),1:l*(i+1))*R(l*i+1:l*(i+1),1:l*(i+1))';
%% Solve the Riccati equation (8)
iLam = inv(Lam0);
[W,d] = eig([ A'-C'*iLam*G' zeros(n,n);...
-G*iLam*G' eye(n) ], ...
[eye(n) -C'*iLam*C; ...
zeros(n,n) A-G*iLam*C ]);
d = diag(d);
[e,index] = sort(abs(d)); % sort on magnitude of eigenvalues
if (~((e(n) < 1) & (e(n+1)>1)))
disp('Can''t order eigenvalues');
end
% select vectors with eigenvalues inside unit circle
WW = W(:,index(1:n));
%%% State covariance:
P = WW(n+1:2*n,:)*inv(WW(1:n,:));
Q = P-A*P*A';
S = G-A*P*C';
R = Lam0-C*P*C';
oe='not_oe';
B=[];
D=[];
else % i.e if not timeseries
invma= R((1:2*i*m+i*l),(1:2*i*m+i*l));
if rank(invma)<i
disp('WARNING: Your input signal is apparently not persistently exciting.')
disp(' Use other input or lower model order.')
disp(' No model returned.')
TH=[];failflag=1;return
end
L = (R(2*m*i+l*i+1:2*m*i+2*l*i,1:2*i*m+i*l)/invma);
Li1 = L(:,1:m*i);
Li3 = L(:,2*m*i+1:2*m*i+l*i);
%%% Gamma_i:
[UU,SS,VV] = svd([Li1 zeros(l*i,m*i) Li3]*R((1:2*i*m+i*l),(1:2*i*m+i*l)));
if ordchoice
dS=diag(SS);
testo=log(diag(SS));
ndef=min(n,max(find(testo>(max(testo)+min(testo))/2)));
if ~def_order
[dum,dum,dum,xx,yy]=makebars(1:length(dS),log(dS)');
mdS=floor(min(log(dS)));
zer=find(yy==0);yy(zer)=mdS*ones(size(zer));
if ~strcmp(arg,'gui')
handax=gca;
else
hfig=XIDplotw(10,1);
end
axes(handax)
line(xx,yy,'color','y');%%
ylabel('Log of Singular values');xlabel('Model order')
title('Model singular values vs order')
text(0.97,0.95,'Red: Default Choice','units','norm','fontsize',10,...
'HorizontalAlignment','right');
set(hfig,'vis','on')
patch(xx(1:ndef*5-3),yy(1:ndef*5-3),'y');%%
patch(xx(5*ndef-3:5*ndef),yy(5*ndef-3:5*ndef),'r');%%%
patch(xx(5*ndef:n*5+1),yy(5*ndef:n*5+1),'y');%%
set(handax,'vis','on')
if strcmp(arg,'gui')
set(XIDplotw(10,3),'userdata',[[ndef,1:length(dS)];[dS(ndef),dS']]);
return
end
title('Select model order in command window.')
n=input('Select model order:(''Return'' gives default) ');
if isempty(n)
%testo=log(diag(SS));
n=ndef;
disp(['Order chosen to ',int2str(n)]);
end
else
n=ndef;
end % if def_order
nchoice=n;
end
U1 = UU(:,1:n); Sigma1 = SS(1:n,1:n);
%%% LS problem:
BB0 = [U1(1:l*(i-1),:)\R(2*m*i+l*i+l+1:2*m*i+2*l*i,1:2*i*m+(i+1)*l);
R(2*m*i+l*i+1:2*m*i+l*i+l,1:2*m*i+l*(i+1))];
BB = BB0(:,1:2*i*m+i*l);
AA = [U1'*R(2*m*i+l*i+1:2*m*i+2*l*i,1:2*i*m+i*l);
R(m*i+1:2*m*i,1:2*m*i+l*i)];
K = BB/AA;
%%% Extraction of system matrices:
A = K(1:n,1:n);
C = K(n+1:n+l,1:n);
G1 = U1(1:l,:)'; G2 = U1(l+1:l*i,:)'; U2 = U1(1:l*(i-1),:);
% Gamma = obs_mat(A,C,i); G = pinv(Gamma); G1 = G(:,1:l);
% G2 = G(:,l+1:l*i); U2 = Gamma(1:l*(i-1),:)
S1 = [-A*G1 eye(n)-A*G2*U2 ; eye(l)-C*G1 -C*G2*U2];
S2 = [ssssaux('idspech',pinv(U2)-A*G2,l) ; ssssaux('idspech',-C*G2,l)];
S2 = S2*[eye(l) zeros(l,n); zeros((i-2)*l,l) U2(1:l*(i-2),:)];
k = [K(:,n+1:n+m); ssssaux('idblockc',K(1:n,n+m+1:n+m*i),m);...
ssssaux('idblockc',K(n+1:n+l,n+m+1:n+m*i),m)];
DB = [S1 ; S2]\k;
D = DB(1:l,:);
B = DB(l+1:l+n,:);
end
if ~timeseries
res = BB - K*AA;
BBrest=BB0(:,2*i*m+i*l+1:2*i*m+(i+1)*l);
QRS=res*res'+BBrest*BBrest';
Q = QRS(1:n,1:n);
S = QRS(1:n,n+1:n+l);
R = QRS(n+1:n+l,n+1:n+l);
end
errflag=0;
eval('K1=ssssaux(''kric'',A,C,Q,R,S);','errflag=1;')
if errflag,
disp(' No model returned.')
TH=[];failflag=1;return,
end
if ~strcmp(lower(oe),'oe')
K=K1;
else
K=zeros(n,l);
end
THprel = ms2th(real(modstruc(A,B,C,D,K)));Kprel=K1;
if tryaux
e=pe(z,THprel);
Vloss=det(e'*e/Ncap);
if trace,
mess=['auxord = ',int2str(i),' gives a fit of ',num2str(Vloss)];
disp(mess)
end
if strcmp(arg,'gui')
[dum,dum,dum,xx,yy]=makebars(i,Vloss);
patch(xx,yy,'y');drawnow
end %if 'gui'
if isempty(bestloss),
update=1;
elseif Vloss<bestloss
update=1;
else
update=0;
end
if update
TH=THprel;bestloss=Vloss;bestchoice=i;KK=Kprel;
end
else
TH=THprel;bestchoice=i;KK=K1;
end % if tryaux
end % for icount=auxord
if tryaux&strcmp(arg,'gui')
[dum,dum,dum,xx,yy]=makebars(bestchoice,bestloss);
patch(xx,yy,'red')
end
e=pe(z,TH);
if (refine&~timeseries)|dkx(3) % refine the estimate of B and D and x0
K1=KK;
[A,B,C,D,K,X0]=th2ss(TH);
lam=e'*e/Ncap;sqrlam=inv(sqrtm(lam));
[ny,nx]=size(C);
[nx,nu]=size(B);nz=ny+nu;
n=nu*nx+nx;if strcmp(dmat,'d'),n=n+ny*length(find(nk==0));end
% *** Compute the gradient PSI. If N>M do it in portions ***
rowmax=max(n*ny,nx+nz);
M=floor(maxsize/rowmax);
R=zeros(n,n);Fcap=zeros(n,1);
for kc=1:M:Ncap-1
jj=(kc:min(Ncap,kc-1+M));
if jj(length(jj))<Ncap,jjz=[jj,jj(length(jj))+1];else jjz=jj;end
psitemp=zeros(length(jj),ny);
psi=zeros(ny*length(jj),n);
x=ltitr(A-K1*C,K1,z(jjz,1:ny),X0);
%We use the good K even for an OE model
yh=(C*x(1:length(jj),:)')';
e=(z(jj,1:ny)-yh)*sqrlam;
[nxr,nxc]=size(x);X0=x(nxr,:)';
evec=e(:);
kl=1;
for kx=1:nx,for ku=1:nu
dB=zeros(nx,1);dB(kx,1)=1;
if kc==1,dX=zeros(nx,1);else dX=dXk(:,kl);end
psix=ltitr(A-K1*C,dB,z(jjz,ny+ku),dX);[rp,cp]=size(psix);
dXk(:,kl)=psix(rp,:)';
psitemp=(C*psix(1:length(jj),:)')'*sqrlam;
psi(:,kl)=psitemp(:);kl=kl+1;
end,end
if strcmp(dmat,'d'),
for ky=1:ny,for ku=find(nk==0);
psitemp=...
[zeros(length(jjz),ky-1),z(jjz,ny+ku),zeros(length(jjz),ny-ky)]...
*sqrlam;psitemp=psitemp(1:length(jj),:);
psi(:,kl)=psitemp(:);kl=kl+1;
end,end
end
%% x0
for kx=1:nx
if kc==1
x0dum=zeros(nx,1);x0dum(kx,1)=1;
else
x0dum=X00(:,kl);
end
psix=ltitr(A-K1*C,zeros(nx,1),zeros(length(jjz),1),x0dum);
[rp,cp]=size(psix);
X00(:,kl)=psix(rp,:)';
psitemp=(C*psix(1:length(jj),:)')'*sqrlam;
psi(:,kl)=psitemp(:);kl=kl+1;
end
R=R+psi'*psi; Fcap=Fcap+ psi'*evec;
end
% *** Compute the estimate of B and D ***
if Ncap>M, g=R\Fcap; else g=psi\evec;end
B1=reshape(g(1:nx*nu),nu,nx)';
D1=zeros(ny,nu);
if strcmp(dmat,'d'),
nud=length(find(nk==0));
Dtemp=reshape(g(nx*nu+1:n-nx),nud,ny)';
D1(:,find(nk==0))=Dtemp;
B1=B1+K1*D1;
end
if dkx(3),x0=g(n-nx+1:n);else x0=zeros(nx,1);end
if refine, B=B1;D=D1;end
if max(abs(eig(A-K*C)))>1 % Safety check. Theoretically this cannot
% happen
eval('K=ssssaux(''kric'',A,C,K*K'',eye(l),K);','');
end
if ~isempty(delayu)
[A,B,C,D,K,x0]=ssssaux('expand',A,B,C,D,K,x0,nk);
end
TH = ms2th(real(modstruc(A,B,C,D,K,x0)));
e=pe(z,TH);
end % if refine
lambda=e'*e/Ncap;Vloss=det(lambda);
npar=l*n+m*n;if ~strcmp(lower(oe),'oe'),npar=npar+l*n;end
Vfpe=Vloss*(1+npar/Ncap)/(1-npar/Ncap);
TH(1,1)=Vloss;TH(2,1)=Vfpe;
TH(2,7)=50;TH=sett(TH,Tsamp);
rarg=TH(1,6);
TH(4+rarg:3+l+rarg,1:l)=lambda;