← Back to branch summary

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

« back to all changes in this revision

Viewing changes to macros/auto/canon.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/auto/canon.sci
 
1
function [ac,bc,u,ind]=canon(a,b)
 
2
//[ac,bc,u,ind]=canon(a,b)  gives the canonical controllable form
 
3
//of the pair (a,b).
 
4
//
 
5
//ind    controllability indices,
 
6
//ac,bc  canonical form
 
7
//u      current basis i.e. ac=inv(u)*a*u,bc=inv(u)*b
 
8
//
 
9
//See also : obsv_mat, cont_mat, ctr_gram, contrss
 
10
//Author:F. D. (Inria)
 
11
//!
 
12
//1: block-hessenberg form
 
13
// Copyright INRIA
 
14
[ac,bc,u,ro]=contr(a,b,[1.d-10*norm([a,b],1),1.d-10]);
 
15
//2:zeroing what is th the right of under-diagonal blocks
 
16
[na,ni]=size(b);[l,k]=size(ro);k0=na+1;k1=k0-ro(k);
 
17
for kk=k:-1:2,
 
18
  k2=k1-ro(kk-1);rows=k1:k0-1;cols=k1:na;intsc=k2:k1-1;
 
19
  i=eye(na,na);i(intsc,cols)=-ac(rows,intsc)\ac(rows,cols);
 
20
  im1=2*eye()-i;ac=im1*ac*i,bc=im1*bc;u=u*i;k0=k1;k1=k2;
 
21
end;
 
22
//3: compression of under-the-diagonal blocks
 
23
i=eye(na,na);n=1;m=ro(1)+1;
 
24
for kk=1:k-1,
 
25
   c=n:n+ro(kk)-1;z=ac(m:m+ro(kk+1)-1,c);
 
26
   [x,s,v]=svd(z);i(c,c)=v;n=n+ro(kk);m=m+ro(kk+1);
 
27
end;
 
28
ac=i'*ac*i,bc=i'*bc;u=u*i;
 
29
//4. normalization of blocks
 
30
j=eye(na,na);i=eye(na,na);k0=na+1;k1=k0-ro(k);
 
31
for kk=k:-1:2,
 
32
   k2=k1-ro(kk-1);rows=k1:k0-1;long=k2:k2+k0-k1-1;
 
33
   i=eye(na,na);j=eye(na,na);
 
34
   z=ac(rows,long);j(long,long)=z;
 
35
   i(long,long)=inv(z);
 
36
   ac=j*ac*i,bc=j*bc;u=u*i;k0=k1;k1=k2;
 
37
end;
 
38
// controllability indices
 
39
ind=0*ones(1,ni);[xx,mi]=size(ro);
 
40
for k=1:ni,for kk=1:mi,
 
41
   if ro(kk)>=k then ind(k)=ind(k)+1;end;
 
42
end;end
 
43
//final permutation:
 
44
v=ones(1,na);
 
45
for k=1:ind(1)-1,
 
46
       index=sum(ro(1:k))+1;v(k+1)=index;
 
47
end;
 
48
k0=1;kmin=ind(1)+1;
 
49
for kk=2:ni,
 
50
   numb=ind(kk);
 
51
   kmax=kmin+numb-1;
 
52
   v(kmin:kmax)=v(k0:k0+numb-1)+ones(1,ind(kk));
 
53
   k0=kmin;kmin=kmax+1;
 
54
end;
 
55
ac=ac(v,v),bc=bc(v,:),u=u(:,v);