TP Optimisation : Étude de la fonction de Rosenbrock

22 juin 2019 0 Par Edouard

Enoncé :

1 TP1
1.1 Exercice 1
Source d’aide : http://kmdb.pagesperso-orange.fr/_src/_python/_formation_2010/python_formation_exemples_optimisation.html
On définit la fonction de Rosenbrock, également appelée Rosenbrock banana, par

1.1.1 Etude théorique
(1) Trouver les points critiques de f et démontrer que f admet un unique minimum global qu’elle atteint en .
(2) Déterminer la matrice .
1.1.2 Etude numérique
(1) Représenter la fonction de Rosenbrock sur le pavé .
(2) On souhaite comparer deux méthodes bien connues de minimisation. La méthode de gradient à pas fixe, puis à pas optimal. Pour comparer les deux méthodes, on choisit comme initialisation des algorithmes le point .
a) Méthode du gradient à pas fixe.
Écrire un programme « PasFixe » mettant en œuvre la méthode de gradient à pas fixe. Ce programme prendra comme arguments et , où désigne l’initialisation de la méthode et , le pas de la méthode.
b) Méthode du gradient à pas optimal.
Écrire une fonction « pasOptimal », des deux variables qui renvoie
la quantité .
Écrire un programme mettant en œuvre la méthode du gradient à pas optimal pour
Minimiser . On pourra initialement effectuer la recherche du pas optimal. On utilisera un double critère d’arrêt : un critère portant sur le nombre total d’itérations de la méthode, et un autre critère relatif à la précision de la méthode
c) Comparaison des méthodes.
i. Sur la même figure que précédemment, représenter les itérés obtenus par la méthode du gradient à pas fixe.
ii. Proposer une explication du phénomène observé pour la méthode du gradient à pas optimal.
(3) À votre avis, comment évoluent les graphes de convergence si l’on remplace la fonction ci-dessus par la fonction et que l’on fait grandir le paramètre positif ?

 

1.2 Exercice 2
L’énergie rayonnante d’un corps noir dans l’intervalle d’émission , par unité de surface et de temps, est appelée émittance monochromatique maximale du corps noir et est notée . Sa valeur est donnée par la loi de Planck

Les constantes dans cette relation sont identifiées par :
• : vitesse de la lumière dans le vide.
• : constante de Planck.
• : constante de Boltzmann.
• : longueur d’onde (m)
• : température absolue de la surface du corps noir (K).
• : indice de réfraction du milieu (ici le vide).
(1) Tracer sur un même graphique la fonction pour les valeurs de : 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800. Associer chaque courbé tracée à la valeur de T correspondante. On prendra .
(2) On souhaite trouver la valeur de qui maximise l’émittance monochromatique pour une température de surface T donnée. À quelle contrainte est-on soumis si l’on souhaite utiliser la méthode de la section dorée ?
Programmer alors cette méthode pour déterminer suivant les différentes valeurs de T.

Solution :

TP1_Roge_Poupa

 

Code :

import math as mp
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#Etude numerique

#Question 1

def Rosenbrock(X):
return ((X[0]-1)**2 + 10*(X[0]**2-X[1])**2)

#Creation du graphique en 3D de la fonction de Rosenbrock

def graph_Rosenbrock():

ax = Axes3D(plt.figure())
X = np.linspace(-10,10,50)
Y = np.linspace(-10,10,50)
X, Y = np.meshgrid(X, Y)
Z = Rosenbrock([X,Y])
ax.plot_surface(X, Y, Z)
plt.title(« Fonction de Rosenbrock »)
plt.xlabel(« X »)
plt.ylabel(« Y »)
plt.legend()
plt.show()

#Question 2
#Question 2a) : Methode du gradient a pas fixe

def gradient_f(X):
return([2*X[0]-2+40*X[0]**3-40*X[0]*X[1],-20*X[0]**2+20*X[1]])

def norme(u):
x=u[0]
y=u[1]
return mp.sqrt(x**2+y**2)

# l enonce propose la fonction PasFixe avec seulement pour argument x (l initialisation) et beta le pas
#Par consequent kmax et epsilon (l erreur acceptee) seront donnes dans le programme

def PasFixe(x,beta):
k=0
x= [x]
kmax=10000 #tres eleve juste pour ne pas bloquer la resolution mais eviter une boucle infinie ou tres longue
epsilon=10**(-3)

while norme(gradient_f(x[k]))>=epsilon and k<=kmax:
a=x[k][0]+beta*(-gradient_f(x[k])[0])
b=x[k][1]+beta*(-gradient_f(x[k])[1])
A=(a,b)
x.append(A)
k+=1

return(k,x[-1],norme(gradient_f(x[-1])),x)

PF = PasFixe([0,1],0.01)
print (« PasFixe : [x,y] = « ,PF[0],’ et k = ‘,PF[1], ‘ et Epsilon = ‘,PF[2])

#Question 2 b : Methode du gradient a pas optimal

#Fonction a minimiser

def Omega (X,alpha):
return ([X[0]-alpha*gradient_f(X)[0],X[1]-alpha*gradient_f(X)[1]])

def SectionDoree(X0,a,b,tolerance): #Recherche du pas optimal
tau = (1+np.sqrt(5))/2
it = 0
err = b-a
while np.abs(err)>tolerance :
aprime = a + (b-a)/(tau*tau)
bprime = a + (b-a)/tau
c , d = Omega(X0,aprime) , Omega(X0,bprime)
if c > d :
a = aprime
elif c < d :
b = bprime
else :
a , b = aprime , bprime
err = b – a
it = it + 1
return (a + b)/2

alpha = SectionDoree([0,1],0,1,0.01)

def PasOptimal (X0,alpha):

x=X0[0]
y=X0[1]
listeXY = [[x,y]]
k=0
kmax=10000

gradx,grady = gradient_f([x,y])

while (norme(gradient_f([x,y]))>0.01) and (k<=kmax):

X = [x,y]
x,y = Omega(X,alpha)
k+=1
gradx,grady = gradient_f([x,y])
listeXY.append([x,y])
X = [x,y]

a=X0[0]-alpha*gradient_f(X0)[0]
b= X0[1]-alpha*gradient_f(X0)[1]

Resp=Rosenbrock([a,b])

return (X,k,norme(gradient_f([x,y])),listeXY,Resp)

PO = PasOptimal([0,1],alpha)
print (« PasOptimal : [x,y] = « ,PO[0],’ et k = ‘,PO[1], ‘ et Epsilon = ‘,PO[2])

#Question 2c) Comparaison des methodes
#i. Affichage des courbes des iteres pour la methode du pas fixe et du pas optimal

def CreerlistepointsX(liste):
A = []
for i in range (len(liste)):
A.append(liste[i][0])

return (A)

def CreerlistepointsY(liste):
B= []
for i in range (len(liste)):
B.append(liste[i][1])
return (B)

A = CreerlistepointsX(PasFixe([0,1],0.01)[3]) # Création de la liste contenant les valeurs de X
B = CreerlistepointsY(PasFixe([0,1],0.01)[3]) # Création de la liste contenant les valeurs de Y

C = CreerlistepointsX(PasOptimal([0,1],alpha)[3]) # Création de la liste contenant les valeurs de X
D = CreerlistepointsY(PasOptimal([0,1],alpha)[3]) # Création de la liste contenant les valeurs de Y

plt.plot(C,D,’o’,color=’y’, label = « Gradient à pas optimal »)
plt.plot(A,B,’o’,color=’r’, label = « Gradient à pas fixe »)
plt.legend()
plt.title(« Convergence des solutions »)
plt.xlabel(« X »)
plt.ylabel(« Y »)

ax = Axes3D(plt.figure())
X = np.linspace(-10,10,50)
Y = np.linspace(-10,10,50)
X, Y = np.meshgrid(X, Y)
Z = Rosenbrock([X,Y])
ax.plot_surface(X, Y, Z)
plt.title(« Fonction de Rosenbrock »)
plt.xlabel(« X »)
plt.ylabel(« Y »)
ax.plot(C,D,’o’,color=’y’, label = « Gradient à pas optimal »)
ax.plot(A,B,’o’,color=’r’, label = « Gradient à pas fixe »)
plt.legend()
plt.show() # affichage de la nappe

### Comparaison des méthodes
### Pour atteindre une précision aussi fine que celle obtenue par la méthode du gradient à pas fixe,
### il faut beaucoup plus d’itérations avec la méthode du gradient à pas optimal

### Évolution des graphes de convergence quand p varie

### Quand le paramètre P augmente, le minimum du graphe de convergence diminue et inversement (minimum augmente si P diminue)

 

 

#Exercice 2

#Question 1

def emittance(L,T) : # On définie la définition de la fonction émittance
C = 2.997*10**8
h = 6.625 * 10**-34
k = 1.38*10**-23
n = 1
return (((2 * np.pi* h * C**2)/(n**2 * L**5))/(np.exp((h*C)/(n*k*T*L))-1))

ListeLambda = np.linspace(10**-7,2*10**-5,1000) #Permet de créer une liste de lambda de 10**-7 à 2*10**-5

def Liste_emittance(Liste,T): #Création des listes de valeur de la fonction emittance pour différents lambda
liste = np.zeros(len(Liste))
for i in range (len(liste)):
liste[i]=emittance(Liste[i],T)
return (liste)

def TracerCourbes(Tinit,Tfin): #Permet de tracer les courbes pour tous T
T = Tinit
while T !=Tfin:
plt.plot(ListeLambda,Liste_emittance(ListeLambda,T),label = str(T) +  » K » )
T+=50
plt.legend()
plt.xlabel(« Lambda (en m) »)
plt.ylabel(« Emittance »)
plt.title(« Evolution de M en fonction de Lambda pour différents T »)
plt.show()

TracerCourbes(300,800) #Lance le tracer des courbes pour différents T

##### Q2 ####

### Si on utilise la méthode de la section dorée, il faut connaitre un intervalle sur lequel l’emittance admet un unique maximum
### D’après l’affichage de l’emittance pour différents T, on prendra comme intervalle [10**-7,2*10**-5]

# La méthode de la section dorée cherche le minimum d’une fonction or dans notre cas, nous recherchons le maximum.
# Par conséquent, pour utiliser la section dorée nous avons changé le signe de l’inégalité.

def SectionDoree_Max(a,b,tolerance,T): #Section dorée pour la recherche d’un maximum
tau = (1+np.sqrt(5))/2
c = b-(b-a)/tau
d = a+(b-a)/tau
k=0
kmax = 10000

while np.abs(c-d)>tolerance and k<=kmax:
if emittance(c,T)>emittance(d,T):
b=d
else:
a=c

c = b-(b-a)/tau
d = a+(b-a)/tau
return ((b+a)/2)

def RechercheMax(T,tolerance): #Cette fonction cherche l’abscisse correspondant au maximum de la fonction M pour un T donné.
intervalle = [10**-7,2*10**-5]
min = SectionDoree_Max(intervalle[0],intervalle[1],tolerance,T)
return (min)

print (« Lambda* (T = 300 K) = « ,RechercheMax(300,10**-8))
print (« Lambda* (T = 400 K) = « ,RechercheMax(400,10**-8))
print (« Lambda* (T = 450 K) = « ,RechercheMax(450,10**-8))
print (« Lambda* (T = 500 K) = « ,RechercheMax(500,10**-8))
print (« Lambda* (T = 600 K) = « ,RechercheMax(600,10**-8))
print (« Lambda* (T = 800 K) = « ,RechercheMax(800,10**-8))