Coverage for src/ensae_teaching_cs/td_1a/optimisation_contrainte.py: 86%

63 statements  

« prev     ^ index     » next       coverage.py v7.1.0, created at 2023-01-27 05:44 +0100

1# -*- coding: utf-8 -*- 

2""" 

3@file 

4@brief quelques fonctions sur l'optimisation quadratique avec contraintes 

5""" 

6 

7import random 

8 

9 

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>`_. 

14 

15 Elle répond aux trois cas : 

16 

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``. 

23 

24 Voir @see fn exercice_particulier1. 

25 Le problème d'optimisation est le suivant : 

26 

27 .. math:: 

28 

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. 

31 

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 

49 

50 

51def Arrow_Hurwicz(F, C, X0, p0, epsilon=0.1, rho=0.1, 

52 do_print=True): 

53 """ 

54 

55 .. _code_Arrow_Hurwicz: 

56 

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 : 

59 

60 .. math:: 

61 

62 \\left\\{ \\begin{array}{l} \\min_{X} f(X) \\\\ sous \\; contrainte \\; C(x) = 0 \\end{array}\\right. 

63 

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))] 

79 

80 th, dt = C(Xt) 

81 pt = p0 + rho * th 

82 

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(f"i {itern} diff {diff:0.000}", 

89 ":", f, X0, p0, th) 

90 

91 return {'x': X0, 'iteration': itern} 

92 

93 

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 

102 

103 

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 

112 

113 

114def exercice_particulier1(): 

115 """ 

116 .. exref:: 

117 :title: solver.cp de cvxopt 

118 :tag: Computer Science 

119 

120 On résoud le problème suivant avec `cvxopt <http://cvxopt.org/userguide/index.html>`_ : 

121 

122 .. math:: 

123 

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. 

126 

127 Qui s'implémente à l'aide de la fonction suivante : 

128 

129 :: 

130 

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 

145 

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) 

150 

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 

160 

161 

162def exercice_particulier2(): 

163 """ 

164 .. exref:: 

165 :title: algorithme de Arrow-Hurwicz 

166 :tag: Computer Science 

167 

168 On résoud le problème suivant avec l'algorithme de 

169 `Arrow-Hurwicz <https://hal.archives-ouvertes.fr/hal-00490826/document>`_. 

170 

171 .. math:: 

172 

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. 

175 

176 Qui s'implémente à l'aide de la fonction suivante : 

177 

178 :: 

179 

180 import random 

181 from ensae_teaching_cs.td_1a.optimisation_contrainte import Arrow_Hurwicz 

182 

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 

188 

189 def contrainte(X) : 

190 x,y = X 

191 f = x+2*y-1 

192 d = [ 1,2] 

193 return f, d 

194 

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 

203 

204 

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'])