function q=car(P,A,B)

% function Q=car(P,A,B) : Q has maximal trace and satisfies P-APA'=BQB'.
%
% P is a state covariance for x_k=Ax_{k-1}+Bu_k
%
% Q     is an input variance which is consistent with the state covariance P,
%       corresponds to a white noise process, and has maximal-trace
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n,m]=size(B); Neps=1e-15;
setlmis([])
Q=lmivar(1,[m 1]);
P0=lmivar(1,[n 1]);

%1st Q>0
lmiterm([-1 1 1 Q],1,1);

%2nd P0-A P0 A' - B Q B' < eps
lmiterm([2 1 1 P0],1,1);
lmiterm([2 1 1 P0],-A,A');
lmiterm([2 1 1 Q],-B,B');
lmiterm([-2 1 1 0],Neps);

%3nd P0-A P0 A' - B Q B' > -eps
lmiterm([-3 1 1 P0],1,1);
lmiterm([-3 1 1 P0],-A,A');
lmiterm([-3 1 1 Q],-B,B');
lmiterm([3 1 1 0],-Neps);

%4th P0<P
lmiterm([-4 1 1 0],P);
lmiterm([4 1 1 P0],1,1);

lmisys=getlmis;
c=mat2dec(lmisys,eye(m,m),zeros(n,n));
options=[1e-13,0,0,0,0];
[copt,xopt]=mincx(lmisys,-c,options);
trcq=-copt;

%qnumbers=mat2dec(lmisys,ones(m,m),zeros(n,n));
q=dec2mat(lmisys,xopt,Q);