summaryrefslogtreecommitdiffstats
path: root/scilab/modules
diff options
context:
space:
mode:
authorSerge Steer <serge.steer@scilab.org>2010-05-18 17:42:24 +0200
committerYann Collette <yann.collette@scilab.org>2010-05-19 13:06:44 +0200
commit6c6315bec3dce92505f7fde38e9d8d597996a83b (patch)
tree67bc1ee6a6997f477876e6056bf02f06ebdb5334 /scilab/modules
parent5d21d422c4f71f9ac9c3886075d3dc61ad95defe (diff)
downloadscilab-6c6315bec3dce92505f7fde38e9d8d597996a83b.zip
scilab-6c6315bec3dce92505f7fde38e9d8d597996a83b.tar.gz
bug 3796 fix: pb with tf2ss
Change-Id: I4f63f4ddb56bf300619f5ce1dca60ff93fe4f1f2
Diffstat (limited to 'scilab/modules')
-rw-r--r--scilab/modules/cacsd/macros/tf2ss.sci127
-rw-r--r--scilab/modules/cacsd/tests/nonreg_tests/bug_3796.dia.ref43
-rw-r--r--scilab/modules/cacsd/tests/nonreg_tests/bug_3796.tst53
3 files changed, 159 insertions, 64 deletions
diff --git a/scilab/modules/cacsd/macros/tf2ss.sci b/scilab/modules/cacsd/macros/tf2ss.sci
index c6e79fa..d33e554 100644
--- a/scilab/modules/cacsd/macros/tf2ss.sci
+++ b/scilab/modules/cacsd/macros/tf2ss.sci
@@ -1,85 +1,84 @@
1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2// Copyright (C) INRIA - 2// Copyright (C) 2010 - INRIA - Serge Steer
3// 3//
4// This file must be used under the terms of the CeCILL. 4// This file must be used under the terms of the CeCILL.
5// This source file is licensed as described in the file COPYING, which 5// This source file is licensed as described in the file COPYING, which
6// you should have received as part of this distribution. The terms 6// you should have received as part of this distribution. The terms
7// are also available at 7// are also available at
8// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt 8// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9 9
10function [sl]=tf2ss(h,tol) 10function [sl]=tf2ss(h,tol)
11// Transfer function to state-space. 11// Transfer function to state-space.
12//Syntax : sl=tf2ss(h [,tol]) 12//Syntax : sl=tf2ss(h [,tol])
13// h = transfer matrix 13// h = transfer matrix
14// sl = linear system in state-space representation (syslin list) 14// sl = linear system in state-space representation (syslin list)
15//! 15//!
16 16
17[lhs,rhs]=argn(0) 17 [lhs,rhs]=argn(0)
18 18
19if type(h)<=2 then 19 if type(h)<=2 then
20 sl=syslin([],[],[],[],h);return; 20 sl=syslin([],[],[],[],h);
21end 21 return;
22[num,den]=h(['num','den']);
23//
24[nd,md]=size(den)
25a=[];b=[];c=[];d=[];n1=0; // s=[]
26for k=1:md
27 for l=1:nd
28 [r,q]=pdiv(num(l,k),den(l,k));
29 dk(l)=q;
30 num(l,k)=r;
31 end 22 end
32 if nd<>1 then 23 [num,den]=h(['num','den']);
33 [nk,pp]=cmndred(num(:,k),den(:,k)); 24 //
34 else 25 [nd,md]=size(den)
35 pp=den(k);nk=num(k); 26 a=[];b=[];c=[];d=[];n1=0; // s=[]
36 end; 27 for k=1:md
37 28 for l=1:nd
38 slk=cont_frm(nk,pp);// 29 [r,q]=pdiv(num(l,k),den(l,k));
39 [ak,bk,ck,dk1]=slk(2:5); 30 dk(l)=q;
40 // [s sk] 31 num(l,k)=r;
41 [n2,m2]=size(bk); 32 end
42 if n2<>0 then 33 if nd<>1 then
43 a(n1+1:n1+n2,n1+1:n1+n2)=ak; 34 [nk,pp]=cmndred(num(:,k),den(:,k));
44 b(n1+1:n1+n2,k)=bk; 35 else
45 c=[c ck]; 36 pp=den(k);nk=num(k);
46 n1=n1+n2; 37 end;
47 else 38
48 if n1<>0 then b(n1,k)=0;end 39 slk=cont_frm(nk,pp);
40 [ak,bk,ck,dk1]=slk(2:5);
41 // [s sk]
42 [n2,m2]=size(bk);
43 if n2<>0 then
44 a(n1+1:n1+n2,n1+1:n1+n2)=ak;
45 b(n1+1:n1+n2,k)=bk;
46 c=[c ck];
47 n1=n1+n2;
48 else
49 if n1<>0 then b(n1,k)=0;end
50 end;
51 d=[d dk];
49 end; 52 end;
50 d=[d dk];
51end;
52 53
53if degree(d)==0 then d=coeff(d),end 54 if degree(d)==0 then d=coeff(d),end
54 55
55if n1<>0 then 56 if n1<>0 then
56 nrmb=norm(b,1);nrmc=norm(c,1);fact=sqrt(nrmc*nrmb); 57 nrmb=norm(b,1);
57 b=b*fact/nrmb;c=c*fact/nrmc;atmp=a; 58 nrmc=norm(c,1);
58 [a,u]=balanc(a); 59 fact=sqrt(nrmc*nrmb);
59 if rcond(u)< %eps*n1*n1*100 60 b=b*fact/nrmb;
60 nn=size(a,1);u=eye(nn,nn);a=atmp; 61 c=c*fact/nrmc;
61 end 62 atmp=a;
62 //apply transformation u without matrix inversion 63 [a,u]=balanc(a);
63 [k,l]=find(u<>0) //get the permutation 64 //next lines commented out to fix bug 3796
64 u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:); 65 // if rcond(u)< %eps*n1*n1*100
65 //c=c*u;b=u^(-1)*b; 66 // nn=size(a,1);u=eye(nn,nn);a=atmp;
67 // end
68
69 //apply transformation u without matrix inversion
70 [k,l]=find(u<>0) //get the permutation
71 u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:);
66 72
67 if rhs<2 then 73 if rhs<2 then
68 [no,u]=contr(a',c',%eps); 74 [no,u]=contr(a',c',%eps);
75 else
76 [no,u]=contr(a',c',tol);
77 end
78 u=u(:,1:no);
79 a=u'*a*u;b=u'*b;c=c*u;
80 sl=syslin(h('dt'),a,b,c,d);
69 else 81 else
70 [no,u]=contr(a',c',tol); 82 sl=syslin(h('dt'),[],[],[],d)
71 end 83 end
72 u=u(:,1:no);
73 a=u'*a*u;b=u'*b;c=c*u;
74 //[a,u]=balanc(a);
75
76 //apply transformation u without matrix inversion
77 //[k,l]=find(u<>0) //get the permutation
78 //u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:);
79 //c=c*u;b=u^(-1)*b;
80
81 sl=syslin(h('dt'),a,b,c,d);
82else
83 sl=syslin(h('dt'),[],[],[],d)
84end
85endfunction 84endfunction
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_3796.dia.ref b/scilab/modules/cacsd/tests/nonreg_tests/bug_3796.dia.ref
new file mode 100644
index 0000000..38071e1
--- /dev/null
+++ b/scilab/modules/cacsd/tests/nonreg_tests/bug_3796.dia.ref
@@ -0,0 +1,43 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2010 - INRIA - Serge Steer
4//
5// This file is distributed under the same license as the Scilab package.
6// =============================================================================
7// <-- JVM NOT MANDATORY -->
8// <-- Non-regression test for bug 3796 -->
9//
10// <-- Bugzilla URL -->
11// http://bugzilla.scilab.org/show_bug.cgi?id=3796
12//
13// <-- Short Description -->
14//In some situation the "tf2ss()" function is not capable to compute
15//correctly the state space representation of a transfer function.
16s = poly(0,'s');
17E = 15.0 ;
18Vo = 5.0 ;
19D = Vo/E ;
20Delta_Vo = 0.1;
21Ts = 1*10^(-3);
22t_on = D*Ts ;
23t_of = Ts-t_on;
24Ro = 5.1 ;
25L = 2.28*10^(-3);
26Error=0.1; // ramp error desired
27Delta_i_L=(E*(1-D)*D*Ts)/L;
28C = (Delta_i_L*Ts)/(8*Delta_Vo);
29fdt_Vo_d = syslin('c',E,((L*C)*s^(2)+(L/Ro)*s+1) );
30fa0 = 50; // [Hertz]
31wa0 = 2*%pi*fa0; // [rad/sec] pulsation
32Kc=1/(E*Error);
33Ctipo=syslin('c',Kc,s);
34K=10^(35/20);
35ma=20;
36wta=0.8;
37Ca=syslin('c',1+(wta/wa0)*s,1+s*(wta/wa0)/ma);
38Cds = K*Ctipo*Ca*Ca;
39Lds = Cds*fdt_Vo_d;
40uno = syslin('c',1,1);
41Wds =Lds/.uno;
42Wds_ss = tf2ss(Wds);
43if size(Wds_ss.A)<>5 then bugmes();quit;end
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_3796.tst b/scilab/modules/cacsd/tests/nonreg_tests/bug_3796.tst
new file mode 100644
index 0000000..1b4a434
--- /dev/null
+++ b/scilab/modules/cacsd/tests/nonreg_tests/bug_3796.tst
@@ -0,0 +1,53 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2010 - INRIA - Serge Steer
4//
5// This file is distributed under the same license as the Scilab package.
6// =============================================================================
7
8// <-- JVM NOT MANDATORY -->
9
10// <-- Non-regression test for bug 3796 -->
11//
12// <-- Bugzilla URL -->
13// http://bugzilla.scilab.org/show_bug.cgi?id=3796
14//
15// <-- Short Description -->
16
17//In some situation the "tf2ss()" function is not capable to compute
18//correctly the state space representation of a transfer function.
19
20s = poly(0,'s');
21
22E = 15.0 ;
23Vo = 5.0 ;
24
25D = Vo/E ;
26Delta_Vo = 0.1;
27Ts = 1*10^(-3);
28t_on = D*Ts ;
29t_of = Ts-t_on;
30Ro = 5.1 ;
31L = 2.28*10^(-3);
32
33Error=0.1; // ramp error desired
34Delta_i_L=(E*(1-D)*D*Ts)/L;
35C = (Delta_i_L*Ts)/(8*Delta_Vo);
36fdt_Vo_d = syslin('c',E,((L*C)*s^(2)+(L/Ro)*s+1) );
37fa0 = 50; // [Hertz]
38wa0 = 2*%pi*fa0; // [rad/sec] pulsation
39
40Kc=1/(E*Error);
41Ctipo=syslin('c',Kc,s);
42
43K=10^(35/20);
44ma=20;
45wta=0.8;
46Ca=syslin('c',1+(wta/wa0)*s,1+s*(wta/wa0)/ma);
47Cds = K*Ctipo*Ca*Ca;
48Lds = Cds*fdt_Vo_d;
49uno = syslin('c',1,1);
50Wds =Lds/.uno;
51
52Wds_ss = tf2ss(Wds);
53if size(Wds_ss.A)<>5 then pause,end