Discussion :

[jordan64011]

Bonjour à tous,
Je tente actuellement de créer un code qui permet d'optimiser une fonction de type :


Les deltaV, V0, a, b sont des constantes connues, u(t) est la fonction de heaviside et t est un vecteur de temps de longueur connue
Les paramètres à optimiser sont donc P1 et P2, des vecteur de longueur inconnue. Ce que je dois optimiser : le temps pendant lequel la pente de la fonction est à plat ou est supérieur à 0, tout en maintenant la valeur en ordonnée entre une valeur minimum Vmin et une maximum Vmax durant tout la durée donnée par le vecteur t. Comme sur cet exemple:



Dans ce cas, P1= [5000 11000] P2=[7000 13000] Vmin = 20 et Vmax = 120 (à peu près).

Je sais que l'optimisation de cette fonction par Matlab n'est pas utile (il suffit de retarder les valeurs des deux paramètres au maximum et le tour est joué) mais je n'ai pas fourni la totalité de l'équation. Si je comprends comment procéder, je pourrai appliquer la méthode à ce qui reste de l'équation.

J'ai déjà codé la fonction ObjFun qui reprend les deux sommations :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
function f_v = ObjFun(p1,p2)
global t conso Ttour Travit Tpit Reservoir
Ess=zeros(1,length(t));
Rem=zeros(1,length(t));
for i=1:length(p1)
    Ess=Ess+conso/Ttour.*(t-p1(i)).*heaviside(t-p1(i))+1/Travit.*(t-(p1(i)+Tpit)).*heaviside(t-(p1(i)+Tpit));
end
for i=1:length(p2)
    Rem=Rem+1/Travit.*(t-p2(i)).*heaviside(t-p2(i))+conso/Ttour.*(t-p2(i)).*heaviside(t-p2(i));
end
f_v=subplus(diff(Reservoir-conso/Ttour.*t+Ess-Rem));
end

Je présente quelques difficultés pour traduire les contraintres Vmin et Vmax et pour le choix de la routine d'optimisation à utiliser (entre fmincon, ga, ...)
Pouvez-vous m'aider ?

Merci d'avance.

[Antoane]

Bonjour,

je n'ai pas tout compris à la demande et ne suis pas spécialiste en optimisation, mais je doute que matlab soit nativement capable de choisir le nombre de paramètres à optimiser. Je pense qu'il faut faire manière détournée :
- en fixant ex ante la longueur des vecteurs P1 et P2 et en réalisant des optimisations pour différentes longueurs ;
- ou en fixant ex ante la longueur des vecteurs P1 et P2 à une "grande valeur" (supérieure à la longueur max acceptable comme résultat d'optimisation) et en utilisant un fonction d'optimisation capable d'optimiser une fonction contenant un paramètre discret qui caractériserait le nombre d'éléments de P1 et P2 à prendre en compte

[jordan64011]

Bonjour, tout d'abord merci pour ta réponse.
En effet mon explication n'est peut être pas la meilleure mais je présente déjà des difficultes pour m'expliquer à l'oral sur ce projet alors à l'écrit...

Il s'agit en réalité d'un code qui servira dans un programme de stratégie de course dans le domaine course endurance automobile.
L'équation que j'ai envoyé permet en réalité de tracer la courbe du volume d'essence en fonction du temps avec Delta(V) qui est la consommation par seconde, a la quantité d'essence remise par seconde et b un temps de parcours de la pitlane avant d'arriver a la pompe à essence. Nous oscillons donc entre la valeur Vmin et Vmax du reservoir d essence avec une courbe de consommation (apparait dès le départ ou apres un temps paramètré P2) et une courbe de ravitaillement (apres un temps P1). Comme je l'ai déjà dit il y aura d'autres paramètres à optimiser par après (changement de pilote,...) mais il me suffit de comprendre comment proceder juste avec le ravitaillement et la consommation pour complexifier par la suite. Voila j'espère avoir été un peu plus clair.

En ce qui concerne ton idée, j'y ai déjà un peu pensé mais me demande comment procéder. Est-ce que je prédefini des vecteur de zeros d'une certaine longueur supérieur au maximum possible en course et essaie d'optimiser chaque élément ? Qu'arrivera t'il alors dans le cas où je suis arrivé en fin de mon vecteur t (temps de course) ? Le programme laissera t-il un zéro tel qu'il est ?

Encore merci de te pencher sur mon post

[Antoane]

Bonjour,

La méthode consistant à fixer ex ante la longueur des vecteurs P1 et P2 et à réaliser des optimisations pour différentes longueurs est la plus couteuse (en temps de calcul) à mettre en œuvre mais la plus simple. Son implémentation est relativement aisée, c'est une optimisation classique. Le choix de l'algo d'optimisation pourrait cependant ne pas être trivial du fait de la non-dérivabilité de la fonction de coût.

Conseil d'ordre général : mieux vaut éviter les variables globales et passer les paramètre depuis la fonction d'optimisation. Par exemple :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
N = 12; % length of P1 and P2
p12_0 = ones(1, 2*N) ; % initial guess 
 
% build the structure containing all parameters
parameters = struct() ; 
parameters.t =t ; 
parameters.conso =conso ;
parameters.Ttour =Ttour ;
parameters.Travit = Travit ;
parameters.Tpit = Tpit ;
parameters.Reservoir = Reservoir ;
 
bf = fminsearch(@(x) ObjFun(parameters, x), p12_0) ;
 
function f_v = ObjFun(parameters, p12)
    % global t conso Ttour Travit Tpit Reservoir
    Ess=zeros(1,length(t));
    Rem=zeros(1,length(t));
    for i=1:length(p1)
        Ess=Ess+conso/Ttour.*(t-p1(i)).*heaviside(t-p1(i))+1/Travit.*(t-(p1(i)+Tpit)).*heaviside(t-(p1(i)+Tpit));
    end
    for i=1:length(p2)
        Rem=Rem+1/Travit.*(t-p2(i)).*heaviside(t-p2(i))+conso/Ttour.*(t-p2(i)).*heaviside(t-p2(i));
    end
    f_v=subplus(diff(Reservoir-conso/Ttour.*t+Ess-Rem));
end

L'autre technique que j'indiquais plus haut consiste à décider a priori qu'il n'y aura pas plus de N_max ravitaillement et donc que length(P1) < N_max + 1 et length(P2) < N_max + 1.
Dans ce cas, le vecteur de paramètres à optimiser serait de la forme : [P1, P2, length_P1, length_P2], initialisé comme un vecteur de longueur 2*N_max + 2.
Le code d'optimisation serait, en se basant sur celui que tu donnes plus haut :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function f_v = ObjFun(parameters, p12)
    p1 = p12(1:parameters.N_max) ;
    p2 = p12(parameters.N_max+1:2*parameters.N_max) ;
    length_p1 = p12(2*parameters.N_max + 1) ;
    length_p2 = p12(2*parameters.N_max + 2) ;
    % global t conso Ttour Travit Tpit Reservoir
    Ess=zeros(1,length(t));
    Rem=zeros(1,length(t));
    for i=1:length_p1
        Ess=Ess+conso/Ttour.*(t-p1(i)).*heaviside(t-p1(i))+1/Travit.*(t-(p1(i)+Tpit)).*heaviside(t-(p1(i)+Tpit));
    end
    for i=1:length_p2
        Rem=Rem+1/Travit.*(t-p2(i)).*heaviside(t-p2(i))+conso/Ttour.*(t-p2(i)).*heaviside(t-p2(i));
    end
    f_v=subplus(diff(Reservoir-conso/Ttour.*t+Ess-Rem));
end

Les mêmes questions quant au choix de l'algo d'optimisation se poseront, avec en plus la difficulté liée au fait que deux éléments du vecteur d'optimisation sont entières... Il existe des algo le permettant nativement, mais je ne saurait pas en dire plus.

[jordan64011]

Bonsoir,

En effet, cette méthode a l'air de fonctionner (j'ai utilisé la deuxième avec Nmax), ce n'est pas encore ce que j'espère comme réponse mais en tout cas ça avance bien. Je ne rencontre plus d'erreurs lors de la routine d'optimisation et arrive à resortir une valeur pour le vecteur "optim".

Je te remercie déjà énormément pour ton aide !

Concernant l'algorithme d'optimisation je pense me tourner vers le "genetic algorithm" (ga) car celui-ci permet de travailler avec des nombres entiers dans les variables d'optimisation et que ce soit pour mes paramètres p1/p2 ou pour les longueurs de ceux-ci (le tout rassemblé en un vecteur) je ne travaille qu'avec des entiers. Je ne sais pas si cette raison est vraiment la meilleure quant au choix du type d'algorithme mais je me dis qu'en supprimant déjà toutes les valeurs décimales de l'optimisation je gagnerai en temps de calcul quelque part (ce qui n'est pas négligeable étant donné que, pour rappel, ce code sert dans le milieu de la course automobile là où tout doit être le plus rapide possible).

Je rencontre désormais un problème avec le système de contrainte.. je ne comprends pas comment contraindre l'image de ma fonction entre les deux valeurs.

Voici ma fonction de contrainte où j'essaie de contraindre entre 6 et Reservoir (=120) ma fonction Vt :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
function [c,ceq]=constraintsFun(Values,Param)
    Reservoir=Values.Reservoir;
    vol=Vt(Values,Param);
    for i=1:length(vol)
        c=[6-vol(i),vol(i)-Reservoir];
        ceq=[];
    end
end

Au cas où voici Vt : (petit changement par rapport à la réponse précédente, p1 et p2 sont, pour le moment, de même longueur, c'est pourquoi le for va de 1 à length_p1 dans les deux cas, ce sera modifié plus tard)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
function fun1 = Vt(Values,Param)
    t=Values.t;
    conso=Values.conso;
    Ttour=Values.Ttour;
    Travit=Values.Travit;
    Tpit=Values.Tpit;
    Reservoir=Values.Reservoir;
    Nmax=Values.Nmax;
    p1=Param(1:Nmax);
    p2=Param(Values.Nmax+1:2*Nmax);
    length_p1=Param(2*Nmax+1);
    Ess=zeros(1,length(t));
    Rem=zeros(1,length(t));
    for i=1:length_p1
        if p1(i)>0
        Ess=Ess+conso/Ttour.*(t-p1(i)).*heaviside(t-p1(i))+1/Travit.*(t-(p1(i)+Tpit)).*heaviside(t-(p1(i)+Tpit));
        end
    end
    for i=1:length_p1
        if p2(i)>0
        Rem=Rem+1/Travit.*(t-p2(i)).*heaviside(t-p2(i))+conso/Ttour.*(t-p2(i)).*heaviside(t-p2(i));
        end
    end
    fun1=Reservoir-conso/Ttour.*t+Ess-Rem;
end

Et j'ai très légèrement modifié ce que je cherche à minimiser, j'ai créé une variable "Losttime" dans ma fonction objectif :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
function Losttime = ObjFun(Values,Param)
    t=Values.t;
    conso=Values.conso;
    Ttour=Values.Ttour;
    Travit=Values.Travit;
    Tpit=Values.Tpit;
    Reservoir=Values.Reservoir;
    Nmax=Values.Nmax;
    p1=Param(1:Nmax);
    p2=Param(Values.Nmax+1:2*Nmax);
    length_p1=Param(2*Nmax+1);
    Ess=zeros(1,length(t));
    Rem=zeros(1,length(t));
    for i=1:length_p1
        if p1(i)>0
            Ess=Ess+conso/Ttour.*(t-p1(i)).*heaviside(t-p1(i))+1/Travit.*(t-(p1(i)+Tpit)).*heaviside(t-(p1(i)+Tpit));
        end
    end
    for i=1:length_p1
        if p2(i)>0
            Rem=Rem+1/Travit.*(t-p2(i)).*heaviside(t-p2(i))+conso/Ttour.*(t-p2(i)).*heaviside(t-p2(i));
        end
    end
    fun1=Reservoir-conso/Ttour.*t+Ess-Rem;
    fun =diff(fun1);
    Losttime=0;
    for i=1:length(fun)
        if fun(i)>(-conso/Ttour)
            Losttime=Losttime+1;
        end
    end
end

Voilà, concrètement ObjFun est la même fonction que Vt ajouté de la petite partie pour le calcul de Losttime qui est le temps perdu (à minimiser) quand la voiture rentre dans les stands pour ravitailler en essence.

Pour l'instant le plot de ma fonction Vt avec les paramètres sortis de l'optimisation ne vérifient pas du tout les contraintes que j'ai posées (Vt descend loin en négatif ou au contraire loin au dessus de 120, ma valeur maximale). Est-ce que ma méthode de contrainte d'inégalité non linéaire semble correcte ?

Merci d'avance.

[Antoane]

Bonjour,

Je pense que tu devrais commencer par un optimisation plus simple pour te familiariser avec la fonction ga et les différentes contraintes applicables. Cela t'éviterait de devoir gérer simultanément les problèmes liés à :
- l'utilisation d'un algo d'optim génétique
- l'optimisation dans plusieurs directions
- la définition de fonctions de coût pas simple et alambiquées.

Comme indiqué plus haut, je ne suis pas spécialiste en optim, ca va donc être de plus en plus difficile de te répondre . Avec un peu de chance quelqu'un de plus compétant passera par là...

Pour l'instant le plot de ma fonction Vt avec les paramètres sortis de l'optimisation ne vérifient pas du tout les contraintes que j'ai posées (Vt descend loin en négatif ou au contraire loin au dessus de 120, ma valeur maximale). Est-ce que ma méthode de contrainte d'inégalité non linéaire semble correcte ?

Je n'ai pas étudié le code en détails ; as-tu vérifié la valeur de "options.ConstraintTolerance" ? Que dit l'optimiseur après avoir terminé ?

Quelques détails additionnels :

Je ne sais pas si cette raison est vraiment la meilleure quant au choix du type d'algorithme mais je me dis qu'en supprimant déjà toutes les valeurs décimales de l'optimisation je gagnerai en temps de calcul quelque part

C'est possible, mais ne me parait pas évident a priori. Le critère de convergence options.FunctionTolerance (https://fr.mathworks.com/help/gads/g...7-d260b7a43e26) permet de chercher une solution plus ou moins "exacte".

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
function [c,ceq]=constraintsFun(Values,Param)
    Reservoir=Values.Reservoir;
    vol=Vt(Values,Param);
    for i=1:length(vol)
        c=[6-vol(i),vol(i)-Reservoir];
        ceq=[];
    end
end

devrait pouvoir se réécrire :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
function [c,ceq]=constraintsFun(Values,Param)
    Reservoir=Values.Reservoir;
    vol=Vt(Values,Param);
 
    c=[6-vol, vol-Reservoir];
    ceq=[];
end