* Bug #8190 fixed - Optimization: Fixed ICSE demos
[scilab.git] / scilab / modules / optimization / demos / icse / navet.sce
1 //
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) ????-2008 - INRIA
4 // Copyright (C) 2010 - DIGITEO
5 //
6 // This file is distributed under the same license as the Scilab package.
7 //
8
9 //                navet.bas : demo de icse
10 //  calcul trajectoire optimale de rentree d'une navette spatiale
11 //  *************************************************************
12
13 function demo_navet()
14
15     cd "SCI/modules/optimization/demos/icse/"
16
17     libn  = ilib_for_link("icsenb", "icsenb.f",[], "f", "", "", "", "-L"+SCI+"/modules/optimization/.libs/ -lscioptimization")
18     nlink = link("./"+libn, ["icsenb", "icsenf"], "f")
19     libn  = ilib_for_link("icsez0", "icsez0.f",[], "f")
20     nlink = link("./"+libn, "icsez0", "f")
21
22     exec("icsecontexte.sce");
23
24     t0   = 0.d0;  // instant initial
25     tf   = 1.d0;  // instant final
26     dtf  = 0;
27     ermx = 1.d-6; // test d'arret absolu sur la valeur du second membre dans la resolution de l'etat
28     iu = [0,0,1]; //  iu   :indications sur la structure du controle
29     //    iu(1)=1 si l'etat initial depend du controle constant,0 sinon
30     //    iu(2)=1 si l'etat initial depend du controle variable,0 sinon
31     //    iu(3)=1 si le second membre depend du controle constant,0 sinon
32     nuc  = 1;     // nombre de parametres independants du temps
33     nuv  = 1;     // nombre de parametres dependants du temps
34     ilin = 0;     // indicateur de linearite :
35     // 0 pour un systeme non affine
36     // 1 pour un systeme affine dont la partie lineaire n'est pas autonome
37     // 2 pour un systeme affine dont la partie lineaire est autonome
38     nti = 150;    //nombre de pas de temps correspondant a dti (premier pas de temps)
39     dti = tf/nti;
40     ntf = 00;     // nombre de pas de temps correspondant a dtf (second pas de temps)
41     // si l'on utilise un seul pas de temps,on doit prendre ntf=0
42     ny   = 4;     // dimension de l'etat a un instant donne
43     nea  = 0;     // nombre d'equations algebriques (eventuellement nul)
44     itmx = 10;    // nombre maximal d'iterations dans la resolution de l'equation d'etat discrete a un pas de temps donne
45     nex = 1;      // nombre d'experiences effectuees
46     nob = 3;      // dimension du vecteur des mesures pour une experience donnee en un instant donne
47     ntob  = 1;    // nombre d'instants de mesure pour une experience donnee
48     ntobi = 1;    // nombre d'instants de mesure correspondant a dti (premier pas de temps)
49
50     nu=nuc+nuv*(nti+ntf+1); // dimension du vecteur des parametres de controle
51
52     //  uc(1,nuc)          :controle constant
53     echtf = 2000;
54     uc    = [2500/echtf];
55
56     //  uv(1,nuv*(nti+ntf)):controle variable
57     //if nuv>0 then uv(1,nuv*(nti+ntf+1))=0; end;
58     alpha0 = .20704/.029244; legu = "alpha initial : ann. cz";     // annulation cz
59     alpha0 = 17.391;         legu = "alpha initial : finesse max"; // finesse maximum
60     legu = " navette americaine ."+legu;
61     if nuv>0 then uv=alpha0*ones(1,nuv*(nti+ntf+1)); end;
62
63     //  itu(1,nitu)        :tableau de travail entier reserve a l'utilisateur
64     itu = [0];
65
66     //  dtu(1,ndtu)        :tableau de travail double precision reserve a l'utilisateur
67     raddeg = %pi/180;
68
69     //dtu=[ 249.9,         ..//s      1
70     //      .078540,       ..//cx0    2
71     //     -.0061592,      ..//cx1    3
72     //      .621408e-3,    ..//cx2    4
73     //     -.207040,       ..//cz0    5
74     //      .029244,       ..//cz1    6
75     //       83388,        ..//zm     7
76     //     3.9860119e14,   ..//zmu    8
77     //     6378166,        ..//rt     9
78     //     1.2,            ..//ro0   10
79     //     6700,           ..//h     11
80     //     raddeg,         ..//      12
81     //     echtf,          ..//echtf 13
82     //     0,0,0,0,0,0,0,  ..// inutilises    14 a 20
83     //     1000,           ..//mise a echelle v    21
84     //     1,              ..//mise a echelle gam  22
85     //     1.e5,           ..//mise a echelle z    23
86     //     1,           ..//mise a echelle l       24
87     //     1.e6      ];      //cpenal              25
88
89     dtu=[ 249.9,         ..
90     .078540,       ..
91     -.0061592,      ..
92     .621408e-3,    ..
93     -.207040,       ..
94     .029244,       ..
95     83388,        ..
96     3.9860119e14,   ..
97     6378166,        ..
98     1.2,            ..
99     6700,           ..
100     raddeg,         ..
101     echtf,          ..
102     0,0,0,0,0,0,0,  ..
103     1000,           ..
104     1,              ..
105     1.e5,           ..
106     1,           ..
107     1.e6      ];      //cpenal              25
108
109     y0=[7803, -1*raddeg, 121920, 0]; // etat initial (valeur arbitraire si iu(1) ou iu(2) est non nul)
110
111     y0=y0./dtu(1,21:24); // mise a l'echelle de y0
112
113     //  tob(1,ntob)        :instants de mesure (compatibilite avec ntob et ntobi)
114     tob       = [1];
115     binf      = -20*ones(1,nu); // borne inf des parametres
116     binf(1,1) = 2500/echtf;
117     bsup      = 40*ones(1,nu); // borne sup des parametres
118     bsup(1,1) = 4000/echtf;
119
120     obs(nob,ny) = 0; // matrice d'observation
121
122     //don=[762,            ..//vfin         1
123     //     -5*raddeg,      ..//gamma final  2
124     //     24384    ];     ..//zfin         3
125     don=[762,            ..
126     -5*raddeg,      ..
127     24384    ];
128     don   = don./dtu(1,21:23); // mise a l'echelle
129     nomf  = "icsenf";          // noms de subroutines de dynamique
130     legfb = " croissant ";
131
132     // changements pour calculer en temps retrograde
133     retro = 1;
134     if retro>0 then
135         legfb     = " retrograde ";
136         don1      = don;
137         don       = y0(1,1:3);
138         y0(1,1:3) = don1;
139         nomf      ="icsenb";
140     end
141
142     nap    = 20;  // nombre d'appels du simulateur
143     imp    = 2;   // niveau de debug pour optim
144     large  = 100; // taille de nu au dela de laquelle on choisit un optimiseur
145
146     // pour les problemes de grande taille (alg='gc' dans l'appel de optim)
147
148     saveFormat = format();
149
150     exec("icseinit.sce");
151
152     exec("icse.sci");
153     exec("icsegen.sci");
154
155     [co, u, g, itv, dtv] = icse(u, nomf, nap, imp);
156
157     disp(u', "Best value:")
158     disp(co, "Final cost:")
159
160     deletefile("libicsenb.so");
161     deletefile("libicsez0.so");
162     deletefile("loader.sce");
163     deletefile("cleaner.sce");
164
165     format(saveFormat(1, $:-1:1));
166
167 endfunction
168
169 demo_navet();
170 clear demo_navet;