1 2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 3// Copyright (C) ????-2008 - INRIA - François DELEBECQUE 4// 5// Copyright (C) 2012 - 2016 - Scilab Enterprises 6// 7// This file is hereby licensed under the terms of the GNU GPL v2.0, 8// pursuant to article 5.3.4 of the CeCILL v.2.1. 9// This file was originally licensed under the terms of the CeCILL v2.1, 10// and continues to be available under such terms. 11// For more information, see the COPYING file which you should have received 12// along with this program. 13 14function [Xp,dima,dimb,dim]=spantwo(A,B) 15 //Given two matrices A and B with same number of rows, 16 //returns a square matrix Xp (not necessarily orthogonal) 17 //such that : 18 // [*,0] (dim-dimb rows) 19 //Xp*[A,B]=[*,*] (dima+dimb-dim rows) 20 // [0,*] (dim-dima rows) 21 // [0,0] 22 //The dima first columns of inv(Xp) span range(A). 23 //Columns dim-dimb+1 to dima of inv(Xp) span the intesection of 24 //range(A) and range(B). 25 //The dim first columns of inv(Xp) span range(A)+range(B). 26 // Ex: A=[1,0,0,4; 27 // 5,6,7,8; 28 // 0,0,11,12; 29 // 0,0,0,16]; 30 //B=[1,2,0,0]';C=[4,0,0,1]; 31 //Sl=ss2ss(syslin('c',A,B,C),rand(A)); 32 //[no,X]=contr(Sl(2),Sl(3));CO=X(:,1:no); //Controllable part 33 //[uo,Y]=unobs(Sl(2),Sl(4));UO=Y(:,1:uo); //Unobservable part 34 //[Xp,dimc,dimu,dim]=spantwo(CO,UO); //Kalman decomposition 35 //Slcan=ss2ss(Sl,inv(Xp)); 36 37 [X1,dim,dima]=spanplus(A,B);Xp=X1'; 38 B1B2B3=Xp*B;B1B2B3=B1B2B3(1:dim,:); 39 [W,dimb]=rowcomp(B1B2B3);W=W(dim:-1:1,:); 40 W11=W(1:dima,1:dima);W21=W(dima+1:dim,1:dima); 41 if rcond(W11)<1.d-10 then 42 // Which is better? 43 B1B2=B1B2B3(1:dima,:);[W,dimb0]=rowcomp(B1B2);W=W(dima:-1:1,:); 44 [n1,n2]=size(A); 45 Xp=blockdiag(W,eye(n1-dima,n1-dima))*Xp; 46 return; 47 end 48 Q=[eye(dima,dima),zeros(dima,dim-dima); 49 -W21*inv(W11),eye(dim-dima,dim-dima)]; 50 Xp(1:dim,:)=Q*W*Xp(1:dim,:); 51endfunction 52