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