Code covered by the BSD License  

Highlights from
Continuous and discrete time optimal reduced order LQG output feedback

Continuous and discrete time optimal reduced order LQG output feedback



17 Jan 2011 (Updated )

Optimal reduced-oder LQG output feedback controllers

% 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