bug 8174 fixed - ss2tf produced bad results if applied to a complex state space... 10/2210/3
Serge Steer [Wed, 6 Oct 2010 09:32:49 +0000 (11:32 +0200)]
Change-Id: I21644f383eebe004379042bd8dd7806c21589124

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

index 64efb2b..89e5087 100644 (file)
@@ -43,10 +43,10 @@ were completely clipped.
 quadrants.
 * Improved H&V positionning of angular labels
 * Polar frame drawn in grey instead of black for a better data visibility.
-* if the basic radius is >= 10^4 or <=10^-4, LaTeX is used to display smart 
+* if the basic radius is >= 10^4 or <=10^-4, LaTeX is used to display smart
 $..10^p$ exponents instead of the D+/-## console style display.
 * Finally, their positioning is still improved in order to avoid their mutual
-overlaying for big values. 
+overlaying for big values.
 
 
 SciNotes:
@@ -104,7 +104,7 @@ Bug Fixes:
 
 * bug 8087 fixed - prettyprint did not handle %inf and %nan.
 
-* bug 8107 fixed - sum(a,1), cumsum(a,1), .. made Scilab 5.0.3-beta-4 crash 
+* bug 8107 fixed - sum(a,1), cumsum(a,1), .. made Scilab 5.0.3-beta-4 crash
                    in some particular context.
 
 * bug 8108 fixed - ATOMS: Modules that start by the letter "n" were not well
@@ -128,6 +128,9 @@ Bug Fixes:
 
 * bug 8164 fixed - Typo in the fileparts help page.
 
+* bug 8174 fixed - ss2tf produced bad results if applied to a complex
+                   state space system.
+
 * bug 8177 fixed - format mode was modified by some macros and not restored.
 
 * bug 8202 fixed - Typo in the the localization.
index 503eb51..4e80e4c 100644 (file)
@@ -1,29 +1,27 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) INRIA - 
-// 
+// Copyright (C)  1985-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 [h,num,den]=ss2tf(sl,rmax)
 // State-space to transfer function.
 // Syntax:
-//   h=ss2tf(sl) 
+//   h=ss2tf(sl)
 //   h=ss2tf(sl,'b')
 //   h=ss2tf(sl,rmax)
 //
 //   sl   : linear system (syslin list)
-//   h    : transfer matrix 
+//   h    : transfer matrix
 //   rmax : optional parameter controlling the conditioning in
 //          block diagonalization method is used (If 'b' is entered
 //          a default value is used)
-//   Method: By default, one uses characteristic polynomial and 
+//   Method: By default, one uses characteristic polynomial and
 //   det(A+Eij)=det(A)+C(i,j) where C is the adjugate matrix of A
-//   Other method used : block-diagonalization (generally 
+//   Other method used : block-diagonalization (generally
 //   this gives a less accurate result).
 //
 //!
@@ -32,98 +30,121 @@ function [h,num,den]=ss2tf(sl,rmax)
     h=sl
     return
   end
-  flag=sl(1);
   if typeof(sl)<>'state-space' then
-    error(msprintf(gettext("%s: Wrong type for input argument #%d: State-space form expected.\n"),"ss2tf",1));
+    error(msprintf(_("%s: Wrong type for input argument #%d: State-space form expected.\n"),"ss2tf",1));
   end
-  if sl(3)==[] then h=sl(5);num=sl(5);den=eye(sl(5));return;end
-  if sl(4)==[] then h=sl(5);num=sl(5);den=eye(sl(5));return;end
-  if size(sl(2),'*')==0 then
-    h=sl(5)
+  //Handle special cases (no input, no output, no state)
+  if sl.B==[] then h=sl.D;num=sl.D;den=eye(sl.D);return;end
+  if sl.C==[] then h=sl.D;num=sl.D;den=eye(sl.D);return;end
+  if size(sl.A,'*')==0 then //no state
+    h=sl.D
     return
-  end 
-  //1
-  domaine=sl(7)
+  end
+
+  //Determine the rational fraction formal variable name
+  domaine=sl.dt
   if type(domaine)==1 then var='z';end
   if domaine=='c' then var='s';end;
   if domaine=='d' then var='z';end;
   if domaine==[] then
-    var='s'; 
-    if type(sl(5))==2 then var=varn(sl(5));end
+    var='s';
+    if type(sl.D)==2 then var=varn(sl.D);end
   end
-  //
+
+  //Determine the algorithm
   [lhs,rhs]=argn(0);
-  meth='p'
+  meth='p';
   if rhs==2 then
-    if type(rmax)==10 then  
-      meth=part(rmax,1),
-      rhs=1
-    else 
-      meth='b'
+    if type(rmax)==10 then
+      meth=part(rmax,1);
+      if and(meth<>['p','b']) then
+        error(msprintf(_( "%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"ss2tf",1,"''p'',''b''"));
+      end
+      rhs=1;
+    else
+      meth='b';
     end
   end
-  //
+
   select meth
-  case 'b' // 
-    a=sl(2)
-    [n1,n1]=size(a)
+  case 'b' // Block diagonalization + Leverrier method
+    a=sl.A;
+    [n1,n1]=size(a);
     z=poly(0,var);
+    //block diagonal decomposition of the state matrix
     if rhs==1 then
-      [a,x,bs]=bdiag(a)
+      [a,x,bs]=bdiag(a);
     else
-      [a,x,bs]=bdiag(a,rmax)
+      [a,x,bs]=bdiag(a,rmax);
     end
     k=1;m=[];v=ones(1,n1);den=1;
-    for n=bs';k1=k:k-1+n;
-      // leverrier
-      h=z*eye(n,n)-a(k1,k1)
-      f=eye(n,n)
+    for n=bs' //loop on blocks
+      k1=k:k-1+n;
+      // Leverrier algorithm
+      h=z*eye(n,n)-a(k1,k1);
+      f=eye(n,n);
       for kl=1:n-1,
-       b=h*f,
-       d=-sum(diag(b))/kl,
-       f=b+eye()*d,
-      end;
-      d=sum(diag(h*f))/n
+        b=h*f;
+        d=-sum(diag(b))/kl;
+        f=b+eye()*d;
+      end
+      d=sum(diag(h*f))/n;
       //
       den=den*d;
-      l=[1:k-1,k+n:n1] ,
-      if l<>[] then v(l)=v(l)*d,end
+      l=[1:k-1,k+n:n1];
+      if l<>[] then v(l)=v(l)*d;end
       m=[m,x(:,k1)*f];
       k=k+n;
-    end;
-    if lhs==3 then h=sl(5),num=sl(4)*m*diag(v)*(x\sl(3));return;end
-    m=sl(4)*m*diag(v)*(x \ sl(3))+sl(5)*den;
-    [m1,n1]=size(m);[m,den]=simp(m,den*ones(m1,n1))
-    h=syslin(domaine,m,den)
+    end
     
-  case 'p' then
-    Den=real(poly(sl(2),var))
+    if lhs==3 then 
+      h=sl.D,
+      num=sl.C*m*diag(v)*(x\sl.B);
+    else
+      m=sl.C*m*diag(v)*(x \ sl.B)+sl.D*den;
+      [m,den]=simp(m,den*ones(m))
+      h=syslin(domaine,m,den)
+    end
+    
+  case 'p' then //Adjugate matrix method
+    Den=poly(sl.A,var) //common denominator
+
     na=degree(Den);den=[];
-    [m,n]=size(sl(5))
-    c=sl(4)
-    for l=1:m
+    [m,n]=size(sl.D)
+    c=sl.C
+    for l=1:m //loop on outputs
       [m,i]=max(abs(c(l,:)));
       if m<>0 then
-       ci=c(l,i)
-       t=eye(na,na)*ci;t(i,:)=[-c(l,1:i-1), 1, -c(l,i+1:na)]
-       al=sl(2)*t;
-       t=eye(na,na)/ci;t(i,:)=[c(l,1:i-1)/ci, 1, c(l,i+1:na)/ci]
-       al=t*al;ai=al(:,i),
-       b=t*sl(3)
-       for k=1:n
-         al(:,i)=ai+b(:,k);
-         [nlk,dlk]=simp(real(poly(al,var)),Den)
-         den(l,k)=dlk;
-         num(l,k)=-(nlk-dlk)*ci
-       end
+        ci=c(l,i)
+        t=eye(na,na)*ci;
+        t(i,:)=[-c(l,1:i-1), 1, -c(l,i+1:na)]
+        al=sl.A*t;
+
+        t=eye(na,na)/ci;
+        t(i,:)=[c(l,1:i-1)/ci, 1, c(l,i+1:na)/ci]
+        al=t*al;ai=al(:,i),
+        b=t*sl.B
+
+        for k=1:n //loop on inputs
+          al(:,i)=ai+b(:,k);
+          [nlk,dlk]=simp(poly(al,var),Den)
+          den(l,k)=dlk;
+          num(l,k)=-(nlk-dlk)*ci;
+        end
+      else
+        num(l,1:n)=0*ones(1,n);
+        den(l,1:n)=ones(1,n);
+      end
+    end
+    if lhs==3 then 
+      h=sl.D;
+    else
+      w=num./den+sl.D;
+      if type(w)==1 then 
+        h=w;   //degenerate case
       else
-       num(l,1:n)=0*ones(1,n);
-       den(l,1:n)=ones(1,n);
+        h=syslin(domaine,w);
       end
     end
-    if lhs==3 then h=sl(5);return;end
-    w=num./den+sl(5);
-    if type(w)==1 then h=w;return;end   //degenerate case
-    h=syslin(domaine,w);
   end
 endfunction
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_8174.dia.ref b/scilab/modules/cacsd/tests/nonreg_tests/bug_8174.dia.ref
new file mode 100644 (file)
index 0000000..85598b6
--- /dev/null
@@ -0,0 +1,25 @@
+// =============================================================================
+// 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 8174 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8174
+//
+// <-- Short Description -->
+// ss2tf produces bad results if applied to a complex state space system
+//A real dynamical system represented in its eigen state space basis
+//(A,B,C) matrices are complex but the system is real
+A = [ 0.75+%i*0.6,0;0, 0.75-%i*0.6];
+B = [ 0.2-%i*0.45; 0.2+%i*0.45];
+C = [-0.04+%i*0.1,-0.04-%i*0.1];
+Sys =syslin('d',A,B,C,0);
+H = ss2tf(Sys);
+if norm(coeff((H.num-(-0.1011+0.074*%z))))>10*%eps then bugmes();quit;end
+if norm(coeff((H.den-(0.9225-1.5*%z+%z^2))))>10*%eps then bugmes();quit;end
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_8174.tst b/scilab/modules/cacsd/tests/nonreg_tests/bug_8174.tst
new file mode 100644 (file)
index 0000000..50d798c
--- /dev/null
@@ -0,0 +1,25 @@
+// =============================================================================
+// 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 8174 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8174
+//
+// <-- Short Description -->
+// ss2tf produces bad results if applied to a complex state space system
+//A real dynamical system represented in its eigen state space basis
+//(A,B,C) matrices are complex but the system is real
+A = [ 0.75+%i*0.6,0;0, 0.75-%i*0.6];
+B = [ 0.2-%i*0.45; 0.2+%i*0.45];
+C = [-0.04+%i*0.1,-0.04-%i*0.1];
+Sys =syslin('d',A,B,C,0);
+H = ss2tf(Sys);
+if norm(coeff((H.num-(-0.1011+0.074*%z))))>10*%eps then pause,end
+if norm(coeff((H.den-(0.9225-1.5*%z+%z^2))))>10*%eps then pause,end