Changeset 6e9d48a in git for Singular/LIB/control.lib

Timestamp:

Apr 3, 2005, 10:31:37 PM (19 years ago)

Author:

Viktor Levandovskyy <levandov@…>

Branches:

(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')

Children:

f52515c83a11c2bf5c2744c7a7775df92a8eb725

Parents:

0582696e2e59bf0b8bc20c14d1424a96c802fec2

Message:

*levandov: iostruct and other improvements by Markus added, smith form corrected


git-svn-id: file:///usr/local/Singular/svn/trunk@7809 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/control.lib (modified) (19 diffs)

Singular/LIB/control.lib

@@ -1 @@
-version="$Id: control.lib,v 1.26 2005-03-17 18:42:48 levandov Exp $";
+version="$Id: control.lib,v 1.27 2005-04-03 20:31:37 plural Exp $";
 category="Applications";
 info="
@@ -13 +13 @@
 
 PROCEDURES:
-control (R);           analysis of controllability-related properties of R (using Ext modules),
-control2 (R);          analysis of controllability-related properties of R (using dimension),
-autonom (R);           analysis of autonomy-related properties of R (using Ext modules),
-autonom2 (R);          analysis of autonomy-related properties of R (using dimension),
-LeftKernel (R);        a left kernel of R,
-RightKernel (R);       a right kernel of R,
-LeftInverse (R);       a left inverse of R,
-RightInverse (R);      a right inverse of R,
-smith (M);             a Smith form of a module M,
-rank (M);              a column rank of M as of matrix,
-genericity (M);        analysis of the genericity of parameters,
-canonize (L);          Groebnerification for modules in the output of control/autonomy procs,
-FindTorsion (R, I);    generators of the submodule of a module R, annihilated by the ideal I.
+control (R);            analysis of controllability-related properties of R (using Ext modules),
+control2 (R);           analysis of controllability-related properties of R (using dimension),
+autonom (R);            analysis of autonomy-related properties of R (using Ext modules),
+autonom2 (R);           analysis of autonomy-related properties of R (using dimension),
+LeftKernel (R);         a left kernel of R,
+RightKernel (R);        a right kernel of R,
+LeftInverse (R);        a left inverse of R,
+RightInverse (R);       a right inverse of R,
+smith (M);              a Smith form of a module M,
+colrank (M);            a column rank of M as of matrix,
+genericity (M);         analysis of the genericity of parameters,
+canonize (L);           Groebnerification for modules in the output of control/autonomy procs,
+iostruct (R);           computes an I/O-structure of behavior given by a module R
+FindTorsion (R, I);     generators of the submodule of a module R, annihilated by the ideal I.
 
 
@@ -60 +61 @@
 // static control_output();      Generating the output for the procedure 'control'
 // static autonom_output();      Generating the output for the procedure 'autonom' and 'autonom2'
+// static extgcd_Our(poly p, poly q)  Computes extgcd of p and q. for versions ealier than 2006 extgcd has a bug and is therefore not used
+// static normalize_Our(matrix N, matrix Q) normalizes the columns of N and divides the columns of Q through the leading coefficients of the columns of N
 
 LIB "homolog.lib"; 
@@ -539 +542 @@
 "
 {
-  if( nrows(R) != rank(transpose(R)) )
+  if( nrows(R) != colrank(transpose(R)) )
   {
     return ("control2 cannot be applied, since R does not have full row rank");
   }  
-  intvec v=option(get);
+  intvec v=Opt_Our();
   module R_std=std(R);
   int d=dim_Our(R_std);
   int NVars=nvars(basering);
   int i=NVars-d;
-  module Ext_1;
-
-  if(i==1)
-    {
-      Ext_1=std(Ext_Our(1,R_std));
-    }
+  module Ext_1=std(Ext_Our(1,R_std));
   matrix T=lift(R,R_std);
   list l=genericity(T);
@@ -572 +570 @@
 };
 //------------------------------------------------------------------------
-proc rank(module M)
-"USAGE: proc rank(M), M a matrix/module
+proc colrank(module M)
+"USAGE: proc colrank(M), M a matrix/module
 PURPOSE: compute the column rank of M as of matrix
 RETURN: int
-NOTE: this procedure uses bareiss-algorithm which might not terminate
+NOTE: this procedure uses bareiss-algorithm which might not terminate in some cases
 "
 {
@@ -592 +590 @@
     [-D(2),D(1),0];
   R=transpose(R);
-  rank(R); 
+  colrank(R); 
 };
   
@@ -602 +600 @@
            NVars:  integer, number of variables in a base ring  
            RC: module, kernel-representation of controllable part of the system
-           R_rank: integer, rank of the representation matrix
+           R_rank: integer, column rank of the representation matrix
 PURPOSE: compute all the autonomy properties of the system which is to be returned in 'autonom' procedure
 RETURN:  list 
@@ -618 +616 @@
                    "kernel representation for controllable part",
                    RC,
-                   "column-rank of the matrix",
+                   "column rank of the matrix",
                    R_rank,
                    DofS,
@@ -692 +690 @@
     {
       RC=LeftKernel(RightKernel(R));
-      R_rank=rank(R);
+      R_rank=colrank(R);
     }
   return( autonom_output(i,NVars,RC,R_rank) );
@@ -731 +729 @@
   {
     RC=LeftKernel(RightKernel(R));
-    R_rank=rank(R);
+    R_rank=colrank(R);
   }
   return(autonom_output(i,NVars,RC,R_rank));
@@ -1143 +1141 @@
   view(CC);
 }
+
+//----------------------------------------------------------------
+
+static proc elementof (int i, intvec v)
+{
+  int b=0;
+  for(int j=1;j<=nrows(v);j++)
+    {
+      if(v[j]==i)
+        {
+          b=1;
+          return (b);
+        }
+    }
+  return (b);
+}
+//-----------------------------------------------------------------
+proc iostruct(module R)
+"USAGE: iostruct( R ); R a module 
+RETURN:  list L with entries: string s, intvec v, module P and module Q
+PURPOSE:  if R is the kernel-representation-matrix of some system, then we output a input-ouput representation Py=Qu of the system, the components that have been chosen as outputs(intvec v) and a comment s
+NOTE:  the procedure uses Bareiss algorithm which might not terminate in some cases
+EXAMPLE:  example iostruct; shows an example
+"
+{
+  list L = bareiss(R);
+  int R_rank = ncols(L[1]);
+  int NCols=ncols(R);
+  intvec v=L[2];
+  int temp;
+  int NRows=nrows(v);
+  int i,j;
+  int b=1;
+  module P;
+  module Q;
+  int n=0;
+  
+  while(b==1)               //sort v through bubblesort
+    {
+      b=0;
+      for(i=1;i<NRows;i++)
+        {
+          if(v[i]>v[i+1])
+          {
+            temp=v[i];
+            v[i]=v[i+1];
+            v[i+1]=temp;
+            b=1;
+          }
+        }
+    }
+  P=R[v];                     //generate P 
+  for(i=1;i<=NCols;i++)       //generate Q
+    {
+      if(elementof(i,v)==1)
+        {
+          i++;
+          continue;
+        }
+      Q=Q,R[i];
+    }
+  Q=simplify(Q,2);
+  string s="The following components have been chosen as outputs: ";
+  return (list(s,v,P,Q));
+}
+example
+{"EXAMPLE:";echo = 2;
+ //Example Antenna
+ ring r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), (c,dp);
+ 
+ module RR;
+ RR = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0],
+   [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta],
+   [0, 0, Dt, -K1, 0, 0, 0, 0, 0],
+   [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta],
+   [0, 0, 0, 0, Dt, -K1, 0, 0, 0],
+   [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta];
+ module R = transpose(RR);
+ view(R);
+ view(iostruct(R));
+};      
+
 //---------------------------------------------------------------
 static proc smdeg(matrix N)  
@@ -1218 +1298 @@
 }
 //---------------------------------------------------------------
-static proc Our_extgcd(poly p, poly q)
+static proc extgcd_Our(poly p, poly q)
 {
   ideal J;   //for extgcd-computations
@@ -1248 +1328 @@
   return(L);
 }
+static proc normalize_Our(matrix N, matrix Q)
+"USAGE: normalize_Our(N,Q), N, Q are two matrices
+PURPOSE: normalizes N and divides the columns of Q through the leading coefficients of the columns of N
+RETURN: normalized matrix N and altered Q(according to the scheme mentioned in purpose). If number of columns of N and Q do not coincide, N and Q are returned unchanged
+NOTE: number of columns of N and Q must coincide. 
+"
+{
+  if(ncols(N) != ncols(Q))
+    {
+      return (N,Q);
+    }   
+  module M = module(N);
+  module S = module(Q);
+  int NCols = ncols(N);
+  number n;
+  for(int i=1;i<=NCols;i++)
+    {
+      n = leadcoef(M[i]);
+      if( n != 0 )
+        {
+          M[i]=M[i]/n;
+          S[i]=S[i]/n;
+        }
+     }
+   N = matrix(M);
+   Q = matrix(S);
+   return (N,Q);
+}               
+        
 //---------------------------------------------------------------
 proc smith( module M )
@@ -1257 +1366 @@
 @*      [3]: a unimodular matrix U,
 @*      [4]: a unimodular matrix V, 
-such that U*M*V=S. An empty list is returned when no Smith Form exists. 
-NOTE: The Smith form only exists over PIDs (principal ideal domains).
+such that U*M*V=S. An warning is returned when no Smith Form exists. 
+NOTE: The Smith form only exists over PIDs (principal ideal domains). Use global ordering for computations!
 "
 {
   if (nvars(basering)>1) //if more than one variable, return empty list
   {
-    list @L;
-    return (@L);
+    string s="The Smith-Form only exists for principal ideal domains";
+    return (s);
   }
   matrix N = matrix(M);         //Typecasting
@@ -1321 +1430 @@
         for(ii=k+1;ii<=n;ii++)
         {
-          N = addrow(N,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
-          N = normalize(N);
+          N = addrow(N,k,-lambda[1,ii]*var(1)^f[1,ii],ii);               
           P = addrow(P,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
+          N,Q=normalize_Our(N,Q);
         }
       }
@@ -1345 +1454 @@
       {
         N=addcol(N,k,-mu[1,ii]*var(1)^g[1,ii],ii);
-        N=normalize(N);
         Q=addcol(Q,k,-mu[1,ii]*var(1)^g[1,ii],ii);
+        N,Q=normalize_Our(N,Q);
       }
     }
     if( (k!=1) && (k<n) && (k<m) )
     {
-      L = Our_extgcd(N[k-1,k-1],N[k,k]);      
+      L = extgcd_Our(N[k-1,k-1],N[k,k]);      
       if ( N[k-1,k-1]!=L[1] )  //means that N[k-1,k-1] is not a divisor of N[k,k]
       {
         N=addrow(N,k-1,L[2],k);
-        normalize(N);
         P=addrow(P,k-1,L[2],k);
+        N,Q=normalize_Our(N,Q); 
         
-        N=addcol(N,k,-L[3],k-1);
-        N=normalize(N);
+        N=addcol(N,k,-L[3],k-1);
         Q=addcol(Q,k,-L[3],k-1);
+        N,Q=normalize_Our(N,Q);
         k=k-2;
       }
@@ -1387 +1496 @@
   matrix N=matrix(M);
   matrix B=L[3]*N*L[4];
-  B=normalize(std(B));
   print(B);
   //------- and yet another example --------------
@@ -1396 +1504 @@
   matrix N2=matrix(M2);
   matrix B2=P[3]*N2*P[4];
-  print(B2);
-  B2=normalize(std(B2));
   print(B2);
 }
@@ -1477 +1583 @@
 "USAGE:  FindTorsion(R, I);   R an ideal/matrix/module, I an ideal
 PURPOSE: computes the Groebner basis of the submodule of R, annihilated by I
-RETURN:  module
+ETURN:  module
 NOTE: especially helpful, when I is the annihilator of the t(R) - the torsion submodule of R. In this case, the result is the explicit presentation of t(R) as
 the submodule of R