Code covered by the BSD License  

Highlights from
Continuous and discrete time optimal reduced order output feedback

Continuous and discrete time optimal reduced order output feedback



17 Jan 2011 (Updated )

Software associated with : International Journal of Control, 83, 12, 2546-2563, 2010

% DELROTIN: Delta Deterministic Parameter Reduced-Order Time-Invariant iNfinite horizon LQG
%           compensation. Optimal compensation of (Pd,Gd,Cd,Vd,Wd) based on (Qd,Rd) and nc<=nx.
%           (Pd,Gd,Cd) deterministic.
%           [f,k,l,sigp,sigs,spr,pt,st,ph,sh,sigc,trps,it,tpst]=delrotin(p,g,c,v,w,q,r,delta,nc,opt,pt,st,ph,sh);
% or
%           [f,k,l,sigp,sigs,spr,pt,st,ph,sh,sigc,trps,it,tpst]=delrotin(func,opt,pt,st,ph,sh);
%           LQG parameters in delta domain: pd,gd,cd,vd,wd,qd,rd
% Input:
%           p,g,c  : p=pd, g=gd, c=cd
%           v,w,q,r: v=delta*vd, w=delta*wd, q=qd/delta, r=rd/delta
%           func  : string containing the function name of the function that
%                   specifies all problem parameters (see e.g. delfu).
%           opt   : optional array containing algorithm parameter settings.
%           opt(1): convergence tolerance, epsl, default 1e-6. 
%           opt(2): iteration limit, maxtps, default 1e20.
%                   iteration stops if trace(pt+st)>maxtps.
%           opt(3): maximum number of iterations, itm, default 5000. 
%           opt(4): plot parameter, itc, default 0.
%                   0    : no plot.
%                   <0   : plot at the end
%                   >0   : plot every itc iterations.
%           opt(5): numerical damping coefficient, 0<b<1, default 0.25.
%                   if opt(5)>1 | opt(5)<0 then b=0
%           opt(6): initial conditions ph, sh, ic=1,2,3, default 1
%                   ic=1: 0.9*eye(nx)+0.1*ones(nx)
%                   ic=2: random
%                   ic=3: eye(nx)
%           opt(7): method to compute projection, m=1,2,3, default 1
%                   m=1: Eigenvalues
%                   m=2: SVD's
%                   m=3: Group inverse
%           opt(8): relative tolerance rank(ph*sh), rtol default 1e-12
%           opt(9): Automatic adjustment numerical damping, cc
%                   cc=0 activated (default) otherwise deactivated
%          opt(10): Convergence of tr(P(i+1)+S(i+1))/tr(P(i)+S(i)) used if non-zero
%          opt(11): Stepsize dt used to iterate the delta projection equations
%                   Default determined by Shannon/De Koning/Van Willigenburg rule
%                   If delta<dt the problem is considered 'almost continuous-time' and instead of delta
%                   dt is used as the 'time-step' to iterate the delta projection equations.
%                   The iterations can then be interpreted as the Euler integration of pt,st,ph,sh.
%                   dt is adjusted downwards if convergence is not obtained.
%          opt(13): If non-zero compensatability messages produced
%          opt(20): stepsize homotopy parameter, default 1
%         opt(20): if 0 then the full-order LQG compensator is computed
%       opt(20:..): individual values homotopy parameter, default none
%      pt,st,ph,sh: initial values pt,st,ph,sh
% Output:
%           f, k, l     : Delta Optimal compensator.
%           sigp, sigs  : Minimum costs (if sigp=sigs).
%           pt,st,ph,sh : Solution optimal projection equations.
%           sigc        : Costs of computed compensator.
%                        (compensator optimal: sigc=sigp=sigs).
%           trps        : Trace(P+S) (array in case of homotopy algorithm)
%           it          : Number of iterations
%           tpst        : Array with trace(P+S) used for convergence detection
% Comments:
%          See also delroeq, gmhfac
%          L.G. Van Willigenburg, W.L. De Koning, 28-11-95.

  function [f,k,l,sigp,sigs,spr,pt,st,ph,sh,sigc,trps,it,tpst]=delrotin(p,g,c,v,w,q,r,delta,nc,opt,pt,st,ph,sh)
% Check number of input arguments and get
% problem data from func=p if p is a string
  if isstr(p);
    if nargin==1; nc=g; no=0;
    elseif nargin==2; opt=g; no=max(size(opt));
    elseif nargin==6; opt=g; no=max(size(opt));
      pt=c; st=v; ph=w; sh=q;
    else error('wrong number of input arguments'); end;
    func=p; [p,g,c,v,w,q,r,delta,nc]=feval(func);
    if nargin==9; no=0;
    elseif nargin==10 ; no=max(size(opt));
    elseif nargin==14; no=max(size(opt));
    else error('wrong number of input arguments'); end;
% Check dimensions
  if nc<1 || nc>nx; error('nc incompatible with p'); end;
  if delta<0; error('delta should be >=0'); end;
% Algorithm parameter settings.
  if no~=0; opth(1:no)=opt(1:no); end;
  if opt(1)>0; epsl=opt(1); else epsl=1e-6; end
  if opt(2)>0; maxtps=opt(2); else maxtps=1e20; end
  if opt(3)>0; itm=opt(3); else itm=5000; end
  if opt(5)==0; b=0.25;
  elseif opt(5)<0 || opt(5)>1; b=0;
  else b=opt(5); end; bb=b;
  if opt(6)>=1 && opt(6)<=3; ic=round(opt(6)); else ic=1; end
  if opt(7)>=1 && opt(7)<=3; m=round(opt(7)); else m=1; end
  if opt(8)>0; rtol=opt(8); else rtol=1e-12; end;
  if opt(9)~=0; cc=1; else cc=0; end;
  if opt(10)~=0; cc=1; b=0; end;
  if opt(11)>0; dt=opt(11);
% Determine stepsize according to Shannon/De Koning/Van Willigenburg
    ep=eig(p); if delta>0; ep=(exp(ep*delta)-1)/delta; end;
    rep=max([abs(real(ep)); 0.1e-12]); iep=max([abs(imag(ep)); 0.03e-12]);
    dt=min([0.1/rep 0.03/iep]);
  if no==20; if opt(20)<0 || opt(20)>1;
               error('0 <= opt(20) <= 1 not satisfied')
             elseif opt(20)~=0; da=1/round(1/opt(20)); av=[0:da:1];
             else av=0; nc=nx;
  elseif no>20; av=sort(opt(20:no));
                av=av-av(1); av=av/av(no-19);
  else av=1; end;
% Initialization recursions.
  it=0; tps=0; tpst=[]; trps=[];
  if nargin~=14 && nargin~=6
    if opt(10)~=0;
      pt=eye(nx); st=pt;
      pt=zeros(nx); st=pt;
  ptc=pt; phc=ph; stc=st; shc=sh;
% pt,ph,st,sh loop
  for a=av;
    conver=0; osc=0; endc=0; iti=0; idt=0; dtr=0; ndt=0;
    while endc~=1;

      it=it+1; iti=iti+1; tpso=tps;



      if nocon==2;
        ptc=pt; phc=ph; stc=st; shc=sh;
    pt=ptc; ph=phc; st=stc; sh=shc;
    trps=[trps tps];

  if nocon~=1;
% Determine compensator and costs
    k=th*ko; l=lo*tg'; f=th*fo*tg';
    sigp=trace(q*pt+(q+lo'*r*lo -2*mc*lo )*ph);
    sigs=trace(v*st+(v+ko *w*ko'-2*me*ko')*sh);
    if delta~=0; sigp=delta*sigp; sigs=delta*sigs; end;
    pa=[p -g*l ; k*c f];
    spr=maxreig(pa,delta); % delta domain stability measure
    %spr=sperad(delta*pa+eye(size(pa))); % discrete-time stability measure

    if nargout>=11 
      qa=[q -mc*l ; -l'*mc' l'*r*l ];
      va=[v  me*k';  k *me' k *w*k'];
      pp=dellyap(pa,va,delta) ; sigc=trace(qa*pp);
%      ss=dellyap(pa',qa,delta); sigc=trace(va*ss);
      if delta~=0; sigc=delta*sigc; end;
    if opt(13)~=0
      disp(' '); disp(['  System nc=' num2str(nc) ' compensatable'])            
      disp(['  Two expressions for the minimal cost'])
      disp(['  Delta operator stability measure = ' num2str(spr)])
    sigc=Inf; sigs=Inf; sigp=Inf; spr=maxreig(p,delta); f=[]; k=[]; l=[];
    if opt(13)~=0
      disp(' '); disp(['  System may not be nc=' num2str(nc) ' compensatable']);

Contact us