Code covered by the BSD License  

Highlights from
Reduced-order inf. horizon time-inv. discr.-time LQG control for systems with white parameters

Reduced-order inf. horizon time-inv. discr.-time LQG control for systems with white parameters



08 Jun 2008 (Updated )

Optimal reduced-order compensation of discrete-time linear systems with white parameters

% SPROTIN1: Stochastic 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) stochastic and independent.
%          [f,k,l,sigp,sigs,spr,pt,st,ph,sh,trps]=...
%           sprotin1(pm,gm,cm,pms,gms,cms,pva,gva,cva,...
%             v,w,q,r,nc,opt,pt,st,ph,sh);
% or
%          [f,k,l,sigp,sigs,spr,pt,st,ph,sh,trps]=...
%          sprotin1(func,opt,pt,st,ph,sh);
% Input:
%          func  : string containing the function name of the function that
%                  specifies all LQG problem parameters (see e.g. cospfu).
%              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 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-9
%          opt(9): Automatic adjustment numerical damping, cc
%                  cc=0 activated (default) otherwise deactivated
%         opt(10): stepsize homotopy parameter, default 1
%      opt(10:..): 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).
%          spr         : Spectral radius closed loop system.
%          pt,st,ph,sh : Solution optimal projection equations.
%          trps        : Trace(pt+st) (array in case of homotopy algorithm)
% Comments:
%          See also sproeq1, gmhfac
%          W.L. De Koning, L.G. Van Willigenburg 28-11-95.
  function [f,k,l,sigp,sigs,spr,pt,st,ph,sh,trps]=...

% Check number of input arguments and get
% problem data from func=p if p is a string
  if isstr(pm);
    if nargin==1; no=0;
    elseif nargin==2; opt=gm; no=max(size(opt));
    elseif nargin==6; opt=gm; no=max(size(opt));
      pt=cm; st=pms; ph=gms; sh=cms;
    else error('wrong number of input arguments'); end;
    func=pm; [pm,gm,cm,pms,gms,cms,...
    if nargin==14; no=0;
    elseif nargin==15; no=max(size(opt));
    elseif nargin==19; no=max(size(opt));
    else error('wrong number of input arguments'); end;
% Check dimensions
  nxs=nx*nx; nus=nu*nu; nys=ny*ny;
  if nc<1 | nc>nx; error('nc incompatible with p'); end;
  if nx1~=1 | nx2~=1; error('nc must be scalar'); 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-9; end;
  if opt(9)~=0; cc=1; else; cc=0; end;
  if no==10; if opt(10)<0 | opt(10)>1;
               error('0 < opt(10) < 1 not satisfied')
             else; da=1/round(1/opt(10)); av=[0:da:1];
  elseif no>10; av=sort(opt(10:no));
                av=av-av(1); av=av/av(no-9);
  else; av=1; end;
% Initialization recursions.
  it=0; tps=0; tpst=[]; trps=[];
  if nargin~=6 & nargin~=19
    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;
    while endc~=1;
      it=it+1; iti=iti+1; tpso=tps; %it, keyboard

      if nocon==2;
        ptc=pt; phc=ph; stc=st; shc=sh;

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

  if ~nocon;
% Determine compensator and costs
    gva=gms-kron(gm,gm); cva=cms-kron(cm,cm);

    ln=r; kn=w; stn=st; ptn=pt;
    ln(:)=gms'*st(:)+r(:)+gva'*sh(:); lo=pinv(ln)*gm'*st*pm;
    kn(:)=cms *pt(:)+w(:)+cva *ph(:); ko=pm*pt*cm'*pinv(kn);

    sigs=trace(v*st+(v+ko *w*ko'-2*me*ko')*sh);
    f=th*fo*tg'; k=th*ko; l=lo*tg';

% Computation and display of spectral radius closed loop system
    pams=[pms            -kron(gm*l,pm)    -kron(pm,gm*l)     gms*kron(l,l)
          kron(k*cm,pm)   kron(f,pm)       -kron(k*cm,gm*l)  -kron(f,gm*l)
          kron(pm,k*cm)  -kron(gm*l,k*cm)   kron(pm,f)       -kron(gm*l,f)
          kron(k,k)*cms   kron(f,k*cm)      kron(k*cm,f)      kron(f,f)    ];
%    disp(['         Spectral radius closed loop system=' num2str(spr)]);
    disp(' '); disp(['  System may not be nc=' num2str(nc) ' compensatable']);
    sigc=Inf; sigs=Inf; sigp=Inf; spr=sqrt(sperad(pms)); f=[]; k=[]; l=[];

Contact us