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