function[indices_image_discretisee]=predec_succes(etat,image)
//etat = vecteur ordonné
//image = vecteur
// indices_image_discretisee = liste
//indices_image_discretisee(1)(i)=j
// SSI etat(j) est le prédécesseur de image(i)
//indices_image_discretisee(2)(i)= j
// SSI etat(j) est le successeur de image(i)
//indices_image_discretisee(3)(i)=
//probabilité d'aller vers le prédécesseur de image(i)

cardinal_etat=prod(size(etat));
ind=dsearch(image,etat);
// image(i) \in [etat(ind(i)), etat(ind(i)+1)[ sauf si
// ind(i)=0, auquel cas image(i)<etat(1) ou image(i)>etat($)
//ind_nd=ind(image_ind);
image_ind_h=find(image>etat($));
ind(image_ind_h)=cardinal_etat;
// ind(i)=0 ssi image(i)<etat(1)
// ind(i)=n ssi image(i)>etat($)
image_ind_m=find(ind>0 & ind<cardinal_etat);
// indices du milieu : image(image_ind_m) \subset [etat(1),etat($)]
// ind(image_ind_m) \subset ]0,n[
image_ind_b=find(ind==0);
// indices du bas : image(image_ind_b) <etat(1)
image_ind_h=find(ind==cardinal_etat);
// indices du haut : image(image_ind_h) >etat($)

ind1=zeros(image);ind2=ind1;
proba=zeros(image);

ind1(image_ind_h)=ind(image_ind_h); 
// envoie les indices pour lesquels l'image de l'etat correspond est
// >etat($) vers cardinal_etat, dernier indice de etat
ind2(image_ind_h)=ind1(image_ind_h);
// prédécesseur = successeur = etat($)
proba(image_ind_h)=ones(image_ind_h);
// probabilité 1 d'aller vers prédécesseur = successeur

ind1(image_ind_m)=ind(image_ind_m) ;
ind2(image_ind_m)=1+ind(image_ind_m) ;
proba(image_ind_m)= ...
(etat(1+ind(image_ind_m))-image(image_ind_m))./ ...
(etat(1+ind(image_ind_m))-etat(ind(image_ind_m)));

ind1(image_ind_b)=ind(image_ind_b)+1; 
// +1 à cause de ind(image_ind_b) composé de 0
// envoie les indices pour lesquels l'image de l'etat correspond est
// <etat(1) vers 1, premier indice de etat
ind2(image_ind_b)=ind1(image_ind_b);
// prédécesseur = successeur = etat(1)
proba(image_ind_b)=ones(image_ind_b);
// probabilité 1 d'aller vers prédécesseur = successeur

indices_image_discretisee=list();
indices_image_discretisee(1)=ind1;
indices_image_discretisee(2)=ind2;
indices_image_discretisee(3)=proba;
endfunction


function[z]=egal(x,y)
z= bool2s(x==y)
endfunction


function[matrice_transition]=discretisation(etat,commande,dynamique_commandee)
matrice_transition=list()
//etat = vecteur ordonné
//commande = vecteur
//dynamique_commandee = fonction(x,u)
//matrice_transition = liste de matrices de transition 
//  sur les indices de etat
cardinal_etat=prod(size(etat));
for l=1:prod(size(commande))
//la liste matrice_transition est indicée par les indices du vecteur commande
image=dynamique_commandee(etat,commande(l));
// vecteur des images du vecteur etat
[indices_image_discretisee]=predec_succes(etat,image)
indices1=indices_image_discretisee(1);
indices2=indices_image_discretisee(2);
probabilites=indices_image_discretisee(3);

M1=zeros(cardinal_etat,cardinal_etat);
M2=zeros(M1);
for i=1:cardinal_etat,
  M1(i,indices1(i))=probabilites(i);
  M2(i,indices2(i))=1-probabilites(i);
end
matrice_transition(l)=M1+M2

//coincidmoins=feval(zmoins(1,:),etat,egal)
// coincidmoins(i,j)=1 ssi zmoins(1,i) == etat(j)
// ssi partant de etat(i), etat(j) est visité
//coincidplus=feval(zplus(1,:),etat,egal)
//coincidmoins=1-sign ( abs ( log( ( exp(-zmoins(1,:)') ) * exp(etat) )));
//coincidplus=1-sign ( abs ( log( ( exp(-zplus(1,:)') ) * exp(etat) ) ) );
//matrice_transition(l)= coincidmoins .* ( ones(etat)' * zmoins(2,:) )' +...
//coincidplus .* ( ones(etat)' * zplus(2,:) )' 
// affectation des probabilités
end
endfunction




function[b]=dyn_plant_control(x,v)
//b = total biomasse
//x = vegetative biomass
//v = control
     b=v*r*x^{power}
endfunction





function [valeur,feedback]=Bell_stoch(matrice_transition,...
cout_instantane,cout_final,minimaxi)
// algorithme de programmation dynamique stochastique
//     pour une chaîne de Markov sur \{1,2...,cardinal_etat\}
// matrice_transition = liste à cardinal_commande éléments
//     composée de matrices de dimension (cardinal_etat,cardinal_etat)
// cout_instantane = liste à cardinal_commande éléments
//     composée de matrices de dimension (cardinal_etat,horizon-1)
// cout_final = vecteur de dimension (cardinal_etat,1)
// valeur = fonction valeur de Bellman,
//     matrice de dimension (cardinal_etat,horizon)
// feedback = commande en feedback sur l'état et le temps
//     matrice de dimension (cardinal_etat,horizon-1)
// MINIMISATION si minimaxi=-1
// MAXIMISATION si minimaxi=1

ee=size(cout_instantane(1)); cardinal_etat=ee(1); 
cardinal_commande=size(cout_instantane);
hh=size(cout_instantane(1)); horizon = 1+hh(2);

valeur=0*ones(cardinal_etat,horizon);
// tableau ou l'on stocke la fonction de Bellman
// valeur(:,t) est un vecteur
// : représente ici le vecteur d'état

valeur(:,horizon)=cout_final;
// La fonction valeur au temps T est cout_final

feedback=0*ones(cardinal_etat,horizon-1);
// tableau ou l'on stocke le feedback u(x,t)

// Calcul rétrograde de la fonction de Bellman et
// de la commande optimale à la date 'temp' par l'équation de Bellman
// On calcule le coût pour la commande choisie connaissant la valeur de la
// fonction coût à la date suivante pour l'état obtenu

for temp=horizon:-1:2,
// équation rétrograde
// (attention, pour des questions d'indices, bien finir en :2)
loc=zeros(cardinal_etat,cardinal_commande);
// variable locale contenant les valeurs de la fonction à minimiser
// loc est une matrice
for j=1:cardinal_commande
loc(:,j)=matrice_transition(j)*valeur(:,temp)+cout_instantane(j)(:,temp-1);
end;

[mmn,jjn]=mini(loc,'c');
[mmx,jjx]=maxi(loc,'c');
// mm est l'extremum atteint
// jj est la commande qui réalise l'extremum

valeur(:,temp-1)=(1-minimaxi)/2 * mmn + (1+minimaxi)/2 * mmx;
// mm;
// coût optimal
feedback(:,temp-1)=(1-minimaxi)/2 * jjn + (1+minimaxi)/2 * jjx;
//jj;
// feedback optimal 
end
endfunction



function [z]=trajopt(matrice_transition,feedback,cout_instantane,...
                                 cout_final,etat_initial)
// z est une liste :
// z(1) est la trajectoire optimale obtenue partant de etat_initial
// z(2) est la trajectoire correspondante des commandes optimales
// z(3) est la trajectoire correspondante des coûts 
// etat_initial = un entier 
// pour le reste, voir la macro Bell_stoch

ee=size(cout_instantane(1)); cardinal_etat=ee(1); 
cardinal_commande=size(cout_instantane);
hh=size(cout_instantane(1)); horizon = 1+hh(2);

z=list();
x=etat_initial;
u=[];
c=[];

for j=1:horizon-1
// remplissage des trajectoires : commande, coût, état
  u=[u,feedback(x($),j)];
  c=[c,c($)+cout_instantane(u($))(x($),j)];
  mat_trans=full(matrice_transition(u($)));
  x=[x,grand(1,'markov',mat_trans,x($))];
end
  c=[c,c($)+cout_final(x($))];

z(1)=x;z(2)=u;z(3)=c;
endfunction