www.pudn.com > Closid30.rar > CLDARE.M, change:1998-05-19,size:3259b


function X = cldare(F,G,R,H,option) 
 
% X = cldare(F,G,R,H,option) 
% 
% Returns the stabilizing solution to the discrete-time Riccati equation 
% that is given by 
%                                 -1 
%   X - F'*X*F + F'*X*G*(R+G'*X*G)  *G'*X*F - H = 0 
% 
% where (F,G) is stabilizable, (C,F) is detectable where C'C = H,  
% corresponding to the state space system 
% 
%  x(k+1) = Fx(k) + Gu(k) 
%    y(k) = Cx(k) 
% 
% METHOD: 
% 'option' can be used to control how the computation is performed: 
% option = 'SCH' Schur method, NO SINGULAR F MATRIX. 
%            See: Laub,A.J., A Schur method for solving Algebraic Ricatti 
%            equations, IEEE AC-24, no.6, pp 913-921, 1979. 
%          'GEV' Generalized eigen vector SINGULAR F MATRIX SIMPLE POLES 
%            See: Pappas,T, A.J.Laub, N.R.Sandell, On the numerical solution 
%            of the discrete-time algebraic riccati equation, IEEE 
%            AC-25, no.4, pp 631-641, 1980. 
 
%************************************************************************** 
% 
% (c) P.M.J. Van den Hof, november 1996. 
% Mechanical Engineering Systems and Control Group 
% Delft University of Technology 
 
% file         : cldare.m 
% info         : This file is based on the routines dare.m and dare2.m  
%                written by P.M.M. Bongers (april 1989), and modified by  
%                R.A. de Callafon, june 1996.  
% functions    : schur, rsf2csf, schord 
% last update  : 19-05-1998, Paul Van den Hof 
 
%**************************************************************************** 
% INPUT CHECKS 
if nargin<5, 
    error('Number of input arguments is not correct'); 
end; 
G=G*inv(R)*G'; 
ind=[]; 
 
[nr,nc] = size(F); 
n = nr; 
if (nr ~= nc), 
    error('Matrix F must be square.') 
end; 
[nr,nc] = size(G); 
if (nr~=n | nc~=n), 
    error('Incorrectly dimensioned G matrix') 
end; 
[nr,nc] = size(H); 
if (nr~=n | nc~=n), 
    error('Incorrectly dimensioned H matrix') 
end; 
 
if (option=='SCH')|(option=='sch'), 
    if (rank(F)~=n), 
        error('Matrix F is not invertible, use option GEV'); 
    end; 
    [q,t] = schur([F+G/F'*H  -G/F'; -F'\H  inv(F)']); 
    [q,t] = rsf2csf(q,t); 
    ns = 0; 
    for i = 1:2*n, 
        if (abs(t(i,i)) < 1-100*eps), 
            index = [ index -1 ]; 
            ns = ns + 1; 
        elseif (abs(t(i,i)) > 1+100*eps), 
            index = [ index 1 ]; 
        else, 
            index = [ index 0 ]; 
        end 
    end; 
    if (ns ~= n), 
        error('Poles on unit circle: (F,G) may be uncontrollable'); 
    end; 
    [q,t] = schord(q,t,index); 
    X = real(q(n+1:n+n,1:n)/q(1:n,1:n)); 
elseif (option=='GEV')|(option=='gev'), 
    L=[eye(n) G;zeros(n) F']; 
    M=[F zeros(n);-H eye(n)]; 
    [ev,ew]=eig(M,L); 
    ew = diag(ew); 
    ind=find(isnan(ew)==1 | ew == inf); 
    if length(ind)~=0, 
        ew(ind)=100; 
    end 
    [ew,index] = sort(abs(ew));    % sort on magnitude of eigenvalues 
    if (~( (ew(n)<1) & (ew(n+1)>1) ) ) 
      error('Sorry: I can not order eigenvalues, (F,G) may be uncontrollable.') %(*) 
    end 
    % select vectors with eigenvalues inside unit circle 
    X=real(ev(n+1:2*n,index(1:n))/ev(1:n,index(1:n))); 
else 
    error('Incorrect option: can only be GEV or SCH'); 
end