luenberger.m

控制系统计算机辅助设计——MATLAB语言与应用(源代码）

function T=luenberger(A,B)
n=size(A,1); p=size(B,2); S=[]; sigmas=[]; k=1;
for i=1:p
   for j=0:n-1
      S=[S,A^j*B(:,i)];
      if rank(S)==k, k=k+1;
      else, sigmas(i)=j-1; S=S(:,1:end-1); break; end,
   end
   if k>n, break; end
end
k=k-1; % 如果不是完全可控，则用随机数补足满秩矩阵
if k<n
    while rank(S)~=n, S(:,k+1:n)=floor(10*rand(n,n-k)); end
end
L=inv(S); iT=[];
for i=1:p
    for j=0:sigmas(i)
        iT=[iT; L(i+sum(sigmas(1:i)),:)*A^j];
    end,
end
if k<n, iT(k+1:n,:)=L(k+1:end,:); end  % 不可控时补足满秩矩阵
T=inv(iT);  % 构造变换矩阵