UNIVERSITÉ ABOUBEKR BELKAID TLEMCEN
Faculté de Technologie
Département Génie Électrique et Électronique
Commande Optimale
Problème du Brachistochrone
On se propose de trouver la courbe y(x), xA ≤ x ≤ xB dans le plan vertical (x, y) joignant les points
A = (xA , xB ) et B = (xB , yB ), xA < xB tel que un point matériel glissant le long de la courbe y(x) sans
frottement de A vers B sous l’effet de la gravité et avec une vitesse initiale vA ≥ 0 atteint B en temps minimum.
Considérer le cas où xA = yA = 0. Montrer que le problème peut être formulé mathématiquement comme
le problème de calcul de variations suivant :
Rx q ẏ(x)2
minimum de J(y) = 0 B v21+ −2gy(x)
dx (1)
A
Sous contraintes y ∈ D := {y ∈ C 1 [0, xB ] : y(0) = 0, y(xB ) = yB }
g représente l’accélération de la pesanteur. On veut calculer une solution approximative à ce problème en
paramétrisant la courbe y(x) en utilisant des polynômes de Lagrange d’ordre N ≥ 1 :
N N
X Y x − Xj
Y (ỹ; x) = ỹk
j=0
Xk − Xj
k=0
j6=k
Où ỹ = (ỹ0, · · · , ỹN )T ; et X0 , · · · , XN sont N + 1 points dans [0, XB ]. On considŕera des points ègalement
espacés :
xB
X0 = 0, Xk := Xk−1 + .
N
En utilisant cette paramétrisation, le problème devient
Rx q 2
minimum de f (ỹ) := 0 B v21+Z(ỹ;x)
−2gY (ỹ;x)
dx (2)
A
sous contraintes Y (ỹ; 0) = 0
Y (ỹ; xB ) = yB
2
vA
ỹk ≤ 2g , k = 0, · · · , N
Avec
k
Xk := xB , k = 0, · · · , N (3)
N
N N
X Y x − Xj
Y (ỹ; x) := ỹk
j=0
Xk − Xj
k=0
j6=k
N N N
X X 1 Y x − Xj
Z(ỹ; x) := ỹk
i=0
Xk − Xi j=0 Xk − Xj
k=0
i6=k j6=i,k
(4)
a) Trouver une solution analytique pour ce problème
b) Trouver une solution optimale à ce problème en utilisant la fonction fmincom de MATLAB.
– Considérer N=10, et prendre xB = 1,yB = −0.5 et vA = 1.
– Prendre comme condition initiale ỹk0 = k
N
yB , pour k = 0, · · · N.
– Dans la fonction fmincon utiliser l’algorithme SQP, avec une mise à jour Quasi Newton et recherche linéaire.
– Choisir une tolérance pour la solution, la fonction et les contraintes égale à 10−7
– utiliser la fonction quadl pour calculer la valeur du coût avec une tolérance égale à 10−7
c) Representer sur un même graphe la courbe analytique et la courbe approchée.
1
• Donner le vecteur solution et le temps mis par le mobile pour aller de A à B.
• Trouver le temps que mettra le mobile pour aller de A à B en suivant une ligne droite.
Une solution possible pour le calcul de la solution approchée est la suivante :
%fonction brachit_main
clear all
%Nombre de points
N=11;
%parametres
xB=1;
yB=-.5;
vA=1;
g=10;
for k=1:N
y0(k)=(k-1.)/(N-1.)*yB;
yu(k)=vA^2/(2*g);
xk(k)=(k-1.)/(N-1.)*xB;
end
% Options d’Optimisation
options=optimset(’GradObj’,’off’,’GradConstr’,’off’,’Display’,’iter’,...
’LargeScale’,’off’,’HessUpdate’,’bfgs’,’Diagnostics’,’on’,...
’TolX’,1e-7,’TolFun’,1e-7,’TolCon’,1e-7,’MaxIter’,1000,...
’MaxFunEval’,1000,’DiffMinChange’,1e-6);
%Résolution du problème Nl
[yopt,fopt,iout]=fmincon( @brachistFun,y0,[],[],...
[],[],[],yu, @brachistCtr,...
options);
%Representation des resultats
clf
figure(1)
x=[0:.01:xB];
[dum P]=size(x);
for k=1:P
y(k)=lagrange(yopt,xk,x(k));
end
plot(x,y,’r’);
//////////////////////////////////////////////////////////////////////////////////////////
% Fonction à minimiser
% Cout Brachistochrone
%f=int_0^xB sqrt((1-dy^2)/(vA^2/(2*g)-y)dx
function [f]=brachistFun(yk)
%parameters
f=0.;
xB=1;
qtol=1e-7;
disp=0;
f=f+quadl(@brachistInt,0,xB,qtol,disp,yk);
end
/////////////////////////////////////////////////////////////////////////////////////////
% Integrand Brachistochrone
% l=sqrt((1_dy^2)/(vA^2/(2*g)-y)
function l=brachistInt(x,yk)
N=11;
%parametres
xB=1;
yB=-.5;
vA=1;
g=10;
for k=1:N
2
xk(k)=(k-1.)/(N-1.)*xB;
end
[dum P]=size(x);
for i=1:P
l(i)=sqrt((1+dlagrange(yk,xk,x(i))^2) ...
/(vA^2-2.*g*lagrange(yk,xk,x(i))));
end
end
///////////////////////////////////////////////////////////////////////////////////////////
y=lagrange(wk,xk,x)
% arguments
% wk : vecteur ligne de dimension N+1 contenant les coefficients
% d’interpolation
%xk : vecteur ligne de dimension N+1 contenant lees points d’interpolation
%xk(1)<....<xk(N+1)
% x les points en lesquels le polynome de lagrange est évalué
% Sortie
% Valeur du polynome de Lagrange en x
function y=lagrange(wk,xk,x)
N=11;
%parametres
xB=1;
yB=-.5;
vA=1;
g=10;
for k=1:N
y0(k)=(k-1.)/(N-1.)*yB;
yu(k)=vA^2/(2*g);
end
[dum,N]=size(wk);
y=0.;
for k=1:N
produit=wk(k);
for j=1:N
if j ~=k
produit=produit*(x-xk(j))/(xk(k)-xk(j));
end
end
y=y+produit;
end
end
/////////////////////////////////////////////////////////////////////////////////////////////
%dy=dlagrange(wk,xk,x)
%arguments
%wk : vecteur ligne de dimension N+1 contenant les coefficients
%d’interpolation
%xk : vecteur ligne de dimension N+1 contenant lees points d’interpolation
%xk(1)<....<xk(N+1)
% x les points en lesquels le polynome de lagrange est évalué
% Sortie
% y :Valeur du polynome de Lagrange en x
function dy=dlagrange(wk,xk,x)
N=11;
%parametres
xB=1;
yB=-.5;
vA=1;
g=10;
for k=1:N
3
y0(k)=(k-1.)/(N-1.)*yB;
yu(k)=vA^2/(2*g);
end
[dum,N]=size(wk);
dy=0.;
for k=1:N
somme=0. ;
for i=1:N
if i ~=k
produit=1.;
for j=1:N
if j ~=k && j ~=i
produit=produit*(x-xk(j))/(xk(k)-xk(j));
end
end
somme=somme +produit/(xk(k)-xk(i));
end
end
dy=dy+wk(k)*somme;
end
end
//////////////////////////////////////////////////////////////////////////////////////////
%points terminaux du brachistochrone
% geq(1)=y(0)-0;
% geq(2)=y(xB)-yB
function [gin,geq]=brachistCtr(yk)
N=11;
%parametres
xB=1;
yB=-.5;
vA=1;
g=10;
for k=1:N
xk(k)=(k-1.)/(N-1.)*xB;
end
gin=[];
geq(1)=lagrange(yk,xk,0.);
geq(2)=lagrange(yk,xk,xB)-yB;
end