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