CMSCOPE changed

18 files changed, 1185 insertions, 0 deletions

diff --git a/modules/optimization/demos/icse/README b/modules/optimization/demos/icse/README
new file mode 100755
index 000000000..362acf94e
--- /dev/null
+++ b/modules/optimization/demos/icse/README
@@ -0,0 +1,45 @@
+              Use of ICSE within SCILAB
+
+-1- The user must provide a ``simulator'' program for his problem:
+ 
+      this simulator is made of a main routine :
+
+      subroutine mysim(ind,nu,u,co,g,itv,rtv,dtv)
+      external secmbr, cost, inist
+      call icse(ind,nu,u,co,g,itv,rtv,dtv,secmbr,cost,inist)
+      end
+
+      and three additional subroutines secmbr, cost and inist
+
+     where secmbr, cost, inist are symbolic name of particular routines
+       secmbr : computes the RHS 
+       cost   : computes the cost function
+       inist  : computes the initial state
+
+       cost may be replaced by the predefined routine icsec2 for quadratic cost.
+
+       inist may be replaced by the predefined routine icsei for standard
+            initialization.
+
+ It is better to put the four routines in a single file called mysim.f
+
+-2- compile (f77) the file mysim.f 
+
+-3- Link mysim with scilab :
+    link('mysim.o',mysim)
+
+    Otherwise you have to link your program as follows:
+ .modify the routines/default/foptim.f routine to add the call to mysim 
+ .put your mysim.f file in the routines/default directory, 
+ .append mysim.o in ICSE variable definition within routines/default/Makefile
+ .re-generate a new scilab with the main Makefile
+     
+
+
+
+
+
+
+
+
+
diff --git a/modules/optimization/demos/icse/icob.sci b/modules/optimization/demos/icse/icob.sci
new file mode 100755
index 000000000..40b6af68a
--- /dev/null
+++ b/modules/optimization/demos/icse/icob.sci
@@ -0,0 +1,26 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - DIGITEO - Yann COLLETTE
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function ytob=icob(dtv)
+    // <ytob>=icob(dtv)
+    // extraction of the state at measure instant
+    // input variables:
+    // dtv(ndtv) : double precision workarea obtained by the use of
+    //             icse, icsu, icsua or icsuq
+    // output variables :
+    // ytob(ny,ntob) : value of the state at measure instant
+
+    ytob(ny,ntob) = 0.d0;
+    lob = ndtu + ny + ntob + nob*ny + nex*ntob*nob + ...
+    2*nu + ntob*nob + ny*(nti+ntf+1+ny+nuc+nuv+1);
+    ytob = matrix(dtv(lob+1:lob+ny*ntob),ny,ntob)
+endfunction
+
+
+
+
diff --git a/modules/optimization/demos/icse/icot.sci b/modules/optimization/demos/icse/icot.sci
new file mode 100755
index 000000000..9e70b3fc3
--- /dev/null
+++ b/modules/optimization/demos/icse/icot.sci
@@ -0,0 +1,22 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - DIGITEO - Yann COLLETTE
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function [ytot]=icot(dtv)
+    // <ytot>=icot(dtv)
+    // extraction of the total state
+    // input variables :
+    // dtv(ndtv) : double precision workarea obtained by the use of
+    //             icse,icsu,icsua or icsuq
+    // output variables:
+    // ytot(ny,nt) : state of the system
+
+    ytot(ny,nti+ntf) = 0.d0;
+    lot  = ndtu + ny + ntob + nob*ny + nex*ntob*nob + ...
+    2*nu + ntob*nob + ny*(1+ny+nuc+nuv);
+    ytot = matrix(dtv(lot+1:lot+ny*(nti+ntf)),ny,nti+ntf)
+endfunction
diff --git a/modules/optimization/demos/icse/icse.contexte b/modules/optimization/demos/icse/icse.contexte
new file mode 100755
index 000000000..d237f9483
--- /dev/null
+++ b/modules/optimization/demos/icse/icse.contexte
@@ -0,0 +1 @@
+clear nu u uc uv itu dtu y0 tob binf bsup b fy fu obs don;
diff --git a/modules/optimization/demos/icse/icse.dem.gateway.sce b/modules/optimization/demos/icse/icse.dem.gateway.sce
new file mode 100755
index 000000000..319b42c16
--- /dev/null
+++ b/modules/optimization/demos/icse/icse.dem.gateway.sce
@@ -0,0 +1,14 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA
+// Copyright (C) 2010 - DIGITEO - Yann COLLETTE
+// Copyright (C) 2010 - DIGITEO - Allan CORNET
+//
+// This file is released under the 3-clause BSD license. See COPYING-BSD.
+
+
+subdemolist = [_("LQV"),                               "lqv.sce"; ..
+_("Spaceship landing trajectory"),      "navet.sce"; ..
+_("Computation of optimal parameters"), "seros.sce";];
+
+subdemolist(:,2) = SCI + "/modules/optimization/demos/icse/" + subdemolist(:,2);
+
diff --git a/modules/optimization/demos/icse/icse.sci b/modules/optimization/demos/icse/icse.sci
new file mode 100755
index 000000000..1198e0401
--- /dev/null
+++ b/modules/optimization/demos/icse/icse.sci
@@ -0,0 +1,30 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) ????-2010 - DIGITEO - Yann COLLETTE
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function [co,u,g,itv,dtv]=icse(u,simu,nap,imp)
+    // Computation of the optimal control without scaling of the control
+    // and equal weighting of the observations
+    // input variables :
+    // 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
+    // imp       : debug value during optimization
+    // output variables :
+    // co        : final cost
+    // u(nu)     : final parameters
+    // g(nu)     : final gradient
+    // itv(nitv) : work area (fortran integers)
+    // dtv(ndtv) : work area (fortran double precision)
+    // Use the macros icot and icob to extract the state
+    df0 = 1;
+    nu  = prod(size(u))
+    ech = ones(1,nu);
+    cof = ones(1,nob*ntob);
+    [co,u,g,itv,dtv] = icsegen(u,simu,nap,imp)
+endfunction
diff --git a/modules/optimization/demos/icse/icsegen.sci b/modules/optimization/demos/icse/icsegen.sci
new file mode 100755
index 000000000..9bc58fed2
--- /dev/null
+++ b/modules/optimization/demos/icse/icsegen.sci
@@ -0,0 +1,45 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - DIGITEO
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function [co,u,g,itv,dtv]=icsegen(u,simu,nap,imp,ech,cof)
+    // Computation of the optimal control with scaling of the control and
+    // weighting of the observations
+    // Syntax
+    // [co,u,g,itv,dtv]=icsegen(u,simu,nap,imp,ech,cof)
+    //
+    // input variables :
+    // 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
+    // imp       : debug value during optimization
+    // ech(1,nu) : scaling coeff of the control
+    // cof(1,ntob*nob) : weighting coeff of the observations
+    // output variables :
+    // co        : final cost
+    // u(nu)     : final parameters
+    // g(nu)     : final gradient
+    // itv(nitv) : work area (fortran integers)
+    // dtv(ndtv) : work area (fortran double precision)
+    // Use the icot and icob macro to extract the state
+
+    if nu<large then
+        alg="qn";
+    else
+        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, ...
+    "ar",nap, "ti",itv,"td",dtv,"si","sd");
+    debug(0);
+endfunction
diff --git a/modules/optimization/demos/icse/icseinit.sce b/modules/optimization/demos/icse/icseinit.sce
new file mode 100755
index 000000000..63c1210ae
--- /dev/null
+++ b/modules/optimization/demos/icse/icseinit.sce
@@ -0,0 +1,56 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - DIGITEO - Yann COLLETTE
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+// icse.init.bas : initialization and tests for icse
+//**************************************************
+//
+
+// creation of y0 (initial state)
+if exists("y0")==0 then
+    y0 = ones(1,ny);
+end
+// creation of b,fy and fu and transformation to line vectors
+if exists("b")==0 then
+    b = ones(1,ny);
+end
+if exists("fy")==0 then
+    fy1 = ones(1,ny*ny);
+else
+    fy1 = matrix(fy,1,ny*ny);
+end
+if exists("fu")==0 then
+    fu1 = ones(1,ny*(nuc+nuv));
+else
+    fu1 = matrix(fu,1,ny*(nuc+nuv));
+end
+
+format("e");
+iu(5) = 0;
+[xx,nitu] = size(itu);
+[yy,ndtu] = size(dtu);
+
+if xx+yy>2 then
+    error("itu and / or dtu is not a line vector");
+end
+u = [];
+if nuc>0 then u = uc; end;
+if nuv>0 then u = [u,uv]; end;
+nu = nuc+nuv*(nti+ntf+1);
+
+if size(u)<>[1,nu] then
+    error("dimensions of the control are incompatible");
+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",..
+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",..
+[1,1],22,"i",[1,1],23,"i",[1,1],24,"i");
+
diff --git a/modules/optimization/demos/icse/icsenb.f b/modules/optimization/demos/icse/icsenb.f
new file mode 100755
index 000000000..87506e3d6
--- /dev/null
+++ b/modules/optimization/demos/icse/icsenb.f
@@ -0,0 +1,291 @@
+      subroutine icsenb(ind,nu,u,co,g,itv,rtv,dtv)
+c     Copyright INRIA
+      external nbsec,nvcout,icsei
+c     controle de la rentree de la navette en backward
+      call icse(ind,nu,u,co,g,itv,rtv,dtv,nbsec,nvcout,icsei)
+      end
+      subroutine icsenf(ind,nu,u,co,g,itv,rtv,dtv)
+      external nfsec,nvcout,icsei
+c     controle de la rentree de la navette en forward
+      call icse(ind,nu,u,co,g,itv,rtv,dtv,nfsec,nvcout,icsei)
+      end
+
+      subroutine nbsec(indf,t,y,uc,uv,f,fy,fu,b,itu,dtu,
+     & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
+     & itmx,nex,nob,ntob,ntobi,nitu,ndtu)
+c
+c      Second membre de l'equation d'etat :
+c       Parametres d'entree :
+c        indf     : vaut 1,2,3 suivant qu'on veut calculer f,fy,fu
+c        t        : instant courant
+c        y(ny)    : etat a un instant donne
+c        uc(nuc)  : controle independant du temps
+c        uv(nuv)  : controle dependant du temps, a l'instant t
+c        b(ny)    : terme constant dans le cas lineaire quadratique
+c       Parametres de sortie :
+c         indf    : >0 si  le calcul s'est  correctement effectue,<=0
+c                   sinon
+c        f(ny)    : second membre
+c        fy(ny,ny): jacobien de f par rapport a y
+c        fu(ny,nuc+nuv) : derivee de f par rapport au controle
+c       Tableaux de travail reserves a l'utilisateur :
+c        itu(nitu): tableau entier
+c        dtu(ndtu): tableau double precision
+c       (nitu et ndtu sont initialises par le common icsez).
+c!
+      implicit double precision (a-h,o-z)
+      dimension y(ny),uc(*),uv(*),f(ny),fy(ny,ny),fu(ny,*),
+     &     b(ny),itu(*),dtu(*),iu(5)
+      call navetb(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
+     &  ny,nuc,nuv,nitu,ndtu)
+      end
+
+      subroutine nfsec(indf,t,y,uc,uv,f,fy,fu,b,itu,dtu,
+     & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
+     & itmx,nex,nob,ntob,ntobi,nitu,ndtu)
+c
+c      Second membre de l'equation d'etat :
+c       Parametres d'entree :
+c        indf     : vaut 1,2,3 suivant qu'on veut calculer f,fy,fu
+c        t        : instant courant
+c        y(ny)    : etat a un instant donne
+c        uc(nuc)  : controle independant du temps
+c        uv(nuv)  : controle dependant du temps, a l'instant t
+c        b(ny)    : terme constant dans le cas lineaire quadratique
+c       Parametres de sortie :
+c         indf    : >0 si  le calcul s'est  correctement effectue,<=0
+c                   sinon
+c        f(ny)    : second membre
+c        fy(ny,ny): jacobien de f par rapport a y
+c        fu(ny,nuc+nuv) : derivee de f par rapport au controle
+c       Tableaux de travail reserves a l'utilisateur :
+c        itu(nitu): tableau entier
+c        dtu(ndtu): tableau double precision
+c       (nitu et ndtu sont initialises par le common icsez).
+c!
+      implicit double precision (a-h,o-z)
+      dimension y(ny),uc(*),uv(*),f(ny),fy(ny,ny),fu(ny,*),
+     &     b(ny),itu(*),dtu(*),iu(5)
+      call navetf(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
+     &  ny,nuc,nuv,nitu,ndtu)
+      end
+
+      subroutine nvcout(indc,nu,tob,obs,cof,ytob,ob,u,
+     & c,cy,g,yob,d,itu,dtu,
+     & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
+     & itmx,nex,nob,ntob,ntobi,nitu,ndtu)
+c
+c
+c     critere standard des moindres carres
+c
+c
+c      Cout ponctuel :
+c       Parametres d'entree :
+c        indc     : 1 si on desire calculer c2,2 si on desire
+c                   calculer c2y,c2u
+c        tob      : instants de mesure
+c        obs      : matrice d'observation
+c        cof      : coefficients de ponderation du cout
+c        ytob     : valeur de l'etat aux instants d'observation
+c        ob       : mesures
+c        u(nu)    : controle.Le controle variable est stocke a la
+c                   suite du controle suite du constant.
+c       Parametres de sortie :
+c        indc     : comme pour icsec1
+c        c2       : cout
+c        c2y(ny,ntob) : derivee de c2 par rapport a y
+c        g(nu)  : derivee de c2 par rapport a u
+c!
+      implicit double precision (a-h,o-z)
+      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
+      call navetc(indc,nu,tob,obs,cof,ytob,ob,u,
+     &     c,cy,g,yob,d,itu,dtu)
+c
+      end
+      subroutine navetb(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
+     &  ny,nuc,nuv,nitu,ndtu)
+C        test de icse : navette
+C        fb,inria
+C!
+      implicit double precision (a-h,o-z)
+      dimension y(ny),uc(*),uv(*),f(ny),fy(ny,ny),fu(ny,nuc+nuv),itu(nit
+     &u),dtu(ndtu)
+      f(1)=-uc(1)*dtu(13)*(-dtu(8)*sin(y(2)*dtu(22))/(y(3)*dtu(23)+dtu(9
+     &))**2-0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu(2))*dtu
+     &(10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(21)
+      f(2)=-uc(1)*dtu(13)*(y(1)*dtu(21)*cos(y(2)*dtu(22))/(y(3)*dtu(23)+
+     &dtu(9))-dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21)*(y(3)*dtu(23)+dtu(9
+     &))**2)+0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*dtu(21)*exp(-
+     &y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(22)
+      f(3)=-uc(1)*y(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
+      f(4)=-uc(1)*y(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+
+     &dtu(9))*dtu(24))
+      if(indf.eq.2) then
+        fy(1,1)=dtu(1)*uc(1)*y(1)*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu(2))*
+     &  dtu(10)*dtu(13)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/dtu(7)
+        fy(1,2)=uc(1)*dtu(8)*dtu(13)*dtu(22)*cos(y(2)*dtu(22))/(dtu(21)*
+     &  (y(3)*dtu(23)+dtu(9))**2)
+        fy(1,3)=-uc(1)*dtu(13)*dtu(23)*(2*dtu(8)*sin(y(2)*dtu(22))/(y(3)
+     &  *dtu(23)+dtu(9))**3+0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dt
+     &  u(3)+dtu(2))*dtu(10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/(dtu(
+     &  7)*dtu(11)))/dtu(21)
+        fy(1,4)=0
+        fy(2,1)=-uc(1)*dtu(13)*dtu(21)*(cos(y(2)*dtu(22))/(y(3)*dtu(23)+
+     &  dtu(9))+dtu(8)*cos(y(2)*dtu(22))/(y(1)**2*dtu(21)**2*(y(3)*dtu(2
+     &  3)+dtu(9))**2)+0.5*dtu(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*exp(-y(3
+     &  )*dtu(23)/dtu(11))/dtu(7))/dtu(22)
+        fy(2,2)=-uc(1)*dtu(13)*(dtu(8)*sin(y(2)*dtu(22))/(y(1)*dtu(21)*(
+     &  y(3)*dtu(23)+dtu(9))**2)-y(1)*dtu(21)*sin(y(2)*dtu(22))/(y(3)*dt
+     &  u(23)+dtu(9)))
+        fy(2,3)=-uc(1)*dtu(13)*dtu(23)*(-y(1)*dtu(21)*cos(y(2)*dtu(22))/
+     &  (y(3)*dtu(23)+dtu(9))**2+2*dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21
+     &  )*(y(3)*dtu(23)+dtu(9))**3)-0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5)
+     &  )*dtu(10)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/(dtu(7)*dtu(11)))/d
+     &  tu(22)
+        fy(2,4)=0
+        fy(3,1)=-uc(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
+        fy(3,2)=-uc(1)*y(1)*dtu(13)*dtu(21)*dtu(22)*cos(y(2)*dtu(22))/dt
+     &  u(23)
+        fy(3,3)=0
+        fy(3,4)=0
+        fy(4,1)=-uc(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+
+     &  dtu(9))*dtu(24))
+        fy(4,2)=uc(1)*y(1)*dtu(13)*dtu(21)*dtu(22)*sin(y(2)*dtu(22))/((y
+     &  (3)*dtu(23)+dtu(9))*dtu(24))
+        fy(4,3)=uc(1)*y(1)*dtu(13)*dtu(21)*dtu(23)*cos(y(2)*dtu(22))/((y
+     &  (3)*dtu(23)+dtu(9))**2*dtu(24))
+        fy(4,4)=0
+      end if
+      if(indf.eq.3) then
+        fu(1,1)=-dtu(13)*(-dtu(8)*sin(y(2)*dtu(22))/(y(3)*dtu(23)+dtu(9)
+     &  )**2-0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu(2))*dt
+     &  u(10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(21)
+        fu(1,2)=0.5*dtu(1)*uc(1)*y(1)**2*(2*uv(1)*dtu(4)+dtu(3))*dtu(10)
+     &  *dtu(13)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/dtu(7)
+        fu(2,1)=-dtu(13)*(y(1)*dtu(21)*cos(y(2)*dtu(22))/(y(3)*dtu(23)+d
+     &  tu(9))-dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21)*(y(3)*dtu(23)+dtu(
+     &  9))**2)+0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*dtu(21)*ex
+     &  p(-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(22)
+        fu(2,2)=-0.5*dtu(1)*uc(1)*y(1)*dtu(6)*dtu(10)*dtu(13)*dtu(21)*ex
+     &  p(-y(3)*dtu(23)/dtu(11))/(dtu(7)*dtu(22))
+        fu(3,1)=-y(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
+        fu(3,2)=0
+        fu(4,1)=-y(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+d
+     &  tu(9))*dtu(24))
+        fu(4,2)=0
+      end if
+      end
+
+
+      subroutine navetf(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
+     &  ny,nuc,nuv,nitu,ndtu)
+C        test de icse : navette
+C        fb,inria
+C        [" "," "]
+C
+      implicit double precision (a-h,o-z)
+      dimension y(ny),uc(*),uv(*),f(ny),fy(ny,ny),fu(ny,nuc+nuv),itu(nit
+     &u),dtu(ndtu)
+      f(1)=uc(1)*dtu(13)*(-dtu(8)*sin(y(2)*dtu(22))/(y(3)*dtu(23)+dtu(9)
+     &)**2-0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu(2))*dtu(
+     &10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(21)
+      f(2)=uc(1)*dtu(13)*(y(1)*dtu(21)*cos(y(2)*dtu(22))/(y(3)*dtu(23)+d
+     &tu(9))-dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21)*(y(3)*dtu(23)+dtu(9)
+     &)**2)+0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*dtu(21)*exp(-y
+     &(3)*dtu(23)/dtu(11))/dtu(7))/dtu(22)
+      f(3)=uc(1)*y(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
+      f(4)=uc(1)*y(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+d
+     &tu(9))*dtu(24))
+      if(indf.eq.2) then
+        fy(1,1)=-1.0*dtu(1)*uc(1)*y(1)*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu
+     &  (2))*dtu(10)*dtu(13)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/dtu(7)
+        fy(1,2)=-uc(1)*dtu(8)*dtu(13)*dtu(22)*cos(y(2)*dtu(22))/(dtu(21)
+     &  *(y(3)*dtu(23)+dtu(9))**2)
+        fy(1,3)=uc(1)*dtu(13)*dtu(23)*(2*dtu(8)*sin(y(2)*dtu(22))/(y(3)*
+     &  dtu(23)+dtu(9))**3+0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dtu
+     &  (3)+dtu(2))*dtu(10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/(dtu(7
+     &  )*dtu(11)))/dtu(21)
+        fy(1,4)=0
+        fy(2,1)=uc(1)*dtu(13)*dtu(21)*(cos(y(2)*dtu(22))/(y(3)*dtu(23)+d
+     &  tu(9))+dtu(8)*cos(y(2)*dtu(22))/(y(1)**2*dtu(21)**2*(y(3)*dtu(23
+     &  )+dtu(9))**2)+0.5*dtu(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*exp(-y(3)
+     &  *dtu(23)/dtu(11))/dtu(7))/dtu(22)
+        fy(2,2)=uc(1)*dtu(13)*(dtu(8)*sin(y(2)*dtu(22))/(y(1)*dtu(21)*(y
+     &  (3)*dtu(23)+dtu(9))**2)-y(1)*dtu(21)*sin(y(2)*dtu(22))/(y(3)*dtu
+     &  (23)+dtu(9)))
+        fy(2,3)=uc(1)*dtu(13)*dtu(23)*(-y(1)*dtu(21)*cos(y(2)*dtu(22))/(
+     &  y(3)*dtu(23)+dtu(9))**2+2*dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21)
+     &  *(y(3)*dtu(23)+dtu(9))**3)-0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5))
+     &  *dtu(10)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/(dtu(7)*dtu(11)))/dt
+     &  u(22)
+        fy(2,4)=0
+        fy(3,1)=uc(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
+        fy(3,2)=uc(1)*y(1)*dtu(13)*dtu(21)*dtu(22)*cos(y(2)*dtu(22))/dtu
+     &  (23)
+        fy(3,3)=0
+        fy(3,4)=0
+        fy(4,1)=uc(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+d
+     &  tu(9))*dtu(24))
+        fy(4,2)=-uc(1)*y(1)*dtu(13)*dtu(21)*dtu(22)*sin(y(2)*dtu(22))/((
+     &  y(3)*dtu(23)+dtu(9))*dtu(24))
+        fy(4,3)=-uc(1)*y(1)*dtu(13)*dtu(21)*dtu(23)*cos(y(2)*dtu(22))/((
+     &  y(3)*dtu(23)+dtu(9))**2*dtu(24))
+        fy(4,4)=0
+      end if
+      if(indf.eq.3) then
+        fu(1,1)=dtu(13)*(-dtu(8)*sin(y(2)*dtu(22))/(y(3)*dtu(23)+dtu(9))
+     &  **2-0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu(2))*dtu
+     &  (10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(21)
+        fu(1,2)=-0.5*dtu(1)*uc(1)*y(1)**2*(2*uv(1)*dtu(4)+dtu(3))*dtu(10
+     &  )*dtu(13)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/dtu(7)
+        fu(2,1)=dtu(13)*(y(1)*dtu(21)*cos(y(2)*dtu(22))/(y(3)*dtu(23)+dt
+     &  u(9))-dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21)*(y(3)*dtu(23)+dtu(9
+     &  ))**2)+0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*dtu(21)*exp
+     &  (-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(22)
+        fu(2,2)=0.5*dtu(1)*uc(1)*y(1)*dtu(6)*dtu(10)*dtu(13)*dtu(21)*exp
+     &  (-y(3)*dtu(23)/dtu(11))/(dtu(7)*dtu(22))
+        fu(3,1)=y(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
+        fu(3,2)=0
+        fu(4,1)=y(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+dt
+     &  u(9))*dtu(24))
+        fu(4,2)=0
+      end if
+      end
+
+      subroutine navetc(indc,nu,tob,obs,cof,ytob,ob,u,
+     & c,cy,g,yob,d,itu,dtu)
+c
+c     test probleme navette
+c
+c     sous programme appele par icse.f qui donne :
+c     pour indc=1,le cout:c
+c     pour indc=2,la matrice derivee du cout par rapport a l'etat
+c                 calcule aux instants de mesure:cy(ny,ntob)
+c     sorties :
+c
+c     c        double precision
+c              cout
+c     cy       double precision (ny,ntob)
+c              derivee du cout par rapport a l'etat calcule aux
+c              instants de mesure
+c     g        double precision (nu)
+c              le gradient est initialise a la derivee partielle du
+c              cout par rapport au controle
+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(*)
+c
+      cpen=dtu(25)
+      cy(1)=ytob(1)-ob(1)
+      cy(2)=ytob(2)-ob(2)
+      cy(3)=ytob(3)-ob(3)
+      cy(4)=1.0d+0/cpen
+      c=ytob(4)/cpen + ( cy(1)**2 + cy(2)**2 + cy(3)**2 )/2.0d+0
+      do 10 k=1,nu
+10    g(k)=0.0d+0
+      end
diff --git a/modules/optimization/demos/icse/icsest.f b/modules/optimization/demos/icse/icsest.f
new file mode 100755
index 000000000..262ca0b6e
--- /dev/null
+++ b/modules/optimization/demos/icse/icsest.f
@@ -0,0 +1,99 @@
+      subroutine icsest(ind,nu,u,co,g,itv,rtv,dtv)
+c     Copyright INRIA
+      external stsec,icsec2,icsei
+c
+c     identification de parametres d'un systeme lineaire
+c     serotonine G.Launay
+c
+      call icse(ind,nu,u,co,g,itv,rtv,dtv,stsec,icsec2,icsei)
+      end
+
+      subroutine stsec(indf,t,y,uc,uv,f,fy,fu,b,itu,dtu,
+     & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
+     & itmx,nex,nob,ntob,ntobi,nitu,ndtu)
+c
+c      Second membre de l'equation d'etat :
+c       Parametres d'entree :
+c        indf     : vaut 1,2,3 suivant qu'on veut calculer f,fy,fu
+c        t        : instant courant
+c        y(ny)    : etat a un instant donne
+c        uc(nuc)  : controle independant du temps
+c        uv(nuv)  : controle dependant du temps, a l'instant t
+c        b(ny)    : terme constant dans le cas lineaire quadratique
+c       Parametres de sortie :
+c         indf    : >0 si  le calcul s'est  correctement effectue,<=0
+c                   sinon
+c        f(ny)    : second membre
+c        fy(ny,ny): jacobien de f par rapport a y
+c        fu(ny,nuc+nuv) : derivee de f par rapport au controle
+c       Tableaux de travail reserves a l'utilisateur :
+c        itu(nitu): tableau entier
+c        dtu(ndtu): tableau double precision
+c       (nitu et ndtu sont initialises par le common icsez).
+c!
+      implicit double precision (a-h,o-z)
+      dimension y(ny),uc(*),uv(*),f(ny),fy(ny,ny),fu(ny,*),
+     &     b(ny),itu(*),dtu(*),iu(5)
+c
+      call seros(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
+     &  ny,nuc,nuv,nitu,ndtu)
+      end
+
+      subroutine seros(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
+     &  ny,nuc,nuv,nitu,ndtu)
+c     sous programme appele par icse.fortran qui donne :
+c     pour indf=1,les seconds membres du systeme:f(ny)
+c     pour indf=2,la matrice derivee des seconds membres par
+c                 rapport a l'etat:fy(ny,ny)
+c     pour indf=3,la matrice derivee des seconds membres par
+c                 rapport aux parametres:fu(ny,nu)
+c     on prend ici un modele lineaire ,en supposant que toute la 5HT
+c     excretee est degradee avant d'arriver dans le milieu exterieur
+      implicit double precision (a-h,o-z)
+      dimension y(ny),uc(*),uv(*),f(ny),fy(ny,ny),fu(ny,*),
+     &itu(*),dtu(*)
+c     seconds membres:
+      if (indf.eq.1) then
+        f(1)=-(uc(1)+uc(2))*y(1)
+        f(2)=uc(1)*y(1)-(uc(5)+uc(3))*y(2)
+        f(3)=uc(2)*y(1)+uc(3)*y(2)-(uc(4)+uc(6))*y(3)
+        f(4)=uc(4)*y(3)-uc(7)*y(4)
+        return
+      endif
+c     derivee des seconds membres par rapport a l'etat:
+c     pour 1<=i<=ny et 1<=j<=ny,fy(i,j)=d(f(i))/dyj
+      if (indf.eq.2) then
+        do 10 i=1,ny
+        do 10 j=1,ny
+10      fy(i,j)=0.0d+0
+        fy(1,1)=-uc(1)-uc(2)
+        fy(2,1)=uc(1)
+        fy(3,1)=uc(2)
+        fy(2,2)=-uc(3)-uc(5)
+        fy(3,2)=uc(3)
+        fy(3,3)=-uc(4)-uc(6)
+        fy(4,3)=uc(4)
+        fy(4,4)=-uc(7)
+        return
+      endif
+c     derivee des seconds membres par rapport aux parametres:
+c     pour 1<=i<=ny et 1<=j<=nu,fu(i,j)=d(f(i))/duj
+      if (indf.eq.3) then
+        do 20 i=1,ny
+        do 20 j=1,nuc+nuv
+20      fu(i,j)=0.0d+0
+        fu(1,1)=-y(1)
+        fu(2,1)=y(1)
+        fu(1,2)=-y(1)
+        fu(3,2)=y(1)
+        fu(2,3)=-y(2)
+        fu(3,3)=y(2)
+        fu(3,4)=-y(3)
+        fu(4,4)=y(3)
+        fu(2,5)=-y(2)
+        fu(3,6)=-y(3)
+        fu(4,7)=-y(4)
+        return
+      endif
+      end
+
diff --git a/modules/optimization/demos/icse/icsu.sci b/modules/optimization/demos/icse/icsu.sci
new file mode 100755
index 000000000..3fe0553a9
--- /dev/null
+++ b/modules/optimization/demos/icse/icsu.sci
@@ -0,0 +1,47 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - DIGITEO - Yann COLLETTE
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function [co,u,g,itv,dtv]=icsu(u,simu,nap,imp)
+    // Computation of the optimal control with weighting of the initial control.
+    // A scaling brings all the componant of the initial control to 1.
+    // The componant initially null will stay null.
+    // Only works if the lower bound is positive.
+    // Syntax
+    // [co,u,g,itv,dtv] = icsu(u,simu,nap,imp)
+
+    // input variables :
+    // 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
+    // imp       : debug value during optimization
+    // output variables  :
+    // co        : final cost
+    // u(nu)     : final parameters
+    // g(nu)     : final gradient
+    // itv(nitv) : work area (fortran integers)
+    // dtv(ndtv) : work area (fortran double precision)
+    // Use the macros icot and icob to extract the total state
+    // or the state at measure time instants of dtv.
+
+    df0 = 1;
+    if min(binf) <=0 then
+        error("appel de icsu avec binf non strictement positif");
+    end
+    for i=1:nu
+        u(1,i) = max( [binf(1,i),min([u(1,i),bsup(1,i)])] );
+    end
+    ech  = u;
+    binf = binf./u;
+    bsup = bsup./u;
+    u    = ones(1,nu);
+    cof  = ones(1,ntob*nob);
+    [co,u,g,itv,dtv] = icsegen(u,simu,nap,imp);
+    u = ech.*u;
+endfunction
diff --git a/modules/optimization/demos/icse/icsua.sci b/modules/optimization/demos/icse/icsua.sci
new file mode 100755
index 000000000..a5db9dbd6
--- /dev/null
+++ b/modules/optimization/demos/icse/icsua.sci
@@ -0,0 +1,51 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - DIGITEO - Yann COLLETTE
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function [co,u,g,itv,dtv,cof]=icsua(u,simu,nap,imp)
+    // Computation of the optimal control with weighting of the initial control
+    // and arithmetic weighting of the observations.
+    // A scaling brings all the components of the control to 1.
+    // The initially null componants stay null.
+    // Only works if the lower bound is positive.
+    // Weighting formula :
+    // cof(i,j)=nex / (abs(ob(1,j,i)) + ... + abs(ob(nex,j,i)) )
+    // Syntax
+    // [co,u,g,itv,dtv,cof]=icsua(u,simu,nap,imp)
+
+    // input variables :
+    // u(nu)     : initial parameters
+    // nap       : maximum number of call to the simulator
+    // imp       : debug value during optimization
+    // output variables :
+    // co            : final cost
+    // u(nu)         : final parameters
+    // g(nu)         : finale gradient
+    // itv(nitv)     : work area (fortran integers)
+    // dtv(ndtv)     : work area (fortran double precision)
+    // cof(nob,ntob) : weighting coefficient of the cost
+    // Use the macros icot and icob to extract the state
+
+    df0 = 1;
+    if min(binf) <=0 then
+        error("call to icsua with binf not strictly positive");
+    end
+    for i=1:nu
+        u(1,i)=max( [binf(1,i),min([u(1,i),bsup(1,i)])] )
+    end
+    ech  = u;
+    binf = binf./u;
+    bsup = bsup./u;
+    u    = ones(1,nu);
+    ico  = 1;
+    yob  = 0.d0*ones(nob,ntob);
+    ob   = don;
+    [cof] = fort("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;
+endfunction
diff --git a/modules/optimization/demos/icse/icsuq.sci b/modules/optimization/demos/icse/icsuq.sci
new file mode 100755
index 000000000..590420ad2
--- /dev/null
+++ b/modules/optimization/demos/icse/icsuq.sci
@@ -0,0 +1,60 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - DIGITEO - Yann COLLETTE
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function [co,u,g,itv,dtv,cof]=icsuq(u,simu,nap,imp,obs,ytob)
+    // Computation of the optimal control with weighting of the initial control
+    // and geometric weighting of the observation weights.
+    // A scaling bring all the componant of the control to 1.
+    // The componant initially null will stay null.
+    // Only works if the lower bound is positive.
+    // Weighting formula :
+    // cof(i,j) = 1/2*[(ytob(i,j)-ob(1,j,i))**2+..+(ytob(i,j)-ob(nex,j,i)**2]
+    //
+    // Syntax
+    // [co,u,g,itv,dtv,cof]=icsuq(u,nap,imp,obs,ytob)
+    // input variables :
+    // 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
+    // imp           : debug value during optimization
+    // obs(nob,ny)   : observation matrix
+    // ytob(ny,ntob) : initial value of the state at measure time instants
+    //                 obtained by icob after using icse,icsu,icsua or icsuq
+    // output variables :
+    // co            : final cost
+    // u(nu)         : final parameters
+    // g(nu)         : final gradient
+    // itv(nitv)     : work area (fortran integers)
+    // dtv(ndtv)     : work area (fortran double precision)
+    // cof(nob,ntob) : weighting coefficients of the cost
+    // Use the macros icot and icob to extract the state
+
+    df0 = 1;
+    if min(binf) <=0 then
+        error("call to icsuq with binf not strictly positive");
+    end
+    for i=1:nu
+        u(1,i)=max([ binf(1,i),min([u(1,i),bsup(1,i)])] );
+    end
+    ech  = u;
+    binf = binf./u;
+    bsup = bsup./u;
+    u    = ones(1,nu);
+    ico  = 2;
+    yob  = obs*ytob;
+    ob   = don;
+
+    [cof] = fort("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;
+endfunction
diff --git a/modules/optimization/demos/icse/icsvisu.sci b/modules/optimization/demos/icse/icsvisu.sci
new file mode 100755
index 000000000..5fe06da0d
--- /dev/null
+++ b/modules/optimization/demos/icse/icsvisu.sci
@@ -0,0 +1,31 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function [yobt,yob]=icsvisu(vue)
+    // Visualization of the observation of the state
+    // Syntax
+    // [yobt,yob] = icsvisu(vue)
+    // input variables :
+    // vue  : vector of indexes of the componant to be visualized
+    // output variables :
+    // yobt : observation of the state at every instants
+    // yob  : idem, at observation instant
+
+    tob1 = [t0,tob];
+    tobt = t0*ones(1,nti+ntf) + ...
+    [dti*[1:nti],dtf*[(nti*dti/dtf)+1:(nti*dti/dtf)+ntf]];
+    tobt = [t0,tobt];
+
+    [ytob] = icob(dtv);
+    [ytot] = icot(dtv);
+
+    yobt = obs*[y0',ytot];
+    yob  = obs*[y0',ytob];
+    for i=vue
+        plot(tobt,yobt(i,:)); plot(tob1,yob(i,:));
+    end
+endfunction
diff --git a/modules/optimization/demos/icse/lqv.sce b/modules/optimization/demos/icse/lqv.sce
new file mode 100755
index 000000000..3e03c5c59
--- /dev/null
+++ b/modules/optimization/demos/icse/lqv.sce
@@ -0,0 +1,92 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - DIGITEO - Yann COLLETTE
+// Copyright (C) 2010 - DIGITEO - Allan CORNET
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+//                     lqv.bas : demo of icse
+//  *************************************************************
+//
+
+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);
+
+endfunction
+
diff --git a/modules/optimization/demos/icse/navet.sce b/modules/optimization/demos/icse/navet.sce
new file mode 100755
index 000000000..4de0f3822
--- /dev/null
+++ b/modules/optimization/demos/icse/navet.sce
@@ -0,0 +1,152 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - DIGITEO
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+//                navet.bas : demo de icse
+//  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);
diff --git a/modules/optimization/demos/icse/sero.mes b/modules/optimization/demos/icse/sero.mes
new file mode 100755
index 000000000..d949efbf6
--- /dev/null
+++ b/modules/optimization/demos/icse/sero.mes
@@ -0,0 +1,29 @@
+  0.2792531d+00  0.2792531d+00  0.2792531d+00  0.2792531d+00  0.2792531d+00
+  0.2792531d+00  0.2792531d+00  0.2792531d+00  0.3820470d+00  0.3820470d+00
+  0.3820470d+00  0.3820470d+00  0.3820470d+00  0.3820470d+00  0.3820470d+00
+  0.3820470d+00  0.4318366d+00  0.4318366d+00  0.4318366d+00  0.4318366d+00
+  0.4318366d+00  0.4318366d+00  0.4318366d+00  0.4318366d+00  0.4674049d+00
+  0.4674049d+00  0.4674049d+00  0.4674049d+00  0.4674049d+00  0.4674049d+00
+  0.4674049d+00  0.4674049d+00  0.5005941d+00  0.5005941d+00  0.5005941d+00
+  0.5005941d+00  0.5005941d+00  0.5005941d+00  0.5005941d+00  0.5005941d+00
+  0.5163734d+00  0.5163734d+00  0.5163734d+00  0.5163734d+00  0.5163734d+00
+  0.5163734d+00  0.5163734d+00  0.5163734d+00  0.5152169d+00  0.5152169d+00
+  0.5152169d+00  0.5152169d+00  0.5152169d+00  0.5152169d+00  0.5152169d+00
+  0.5152169d+00  0.4957038d+00  0.4957038d+00  0.4957038d+00  0.4957038d+00
+  0.4957038d+00  0.4957038d+00  0.4957038d+00  0.4957038d+00  0.4849529d+00
+  0.4849529d+00  0.4849529d+00  0.4849529d+00  0.4849529d+00  0.4849529d+00
+  0.4849529d+00  0.4849529d+00  0.7521229d-01  0.7521229d-01  0.7521229d-01
+  0.7521229d-01  0.7521229d-01  0.7521229d-01  0.7521229d-01  0.7521229d-01
+  0.1342286d+00  0.1342286d+00  0.1342286d+00  0.1342286d+00  0.1342286d+00
+  0.1342286d+00  0.1342286d+00  0.1342286d+00  0.1750144d+00  0.1750144d+00
+  0.1750144d+00  0.1750144d+00  0.1750144d+00  0.1750144d+00  0.1750144d+00
+  0.1750144d+00  0.2095138d+00  0.2095138d+00  0.2095138d+00  0.2095138d+00
+  0.2095138d+00  0.2095138d+00  0.2095138d+00  0.2095138d+00  0.2454738d+00
+  0.2454738d+00  0.2454738d+00  0.2454738d+00  0.2454738d+00  0.2454738d+00
+  0.2454738d+00  0.2454738d+00  0.2657111d+00  0.2657111d+00  0.2657111d+00
+  0.2657111d+00  0.2657111d+00  0.2657111d+00  0.2657111d+00  0.2657111d+00
+  0.2722920d+00  0.2722920d+00  0.2722920d+00  0.2722920d+00  0.2722920d+00
+  0.2722920d+00  0.2722920d+00  0.2722920d+00  0.2663367d+00  0.2663367d+00
+  0.2663367d+00  0.2663367d+00  0.2663367d+00  0.2663367d+00  0.2663367d+00
+  0.2663367d+00  0.2609468d+00  0.2609468d+00  0.2609468d+00  0.2609468d+00
+  0.2609468d+00  0.2609468d+00  0.2609468d+00  0.2609468d+00
diff --git a/modules/optimization/demos/icse/seros.sce b/modules/optimization/demos/icse/seros.sce
new file mode 100755
index 000000000..3a53973b2
--- /dev/null
+++ b/modules/optimization/demos/icse/seros.sce
@@ -0,0 +1,94 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2010 - DIGITEO - Yann COLLETTE
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+//                sero.bas : demo de icse
+//  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);