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

by

 

17 Jan 2011 (Updated )

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

[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)
% 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);
  else
    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;
  end;
%
% Check dimensions
%
  [nx,ny,nu,q,mc,v,me]=pgcchk(p,g,c,v,w,q,r);
  if nc<1 || nc>nx; error('nc incompatible with p'); end;
  if delta<0; error('delta should be >=0'); end;
%
% Algorithm parameter settings.
%
  opth=zeros(1,19);
  if no~=0; opth(1:no)=opt(1:no); end;
  opt=opth;
  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
  itc=opt(4);
  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);
  else
% 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]);
  end
  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;
             end;
  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;
    else
      pt=zeros(nx); st=pt;
    end
    [ph]=iniphsh(ic,nx,nc);
    [sh]=iniphsh(ic,nx,nc);
  end
  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;

      [ptc,phc,stc,shc,tps,psi1,psi2,ta]=delroeq(ptc,phc,stc,shc,p,g,c,v,w,q,r,mc,me,delta,nc,1,dt,bb,m,rtol);

      [osc,endc,conver,bb,nocon,tpst,b,iti,dt,idt]=rcchk(opt,osc,conver,bb,it,itm,itc,...
      tps,tpst,epsl,maxtps,'DELROTIN','Trace(P+S)',nx,nc,b,a,cc,iti,delta,dt,idt,ndt,dtr);

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

  if nocon~=1;
%
% Determine compensator and costs
%
    lo=pinv(delta*g'*st*g+r)*(g'*st*(delta*p+eye(nx))+mc');
    ko=(me+(delta*p+eye(nx))*pt*c')*pinv(delta*c*pt*c'+w);
    fo=p-g*lo-ko*c;
    [ta,tg,th]=gmhfac(ph,sh,nc,1);
    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;
    end
    if opt(13)~=0
      disp(' '); disp(['  System nc=' num2str(nc) ' compensatable'])            
      disp(['  Two expressions for the minimal cost'])
      disp([sigp,sigs]);
      disp(['  Delta operator stability measure = ' num2str(spr)])
    end
  else
    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']);
    end
  end

Contact us