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);