A=diag(2*rand(1,10)-ones(1,10));
B=[zeros(3,2);rand(7,2)];
C=[rand(2,6),zeros(2,4)];
D=rand(2,2);                                          % generujemy macierze układu
%P=rand(10,10)+eye(10,10);
R=rand(10,10);R=R+R';[P,DD]=eig(R);    % generujemy losową macierz ortogonalną
A=inv(P)*A*P;B=inv(P)*B;C=C*P;           % stosujemy do macierzy układu transformację podobieństwa, żeby ukryć jej strukturę
u=ss(A,B,C,D)                                        % tworzymy układ
%
S=ctrb(u);
rank(S)
O=obsv(u);
rank(O)                                                 % macierze sterowalności i obserwowalności i ich rzędy
%
XS=orth(S);                                          % znajdujemy baze obrazu S = przestrzeni sterowalnej
XNS=null(S');                                        % i bazę przestrzeni niesterowalnej = dopełnienia ortogonalnego przestrzeni sterowalnej
%XS'*XNS
%rank([XS,XNS])
XO=orth(O');                                         % znajdujemy baze jądra O = przestrzeni nieobserwowalnej
XNO=null(O);                                        % i bazę przestrzeni obserwowalnej = dopełnienia ortogonalnego przestrzeni nieobserwowalnej
%
XSNO=null([XNS,XO]');rank(XSNO)             % wyznaczamy bazę przestrzeni sterowalnej nieobserwowalnej
XSO=null([XNS,XSNO]');rank(XSO)             % wyznaczamy bazę przestrzeni sterowalnej obserwowalnej
XNSNO=null([XO,XSNO]');rank(XNSNO)      % wyznaczamy bazę przestrzeni niesterowalnej nieobserwowalnej
XNSO=null([XS,XNO]');rank(XNSO)             % wyznaczamy bazę przestrzeni niesterowalnej obserwowalnej
%
Q=[XSO,XSNO,XNSO,XNSNO];                  % macierz przejścia do bazy sumy prostej powyższych podprzestrzeni
eig(Q)                                                          % sprawdzamy czy wszystko poszło dobrze i macierz nie jest osobliwa
inv(Q)*A*Q
inv(Q)*B
C*Q                                                             % zmieniamy bazę w przestrzeni stanów, macierz A staje się macierzą blokową a w B i C pojawiają się bloki zer