getf SCI/util/testexamples.sci
reinit_for_test()
%U=mopen('SCI/tests/automatic_tests/lqr_data.ref','rb');
A = rand(2, 2);B = rand(2, 1);//two states, one input
Q = diag([2,5]);R = 2;//Usual notations x'Qx + u'Ru
Big = sysdiag(Q, R);//Now we calculate C1 and D12
[w,wp] = fullrf(Big);C1 = wp(:, 1:2);D12 = wp(:, 3:$);//[C1,D12]'*[C1,D12]=Big
P = syslin('c', A, B, C1, D12);//The plant (continuous-time)
[K,X] = lqr(P);
if load_ref('X') then   pause,end,
if load_ref('K') then   pause,end,

%ans = spec(A + B * K);
if load_ref('%ans') then   pause,end,
//check stability
%ans = norm(A' * X + X * A - X * B * inv(R) * B' * X + Q, 1);
if load_ref('%ans') then   pause,end,
//Riccati check
P = syslin('d', A, B, C1, D12);// Discrete time plant
[K,X] = lqr(P);
if load_ref('X') then   pause,end,
if load_ref('K') then   pause,end,

%ans = spec(A + B * K);
if load_ref('%ans') then   pause,end,
//check stability
%ans = norm(A' * X * A - A' * X * B * pinv(B' * X * B + R) * (B' * X * A) + Q - X, 1);
if load_ref('%ans') then   pause,end,
//Riccati check
xdel_run(winsid());

mclose(%U);