1//
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) ????-2008 - INRIA
4// Copyright (C) 2010 - DIGITEO - Yann COLLETTE
5//
6// This file is distributed under the same license as the Scilab package.
7//
8
9//                sero.bas : demo de icse
10//  calcul coefficients optimaux du modele simplifie 5ht-plaquette
11//  **************************************************************
12
13function demo_seros()
14
15    copyfile("SCI/modules/optimization/demos/icse/", "TMPDIR/icse");
16    cd("TMPDIR/icse");
17
18    libn  = ilib_for_link("icsest", "icsest.f", [], "f", "", "", "",  "-L"+SCI+"/modules/optimization/.libs/ -L"+SCI+"/../../lib/scilab/ -lscioptimization");
19    nlink = link("./"+libn, "icsest", "f");
20    libn  = ilib_for_link("icsez0", "icsez0.f", [], "f")
21    nlink = link("./"+libn, "icsez0", "f")
22
23    // contexte : tue les variables de nom reserve
24    exec("icsecontexte.sce");
25
26    t0   = 0.d0;  // instant initial
27    tf   = 18.d1; // instant final
28    dti  = 1;     // premier pas de temps
29    dtf  = 2;     // second pas de temps
30    ermx = 1.d-9; // test d'arret absolu sur la valeur du second membre dans la resolution de l'etat
31    iu = [0,0,1]; //  iu   :indications sur la structure du controle
32    //    iu(1)=1 si l'etat initial depend du controle constant,0 sinon
33    //    iu(2)=1 si l'etat initial depend du controle variable,0 sinon
34    //    iu(3)=1 si le second membre depend du controle constant,0 sinon
35    nuc  = 7;     // nombre de parametres independants du temps
36    nuv  = 0;     // nombre de parametres dependants du temps
37    ilin = 2;     // indicateur de linearite :
38    // 0 pour un systeme non affine
39    // 1 pour un systeme affine dont la partie lineaire n'est pas autonome
40    // 2 pour un systeme affine dont la partie lineaire est autonome
41    nti = 80;     // nombre de pas de temps correspondant a dti (premier pas de temps)
42    ntf = 50;     // nombre de pas de temps correspondant a dtf (second pas de temps)
43    // si l'on utilise un seul pas de temps,on doit prendre ntf=0
44    ny   = 4;     // dimension de l'etat a un instant donne
45    nea  = 0;     // nombre d'equations algebriques (eventuellement nul)
46    itmx = 10;    // nombre maximal d'iterations dans la resolution de l'equation d'etat discrete a un pas de temps donne
47    nex = 8;      // nombre d'experiences effectuees
48    nob = 2;      // dimension du vecteur des mesures pour une experience donnee en un instant donne
49    ntob  = 9;    // nombre d'instants de mesure pour une experience donnee
50    ntobi = 6;    // nombre d'instants de mesure correspondant a dti (premier pas de temps)
51
52    // ne pas modifier l'instruction suivante
53    nu = nuc+nuv*(nti+ntf+1); // dimension du vecteur des parametres de controle
54
55    //  uc(1,nuc)          :controle constant
56    ucref = [2.d-4,1.d-3,1.d-2,5.d-3,2.d-2,1.5d-1,3.d-2];
57    uc    = .1*ucref;
58
59    //  uv(1,nuv*(nti+ntf)):controle variable
60    //if nuv>0, uv(1,nuv*(nti+ntf))=0; end;
61    //  itu(1,nitu)        :tableau de travail entier reserve a l'utilisateur
62    itu = [0];
63
64    //  dtu(1,ndtu)        :tableau de travail double precision reserve a l'utilisateur
65    dtu = [0.d0];
66
67    //  y0(ny)              :etat initial (inutile si iu(1) ou iu(2) est non nul)
68    y0 = [4.d1,0.d0,0.d0,0.d0];
69
70    //  tob(1,ntob)        :instants de mesure (compatibilite avec ntob et ntobi)
71    tob  = [1.d1,2.d1,3.d1,4.d1,6.d1,8.d1,1.1d2,1.6d2,1.8d2];
72    binf = 1.d-17*ones(1,nu); // borne inf des parametres
73    bsup = 1.d1*ones(1,nu);   // borne sup des parametres
74
75    // termes utiles pour une dynamique lineaire ou une observation quadratique
76    // b(1,ny)        = 0; // terme constant d'une dynamique lineaire
77    // fy(ny,ny)      = 0; // derivee de la dynamique par rapport a l'etat
78    // fu(ny,nuc+nuv) = 0; // derivee de la dynamique par rapport au controle
79    obs = [0,1,1,1;0,1,0,1]; // matrice d'observation obs(nob,ny)
80
81    //  don(nex*ntob*nob)  :mesures prealablement entrees dans le fichier sero.mes.
82    //                      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]
83    don = read("sero.mes",1,nex*ntob*nob,"(5d15.7)");
84
85    nap   = 20;   // nombre d'appels du simulateur
86    imp   = 2;    // niveau de debug pour optim
87    large = 100;  // taille de nu au dela de laquelle on choisit un optimiseur
88    // pour les problemes de grande taille (alg='gc' dans l'appel de optim)
89
90    saveFormat = format();
91
92    exec("icseinit.sce");
93
94    exec("icse.sci");
95    exec("icsegen.sci");
96
97    [co, u, g, itv, dtv] = icse(u, "icsest", nap, imp);
98
99    disp("Best value:", u')
100    disp("Final cost:", co)
101
102    deletefile("libicsest.so");
103    deletefile("libicsez0.so");
104    deletefile("loader.sce");
105    deletefile("cleaner.sce");
106
107    format(saveFormat(1, $:-1:1));
108
109endfunction
110
111demo_seros();
112clear demo_seros;
113