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


Gerard Van Willigenburg


17 Jan 2011 (Updated )

Optimal reduced-oder LQG output feedback controllers

% DPROTIN: Discrete-time Deterministic Parameter Reduced-Order Time-Invariant iNfinite horizon LQG compensation.
%          Optimal compensation of (P,G,C,V,W) based on (Q,R) and nc<=nx.
%          (P,G,C) deterministic.
%          [f,k,l,sigp,sigs,spr,pt,st,ph,sh,sigc,trps,it,tpst]=dprotin(p,g,c,v,w,q,r,nc,opt,pt,st,ph,sh);
% or
%          [f,k,l,sigp,sigs,spr,pt,st,ph,sh,sigc,trps,it,tpst]=dprotin(func,opt,pt,st,ph,sh);
% Input:
%          func  : string containing the function name of the function that
%                  specifies all problem parameters (see e.g. cofu).
%              q : q or [q mc] mc=cross term criterion
%              v : v or [v me] me=cross covariance
%          opt   : optional array containing algorithm parameter settings.
%          opt(1): convergence tolerance, epsl, default 1e-6. 
%          opt(2): iteration limit, maxtps, default 1e10.
%                  iteration stops if trace(pt+st)>maxtps.
%          opt(3): maximum number of iterations, itm, default 10000. 
%          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-9
%          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(20): stepsize homotopy parameter
%         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     : Optimal compensator.
%          sigp, sigs  : Minimum costs (if sigp=sigs).
%          pt,st,ph,sh : Solution optimal projection equations.
%          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 dproeq, 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]=dprotin(p,g,c,v,w,q,r,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; 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,nc]=feval(func);
    if nargin==8; no=0;
    elseif nargin==9 ; no=max(size(opt));
    elseif nargin==13; no=max(size(opt));
    else error('wrong number of input arguments'); end;
% Check dimensions
  if nx1~=1||nx2~=1; error('nc must be scalar'); end;
  if nc<1||nc>nx; error('nc incompatible with p'); 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=1e10; end
  if opt(3)>0; itm=opt(3); else itm=10000; 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-9; end;
  if opt(9)~=0; cc=1; else cc=0; end;
  if opt(10)~=0; cc=1; b=0; end;
  delta=1; dt=delta;
  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=1; tpsr=1; tpst=[]; trps=[];
  if nargin~=13 && 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; endc=0; osc=0; iti=0; idt=0; ndt=10; dtr=1;
    while endc~=1;

      it=it+1; iti=iti+1;


      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
    sigp=trace(q*pt+(q+lo'*r*lo -2*mc*lo )*ph);
    sigs=trace(v*st+(v+ko *w*ko'-2*me*ko')*sh);
    f=th*fo*tg'; k=th*ko; l=lo*tg';
    pa=[p -g*l ; k*c f]; spr=sperad(pa);

    if nargout>=11
      qa=[q -mc*l ; -l'*mc' l'*r*l ];
      va=[v  me*k';  k *me' k *w*k'];
      pp=dlyap(pa,va) ; sigc=trace(qa*pp);
%      ss=dlyap(pa',qa); sigc=trace(va*ss);
    disp(' '); disp(['  System may not be nc=' num2str(nc) ' compensatable']);
    sigc=Inf; sigs=Inf; sigp=Inf; spr=sperad(p); f=[]; k=[]; l=[];

Contact us