bug 3796 fix: pb with tf2ss 47/647/3
Serge Steer [Tue, 18 May 2010 15:42:24 +0000 (17:42 +0200)]
Change-Id: I4f63f4ddb56bf300619f5ce1dca60ff93fe4f1f2

scilab/CHANGES_5.3.X
scilab/modules/cacsd/macros/tf2ss.sci
scilab/modules/cacsd/tests/nonreg_tests/bug_3796.dia.ref [new file with mode: 0644]
scilab/modules/cacsd/tests/nonreg_tests/bug_3796.tst [new file with mode: 0644]

index 5593468..2cd3093 100644 (file)
@@ -323,6 +323,10 @@ Bug fixes:
                    file (if filename as no extension or .mat extension) or for an
                    ASCII file if filename has an other extension.
 
+* bug 3796 fixed - In some situation the "tf2ss()" function failed to
+                   compute correctly the state space representation of
+                   a transfer function.
+
 * bug 3811 fixed - Documention about "typeof" and overload prefixes was not
                    up-to-date (help overloading).
 
index c6e79fa..d33e554 100644 (file)
@@ -1,85 +1,84 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) INRIA - 
-// 
+// Copyright (C) 2010 - INRIA - Serge Steer
+//
 // This file must be used under the terms of the CeCILL.
 // This source file is licensed as described in the file COPYING, which
 // you should have received as part of this distribution.  The terms
-// are also available at    
+// are also available at
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 function [sl]=tf2ss(h,tol)
 // Transfer function to state-space.
 //Syntax : sl=tf2ss(h [,tol])
-// h = transfer matrix 
+// h = transfer matrix
 // sl = linear system in state-space representation (syslin list)
 //!
 
-[lhs,rhs]=argn(0)
+  [lhs,rhs]=argn(0)
 
-if type(h)<=2 then 
-  sl=syslin([],[],[],[],h);return;
-end
-[num,den]=h(['num','den']);
-//
-[nd,md]=size(den)
-a=[];b=[];c=[];d=[];n1=0; // s=[]
-for k=1:md
-  for l=1:nd
-    [r,q]=pdiv(num(l,k),den(l,k));
-    dk(l)=q;
-    num(l,k)=r;
+  if type(h)<=2 then
+    sl=syslin([],[],[],[],h);
+    return;
   end
-  if nd<>1 then 
-    [nk,pp]=cmndred(num(:,k),den(:,k));
-  else
-    pp=den(k);nk=num(k);
-  end;
-  
-  slk=cont_frm(nk,pp);//
-  [ak,bk,ck,dk1]=slk(2:5);
-  // [s sk]
-  [n2,m2]=size(bk);
-  if n2<>0 then
-    a(n1+1:n1+n2,n1+1:n1+n2)=ak;
-    b(n1+1:n1+n2,k)=bk;
-    c=[c ck];
-    n1=n1+n2;
-  else
-    if n1<>0 then b(n1,k)=0;end
+  [num,den]=h(['num','den']);
+  //
+  [nd,md]=size(den)
+  a=[];b=[];c=[];d=[];n1=0; // s=[]
+  for k=1:md
+    for l=1:nd
+      [r,q]=pdiv(num(l,k),den(l,k));
+      dk(l)=q;
+      num(l,k)=r;
+    end
+    if nd<>1 then
+      [nk,pp]=cmndred(num(:,k),den(:,k));
+    else
+      pp=den(k);nk=num(k);
+    end;
+
+    slk=cont_frm(nk,pp);
+    [ak,bk,ck,dk1]=slk(2:5);
+    // [s sk]
+    [n2,m2]=size(bk);
+    if n2<>0 then
+      a(n1+1:n1+n2,n1+1:n1+n2)=ak;
+      b(n1+1:n1+n2,k)=bk;
+      c=[c ck];
+      n1=n1+n2;
+    else
+      if n1<>0 then b(n1,k)=0;end
+    end;
+    d=[d dk];
   end;
-  d=[d dk];
-end;
 
-if degree(d)==0 then d=coeff(d),end
+  if degree(d)==0 then d=coeff(d),end
 
-if n1<>0 then
-  nrmb=norm(b,1);nrmc=norm(c,1);fact=sqrt(nrmc*nrmb);
-  b=b*fact/nrmb;c=c*fact/nrmc;atmp=a;
-  [a,u]=balanc(a);
-  if rcond(u)< %eps*n1*n1*100
-  nn=size(a,1);u=eye(nn,nn);a=atmp;
-  end
-  //apply transformation u without matrix inversion
-  [k,l]=find(u<>0) //get the permutation
-  u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:);
-  //c=c*u;b=u^(-1)*b; 
+  if n1<>0 then
+    nrmb=norm(b,1);
+    nrmc=norm(c,1);
+    fact=sqrt(nrmc*nrmb);
+    b=b*fact/nrmb;
+    c=c*fact/nrmc;
+    atmp=a;
+    [a,u]=balanc(a);
+    //next lines commented out to fix bug 3796
+    //   if rcond(u)< %eps*n1*n1*100
+    //     nn=size(a,1);u=eye(nn,nn);a=atmp;
+    //   end
+    
+    //apply transformation u without matrix inversion
+    [k,l]=find(u<>0) //get the permutation
+    u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:);
 
-  if rhs<2 then 
-    [no,u]=contr(a',c',%eps);
+    if rhs<2 then
+      [no,u]=contr(a',c',%eps);
+    else
+      [no,u]=contr(a',c',tol);
+    end
+    u=u(:,1:no);
+    a=u'*a*u;b=u'*b;c=c*u;
+    sl=syslin(h('dt'),a,b,c,d);
   else
-    [no,u]=contr(a',c',tol);
+    sl=syslin(h('dt'),[],[],[],d)
   end
-  u=u(:,1:no);
-  a=u'*a*u;b=u'*b;c=c*u;
-  //[a,u]=balanc(a);
-  
-  //apply transformation u without matrix inversion
-  //[k,l]=find(u<>0) //get the permutation
-  //u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:);
-  //c=c*u;b=u^(-1)*b;
-
-  sl=syslin(h('dt'),a,b,c,d);
-else
-  sl=syslin(h('dt'),[],[],[],d)
-end
 endfunction
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 (file)
index 0000000..38071e1
--- /dev/null
@@ -0,0 +1,43 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+// <-- Non-regression test for bug 3796 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=3796
+//
+// <-- Short Description -->
+//In some situation the "tf2ss()" function is not capable to compute
+//correctly the state space representation of a transfer function.
+s = poly(0,'s');
+E  = 15.0 ;
+Vo =  5.0 ;
+D = Vo/E      ;
+Delta_Vo = 0.1;
+Ts = 1*10^(-3);
+t_on = D*Ts ;
+t_of = Ts-t_on;
+Ro = 5.1 ;
+L  = 2.28*10^(-3);
+Error=0.1;    // ramp error desired
+Delta_i_L=(E*(1-D)*D*Ts)/L;
+C = (Delta_i_L*Ts)/(8*Delta_Vo);
+fdt_Vo_d = syslin('c',E,((L*C)*s^(2)+(L/Ro)*s+1) );
+fa0 = 50;        // [Hertz]
+wa0 = 2*%pi*fa0; // [rad/sec] pulsation
+Kc=1/(E*Error);
+Ctipo=syslin('c',Kc,s);
+K=10^(35/20);
+ma=20;
+wta=0.8;
+Ca=syslin('c',1+(wta/wa0)*s,1+s*(wta/wa0)/ma);
+Cds = K*Ctipo*Ca*Ca;
+Lds = Cds*fdt_Vo_d;
+uno = syslin('c',1,1);
+Wds =Lds/.uno;
+Wds_ss = tf2ss(Wds);
+if 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 (file)
index 0000000..1b4a434
--- /dev/null
@@ -0,0 +1,53 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- JVM NOT MANDATORY -->
+
+// <-- Non-regression test for bug 3796 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=3796
+//
+// <-- Short Description -->
+
+//In some situation the "tf2ss()" function is not capable to compute
+//correctly the state space representation of a transfer function.
+
+s = poly(0,'s');
+
+E  = 15.0 ;
+Vo =  5.0 ;
+
+D = Vo/E      ;
+Delta_Vo = 0.1;
+Ts = 1*10^(-3);
+t_on = D*Ts ;
+t_of = Ts-t_on;
+Ro = 5.1 ;
+L  = 2.28*10^(-3);
+
+Error=0.1;    // ramp error desired
+Delta_i_L=(E*(1-D)*D*Ts)/L;
+C = (Delta_i_L*Ts)/(8*Delta_Vo);
+fdt_Vo_d = syslin('c',E,((L*C)*s^(2)+(L/Ro)*s+1) );
+fa0 = 50;        // [Hertz]
+wa0 = 2*%pi*fa0; // [rad/sec] pulsation
+
+Kc=1/(E*Error);
+Ctipo=syslin('c',Kc,s);
+
+K=10^(35/20);
+ma=20;
+wta=0.8;
+Ca=syslin('c',1+(wta/wa0)*s,1+s*(wta/wa0)/ma);
+Cds = K*Ctipo*Ca*Ca;
+Lds = Cds*fdt_Vo_d;
+uno = syslin('c',1,1);
+Wds =Lds/.uno;
+
+Wds_ss = tf2ss(Wds);
+if size(Wds_ss.A)<>5 then pause,end