* Bug #8190 fixed - Optimization: Fixed ICSE demos 97/14097/4
Paul Bignier [Fri, 21 Mar 2014 11:54:00 +0000 (12:54 +0100)]
Change-Id: Id05ba33f1cc41b66e02d25221d9d17a1204e87ba

13 files changed:
scilab/CHANGES.md
scilab/modules/optimization/demos/icse/icse.sci
scilab/modules/optimization/demos/icse/icsecontexte.sce [moved from scilab/modules/optimization/demos/icse/icse.contexte with 100% similarity]
scilab/modules/optimization/demos/icse/icsegen.sci
scilab/modules/optimization/demos/icse/icseinit.sce
scilab/modules/optimization/demos/icse/icsenb.f
scilab/modules/optimization/demos/icse/icsez0.f [new file with mode: 0644]
scilab/modules/optimization/demos/icse/icsua.sci
scilab/modules/optimization/demos/icse/icsuq.sci
scilab/modules/optimization/demos/icse/lqv.sce
scilab/modules/optimization/demos/icse/navet.sce
scilab/modules/optimization/demos/icse/seros.sce
scilab/modules/optimization/demos/optimization.dem.gateway.sce

index 2791fe8..1ca0cd8 100644 (file)
@@ -324,6 +324,7 @@ Bug Fixes
 * [#7794](http://bugzilla.scilab.org/show_bug.cgi?id=7794): Example of findABCD help page failed.
 * [#7958](http://bugzilla.scilab.org/show_bug.cgi?id=7958): `mrfit`did not allow a fourth parameter as shown in the help page.
 * [#8010](http://bugzilla.scilab.org/show_bug.cgi?id=8010): Permanent variables could be redefined through a syntax like `%i(1,1)=1`
+* [#8190](http://bugzilla.scilab.org/show_bug.cgi?id=8190): Fixed ICSE demos of Optimization module.
 * [#8356](http://bugzilla.scilab.org/show_bug.cgi?id=8356): `sci2exp` applied to lists, tlists or mlists having undefined fields yielded an error or a wrong result.
 * [#8493](http://bugzilla.scilab.org/show_bug.cgi?id=8493): Some trivial simplifications of `p1./p2` with matrices of complex-encoded polynomials were not done.
 * [#8841](http://bugzilla.scilab.org/show_bug.cgi?id=8841): Error in struct extraction, s.a is not equal to s(1).a
index 1198e04..826a960 100644 (file)
@@ -13,7 +13,7 @@ function [co,u,g,itv,dtv]=icse(u,simu,nap,imp)
     // u(nu)     : initial parameters
     // simu      : string containing the name of the sub program which
     //             describes the problem
-    // nap       : maximum number of call to the simulator
+    // nap       : maximum number of calls to the simulator
     // imp       : debug value during optimization
     // output variables :
     // co        : final cost
index 9bc58fe..3667ee3 100644 (file)
@@ -16,7 +16,7 @@ function [co,u,g,itv,dtv]=icsegen(u,simu,nap,imp,ech,cof)
     // u(nu)     : initial parameters
     // simu      : string which contains the name of the sub program which
     //             describes the problem (second member, criterion and initial state)
-    // nap       : maximum number of call to the simulator
+    // nap       : maximum number of calls to the simulator
     // imp       : debug value during optimization
     // ech(1,nu) : scaling coeff of the control
     // cof(1,ntob*nob) : weighting coeff of the observations
@@ -34,12 +34,19 @@ function [co,u,g,itv,dtv]=icsegen(u,simu,nap,imp,ech,cof)
         alg="gc";
     end
 
-    itv       = itu;
-    itv(nitv) = 0;
-    dtv       = [dtu,y0,tob,matrix(obs,1,ny*nob),don,ech,cof,b,fy1,fu1];
-    dtv(ndtv) = 0;
-    debug(imp);
-    [co,u,g,itv,dtv] = optim(simu,"b",binf,bsup,u, alg, df0, ...
+    nitv  = nitu+ny+ntob;
+
+    ldtvt = ny*(nob+nuc+nuv+nti+ntf+ny+3)+ntob*(nob*nex+nob+1)+2*nu+ndtu;
+    ndtv1 = ldtvt+ny*(ny+7);
+    nui   = iu(1)*nuc+iu(2)*nuv*(nti+ntf+1);
+    ndtv2 = ldtvt+ny*(2*ntob+nu+2*ny+nuc+nuv+4)+max(nuc+nuv,nui)+nob*ntob+nob;
+    ndtv  = max(ndtv1,ndtv2);
+
+    itv         = itu;
+    itv($:nitv) = 0;
+    dtv         = [dtu,y0,tob,matrix(obs,1,ny*nob),don,ech,cof,b,fy(1,:),fu(1,:)];
+    dtv($:ndtv) = 0;
+
+    [co,u,g,itv,dtv] = optim(simu,"b",binf,bsup, u, alg, df0, ...
     "ar",nap, "ti",itv,"td",dtv,"si","sd");
-    debug(0);
 endfunction
index 63c1210..4c0275b 100644 (file)
@@ -19,14 +19,14 @@ if exists("b")==0 then
     b = ones(1,ny);
 end
 if exists("fy")==0 then
-    fy1 = ones(1,ny*ny);
+    fy = ones(1,ny*ny);
 else
-    fy1 = matrix(fy,1,ny*ny);
+    fy = fy(:)';
 end
 if exists("fu")==0 then
-    fu1 = ones(1,ny*(nuc+nuv));
+    fu = ones(1,ny*(nuc+nuv));
 else
-    fu1 = matrix(fu,1,ny*(nuc+nuv));
+    fu = fu(:)';
 end
 
 format("e");
@@ -48,7 +48,7 @@ end
 clear xx yy;
 
 //    initialisation du common icsez
-[nitv,nrtv,ndtv]=fort("icse0",nu,1,"i",t0,2,"d",tf,3,"d",dti,4,"d",..
+[nitv,nrtv,ndtv]=call("icsez0",nu,1,"i",t0,2,"d",tf,3,"d",dti,4,"d",..
 dtf,5,"d",ermx,6,"d",iu,7,"i",nuc,8,"i",nuv,9,"i",ilin,10,"i",nti,..
 11,"i",ntf,12,"i",ny,13,"i",nea,14,"i",itmx,15,"i",nex,16,"i",nob,..
 17,"i",ntob,18,"i",ntobi,19,"i",nitu,20,"i",ndtu,21,"i","sort",..
index 87506e3..7aaeaf3 100644 (file)
@@ -102,7 +102,7 @@ c!
      &d(nob),itu(nitu),dtu(ndtu),iu(5)
 c
       call navetc(indc,nu,tob,obs,cof,ytob,ob,u,
-     &     c,cy,g,yob,d,itu,dtu)
+     &     c,cy,g,yob,d,itu,dtu,ntob,nob,ny,nex,nitu)
 c
       end
       subroutine navetb(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
@@ -256,7 +256,7 @@ C
       end
 
       subroutine navetc(indc,nu,tob,obs,cof,ytob,ob,u,
-     & c,cy,g,yob,d,itu,dtu)
+     & c,cy,g,yob,d,itu,dtu,ntob,nob,ny,nex,nitu)
 c
 c     test probleme navette
 c
@@ -278,7 +278,11 @@ c
 c     variables internes:     yob,d
 c
       implicit double precision (a-h,o-z)
-      dimension ytob(4),ob(3),u(nu),cy(4),g(nu),dtu(*)
+      dimension ytob(4),ob(3),u(nu),cy(4),g(nu),dtu(*),tob(ntob),
+     &obs(nob,ny),cof(nob,ntob),yob(nob,ntob),d(nob),itu(nitu)
+      !dimension tob(ntob),obs(nob,ny),cof(nob,ntob),ytob(ny,ntob),
+      !&ob(nex,ntob,nob),u(nu),cy(ny,ntob),g(nu),yob(nob,ntob),
+      !&d(nob),itu(nitu),dtu(ndtu),iu(5)
 c
       cpen=dtu(25)
       cy(1)=ytob(1)-ob(1)
diff --git a/scilab/modules/optimization/demos/icse/icsez0.f b/scilab/modules/optimization/demos/icse/icsez0.f
new file mode 100644 (file)
index 0000000..b1f7231
--- /dev/null
@@ -0,0 +1,58 @@
+c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+c Copyright (C) INRIA
+c
+c This file must be used under the terms of the CeCILL.
+c This source file is licensed as described in the file COPYING, which
+c you should have received as part of this distribution.  The terms
+c are also available at
+c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+c
+      subroutine icsez0(nu,t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,
+     & ny,nea,itmx,nex,nob,ntob,ntobi,nitu,ndtu,nitv,nrtv,ndtv)
+c
+c     programme d'initialisation appele par icse.bas
+c     initialisation des commons icsez icsez0 et nird
+c
+      implicit double precision (a-h,o-z)
+      dimension iu(5),iu0(5)
+      common/icsez/t00,tf0,dti0,dtf0,ermx0,iu0,nuc0,nuv0,ilin0,nti0,
+     & ntf0,ny0,nea0,itmx0,nex0,nob0,ntob0,ntobi0,nitu0,ndtu0
+      common/nird/nitv0,nrtv0,ndtv0
+c
+      t00=t0
+      tf0=tf
+      dti0=dti
+      dtf0=dtf
+      ermx0=ermx
+      do 10 i=1,5
+10    iu0(i)=iu(i)
+      nuc0=nuc
+      nuv0=nuv
+      ilin0=ilin
+      nti0=nti
+      ntf0=ntf
+      ny0=ny
+      nea0=nea
+      itmx0=itmx
+      nex0=nex
+      nob0=nob
+      ntob0=ntob
+      ntobi0=ntobi
+      nitu0=nitu
+      ndtu0=ndtu
+      nitv0=nitu+ny+ntob
+      nrtv0=0
+      ldtvt=ny*(nob+nuc+nuv+nti+ntf+ny+3)+ntob*(nob*nex+nob+1)+2*nu+
+     &ndtu
+      mdtv1=ldtvt+ny*(ny+7)
+      nui=iu(1)*nuc+iu(2)*nuv*(nti+ntf+1)
+      mdtv2=ldtvt+ny*(2*ntob+nu+2*ny+nuc+nuv+4)+max(nuc+nuv,nui)+
+     &nob*ntob+nob
+      ndtv0=max(mdtv1,mdtv2)
+      !ind=0
+      !call icse(ind,nu,zz,zz,zz,zz,zz,zz,zz,zz,zz)
+      nitv=max(1,nitv0)
+      nrtv=max(1,nrtv0)
+      ndtv=max(1,ndtv0)
+      return
+      end
index a5db9db..1f22db3 100644 (file)
@@ -44,7 +44,7 @@ function [co,u,g,itv,dtv,cof]=icsua(u,simu,nap,imp)
     ico  = 1;
     yob  = 0.d0*ones(nob,ntob);
     ob   = don;
-    [cof] = fort("icscof",ico,1,"i",ntob,2,"i",nex,3,"i",...
+    [cof] = call("icscof",ico,1,"i",ntob,2,"i",nex,3,"i",...
     nob,4,"i",yob,5,"d",ob,6,"d","sort",[1,nob*ntob],7,"d");
     [co,u,g,itv,dtv] = icsegen(u,simu,nap,imp);
     u = ech.*u;
index 590420a..31035c7 100644 (file)
@@ -51,7 +51,7 @@ function [co,u,g,itv,dtv,cof]=icsuq(u,simu,nap,imp,obs,ytob)
     yob  = obs*ytob;
     ob   = don;
 
-    [cof] = fort("icscof",ico,1,"i",ntob,2,"i",nex,3,"i",...
+    [cof] = call("icscof",ico,1,"i",ntob,2,"i",nex,3,"i",...
     nob,4,"i",yob,5,"d",ob,6,"d","sort",[1,nob*ntob],7,"d");
 
     [co,u,g,itv,dtv] = icsegen(u,simu,nap,imp);
index 3e03c5c..098ee19 100644 (file)
 //  *************************************************************
 //
 
-exec("SCI/modules/optimization/demos/icse/icse.contexte");  // context
-
-t0   = 0;     // instant initial
-tf   = 20;    // instant final
-dti  = 1;     // premier pas de temps
-dtf  = 1;     // second pas de temps
-ermx = 1.d-6; // test d'arret absolu sur la valeur du second membre dans
-// la resolution de l'etat
-iu = [0,0,1]; //  iu   :indications sur la structure du controle
-//    iu(1)=1 si l'etat initial depend du controle constant,0 sinon
-//    iu(2)=1 si l'etat initial depend du controle variable,0 sinon
-//    iu(3)=1 si le second membre depend du controle constant,0 sinon
-nuc  = 5;     // nombre de parametres independants du temps
-nuv  = 1;     // nombre de parametres dependants du temps
-ilin = 2;     // indicateur de linearite :
-// 0 pour un systeme non affine
-// 1 pour un systeme affine dont la partie lineaire n'est pas autonome
-// ilin=2 pour un systeme affine dont la partie lineaire est autonome
-nti = 10;     // nombre de pas de temps correspondant a dti (premier pas de temps)
-ntf = 10;     // nombre de pas de temps correspondant a dtf (second pas de temps)
-// si l'on utilise un seul pas de temps,on doit prendre ntf=0
-ny   = 4;     // dimension de l'etat a un instant donne
-nea  = 0;     // nombre d'equations algebriques (eventuellement nul)
-itmx = 10;    // nombre maximal d'iterations dans la resolution de
-// l'equation d'etat discrete a un pas de temps donne
-nex = 1;      // nombre d'experiences effectuees
-nob = 2;      // dimension du vecteur des mesures pour une experience donnee
-// en un instant donne
-ntob  = 10;   // nombre d'instants de mesure pour une experience donnee
-ntobi = 5;    // nombre d'instants de mesure correspondant a dti (premier
-// pas de temps)
-
-nu = nuc + nuv * (nti + ntf + 1); // dimension du vecteur des parametres de controle
-
-//  uc(1,nuc)          :controle constant
-uc = 0*ones(1,nuc);
-
-//  uv(1,nuv*(nti+ntf+1)):controle variable
-if nuv>0 then uv(1,nuv*(nti+ntf+1))=0; end;
-
-//  itu(1,nitu)        :tableau de travail entier reserve a
-//                      l'utilisateur
-itu = [0];
-
-//  dtu(1,ndtu)        :tableau de travail double precision reserve
-//                      a l'utilisateur
-dtu = [0];
-
-//  y0(ny)             :etat initial
-//          (valeur arbitraire si iu(1) ou iu(2) est non nul)
-y0 = ones(1,ny);
-
-//  tob(1,ntob)        :instants de mesure (compatibilite avec ntob
-//                      et ntobi)
-tob  = 2*(1:10);
-binf = -10*ones(1,nu); // borne inf des parametres
-bsup =     ones(1,nu); // borne sup des parametres
-
-//  termes utiles pour une dynamique lineaire ou une observation quadratique
-b(1,ny) = 0;                    // terme constant d'une dynamique lineaire
-fy      = 0.1*ones(ny,ny);      // derivee de la dynamique par rapport a l'etat
-fu       =    ones(ny,nuc+nuv); // derivee de la dynamique par rapport au controle
-
-obs(nob,ny) = 0; // matrice d'observation
-obs         = ones(nob,ny);
-
-don=0*ones(1,nex*ntob*nob);
-
-nap   = 20;  // nombre d'appels du simulateur
-imp   = 2;   // niveau de debug pour optim
-large = 100; // taille de nu au dela de laquelle on choisit un optimiseur
-// pour les problemes de grande taille (alg='gc' dans l'appel de optim)
-
-exec("SCI/modules/optimization/demos/icse/icse.contexte");  // context
-
-[co,u,g,itv,dtv] = icse(u,"icsemc",nap,imp);
+function demo_lqv()
+
+    cd "SCI/modules/optimization/demos/icse/"
+
+    libn  = ilib_for_link("icsez0","icsez0.f",[],"f")
+    nlink = link("./"+libn,"icsez0","f")
+
+    exec("icsecontexte.sce");  // context
+
+    t0   = 0;     // instant initial
+    tf   = 20;    // instant final
+    dti  = 1;     // premier pas de temps
+    dtf  = 1;     // second pas de temps
+    ermx = 1.d-6; // test d'arret absolu sur la valeur du second membre dans la resolution de l'etat
+    iu = [0,0,1]; //  iu   :indications sur la structure du controle
+    //    iu(1)=1 si l'etat initial depend du controle constant,0 sinon
+    //    iu(2)=1 si l'etat initial depend du controle variable,0 sinon
+    //    iu(3)=1 si le second membre depend du controle constant,0 sinon
+    nuc  = 5;     // nombre de parametres independants du temps
+    nuv  = 1;     // nombre de parametres dependants du temps
+    ilin = 2;     // indicateur de linearite :
+    // 0 pour un systeme non affine
+    // 1 pour un systeme affine dont la partie lineaire n'est pas autonome
+    // 2 pour un systeme affine dont la partie lineaire est autonome
+    nti = 10;     // nombre de pas de temps correspondant a dti (premier pas de temps)
+    ntf = 10;     // nombre de pas de temps correspondant a dtf (second pas de temps)
+    // si l'on utilise un seul pas de temps, on doit prendre ntf=0
+    ny   = 4;     // dimension de l'etat a un instant donne
+    nea  = 0;     // nombre d'equations algebriques (eventuellement nul)
+    itmx = 10;    // nombre maximal d'iterations dans la resolution de l'equation d'etat discrete a un pas de temps donne
+    nex = 1;      // nombre d'experiences effectuees
+    nob = 2;      // dimension du vecteur des mesures pour une experience donnee en un instant donne
+    ntob  = 10;   // nombre d'instants de mesure pour une experience donnee
+    ntobi = 5;    // nombre d'instants de mesure correspondant a dti (premier pas de temps)
+
+    nu = nuc + nuv * (nti + ntf + 1); // dimension du vecteur des parametres de controle
+
+    //  uc(1,nuc)          :controle constant
+    uc = 100*ones(1,nuc);
+
+    //  uv(1,nuv*(nti+ntf+1)):controle variable
+    if nuv>0 then uv(1,nuv*(nti+ntf+1))=0; end;
+
+    //  itu(1,nitu)        :tableau de travail entier reserve a l'utilisateur
+    itu = [0];
+
+    //  dtu(1,ndtu)        :tableau de travail double precision reserve a l'utilisateur
+    dtu = [0];
+
+    //  y0(ny)             :etat initial (valeur arbitraire si iu(1) ou iu(2) est non nul)
+    y0 = ones(1,ny);
+
+    //  tob(1,ntob)        :instants de mesure (compatibilite avec ntob et ntobi)
+    tob  = 2*(1:10);
+    binf = -10*ones(1,nu); // borne inf des parametres
+    bsup =     ones(1,nu); // borne sup des parametres
+
+    //  termes utiles pour une dynamique lineaire ou une observation quadratique
+    b(1,ny) = 0;                    // terme constant d'une dynamique lineaire
+    fy      = 0.1*ones(ny,ny);      // derivee de la dynamique par rapport a l'etat
+    fu       =    ones(ny,nuc+nuv); // derivee de la dynamique par rapport au controle
+
+    obs(nob,ny) = 0; // matrice d'observation
+    obs         = ones(nob,ny);
+
+    don=0*ones(1,nex*ntob*nob);
+
+    nap   = 20;  // nombre d'appels du simulateur
+    imp   = 2;   // niveau de debug pour optim
+    large = 100; // taille de nu au dela de laquelle on choisit un optimiseur
+    // pour les problemes de grande taille (alg='gc' dans l'appel de optim)
+
+    saveFormat = format();
+
+    exec("icseinit.sce");
+
+    exec("icse.sci");
+    exec("icsegen.sci");
+
+    [co, u, g, itv, dtv] = icse(u, "icsemc", nap, imp);
+
+    disp(u', "Best value:")
+    disp(co, "Final cost:")
+
+    deletefile("libicsez0.so");
+    deletefile("loader.sce");
+    deletefile("cleaner.sce");
+
+    format(saveFormat(1, $:-1:1));
 
 endfunction
 
+demo_lqv();
+clear demo_lqv;
+
index 4de0f38..b22e60a 100644 (file)
 //  calcul trajectoire optimale de rentree d'une navette spatiale
 //  *************************************************************
 
-libn  = ilib_for_link("icsenb","icsenb.o",[],"f")
-nlink = link("./"+libn,["icsenb","icsenf"],"f")
-
-exec("icse.contexte");
-
-t0   = 0.d0;  // instant initial
-tf   = 1.d0;  // instant final
-dtf  = 0;
-ermx = 1.d-6; // test d'arret absolu sur la valeur du second membre dans
-// la resolution de l'etat
-iu = [0,0,1]; //  iu   :indications sur la structure du controle
-//    iu(1)=1 si l'etat initial depend du controle constant,0 sinon
-//    iu(2)=1 si l'etat initial depend du controle variable,0 sinon
-//    iu(3)=1 si le second membre depend du controle constant,0 sinon
-nuc  = 1;     // nombre de parametres independants du temps
-nuv  = 1;     // nombre de parametres dependants du temps
-ilin = 0;     // indicateur de linearite :
-// 0 pour un systeme non affine
-// 1 pour un systeme affine dont la partie lineaire n'est pas autonome
-// ilin=2 pour un systeme affine dont la partie lineaire est autonome
-nti = 150;    //nombre de pas de temps correspondant a dti (premier pas de temps)
-dti = tf/nti;
-ntf = 00;     // nombre de pas de temps correspondant a dtf (second pas de temps)
-// si l'on utilise un seul pas de temps,on doit prendre ntf=0
-ny   = 4;     // dimension de l'etat a un instant donne
-nea  = 0;     // nombre d'equations algebriques (eventuellement nul)
-itmx = 10;    // nombre maximal d'iterations dans la resolution de
-// l'equation d'etat discrete a un pas de temps donne
-nex = 1;      // nombre d'experiences effectuees
-nob = 3;      // dimension du vecteur des mesures pour une experience donnee
-// en un instant donne
-ntob  = 1;    // nombre d'instants de mesure pour une experience donnee
-ntobi = 1;    // nombre d'instants de mesure correspondant a dti (premier
-// pas de temps)
-
-nu=nuc+nuv*(nti+ntf+1); // dimension du vecteur des parametres de controle
-
-//  uc(1,nuc)          :controle constant
-echtf = 2000;
-uc    = [2500/echtf];
-
-//  uv(1,nuv*(nti+ntf)):controle variable
-//if nuv>0 then uv(1,nuv*(nti+ntf+1))=0; end;
-alpha0 = .20704/.029244; legu = "alpha initial : ann. cz";     // annulation cz
-alpha0 = 17.391;         legu = "alpha initial : finesse max"; // finesse maximum
-legu = " navette americaine ."+legu;
-if nuv>0 then uv=alpha0*ones(1,nuv*(nti+ntf+1)); end;
-
-//  itu(1,nitu)        :tableau de travail entier reserve a
-//                      l'utilisateur
-itu = [0];
-
-//  dtu(1,ndtu)        :tableau de travail double precision reserve
-//                      a l'utilisateur
-raddeg = %pi/180;
-
-//dtu=[ 249.9,         ..//s      1
-//      .078540,       ..//cx0    2
-//     -.0061592,      ..//cx1    3
-//      .621408e-3,    ..//cx2    4
-//     -.207040,       ..//cz0    5
-//      .029244,       ..//cz1    6
-//       83388,        ..//zm     7
-//     3.9860119e14,   ..//zmu    8
-//     6378166,        ..//rt     9
-//     1.2,            ..//ro0   10
-//     6700,           ..//h     11
-//     raddeg,         ..//      12
-//     echtf,          ..//echtf 13
-//     0,0,0,0,0,0,0,  ..// inutilises    14 a 20
-//     1000,           ..//mise a echelle v    21
-//     1,              ..//mise a echelle gam  22
-//     1.e5,           ..//mise a echelle z    23
-//     1,           ..//mise a echelle l       24
-//     1.e6      ];      //cpenal              25
-
-dtu=[ 249.9,         ..
-.078540,       ..
--.0061592,      ..
-.621408e-3,    ..
--.207040,       ..
-.029244,       ..
-83388,        ..
-3.9860119e14,   ..
-6378166,        ..
-1.2,            ..
-6700,           ..
-raddeg,         ..
-echtf,          ..
-0,0,0,0,0,0,0,  ..
-1000,           ..
-1,              ..
-1.e5,           ..
-1,           ..
-1.e6      ];      //cpenal              25
-
-y0=[7803, -1*raddeg, 121920, 0]; // etat initial
-//          (valeur arbitraire si iu(1) ou iu(2) est non nul)
-
-y0=y0./dtu(1,21:24); // mise a l'echelle de y0
-
-//  tob(1,ntob)        :instants de mesure (compatibilite avec ntob
-//                      et ntobi)
-tob       = [1];
-binf      = -20*ones(1,nu); // borne inf des parametres
-binf(1,1) = 2500/echtf;
-bsup      = 40*ones(1,nu); // borne sup des parametres
-bsup(1,1) = 4000/echtf;
-
-obs(nob,ny) = 0; // matrice d'observation
-
-//don=[762,            ..//vfin         1
-//     -5*raddeg,      ..//gamma final  2
-//     24384    ];     ..//zfin         3
-don=[762,            ..
--5*raddeg,      ..
-24384    ];
-don   = don./dtu(1,21:23); // mise a l'echelle
-nomf  = "icsenf";          // noms de subroutines de dynamique
-legfb = " croissant ";
-
-// changements pour calculer en temps retrograde
-retro = 1;
-if retro>0 then
-    legfb     = " retrograde ";
-    don1      = don;
-    don       = y0(1,1:3);
-    y0(1,1:3) = don1;
-    nomf      ="icsenb";
-end
-
-nap    = 20;  // nombre d'appels du simulateur
-imp    = 2;   // niveau de debug pour optim
-large  = 100; // taille de nu au dela de laquelle on choisit un optimiseur
-
-// pour les problemes de grande taille (alg='gc' dans l'appel de optim)
-
-exec("icseinit.sce");
-
-[co,u,g,itv,dtv]=icse(u,nomf,nap,imp);
+function demo_navet()
+
+    cd "SCI/modules/optimization/demos/icse/"
+
+    libn  = ilib_for_link("icsenb", "icsenb.f",[], "f", "", "", "", "-L"+SCI+"/modules/optimization/.libs/ -lscioptimization")
+    nlink = link("./"+libn, ["icsenb", "icsenf"], "f")
+    libn  = ilib_for_link("icsez0", "icsez0.f",[], "f")
+    nlink = link("./"+libn, "icsez0", "f")
+
+    exec("icsecontexte.sce");
+
+    t0   = 0.d0;  // instant initial
+    tf   = 1.d0;  // instant final
+    dtf  = 0;
+    ermx = 1.d-6; // test d'arret absolu sur la valeur du second membre dans la resolution de l'etat
+    iu = [0,0,1]; //  iu   :indications sur la structure du controle
+    //    iu(1)=1 si l'etat initial depend du controle constant,0 sinon
+    //    iu(2)=1 si l'etat initial depend du controle variable,0 sinon
+    //    iu(3)=1 si le second membre depend du controle constant,0 sinon
+    nuc  = 1;     // nombre de parametres independants du temps
+    nuv  = 1;     // nombre de parametres dependants du temps
+    ilin = 0;     // indicateur de linearite :
+    // 0 pour un systeme non affine
+    // 1 pour un systeme affine dont la partie lineaire n'est pas autonome
+    // 2 pour un systeme affine dont la partie lineaire est autonome
+    nti = 150;    //nombre de pas de temps correspondant a dti (premier pas de temps)
+    dti = tf/nti;
+    ntf = 00;     // nombre de pas de temps correspondant a dtf (second pas de temps)
+    // si l'on utilise un seul pas de temps,on doit prendre ntf=0
+    ny   = 4;     // dimension de l'etat a un instant donne
+    nea  = 0;     // nombre d'equations algebriques (eventuellement nul)
+    itmx = 10;    // nombre maximal d'iterations dans la resolution de l'equation d'etat discrete a un pas de temps donne
+    nex = 1;      // nombre d'experiences effectuees
+    nob = 3;      // dimension du vecteur des mesures pour une experience donnee en un instant donne
+    ntob  = 1;    // nombre d'instants de mesure pour une experience donnee
+    ntobi = 1;    // nombre d'instants de mesure correspondant a dti (premier pas de temps)
+
+    nu=nuc+nuv*(nti+ntf+1); // dimension du vecteur des parametres de controle
+
+    //  uc(1,nuc)          :controle constant
+    echtf = 2000;
+    uc    = [2500/echtf];
+
+    //  uv(1,nuv*(nti+ntf)):controle variable
+    //if nuv>0 then uv(1,nuv*(nti+ntf+1))=0; end;
+    alpha0 = .20704/.029244; legu = "alpha initial : ann. cz";     // annulation cz
+    alpha0 = 17.391;         legu = "alpha initial : finesse max"; // finesse maximum
+    legu = " navette americaine ."+legu;
+    if nuv>0 then uv=alpha0*ones(1,nuv*(nti+ntf+1)); end;
+
+    //  itu(1,nitu)        :tableau de travail entier reserve a l'utilisateur
+    itu = [0];
+
+    //  dtu(1,ndtu)        :tableau de travail double precision reserve a l'utilisateur
+    raddeg = %pi/180;
+
+    //dtu=[ 249.9,         ..//s      1
+    //      .078540,       ..//cx0    2
+    //     -.0061592,      ..//cx1    3
+    //      .621408e-3,    ..//cx2    4
+    //     -.207040,       ..//cz0    5
+    //      .029244,       ..//cz1    6
+    //       83388,        ..//zm     7
+    //     3.9860119e14,   ..//zmu    8
+    //     6378166,        ..//rt     9
+    //     1.2,            ..//ro0   10
+    //     6700,           ..//h     11
+    //     raddeg,         ..//      12
+    //     echtf,          ..//echtf 13
+    //     0,0,0,0,0,0,0,  ..// inutilises    14 a 20
+    //     1000,           ..//mise a echelle v    21
+    //     1,              ..//mise a echelle gam  22
+    //     1.e5,           ..//mise a echelle z    23
+    //     1,           ..//mise a echelle l       24
+    //     1.e6      ];      //cpenal              25
+
+    dtu=[ 249.9,         ..
+    .078540,       ..
+    -.0061592,      ..
+    .621408e-3,    ..
+    -.207040,       ..
+    .029244,       ..
+    83388,        ..
+    3.9860119e14,   ..
+    6378166,        ..
+    1.2,            ..
+    6700,           ..
+    raddeg,         ..
+    echtf,          ..
+    0,0,0,0,0,0,0,  ..
+    1000,           ..
+    1,              ..
+    1.e5,           ..
+    1,           ..
+    1.e6      ];      //cpenal              25
+
+    y0=[7803, -1*raddeg, 121920, 0]; // etat initial (valeur arbitraire si iu(1) ou iu(2) est non nul)
+
+    y0=y0./dtu(1,21:24); // mise a l'echelle de y0
+
+    //  tob(1,ntob)        :instants de mesure (compatibilite avec ntob et ntobi)
+    tob       = [1];
+    binf      = -20*ones(1,nu); // borne inf des parametres
+    binf(1,1) = 2500/echtf;
+    bsup      = 40*ones(1,nu); // borne sup des parametres
+    bsup(1,1) = 4000/echtf;
+
+    obs(nob,ny) = 0; // matrice d'observation
+
+    //don=[762,            ..//vfin         1
+    //     -5*raddeg,      ..//gamma final  2
+    //     24384    ];     ..//zfin         3
+    don=[762,            ..
+    -5*raddeg,      ..
+    24384    ];
+    don   = don./dtu(1,21:23); // mise a l'echelle
+    nomf  = "icsenf";          // noms de subroutines de dynamique
+    legfb = " croissant ";
+
+    // changements pour calculer en temps retrograde
+    retro = 1;
+    if retro>0 then
+        legfb     = " retrograde ";
+        don1      = don;
+        don       = y0(1,1:3);
+        y0(1,1:3) = don1;
+        nomf      ="icsenb";
+    end
+
+    nap    = 20;  // nombre d'appels du simulateur
+    imp    = 2;   // niveau de debug pour optim
+    large  = 100; // taille de nu au dela de laquelle on choisit un optimiseur
+
+    // pour les problemes de grande taille (alg='gc' dans l'appel de optim)
+
+    saveFormat = format();
+
+    exec("icseinit.sce");
+
+    exec("icse.sci");
+    exec("icsegen.sci");
+
+    [co, u, g, itv, dtv] = icse(u, nomf, nap, imp);
+
+    disp(u', "Best value:")
+    disp(co, "Final cost:")
+
+    deletefile("libicsenb.so");
+    deletefile("libicsez0.so");
+    deletefile("loader.sce");
+    deletefile("cleaner.sce");
+
+    format(saveFormat(1, $:-1:1));
+
+endfunction
+
+demo_navet();
+clear demo_navet;
index 3a53973..c4391cb 100644 (file)
 //  calcul coefficients optimaux du modele simplifie 5ht-plaquette
 //  **************************************************************
 
-libn  = ilib_for_link("icsest","icsest.f",[],"f");
-nlink = link("./"+libn,"icsest","f");
-
-// contexte : tue les variables de nom reserve
-exec("icse.contexte");
-
-t0   = 0.d0;  // instant initial
-tf   = 18.d1; // instant final
-dti  = 1;     // premier pas de temps
-dtf  = 2;     // second pas de temps
-ermx = 1.d-9; // test d'arret absolu sur la valeur du second membre dans
-// la resolution de l'etat
-iu = [0,0,1]; //  iu   :indications sur la structure du controle
-//    iu(1)=1 si l'etat initial depend du controle constant,0 sinon
-//    iu(2)=1 si l'etat initial depend du controle variable,0 sinon
-//    iu(3)=1 si le second membre depend du controle constant,0 sinon
-nuc  = 7;     // nombre de parametres independants du temps
-nuv  = 0;     // nombre de parametres dependants du temps
-ilin = 2;     // indicateur de linearite :
-// 0 pour un systeme non affine
-// 1 pour un systeme affine dont la partie lineaire n'est pas autonome
-// ilin=2 pour un systeme affine dont la partie lineaire est autonome
-nti = 80;     // nombre de pas de temps correspondant a dti (premier pas de temps)
-ntf = 50;     // nombre de pas de temps correspondant a dtf (second pas de temps)
-// si l'on utilise un seul pas de temps,on doit prendre ntf=0
-ny   = 4;     // dimension de l'etat a un instant donne
-nea  = 0;     // nombre d'equations algebriques (eventuellement nul)
-itmx = 10;    // nombre maximal d'iterations dans la resolution de
-// l'equation d'etat discrete a un pas de temps donne
-nex = 8;      // nombre d'experiences effectuees
-nob = 2;      // dimension du vecteur des mesures pour une experience donnee
-// en un instant donne
-ntob  = 9;    // nombre d'instants de mesure pour une experience donnee
-ntobi = 6;    // nombre d'instants de mesure correspondant a dti (premier
-// pas de temps)
-
-// ne pas modifier l'instruction suivante
-nu = nuc+nuv*(nti+ntf+1); // dimension du vecteur des parametres de controle
-
-//  uc(1,nuc)          :controle constant
-ucref = [2.d-4,1.d-3,1.d-2,5.d-3,2.d-2,1.5d-1,3.d-2];
-uc    = .1*ucref;
-
-//  uv(1,nuv*(nti+ntf)):controle variable
-//if nuv>0, uv(1,nuv*(nti+ntf))=0; end;
-//  itu(1,nitu)        :tableau de travail entier reserve a
-//                      l'utilisateur
-itu = [0];
-
-//  dtu(1,ndtu)        :tableau de travail double precision reserve
-//                      a l'utilisateur
-dtu = [0.d0];
-
-//  y0(ny)              :etat initial
-//  (inutile si iu(1) ou iu(2) est non nul)
-y0 = [4.d1,0.d0,0.d0,0.d0];
-
-//  tob(1,ntob)        :instants de mesure (compatibilite avec ntob
-//                      et ntobi)
-tob  = [1.d1,2.d1,3.d1,4.d1,6.d1,8.d1,1.1d2,1.6d2,1.8d2];
-binf = 1.d-17*ones(1,nu); // borne inf des parametres
-bsup = 1.d1*ones(1,nu);   // borne sup des parametres
-
-// termes utiles pour une dynamique lineaire ou une observation quadratique
-// b(1,ny)        = 0; // terme constant d'une dynamique lineaire
-// fy(ny,ny)      = 0; // derivee de la dynamique par rapport a l'etat
-// fu(ny,nuc+nuv) = 0; // derivee de la dynamique par rapport au controle
-obs = [0,1,1,1;0,1,0,1]; // matrice d'observation obs(nob,ny)
-
-//  don(nex*ntob*nob)  :mesures prealablement entrees dans le fichier
-//                      sero.mes.Il s'agit de donnees simulees avec
-//                      uc=[2.d-4,1.d-3,1.d-2,1.d-7,1.d-6,1.d-9,1.d-7]
-don = read("sero.mes",1,nex*ntob*nob,"(5d15.7)");
-
-nap   = 20;   // nombre d'appels du simulateur
-imp   = 2;    // niveau de debug pour optim
-large = 100;  // taille de nu au dela de laquelle on choisit un optimiseur
-// pour les problemes de grande taille (alg='gc' dans l'appel de optim)
-
-exec("icseinit.sce");
-
-[co,u,g,itv,dtv]=icse(u,"icsest",nap,imp);
+function demo_seros()
+
+    cd "SCI/modules/optimization/demos/icse/"
+
+    libn  = ilib_for_link("icsest", "icsest.f", [], "f", "", "", "", "-L"+SCI+"/modules/optimization/.libs/ -lscioptimization");
+    nlink = link("./"+libn, "icsest", "f");
+    libn  = ilib_for_link("icsez0", "icsez0.f", [], "f")
+    nlink = link("./"+libn, "icsez0", "f")
+
+    // contexte : tue les variables de nom reserve
+    exec("icsecontexte.sce");
+
+    t0   = 0.d0;  // instant initial
+    tf   = 18.d1; // instant final
+    dti  = 1;     // premier pas de temps
+    dtf  = 2;     // second pas de temps
+    ermx = 1.d-9; // test d'arret absolu sur la valeur du second membre dans la resolution de l'etat
+    iu = [0,0,1]; //  iu   :indications sur la structure du controle
+    //    iu(1)=1 si l'etat initial depend du controle constant,0 sinon
+    //    iu(2)=1 si l'etat initial depend du controle variable,0 sinon
+    //    iu(3)=1 si le second membre depend du controle constant,0 sinon
+    nuc  = 7;     // nombre de parametres independants du temps
+    nuv  = 0;     // nombre de parametres dependants du temps
+    ilin = 2;     // indicateur de linearite :
+    // 0 pour un systeme non affine
+    // 1 pour un systeme affine dont la partie lineaire n'est pas autonome
+    // 2 pour un systeme affine dont la partie lineaire est autonome
+    nti = 80;     // nombre de pas de temps correspondant a dti (premier pas de temps)
+    ntf = 50;     // nombre de pas de temps correspondant a dtf (second pas de temps)
+    // si l'on utilise un seul pas de temps,on doit prendre ntf=0
+    ny   = 4;     // dimension de l'etat a un instant donne
+    nea  = 0;     // nombre d'equations algebriques (eventuellement nul)
+    itmx = 10;    // nombre maximal d'iterations dans la resolution de l'equation d'etat discrete a un pas de temps donne
+    nex = 8;      // nombre d'experiences effectuees
+    nob = 2;      // dimension du vecteur des mesures pour une experience donnee en un instant donne
+    ntob  = 9;    // nombre d'instants de mesure pour une experience donnee
+    ntobi = 6;    // nombre d'instants de mesure correspondant a dti (premier pas de temps)
+
+    // ne pas modifier l'instruction suivante
+    nu = nuc+nuv*(nti+ntf+1); // dimension du vecteur des parametres de controle
+
+    //  uc(1,nuc)          :controle constant
+    ucref = [2.d-4,1.d-3,1.d-2,5.d-3,2.d-2,1.5d-1,3.d-2];
+    uc    = .1*ucref;
+
+    //  uv(1,nuv*(nti+ntf)):controle variable
+    //if nuv>0, uv(1,nuv*(nti+ntf))=0; end;
+    //  itu(1,nitu)        :tableau de travail entier reserve a l'utilisateur
+    itu = [0];
+
+    //  dtu(1,ndtu)        :tableau de travail double precision reserve a l'utilisateur
+    dtu = [0.d0];
+
+    //  y0(ny)              :etat initial (inutile si iu(1) ou iu(2) est non nul)
+    y0 = [4.d1,0.d0,0.d0,0.d0];
+
+    //  tob(1,ntob)        :instants de mesure (compatibilite avec ntob et ntobi)
+    tob  = [1.d1,2.d1,3.d1,4.d1,6.d1,8.d1,1.1d2,1.6d2,1.8d2];
+    binf = 1.d-17*ones(1,nu); // borne inf des parametres
+    bsup = 1.d1*ones(1,nu);   // borne sup des parametres
+
+    // termes utiles pour une dynamique lineaire ou une observation quadratique
+    // b(1,ny)        = 0; // terme constant d'une dynamique lineaire
+    // fy(ny,ny)      = 0; // derivee de la dynamique par rapport a l'etat
+    // fu(ny,nuc+nuv) = 0; // derivee de la dynamique par rapport au controle
+    obs = [0,1,1,1;0,1,0,1]; // matrice d'observation obs(nob,ny)
+
+    //  don(nex*ntob*nob)  :mesures prealablement entrees dans le fichier sero.mes.
+    //                      Il s'agit de donnees simulees avec uc=[2.d-4,1.d-3,1.d-2,1.d-7,1.d-6,1.d-9,1.d-7]
+    don = read("sero.mes",1,nex*ntob*nob,"(5d15.7)");
+
+    nap   = 20;   // nombre d'appels du simulateur
+    imp   = 2;    // niveau de debug pour optim
+    large = 100;  // taille de nu au dela de laquelle on choisit un optimiseur
+    // pour les problemes de grande taille (alg='gc' dans l'appel de optim)
+
+    saveFormat = format();
+
+    exec("icseinit.sce");
+
+    exec("icse.sci");
+    exec("icsegen.sci");
+
+    [co, u, g, itv, dtv] = icse(u, "icsest", nap, imp);
+
+    disp(u', "Best value:")
+    disp(co, "Final cost:")
+
+    deletefile("libicsest.so");
+    deletefile("libicsez0.so");
+    deletefile("loader.sce");
+    deletefile("cleaner.sce");
+
+    format(saveFormat(1, $:-1:1));
+
+endfunction
+
+demo_seros();
+clear demo_seros;
index fc86770..a9da31f 100644 (file)
@@ -11,7 +11,8 @@ function subdemolist = demo_gateway()
 
     subdemolist = [_("Non linear data fitting"), "datafit/datafit.dem.gateway.sce"; ..
     _("Optimisation"),            "optim/optim.dem.gateway.sce"; ..
-    _("fminsearch"),              "neldermead/neldermead.dem.gateway.sce"];
+    _("fminsearch"),              "neldermead/neldermead.dem.gateway.sce"; ..
+    _("ICSE"),              "icse/icse.dem.gateway.sce"];
 
     if with_module("genetic_algorithms") then
         subdemolist = [subdemolist; ..