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]=abcchk(a,b,c,v,w,q,r);
```% ABCCHK: Check and generate dimensions
%         system (A,B,C) or system (A,B,C,V,W)
%         or LQG problem (A,B,C,V,W,ME,Q,R,MC).
%
%         [nx,ny,nu,q,mc,v,me]=abcchk(a,b,c,v,w,q,r)
%
%         System (A,B,C,V,W,ME): Criterion (Q,R,MC)
%
%         xdot(t) = A x(t) + B u(t) + v(t)
%         y(t)    = C x(t) + w(t)
%
%         E{v(t)}=E{w(t)}=0
%         E{v(t)'v(t)}=V*delta(t), E{w(t)'w(t)}=W*delta(t)
%
%         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]=abcchk(a,b,c,v,w,q,r);

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

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

if nargin>=5

[nx1,nx2]=size(v);
if nx1~=nx | (nx2~=nx & nx2~=nx+ny); error('  v incompatible with a'); end;
if nx2==nx+ny;
if nargout<7; error('  v incompatible with a');
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 a'); end;
if nx2==nx+nu;
if nargout<5; error('  q incompatible with a');
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 b'); 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;
```