getf SCI/util/testexamples.sci
reinit_for_test()
%U=mopen('SCI/tests/automatic_tests/lqe_data.ref','rb');
//Assume the equations
//.
//x = Ax + Ge
//y = Cx + v
//with
//E ee' = Q_e,    Evv' = R,    Eev' = N
//
//This is equivalent to
//.
//x = Ax  + B1 w
//y = C2x + D21 w
//with E { [Ge ]  [Ge v]' } = E { [B1w ] [B1w D21w]' } = bigR =
//         [ v ]                  [D21w]
//
//[B1*B1'  B1*D21';
// D21*B1'  D21*D21']
//=
//[G*Q_e*G' G*N;
// N*G' R]

//To find (B1,D21) given (G,Q_e,R,N) form bigR =[G*Q_e*G' G*N;N'*G' R].
//Then [W,Wt]=fullrf(bigR);  B1=W(1:size(G,1),:);
//D21=W(($+1-size(C2,1)):$,:)
//
//P21=syslin('c',A,B1,C2,D21);
//[K,X]=lqe(P21);

//Example:
nx = 5;ne = 2;ny = 3;
A = -diag(1:nx);G = ones(nx, ne);
C = ones(ny, nx);Q_e(ne, ne) = 1;R = diag(1:ny);N = zeros(ne, ny);
bigR = [G * Q_e * G',G * N;N' * G',R];
[W,Wt] = fullrf(bigR);B1 = W(1:size(G, 1), :);
D21 = W($ + 1 - size(C, 1):$, :);
C2 = C;
P21 = syslin('c', A, B1, C2, D21);
[K,X] = lqe(P21);
//Riccati check:
S = G * N;Q = B1 * B1';
%ans = (A - S * inv(R) * C2) * X + X * ((A - S * inv(R) * C2)') - X * C2' * inv(R) * C2 * X + Q - S * inv(R) * S';
if load_ref('%ans') then   pause,end,


//Stability check:
%ans = spec(A + K * C);
if load_ref('%ans') then   pause,end,

xdel_run(winsid());

mclose(%U);