4de0f38226a5228b0145dffaa2eef03d62cc74bf
[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 libn  = ilib_for_link("icsenb","icsenb.o",[],"f")
14 nlink = link("./"+libn,["icsenb","icsenf"],"f")
15
16 exec("icse.contexte");
17
18 t0   = 0.d0;  // instant initial
19 tf   = 1.d0;  // instant final
20 dtf  = 0;
21 ermx = 1.d-6; // test d'arret absolu sur la valeur du second membre dans
22 // la resolution de l'etat
23 iu = [0,0,1]; //  iu   :indications sur la structure du controle
24 //    iu(1)=1 si l'etat initial depend du controle constant,0 sinon
25 //    iu(2)=1 si l'etat initial depend du controle variable,0 sinon
26 //    iu(3)=1 si le second membre depend du controle constant,0 sinon
27 nuc  = 1;     // nombre de parametres independants du temps
28 nuv  = 1;     // nombre de parametres dependants du temps
29 ilin = 0;     // indicateur de linearite :
30 // 0 pour un systeme non affine
31 // 1 pour un systeme affine dont la partie lineaire n'est pas autonome
32 // ilin=2 pour un systeme affine dont la partie lineaire est autonome
33 nti = 150;    //nombre de pas de temps correspondant a dti (premier pas de temps)
34 dti = tf/nti;
35 ntf = 00;     // nombre de pas de temps correspondant a dtf (second pas de temps)
36 // si l'on utilise un seul pas de temps,on doit prendre ntf=0
37 ny   = 4;     // dimension de l'etat a un instant donne
38 nea  = 0;     // nombre d'equations algebriques (eventuellement nul)
39 itmx = 10;    // nombre maximal d'iterations dans la resolution de
40 // l'equation d'etat discrete a un pas de temps donne
41 nex = 1;      // nombre d'experiences effectuees
42 nob = 3;      // dimension du vecteur des mesures pour une experience donnee
43 // en un instant donne
44 ntob  = 1;    // nombre d'instants de mesure pour une experience donnee
45 ntobi = 1;    // nombre d'instants de mesure correspondant a dti (premier
46 // pas de temps)
47
48 nu=nuc+nuv*(nti+ntf+1); // dimension du vecteur des parametres de controle
49
50 //  uc(1,nuc)          :controle constant
51 echtf = 2000;
52 uc    = [2500/echtf];
53
54 //  uv(1,nuv*(nti+ntf)):controle variable
55 //if nuv>0 then uv(1,nuv*(nti+ntf+1))=0; end;
56 alpha0 = .20704/.029244; legu = "alpha initial : ann. cz";     // annulation cz
57 alpha0 = 17.391;         legu = "alpha initial : finesse max"; // finesse maximum
58 legu = " navette americaine ."+legu;
59 if nuv>0 then uv=alpha0*ones(1,nuv*(nti+ntf+1)); end;
60
61 //  itu(1,nitu)        :tableau de travail entier reserve a
62 //                      l'utilisateur
63 itu = [0];
64
65 //  dtu(1,ndtu)        :tableau de travail double precision reserve
66 //                      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
110 //          (valeur arbitraire si iu(1) ou iu(2) est non nul)
111
112 y0=y0./dtu(1,21:24); // mise a l'echelle de y0
113
114 //  tob(1,ntob)        :instants de mesure (compatibilite avec ntob
115 //                      et ntobi)
116 tob       = [1];
117 binf      = -20*ones(1,nu); // borne inf des parametres
118 binf(1,1) = 2500/echtf;
119 bsup      = 40*ones(1,nu); // borne sup des parametres
120 bsup(1,1) = 4000/echtf;
121
122 obs(nob,ny) = 0; // matrice d'observation
123
124 //don=[762,            ..//vfin         1
125 //     -5*raddeg,      ..//gamma final  2
126 //     24384    ];     ..//zfin         3
127 don=[762,            ..
128 -5*raddeg,      ..
129 24384    ];
130 don   = don./dtu(1,21:23); // mise a l'echelle
131 nomf  = "icsenf";          // noms de subroutines de dynamique
132 legfb = " croissant ";
133
134 // changements pour calculer en temps retrograde
135 retro = 1;
136 if retro>0 then
137     legfb     = " retrograde ";
138     don1      = don;
139     don       = y0(1,1:3);
140     y0(1,1:3) = don1;
141     nomf      ="icsenb";
142 end
143
144 nap    = 20;  // nombre d'appels du simulateur
145 imp    = 2;   // niveau de debug pour optim
146 large  = 100; // taille de nu au dela de laquelle on choisit un optimiseur
147
148 // pour les problemes de grande taille (alg='gc' dans l'appel de optim)
149
150 exec("icseinit.sce");
151
152 [co,u,g,itv,dtv]=icse(u,nomf,nap,imp);