from math import *
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (10, 8)
Minimisation de la fonctionnelle $f(x) = x^4-7x+8$ pour $x\in\mathbb{R}$.
def f(x):
return x**4 - 7*x + 8
def gradf(x):
return 4*x**3 -7
delta = 0.025
X = np.arange(0.8, 1.5, delta)
Y = f(X)
plt.plot(X,Y,'r')
ListRho = [0.125,0.1,0.01]
for rho in ListRho:
U = 1
N = 30
ListU = [U]
for j in range(N):
U = U - rho*gradf(U)
ListU += [U]
ListU = np.array(ListU)
plt.plot(ListU, f(ListU),'o-')
plt.legend(['f']+['rho='+str(rho) for rho in ListRho])
plt.show()
Minimisation de fonctionnelles quadratiques (ou résolution de $Ax=b$ pour $A$ symétrique définie positive).
def f1(x,y):
return 2*x**2 + 3*x + y**2 - 2
def f2(x,y):
return y**2 - 2*y + x**2 + 1
def gradf1(x,y):
return np.array([4*x + 3 , 2*y])
def gradf2(x,y):
return np.array([2*x , 2*y - 2])
# Methode du gradient à pas fixe, d'abord sur un exemple
rho = 0.1
U = np.array([1,1])
N = 30
print(U,'\n\t',f1(U[0],U[1]))
for j in range(N):
U = U - rho*gradf1(U[0],U[1])
print(U,'\n\t',f1(U[0],U[1]))
[1 1] 4 [0.3 0.8] -0.28000000000000025 [-0.12 0.64] -1.9216 [-0.372 0.512] -2.577088 [-0.5232 0.4096] -2.85435136 [-0.61392 0.32768] -2.9805902848 [-0.668352 0.262144] -3.042947731456 [-0.7010112 0.2097152] -3.07621972983808 [-0.72060672 0.16777216] -3.095124572510618 [-0.73236403 0.13421773] -3.106363546755924 [-0.73941842 0.10737418] -3.1132468452494777 [-0.74365105 0.08589935] -3.1175406840769124 [-0.74619063 0.06871948] -3.120248610931433 [-0.74771438 0.05497558] -3.1219672373201126 [-0.74862863 0.04398047] -3.12306195736151 [-0.74917718 0.03518437] -3.1237607058829564 [-0.74950631 0.0281475 ] -3.1242072309068645 [-0.74970378 0.022518 ] -3.124492764271431 [-0.74982227 0.0180144 ] -3.1246754182704897 [-0.74989336 0.01441152] -3.1247922853823518 [-0.74993602 0.01152922] -3.124867069012831 [-0.74996161 0.00922337] -3.124914926460737 [-0.74997697 0.0073787 ] -3.124945553760181 [-0.74998618 0.00590296] -3.1249651547036272 [-0.74999171 0.00472237] -3.1249776991172817 [-0.74999502 0.00377789] -3.124985727473566 [-0.74999701 0.00302231] -3.124990865596944 [-0.74999821 0.00241785] -3.1249941539870347 [-0.74999893 0.00193428] -3.124996258553499 [-0.74999936 0.00154743] -3.124997605474886 [-0.74999961 0.00123794] -3.12499846750416
# La même chose, avec stockage et un tracé graphique en sortie
rho = 0.1
U = np.array([1,1])
N = 30
ListU = [U]
for j in range(N):
U = U - rho*gradf1(U[0],U[1])
ListU += [U]
plt.plot([ListU[j][0] for j in range(N+1)],[ListU[j][1] for j in range(N+1)],'o-')
plt.title('minimisation de f1 / itérés dans le plan R^2')
plt.show()
# La même chose, avec stockage et un tracé graphique en sortie
rho = 0.1
U = np.array([3,-2])
N = 30
ListU = [U]
for j in range(N):
U = U - rho*gradf2(U[0],U[1])
ListU += [U]
plt.plot([ListU[j][0] for j in range(N+1)],[ListU[j][1] for j in range(N+1)],'o-')
plt.title('minimisation de f2 / itérés dans le plan R^2')
plt.show()
# Utilisation de contour / exemple extrait de la documentation
# https://matplotlib.org/stable/gallery/images_contours_and_fields/contour_demo.html
delta = 0.025
x = np.arange(-3.0, 3.0, delta)
y = np.arange(-2.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z1 = np.exp(-X**2 - Y**2)
Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2)
Z = (Z1 - Z2) * 2
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z)
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('Exemple d\'utilisation de contour')
plt.show()
delta = 0.025
x = np.arange(-1.0, 2.0, delta)
y = np.arange(-1.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z = f1(X,Y)
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z)
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('Min. de f1 / itérés dans R^2 et lignes de niveau')
rho = 0.1
U = np.array([1,2])
N = 30
ListU = [U]
for j in range(N):
U = U - rho*gradf1(U[0],U[1])
ListU += [U]
plt.plot([ListU[j][0] for j in range(N+1)],[ListU[j][1] for j in range(N+1)],'o-')
plt.axis('square')
plt.show()
delta = 0.025
x = np.arange(-1.0, 1.0, delta)
y = np.arange(-0.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z = f2(X,Y)
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z)
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('Min. de f2 / itérés dans R^2 et lignes de niveau')
rho = 0.1
U = np.array([1,0])
N = 30
ListU = [U]
for j in range(N):
U = U - rho*gradf2(U[0],U[1])
ListU += [U]
plt.plot([ListU[j][0] for j in range(N+1)],[ListU[j][1] for j in range(N+1)],'o-')
plt.axis('square')
plt.show()
delta = 0.025
x = np.arange(-2.5, 1.1, delta)
y = np.arange(-1.5, 2, delta)
X, Y = np.meshgrid(x, y)
Z = f1(X,Y)
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z)
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('Min. de f1 / différentes valeurs de rho')
ListRho = np.linspace(0.001,0.500,5)
for rho in ListRho:
U = np.array([1,2])
N = 30
ListU = [U]
for j in range(N):
U = U - rho*gradf1(U[0],U[1])
ListU += [U]
plt.plot([ListU[j][0] for j in range(N+1)],[ListU[j][1] for j in range(N+1)],'o-')
plt.legend(['rho='+str(rho) for rho in ListRho])
plt.axis('square')
plt.show()
rho = 0.1
U = np.array([1,1])
N = 30
ListU = [U]
for j in range(N):
U = U - rho*gradf1(U[0],U[1])
ListU += [U]
delta = 0.025
x = np.arange(-1.5, 1.0, delta)
y = np.arange(-1.0, 1.5, delta)
X, Y = np.meshgrid(x, y)
Z = f1(X,Y)
fig, ax = plt.subplots()
levels = [f1(U[0],U[1]) for U in ListU]
levels = np.sort(levels)
CS = ax.contour(X, Y, Z, levels)
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('Min. de f1 / lignes de niveau parcourues')
plt.plot([ListU[j][0] for j in range(N+1)],[ListU[j][1] for j in range(N+1)],'o-')
plt.axis('square')
plt.show()
rho = 0.1
U = np.array([1,1])
N = 30
ListU = [U]
for j in range(N):
U = U - rho*gradf2(U[0],U[1])
ListU += [U]
delta = 0.025
x = np.arange(-1.0, 1.0, delta)
y = np.arange(0.0, 2, delta)
X, Y = np.meshgrid(x, y)
Z = f2(X,Y)
fig, ax = plt.subplots()
levels = [f2(U[0],U[1]) for U in ListU]
levels = np.sort(levels)
CS = ax.contour(X, Y, Z, levels)
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('Min de f2 / lignes de niveau parcourues')
plt.plot([ListU[j][0] for j in range(N+1)],[ListU[j][1] for j in range(N+1)],'o-')
plt.axis('square')
plt.show()
U = np.array([0,1])
N = 20
n = 0
ListU = [U]
G = gradf1(U[0],U[1])
A = 2*np.array([[2,0],[0,1]])
print(A)
while n<N and np.linalg.norm(G)>1e-12 :
rho = np.linalg.norm(G)**2/np.vdot(np.dot(A,G),G)
U = U - rho*G
G = gradf1(U[0],U[1])
print(U,'\n n=',n, '\t/\t',f1(U[0],U[1]))
ListU += [U]
n += 1
N = n
delta = 0.025
x = np.arange(-1.5, 1.0, delta)
y = np.arange(-1.0, 1.5, delta)
X, Y = np.meshgrid(x, y)
Z = f1(X,Y)
fig, ax = plt.subplots()
levels = [f1(U[0],U[1]) for U in ListU[:8]]
levels = np.sort(levels)
CS = ax.contour(X, Y, Z, levels)
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('Min de f1 / méthode du gradient à pas optimal')
plt.plot([ListU[j][0] for j in range(N+1)],[ListU[j][1] for j in range(N+1)],'o-')
plt.axis('square')
plt.show()
[[4 0] [0 2]] [-0.88636364 0.40909091] n= 0 / -2.9204545454545454 [-0.67780749 0.09625668] n= 1 / -3.10531113271755 [-0.76312591 0.03937773] n= 2 / -3.1231048149139893 [-0.74305099 0.00926535] n= 3 / -3.1248175757671217 [-0.75126346 0.00379037] n= 4 / -3.124982440448172 [-0.74933111 0.00089185] n= 5 / -3.12499830977576 [-7.50121616e-01 3.64848466e-04] n= 6 / -3.1249998373046184 [-7.49935615e-01 8.58466978e-05] n= 7 / -3.124999984339482 [-7.50011706e-01 3.51191037e-05] n= 8 / -3.1249999984925703 [-7.49993803e-01 8.26331851e-06] n= 9 / -3.1249999998548996 [-7.50001127e-01 3.38044848e-06] n= 10 / -3.1249999999860334 [-7.49999403e-01 7.95399642e-07] n= 11 / -3.1249999999986557 [-7.50000108e-01 3.25390763e-07] n= 12 / -3.1249999999998708 [-7.49999943e-01 7.65625324e-08] n= 13 / -3.1249999999999876 [-7.5000001e-01 3.1321036e-08] n= 14 / -3.124999999999999 [-7.49999994e-01 7.36965555e-09] n= 15 / -3.125 [-7.50000001e-01 3.01485909e-09] n= 16 / -3.125 [-7.49999999e-01 7.09378594e-10] n= 17 / -3.125 [-7.50000000e-01 2.90200339e-10] n= 18 / -3.125 [-7.50000000e-01 6.82823974e-11] n= 19 / -3.125
U = np.array([1,1])
N = 30
n = 0
ListU = [U]
G = gradf2(U[0],U[1])
A = 2*np.array([[1,0],[0,1]])
print(A)
while n<N and np.linalg.norm(G)>1e-12 :
rho = np.linalg.norm(G)**2/np.vdot(np.dot(A,G),G)
U = U - rho*G
G = gradf2(U[0],U[1])
print(U,'\n n=',n, '\t/\t',f2(U[0],U[1]))
ListU += [U]
n += 1
N = n
delta = 0.025
x = np.arange(-1.5, 1.0, delta)
y = np.arange(-1.0, 1.5, delta)
X, Y = np.meshgrid(x, y)
Z = f2(X,Y)
fig, ax = plt.subplots()
levels = [f2(U[0],U[1]) for U in ListU[:]]
levels = np.sort(levels)
CS = ax.contour(X, Y, Z, levels)
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('Min. de f2 / méthode du gradient à pas optimal')
plt.plot([ListU[j][0] for j in range(N+1)],[ListU[j][1] for j in range(N+1)],'o-')
plt.axis('square')
plt.show()
[[2 0] [0 2]] [0. 1.] n= 0 / 0.0
def ProjK(X):
x,y = X[0],X[1]
px,py = x,y
z1 = x+y-1
z2 = x**2+y**2-1
if x>0 and y>0 and z2>0:
px = x/sqrt(x**2+y**2)
py = y/sqrt(x**2+y**2)
if abs(x-y)<=1 and x+y<1:
px = x + 0.5*(1-x-y)
py = y + 0.5*(1-x-y)
if x-y>1 and y<0 :
px,py = 1,0
if x-y<-1 and x<0:
px,py = 0,1
return np.array([px,py])
for i in range(500):
X = -2+4*np.random.rand(2)
Y = ProjK(X)
plt.plot([X[0],Y[0]],[X[1],Y[1]],'-.k',linewidth=0.5)
plt.plot([X[0]],[X[1]],'xb',markersize=3)
plt.plot([Y[0]],[Y[1]],'.r',markersize=2)
plt.axis('square')
plt.title('Fonction de projection sur le convexe fermé non-vide K')
plt.show()
def h(x,y):
return x**2-y
def gradh(x,y):
return np.array([2*x , -1])
rho = 0.05
U = np.array([1,1])
N = 30
ListU = [U]
U = ProjK(U)
ListPU = [U]
for j in range(N):
U = U - rho*gradh(U[0],U[1])
ListU += [U]
U = ProjK(U)
ListPU += [U]
fig, ax = plt.subplots()
x = np.arange(0,1,0.01)
plt.plot(np.cos(2*pi*x),np.sin(2*pi*x),'--k')
plt.plot(1-4*(x-0.5),4*(x-0.5),'--g')
for j in range(len(ListU)-1):
plt.plot([ListU[j][0],ListPU[j][0]],[ListU[j][1],ListPU[j][1]],'o-c')
plt.plot([ListPU[j][0],ListU[j+1][0]],[ListPU[j][1],ListU[j+1][1]],'o-b')
plt.legend(['g1','g2','uk','P(uk)'])
plt.title('Gradient avec projection pour minimisation sous contraintes')
delta = 0.025
x = np.arange(-1.5, 1.0, delta)
y = np.arange(-1.0, 1.5, delta)
X, Y = np.meshgrid(x, y)
Z = h(X,Y)
levels = [h(U[0],U[1]) for U in ListPU[:]]
levels = np.sort(levels)
CS = ax.contour(X, Y, Z, levels)
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('Min. de f2 / méthode du gradient à pas optimal')
plt.axis('square')
plt.xlim((-0.2,1.2))
plt.ylim((-0.2,1.2))
plt.show()
def h(x,y):
return x**2-y
def gradh(x,y):
return np.array([2*x , -1])
def ga2(x,y):
if x**2+y**2-1>=0:
return x**2+y**2-1
else:
return 0
def gradga2(x,y):
if x**2+y**2-1>=0:
return np.array([2*x , 2*y])
else:
return np.array([0 , 0])
def gb2(x,y):
if 1-x-y>=0:
return (1-x-y)**2
else:
return 0
def gradgb2(x,y):
if 1-x-y>=0:
return np.array([-2*x*(1-x-y) , -2*y*(1-x-y)])
else:
return np.array([0 , 0])
epsilon = 1
rho = 0.01
U = np.array([1,1])
N = 200
delta = 0.005
x = np.arange(-0.2, 1.2, delta)
y = np.arange(-0.2, 1.2, delta)
X, Y = np.meshgrid(x, y)
Z = h(X,Y)
# Pour voir les lignes de niveau de pénalisation :
#Z = h(X,Y) + 1/epsilon*((X**2+Y**2-1)*(X**2+Y**2-1>=0) + (1-X-Y)**2*(1-X-Y>=0))
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z, levels=np.arange(-1,1,0.1))
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('f2')
ListU = [U]
for j in range(N):
U = U - rho*(gradh(U[0],U[1]) + gradga2(U[0],U[1])/epsilon + gradgb2(U[0],U[1])/epsilon)
ListU += [U]
plt.plot([U[0] for U in ListU],[U[1] for U in ListU],'.-')
x = np.arange(0,1,0.01)
plt.plot(np.cos(2*pi*x),np.sin(2*pi*x),'--k')
plt.plot(1-4*(x-0.5),4*(x-0.5),'--g')
plt.legend(['uk','g1','g2'])
plt.title('Gradient avec pénalisation pour minimisation sous contraintes')
plt.axis('square')
plt.xlim((-0.2,1.2))
plt.ylim((-0.2,1.2))
plt.show()
epsilon = 1
rho = 0.01
U = np.array([1,1])
N = 200
delta = 0.005
x = np.arange(-0.2, 1.2, delta)
y = np.arange(-0.2, 1.2, delta)
X, Y = np.meshgrid(x, y)
#Z = h(X,Y)
# Pour voir les lignes de niveau de pénalisation :
Z = h(X,Y) + 1/epsilon*((X**2+Y**2-1)*(X**2+Y**2-1>=0) + (1-X-Y)**2*(1-X-Y>=0))
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z, levels=np.arange(-1,1,0.1))
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title('f2')
ListU = [U]
for j in range(N):
U = U - rho*(gradh(U[0],U[1]) + gradga2(U[0],U[1])/epsilon + gradgb2(U[0],U[1])/epsilon)
ListU += [U]
plt.plot([U[0] for U in ListU],[U[1] for U in ListU],'.-')
x = np.arange(0,1,0.01)
plt.plot(np.cos(2*pi*x),np.sin(2*pi*x),'--k')
plt.plot(1-4*(x-0.5),4*(x-0.5),'--g')
plt.legend(['uk','g1','g2'])
plt.title('Meme figure avec les lignes de niveau de pénalisation')
plt.axis('square')
plt.xlim((-0.2,1.2))
plt.ylim((-0.2,1.2))
plt.show()
def f(x,y):
return 0.5*( 2*(3-sqrt(x**2+y**2))**2
+ 1*(2-sqrt((x-3)**2+y**2))**2 )
def gradf(x,y):
L1 = sqrt(x**2+y**2)
L2 = sqrt((x-3)**2+y**2)
return np.array([-0.5*( 2*2*(3-L1)*2*x/L1 + 1*2*(2-L2)*2*(x-3)/L2 ) ,
-0.5*( 2*2*(3-L1)*2*y/L1 + 1*2*(2-L2)*2*y/L2 )])
rho = 0.1
U = np.array([1,0.2])
N = 200
ListU = [U]
for j in range(N):
U = U - rho*gradf(U[0],U[1])
ListU += [U]
plt.plot([U[0] for U in ListU],[U[1] for U in ListU],'o-')
x = np.arange(0,0.5,0.001)
plt.plot(1.5+1.5*np.cos(2*pi*x),2*np.sin(2*pi*x),'-k',linewidth=2)
plt.title('Methode de gradient à pas fixe, en l\'absence de la contrainte')
plt.axis('square')
plt.show()
delta = 0.005
x = np.arange(0., 3., delta)
y = np.arange(0., 3., delta)
X, Y = np.meshgrid(x, y)
Z = 0.5*( 2*(3-np.sqrt(X**2+Y**2))**2
+ 1*(2-np.sqrt((3-X)**2+Y**2))**2 )
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z, levels=np.arange(0,10,0.1))
ax.clabel(CS, inline=True, fontsize=10)
rho = 0.1
U = np.array([1,0.2])
N = 200
ListU = [U]
for j in range(N):
U = U - rho*gradf(U[0],U[1])
ListU += [U]
plt.plot([U[0] for U in ListU],[U[1] for U in ListU],'o-')
x = np.arange(0,0.5,0.001)
plt.plot(1.5+1.5*np.cos(2*pi*x),2*np.sin(2*pi*x),'-k',linewidth=2)
plt.title('Methode de gradient à pas fixe, en l\'absence de la contrainte')
plt.axis('square')
plt.show()
def g2(x,y):
if (x-1.5)**2/2.25+y**2/4>=1:
return (x-1.5)**2/2.25 + y**2/4
else:
return 0
def gradg2(x,y):
if (x-1.5)**2/2.25+y**2/4>=1:
return np.array([(x-1.5)/1.125 , y/2])
else:
return np.array([0 , 0])
epsilon = 1
rho = 0.05
U = np.array([1,0.2])
N = 50
ListU = [U]
for j in range(N):
U = U - rho*(gradf(U[0],U[1]) + gradg2(U[0],U[1])/epsilon)
ListU += [U]
plt.plot([U[0] for U in ListU],[U[1] for U in ListU],'o-')
x = np.arange(0,0.5,0.001)
plt.plot(1.5+1.5*np.cos(2*pi*x),2*np.sin(2*pi*x),'-k',linewidth=2)
plt.title('Methode de minimisation sous contrainte par pénalisation')
plt.axis('square')
plt.show()
delta = 0.005
x = np.arange(0., 3., delta)
y = np.arange(0., 3., delta)
X, Y = np.meshgrid(x, y)
Z = 0.5*( 2*(3-np.sqrt(X**2+Y**2))**2
+ 1*(2-np.sqrt((3-X)**2+Y**2))**2 )
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z, levels=np.arange(0,10,0.2))
ax.clabel(CS, inline=True, fontsize=10)
epsilon = 1
rho = 0.05
U = np.array([1,0.2])
N = 50
ListU = [U]
for j in range(N):
U = U - rho*(gradf(U[0],U[1]) + gradg2(U[0],U[1])/epsilon)
ListU += [U]
plt.plot([U[0] for U in ListU],[U[1] for U in ListU],'o-')
x = np.arange(0,0.5,0.001)
plt.plot(1.5+1.5*np.cos(2*pi*x),2*np.sin(2*pi*x),'-k',linewidth=2)
CS = ax.contour(X, Y, Z, levels=[f(U[0],U[1])])
plt.title('Methode de minimisation sous contrainte par pénalisation')
plt.axis('square')
plt.show()
delta = 0.005
x = np.arange(0., 3., delta)
y = np.arange(0., 3., delta)
X, Y = np.meshgrid(x, y)
Z = 0.5*( 2*(3-np.sqrt(X**2+Y**2))**2
+ 1*(2-np.sqrt((3-X)**2+Y**2))**2 )
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z, levels=np.arange(0,0.2,0.02))
ax.clabel(CS, inline=True, fontsize=10)
epsilon = 1
rho = 0.05
U = np.array([1,0.2])
N = 50
ListU = [U]
for j in range(N):
U = U - rho*(gradf(U[0],U[1]) + gradg2(U[0],U[1])/epsilon)
ListU += [U]
plt.plot([U[0] for U in ListU],[U[1] for U in ListU],'o-')
x = np.arange(0,0.5,0.001)
plt.plot(1.5+1.5*np.cos(2*pi*x),2*np.sin(2*pi*x),'-k',linewidth=2)
plt.title('Methode de minimisation sous contrainte par pénalisation / ZOOM')
plt.axis('square')
plt.xlim((2.0,2.5))
plt.ylim((1.5,2.0))
plt.show()
delta = 0.005
x = np.arange(0., 3., delta)
y = np.arange(0., 3., delta)
X, Y = np.meshgrid(x, y)
Z = 0.5*( 2*(3-np.sqrt(X**2+Y**2))**2
+ 1*(2-np.sqrt((3-X)**2+Y**2))**2 ) + ((X-1.5)**2/2.25 + Y**2/4)*((X-1.5)**2/2.25+Y**2/4>=1)
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z, levels=np.arange(0,10,0.2))
ax.clabel(CS, inline=True, fontsize=10)
epsilon = 1
rho = 0.05
U = np.array([1,0.2])
N = 50
ListU = [U]
for j in range(N):
U = U - rho*(gradf(U[0],U[1]) + gradg2(U[0],U[1])/epsilon)
ListU += [U]
plt.plot([U[0] for U in ListU],[U[1] for U in ListU],'o-')
x = np.arange(0,0.5,0.001)
plt.plot(1.5+1.5*np.cos(2*pi*x),2*np.sin(2*pi*x),'-k',linewidth=2)
plt.title('Methode de minimisation sous contrainte par pénalisation')
plt.axis('square')
plt.show()
