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

% DPROTIN: Deterministic Parameter Discrete-time 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