← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to macros/robust/gamatest.sci

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
macros/robust/gamatest.sci
 
1
function [ok,sk,stint]=gamatest(P12,P21,gama)
 
2
// Test if gama is ok. If ok, sk=controller. stint is true (%T) if
 
3
// internal stability is achieved.
 
4
// Copyright INRIA
 
5
write(%io(2),gama);
 
6
ok=%F;sk=[];stint=%F;
 
7
[A,B2,C1,D12]=P12(2:5);
 
8
[A,B1,C2,D21]=P21(2:5);
 
9
[nx,nx]=size(A);
 
10
[nx,nw]=size(B1);[ny,nx]=size(C2);
 
11
Q=C1'*C1;R=D12'*D12;S=C1'*D12;
 
12
V=B1*B1';L=D21*B1';N=D21*D21';
 
13
// Two Riccati
 
14
[K,X,errx]=leqr(P12,V/gama/gama);
 
15
[H,Y,erry]=leqe(P21,Q/gama/gama);
 
16
//tests
 
17
notok0=or([h==[],k==[]]);
 
18
if notok0 then write(%io(2),'fail');return;end
 
19
dom=p12(7);
 
20
// ------------------>    Discrete time case:
 
21
if dom=='d' then
 
22
ok1=and(abs(spec(A+B2*K)) < ones(nx,1));
 
23
ok2=and(abs(spec(A+H*C2)) < ones(nx,1));
 
24
ok3=and(real(spec(gama*gama*eye()-B1'*X*B1)) > zeros(nw,1));
 
25
ok4=and(real(spec(gama*gama*eye()-C1*Y*C1')) > zeros(ny,1));
 
26
ok5=and(real(spec(inv(X)+B2*inv(R)*B2'-V/gama/gama)) > zeros(nx,1));
 
27
ok6=and(real(spec(inv(Y)+C2'*inv(N)*C2-Q/gama/gama)) > zeros(nx,1));
 
28
ok7=mini(real(spec(eye()-Y*X/gama/gama))) > 100*%eps;
 
29
ok=and([ok3,ok4,ok7]);
 
30
if ~ok then write(%io(2),'fail');return;end
 
31
E=eye()-Y*X/gama/gama;
 
32
W=(A-L'*inv(N)*C2)*inv((inv(X)+B2*inv(R)*B2'-V));
 
33
if rcond(E) > 1.d-4 then
 
34
          Z=inv(E);
 
35
          Ak=A+B2*K*Z+W*(S*K*Z+Q)/gama/gama+H*C2;Bk=H;Ck=-K*Z;
 
36
          Sk=syslin('d',Ak,Bk,Ck);
 
37
          stint=and(abs(spec(h_cl(P,size(C2*B2),sk))) < ones(2*nx,1))
 
38
//Sk1=lqg(P,size(c2*B2));
 
39
else
 
40
          Ak=A*E+B2*K+W*(S*K+Q*E)/gama/gama+H*C2*E;Bk=H;Ck=-K;
 
41
          Sk=des2ss(Ak,Bk,Ck,0*Ck*Bk,E);
 
42
          stint=and(abs(spec(h_cl(P,size(C2*B2),sk))) < ones(2*nx,1))
 
43
end
 
44
end
 
45
// ----------------->       Continuous time case
 
46
if dom=='c' then
 
47
ok1=and(real(spec(A+B2*K)) < zeros(nx,1));
 
48
ok2=and(real(spec(A+H*C2)) < zeros(nx,1));
 
49
ok3=mini(real(spec(eye()-Y*X/gama/gama))) > 100*%eps;
 
50
ok=and([ok1,ok2,ok3]);
 
51
if ~ok then write(%io(2),'fail');return;end
 
52
Z=inv(eye()-Y*X/gama/gama);;
 
53
Ak=A+B2*K*Z+Y*(S*K*Z+Q)/gama/gama+H*C2;Bk=H;Ck=-K*Z;
 
54
Sk=syslin('c',Ak,Bk,Ck);
 
55
stint=and(real(spec(h_cl(P,size(C2*B2),sk))) < zeros(2*nx,1))
 
56
//Sk1=lqg(P,size(c2*B2));
 
57
end
 
58
write(%io(2),'ok')
 
59
 
 
60