← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to macros/tdcs/mineInit.sci

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
macros/tdcs/mineInit.sci
 
1
function []=mineInit()
 
2
// Macro qui initialise les donnees du
 
3
// td3
 
4
//!
 
5
// Copyright INRIA
 
6
n1=30
 
7
n2=40;
 
8
te=n2/2;
 
9
xe=n1/2;
 
10
k0=1;
 
11
kc=0.01;
 
12
[n1,n2,te,xe,k0,kc]=resume(...
 
13
        n1,n2,te,xe,k0,kc);
 
14
 
 
15
function [cout,feed]=mine(n1,n2,uvect)
 
16
//[cout,feed]=mine(n1,n2,uvect)
 
17
// extraction optimale d'un minerai d'une mine a ciel ouvert
 
18
// en avancant progressivement (k=1,2,...,n2) et en
 
19
// prelevant a l'abcisse k+1 la tranche de profondeur x(k+1)
 
20
// par une commande u a partir de la profondeur x(k)
 
21
// a l'abcisse k.
 
22
//
 
23
// la resolution du probleme se fait par programmation dynamique
 
24
// en calculant la commande optimale maximisant le critere :
 
25
//
 
26
//   -- n2-1
 
27
//   \
 
28
//   /        f(x(k),k) + V_F(x,n2)
 
29
//   -- k=1
 
30
//  avec : x(k+1)=x(k) + u
 
31
//  x(k) est la profondeur de la mine a l'abcisse k (x=1
 
32
//      est le niveau du sol)
 
33
// la fonction gain instantane f(i,k) represente le benefice
 
34
//      si on creuse la profondeur i a l'abcisse k
 
35
//  V_F(i,n2) est un cout final destine a forcer l'etat final
 
36
//       a valoir 1 (pour sortir de la mine ...)
 
37
//       V_F(i,n2) vaut -10000 en tout les points i\ne1 et 0
 
38
//       pour i=1
 
39
//
 
40
// le programme mine necessite de connaitre
 
41
//    n1    : l'etat est discretise en n1 points
 
42
//            1 c'est l'air 
 
43
//            2:n1+1 du sol au fond 
 
44
//            n1+2 zone interdite
 
45
//    n2    : nombre d'etapes 
 
46
//    uvect : vecteur ligne des valeurs discretes que peut
 
47
//          prendre u (3 valeurs, on monte, on descend ou on avance)
 
48
// et le programme mine retourne alors deux matrices
 
49
//    cout(n1,n2) : valeurs de la fonction de Bellman
 
50
//    feed(n1,n2) : valeurs de la commande a appliquer
 
51
//          lorsque l'etat varie de 1 a n1 et
 
52
//          l'etape de 1 a n2.
 
53
//    n1    : l'etat est discretise en n1 points
 
54
//            on rajoute dans le calcul deux couches fictive 
 
55
//            en 1 l'air en n1+2 le fond de la mine 
 
56
//            1 c'est l'air 
 
57
//            2:n1+1 du sol au fond 
 
58
//            n1+2 zone interdite
 
59
//!
 
60
usize=prod(size(uvect))
 
61
xgr=1:(n1+2)
 
62
// tableau ou l'on stocke la fonction de Bellman
 
63
cout=0*ones(n1+2,n2);
 
64
// tableau ou l'on stocke le feedback u(x,t)
 
65
feed=0*ones(n1,n2);
 
66
// calul de la fonction de Bellman au  temps final
 
67
penal=10000;
 
68
// Le cout final est le cout au point de sortie de la zone d'extraction 
 
69
// on veut evidement que ce point de sortie soit le sol 
 
70
// le cout est donc 0 au niveau du sol et -inf ailleurs 
 
71
cout(:,n2)=-penal*ones(n1+2,1);
 
72
cout(2,n2)=0;
 
73
// calcul retrograde de la fonction de Bellman et
 
74
// du controle optimal au temps temp par l'equation de Bellman
 
75
for temp=n2:-1:2,
 
76
  loc=list();
 
77
  for i=1:usize,
 
78
    newx=mini(maxi(xgr+uvect(i)*ones(xgr),1*ones(xgr)),(n1+2)*ones(xgr)),
 
79
    loc(i)=cout(newx,temp)+ff_o(xgr,temp-1),
 
80
  end;
 
81
  [mm,kk]=maxi(loc),
 
82
  cout(xgr,temp-1)=mm;
 
83
  cout(1,temp-1)=-penal;
 
84
  cout(n1+2,temp-1)=-penal;
 
85
  feed(xgr,temp-1)=uvect(kk)';
 
86
end
 
87
 
 
88
function [y]=ff_o(x,t)
 
89
//[y]=ff_o(x,t)
 
90
// gain instantane apparaissant dans le critere du
 
91
// programme mine.
 
92
//
 
93
// pour des raisons techniques, l'argument x doit
 
94
// etre un vecteur colonne et donc egalement la sortie y.
 
95
// en sortie y=[ ff_0(x(1),t),...ff_o(x(n),t)];
 
96
//!
 
97
xxloc=ones(x);
 
98
y=k0*(1-(t-te)**2/(n2**2))*xxloc - (x-xe*xxloc)**2/(n1**2) -kc*(x-1*xxloc)
 
99
y=y';
 
100
y(1:2)=[0;0]
 
101
 
 
102
function []=showcost(n1,n2)
 
103
//[]=showcost(n1,n2)
 
104
// Montre en 3d la fonction de gain instantanee (ff)
 
105
// x: profondeur (n1)
 
106
// y: abscisse  (n2)
 
107
// en z : valeur de ff_o(x,t)
 
108
//!
 
109
[lhs,rhs]=argn(0)
 
110
m=[];
 
111
for i=1:n2,m=[m,ff_o(1:n1,i)],end
 
112
contour(1:n2,1:n1,m',10,0,0,' ',[2,2,0])
 
113
 
 
114
function []=trajopt(feed)
 
115
//[]=trajopt(feed)
 
116
// feed est la matrice de feedback calculee par mine
 
117
// trajopt calcule et dessine la trajectoire et le controle
 
118
// optimaux pour un point de depart (1,1)
 
119
//!
 
120
[n1,n2]=size(feed)
 
121
xopt=0*ones(1,n2)
 
122
uopt=0*ones(1,n2+1)
 
123
xopt(1)=2;
 
124
for i=2:(n2+1),xopt(i)=feed(xopt(i-1),i-1)+xopt(i-1),
 
125
             uopt(i-1)=feed(xopt(i-1),i-1),end
 
126
plot2d2("gnn",[1:(n2+1)]',[-xopt]',[2],"111",...
 
127
       "trajectoire optimale",...
 
128
       [1,-n1,n2+1,2]);
 
129
plot2d2("gnn",[1:(n2)]',uopt(1:n2)',[1,-4],"111",...
 
130
       "commande optimale",...
 
131
       [1,-n1,n2+1,2]);
 
132
 
 
133
 
 
134