Coverage for src/ensae_teaching_cs/td_1a/optimisation_contrainte.py : 86%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
2"""
3@file
4@brief quelques fonctions sur l'optimisation quadratique avec contraintes
5"""
7import random
10def f_df_H(x=None, z=None):
11 """
12 Fonction demandée par la fonction
13 `solvers.cp <http://cvxopt.org/userguide/solvers.html#problems-with-nonlinear-objectives>`_.
15 Elle répond aux trois cas :
17 * ``F()`` : la fonction doit retourne le nombre de contraintes
18 non linéaires (:math:`f_k`) et un premier point :math:`X_0`,
19 * ``F(x)`` : la fonction retourne l'évaluation de :math:`f_0`
20 et de son gradient au point ``x``,
21 * ``F(x,z)`` : la fonction retourne l'évaluation de :math:`f_0`,
22 son gradient et de la dérivée seconde au point ``x``.
24 Voir @see fn exercice_particulier1.
25 Le problème d'optimisation est le suivant :
27 .. math::
29 \\left\\{ \\begin{array}{l} \\min_{x,y} \\left\\{ x^2 + y^2 - xy + y \\right\\}
30 \\\\ sous \\; contrainte \\; x + 2y = 1 \\end{array}\\right.
32 """
33 if x is None:
34 from cvxopt import matrix
35 # cas 1
36 x0 = matrix([[random.random(), random.random()]])
37 return 0, x0
38 else:
39 matrix = x.__class__
40 f = x[0] ** 2 + x[1] ** 2 - x[0] * x[1] + x[1]
41 d = matrix([x[0] * 2 - x[1], x[1] * 2 - x[0] + 1]).T
42 h = matrix([[2.0, -1.0], [-1.0, 2.0]])
43 if z is None:
44 # cas 2
45 return f, d
46 else:
47 # cas 3
48 return f, d, h
51def Arrow_Hurwicz(F, C, X0, p0, epsilon=0.1, rho=0.1,
52 do_print=True):
53 """
55 .. _code_Arrow_Hurwicz:
57 On implémente l'algorithme de `Arrow-Hurwicz <https://hal.archives-ouvertes.fr/hal-00490826/document>`_
58 d'une façon générique. Cela correspond au problème d'optimisation :
60 .. math::
62 \\left\\{ \\begin{array}{l} \\min_{X} f(X) \\\\ sous \\; contrainte \\; C(x) = 0 \\end{array}\\right.
64 @param F fonction qui retourne :math:`f(X)` et :math:`\\nabla f(X)`
65 @param C fonction qui retourne :math:`C(X)` et :math:`\\nabla C(X)`
66 @param X0 premier :math:`X`
67 @param p0 premier :math:`\\rho`
68 @param epsilon paramètre multiplicatif devant le gradient
69 @param rho paramètre multiplicatif durant la mise à jour prenant en compte la contrainte
70 @param do_print pour écrire des résultats intermédiaire en fonction des itérations
71 @return un dictionnaire ``{'x': solution, 'iteration': nombre d'itérations }``
72 """
73 diff = 1
74 itern = 0
75 while diff > 1e-10:
76 f, d = F(X0)
77 th, dt = C(X0)
78 Xt = [X0[i] - epsilon * (d[i] + dt[i] * p0) for i in range(len(X0))]
80 th, dt = C(Xt)
81 pt = p0 + rho * th
83 itern += 1
84 diff = sum([abs(Xt[i] - X0[i]) for i in range(len(X0))])
85 X0 = Xt
86 p0 = pt
87 if do_print and iter % 100 == 0:
88 print("i {0} diff {1:0.000}".format(itern, diff),
89 ":", f, X0, p0, th)
91 return {'x': X0, 'iteration': itern}
94def f_df(X):
95 """
96 F dans @see fn Arrow_Hurwicz
97 """
98 x, y = X
99 f = x ** 2 + y ** 2 - x * y + y
100 d = [x * 2 - y, y * 2 - x + 1]
101 return f, d
104def contrainte(X):
105 """
106 C dans @see fn Arrow_Hurwicz
107 """
108 x, y = X
109 f = x + 2 * y - 1
110 d = [1, 2]
111 return f, d
114def exercice_particulier1():
115 """
116 .. exref::
117 :title: solver.cp de cvxopt
118 :tag: Computer Science
120 On résoud le problème suivant avec `cvxopt <http://cvxopt.org/userguide/index.html>`_ :
122 .. math::
124 \\left\\{ \\begin{array}{l} \\min_{x,y} \\left \\{ x^2 + y^2 - xy + y \\right \\}
125 \\\\ sous \\; contrainte \\; x + 2y = 1 \\end{array}\\right.
127 Qui s'implémente à l'aide de la fonction suivante :
129 ::
131 def f_df_H(x=None,z=None) :
132 if x is None :
133 # cas 1
134 x0 = matrix ( [[ random.random(), random.random() ]])
135 return 0,x0
136 f = x[0]**2 + x[1]**2 - x[0]*x[1] + x[1]
137 d = matrix ( [ x[0]*2 - x[1], x[1]*2 - x[0] + 1 ] ).T
138 h = matrix ( [ [ 2.0, -1.0], [-1.0, 2.0] ])
139 if z is None:
140 # cas 2
141 return f, d
142 else :
143 # cas 3
144 return f, d, h
146 solvers.options['show_progress'] = False
147 A = matrix([ [ 1.0, 2.0 ] ]).trans()
148 b = matrix ( [[ 1.0] ] )
149 sol = solvers.cp ( f_df_H, A = A, b = b)
151 """
152 from cvxopt import solvers, matrix
153 t = solvers.options.get('show_progress', True)
154 solvers.options['show_progress'] = False
155 A = matrix([[1.0, 2.0]]).trans()
156 b = matrix([[1.0]])
157 sol = solvers.cp(f_df_H, A=A, b=b)
158 solvers.options['show_progress'] = t
159 return sol
162def exercice_particulier2():
163 """
164 .. exref::
165 :title: algorithme de Arrow-Hurwicz
166 :tag: Computer Science
168 On résoud le problème suivant avec l'algorithme de
169 `Arrow-Hurwicz <https://hal.archives-ouvertes.fr/hal-00490826/document>`_.
171 .. math::
173 \\left\\{ \\begin{array}{l} \\min_{x,y} \\left \\{ x^2 + y^2 - xy + y \\right \\} \\\\
174 sous \\; contrainte \\; x + 2y = 1 \\end{array}\\right.
176 Qui s'implémente à l'aide de la fonction suivante :
178 ::
180 import random
181 from ensae_teaching_cs.td_1a.optimisation_contrainte import Arrow_Hurwicz
183 def f_df(X) :
184 x,y = X
185 f = x**2 + y**2 - x*y + y
186 d = [ x*2 - y, y*2 - x + 1 ]
187 return f, d
189 def contrainte(X) :
190 x,y = X
191 f = x+2*y-1
192 d = [ 1,2]
193 return f, d
195 X0 = [ random.random(),random.random() ]
196 p0 = random.random()
197 sol = Arrow_Hurwicz(f_df, contrainte, X0, p0, do_print=False)
198 """
199 X0 = [random.random(), random.random()]
200 p0 = random.random()
201 sol = Arrow_Hurwicz(f_df, contrainte, X0, p0, do_print=False)
202 return sol
205if __name__ == "__main__":
206 sol1 = exercice_particulier1()
207 sol2 = exercice_particulier2()
208 print("cvxopt")
209 print(sol1)
210 print("solution:", sol1['x'].T)
211 print("Arrow_Hurwicz")
212 print(sol2)
213 print("solution:", sol2['x'])