Code covered by the BSD License

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

### Gerard Van Willigenburg (view profile)

08 Jun 2008 (Updated )

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

[nx,ny,nu,q,mc,v,me]=pgcchk(p,g,c,v,w,q,r);
```% PGCCHK: Check and generate dimensions
%         system (P,G,C) or system (P,G,C,V,W)
%         or LQG problem (P,G,C,V,W,ME,Q,R,MC).
%
%         [nx,ny,nu,q,mc,v,me]=pgcchk(p,g,c,v,w,q,r);
%
%         System (P,G,C,V,W,ME): Criterion (Q,R,MC)
%
%         x    = P x + G u  + v
%          k+1    k k   k k    k
%
%         y    = C x  + w
%          k      k k    k
%
%         E{v }=E{w )}=0
%            k     k
%         E{v 'v }=V *delta(k-l), E{w 'w }=W *delta(k-l)
%            k  l   k                k  k   k
%
%         MC = Cross term criterion should be included in v(input)=[v me]
%         ME = Cross covariance should be include in q(input)=[q mc]
%
%         q,mc: output if q(input)=[q mc]
%         v,me: output if v(input)=[v me]
%
%         L.G. Van Willigenburg, W.L. De Koning, 28-11-95.
%
function [nx,ny,nu,q,mc,v,me]=pgcchk(p,g,c,v,w,q,r);

if nargin~=3 & nargin~=5 & nargin~=7;
error('3, 5, or 7 input arguments required');
end

[nx,nx1]=size(p);
if nx1~=nx; error('  p must be square'); end;
[nx1,nu]=size(g);
if nx1>0 & nx1~=nx; error('  g incompatible with p'); end;
[ny,nx1]=size(c);
if nx1~=nx; error('  c incompatible with p'); end;

if nargin>=5

[nx1,nx2]=size(v);
if nx1~=nx | (nx2~=nx & nx2~=nx+ny); error('  v incompatible with p'); end;
if nx2==nx+ny;
if nargout<7; error('  v incompatible with p');
else; h=v(1:nx,1:nx); me=v(1:nx,nx+1:nx2); v=h; end;
else; me=zeros(nx,ny); end;
% Check if v is positive semi-definite and symmetric
if norm(v)~=0; ev=eig(v);
if any(ev < -max(ev)*eps) | norm(v'-v,1) > 100*eps*norm(v,1)
warning('  v is not symmetric or positive semi-definite');
end
end
% Check v-me*pinv(w)*me' is positive semi-definite
h=v-me*pinv(w)*me'; ev=eig(h);
if any(ev < -max(ev)*eps)
error('  v-me*pinv(w)*me'' is not positive semi-definite'); return;
end

[nx1,nx2]=size(w);
if nx1~=ny | nx2~=ny; error('  w incompatible with c'); end;
% Check if w is positive semi-definite and symmetric
if norm(w)~=0; ev=eig(w);
if any(ev < -max(ev)*eps) | norm(w'-w,1) > 100*eps*norm(w,1)
warning('  w is not symmetric and positive semi-definite');
end
end

end;
if nargin>=7

[nx1,nx2]=size(q);
if nx1~=nx | (nx2~=nx & nx2~=nx+nu); error('  q incompatible with p'); end;
if nx2==nx+nu;
if nargout<5; error('  q incompatible with p');
else; h=q(1:nx,1:nx); mc=q(1:nx,nx+1:nx2); q=h; end;
else; mc=zeros(nx,nu); end;
% Check if q is positive semi-definite and symmetric
if norm(q)~=0
ev=eig(q);
if any(ev < -max(ev)*eps) | norm(q'-q,1) > 100*eps*norm(q,1)
warning('  q is not symmetric or positive semi-definite');
end
end;
% Check if q-mc*pinv(r)*mc' is positive semi-definite
h=q-mc*pinv(r)*mc'; ev=eig(h);
if any(ev < -max(ev)*eps)
error('  q-mc*pinv(r)*mc'' is not positive semi-definite'); return;
end

[nx1,nx2]=size(r);
if nx1~=nu | nx2~=nu; error('  r incompatible with g'); end;
% Check if r is positive semi-definite and symmetric
if norm(r)~=0
ev=eig(r);
if any(ev < -max(ev)*eps) | norm(r'-r,1) > 100*eps*norm(r,1)
warning('  r is not symmetric or positive semi-definite');
end
end

end;
```