Source Code for Module pypower.t.t_pips

  1  # Copyright (C) 2010-2011 Power System Engineering Research Center (PSERC) 
  2  # Copyright (C) 2011 Richard Lincoln 
  3  # 
  4  # PYPOWER is free software: you can redistribute it and/or modify 
  5  # it under the terms of the GNU General Public License as published 
  6  # by the Free Software Foundation, either version 3 of the License, 
  7  # or (at your option) any later version. 
  8  # 
  9  # PYPOWER is distributed in the hope that it will be useful, 
 10  # but WITHOUT ANY WARRANTY], without even the implied warranty of 
 11  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 12  # GNU General Public License for more details. 
 13  # 
 14  # You should have received a copy of the GNU General Public License 
 15  # along with PYPOWER. If not, see <http://www.gnu.org/licenses/>. 
 16   
 17  """Tests of PIPS NLP solver. 
 18  """ 
 19   
 20  from math import sqrt 
 21   
 22  from numpy import zeros, ones, array, eye, prod, dot, asscalar, Inf 
 23   
 24  from scipy.sparse import csr_matrix as sparse 
 25  from scipy.sparse import eye as speye 
 26   
 27  from pypower.pips import pips 
 28   
 29  from pypower.t.t_begin import t_begin 
 30  from pypower.t.t_is import t_is 
 31  from pypower.t.t_end import t_end 
 32  from pypower.t.t_ok import t_ok 
 33   
 34   
 35  ## unconstrained banana function 
 36  ## from MATLAB Optimization Toolbox's bandem.m 
 37 -def f2(x, return_hessian=False): 
 38      a = 100 
 39      f = a * (x[1] - x[0]**2)**2 + (1 - x[0])**2 
 40      df = array([ 
 41          4 * a * (x[0]**3 - x[0] * x[1]) + 2 * x[0] - 2, 
 42          2 * a * (x[1] - x[0]**2) 
 43      ]) 
 44   
 45      if not return_hessian: 
 46          return f, df 
 47   
 48      d2f = 4 * a * sparse([ 
 49          [3 * x[0]**2 - x[1] + 1. / (2 * a), -x[0]], 
 50          [                           -x[0],   0.5] 
 51      ]) 
 52      return f, df, d2f 
 53   
 54   
 55  ## unconstrained 3-d quadratic 
 56  ## from http://www.akiti.ca/QuadProgEx0Constr.html 
 57 -def f3(x, return_hessian=False): 
 58      H = sparse([ 
 59          [ 5, -2, -1], 
 60          [-2,  4,  3], 
 61          [-1,  3,  5] 
 62      ], dtype=float) 
 63      c = array([2, -35, -47], float) 
 64      f = 0.5 * dot(x * H, x) + dot(c, x) + 5 
 65      df = H * x + c 
 66      if not return_hessian: 
 67          return f, df 
 68      d2f = H 
 69      return f, df, d2f 
 70   
 71   
 72  ## constrained 4-d QP 
 73  ## from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm 
 74 -def f4(x, return_hessian=False): 
 75      H = sparse([ 
 76          [1003.1, 4.3, 6.3,  5.9], 
 77          [   4.3, 2.2, 2.1,  3.9], 
 78          [   6.3, 2.1, 3.5,  4.8], 
 79          [   5.9, 3.9, 4.8, 10.0] 
 80      ]) 
 81      c = zeros(4) 
 82      f = 0.5 * dot(x * H, x) + dot(c, x) 
 83      df = H * x + c 
 84      if not return_hessian: 
 85          return f, df 
 86      d2f = H 
 87      return f, df, d2f 
 88   
 89   
 90  ## constrained 2-d nonlinear 
 91  ## from http://en.wikipedia.org/wiki/Nonlinear_programming#2-dimensional_example 
 92 -def f5(x, return_hessian=False): 
 93      c = -array([1.0, 1.0]) 
 94      f = dot(c, x) 
 95      df = c 
 96      if not return_hessian: 
 97          return f, df 
 98      d2f = sparse((2, 2)) 
 99      return f, df, d2f 
100   
101 -def gh5(x): 
102      h = dot([[-1.0, -1.0], 
103               [ 1.0,  1.0]], x**2) + [1, -2] 
104      dh = 2 * sparse([[-x[0], x[0]], 
105                       [-x[1], x[1]]]) 
106      g = array([]) 
107      dg = None 
108      return h, g, dh, dg 
109   
110 -def hess5(x, lam, cost_mult): 
111      mu = lam['ineqnonlin'] 
112      Lxx = 2 * dot([-1.0, 1.0], mu) * eye(2) 
113      return Lxx 
114   
115   
116  ## constrained 3-d nonlinear 
117  ## from http://en.wikipedia.org/wiki/Nonlinear_programming#3-dimensional_example 
118 -def f6(x, return_hessian=False): 
119      f = -x[0] * x[1] - x[1] * x[2] 
120      df = -array([x[1], x[0] + x[2], x[1]]) 
121      if not return_hessian: 
122          return f, df 
123      d2f = -sparse([[0, 1, 0], 
124                     [1, 0, 1], 
125                     [0, 1, 0]], dtype=float) 
126      return f, df, d2f 
127   
128 -def gh6(x): 
129      h = dot([[1.0, -1.0, 1.0], 
130               [1.0,  1.0, 1.0]], x**2) + [-2.0, -10.0] 
131      dh = 2 * sparse([[ x[0], x[0]], 
132                       [-x[1], x[1]], 
133                       [ x[2], x[2]]], dtype=float) 
134      g = array([]) 
135      dg = None 
136      return h, g, dh, dg 
137   
138 -def hess6(x, lam, cost_mult=1): 
139      mu = lam['ineqnonlin'] 
140      Lxx = cost_mult * \ 
141          sparse([[ 0, -1,  0], 
142                  [-1,  0, -1], 
143                  [ 0, -1,  0]], dtype=float) + \ 
144          sparse([[2 * dot([1, 1], mu),  0, 0], 
145                  [0, 2 * dot([-1, 1], mu), 0], 
146                  [0, 0,  2 * dot([1, 1], mu)]], dtype=float) 
147      return Lxx 
148   
149   
150  ## constrained 4-d nonlinear 
151  ## Hock & Schittkowski test problem #71 
152 -def f7(x, return_hessian=False): 
153      f = x[0] * x[3] * sum(x[:3]) + x[2] 
154      df = array([ x[0] * x[3] + x[3] * sum(x[:3]), 
155                   x[0] * x[3], 
156                   x[0] * x[3] + 1, 
157                   x[0] * sum(x[:3]) ]) 
158      if not return_hessian: 
159          return f, df 
160      d2f = sparse([ 
161          [2 * x[3],               x[3], x[3], 2 * x[0] + x[1] + x[2]], 
162          [                  x[3],  0.0,  0.0,                   x[0]], 
163          [                  x[3],  0.0,  0.0,                   x[0]], 
164          [2 * x[0] + x[1] + x[2], x[0], x[0],                    0.0] 
165      ]) 
166      return f, df, d2f 
167   
168 -def gh7(x): 
169      g = array( [sum(x**2) - 40.0] ) 
170      h = array( [ -prod(x) + 25.0] ) 
171      dg = sparse( 2 * x ).T 
172      dh = sparse( -prod(x) / x ).T 
173      return h, g, dh, dg 
174   
175 -def hess7(x, lam, sigma=1): 
176      lmbda = asscalar( lam['eqnonlin'] ) 
177      mu    = asscalar( lam['ineqnonlin'] ) 
178      _, _, d2f = f7(x, True) 
179   
180      Lxx = sigma * d2f + lmbda * 2 * speye(4, 4) - \ 
181          mu * sparse([ 
182              [        0.0, x[2] * x[3], x[2] * x[3], x[1] * x[2]], 
183              [x[2] * x[2],         0.0, x[0] * x[3], x[0] * x[2]], 
184              [x[1] * x[3], x[0] * x[3],         0.0, x[0] * x[1]], 
185              [x[1] * x[2], x[0] * x[2], x[0] * x[1],         0.0] 
186          ]) 
187      return Lxx 
188   
189   
190 -def t_pips(quiet=False): 
191      """Tests of pips NLP solver. 
192   
193      @author: Ray Zimmerman (PSERC Cornell) 
194      @author: Richard Lincoln 
195      """ 
196      t_begin(60, quiet) 
197   
198      t = 'unconstrained banana function : ' 
199      ## from MATLAB Optimization Toolbox's bandem.m 
200      f_fcn = f2 
201      x0 = array([-1.9, 2]) 
202      # solution = pips(f_fcn, x0, opt={'verbose': 2}) 
203      solution = pips(f_fcn, x0) 
204      x, f, s, lam, out = solution["x"], solution["f"], solution["eflag"], \ 
205              solution["lmbda"], solution["output"] 
206      t_is(s, 1, 13, [t, 'success']) 
207      t_is(x, [1, 1], 13, [t, 'x']) 
208      t_is(f, 0, 13, [t, 'f']) 
209      t_is(out['hist'][-1]['compcond'], 0, 6, [t, 'compcond']) 
210      t_ok(len(lam['mu_l']) == 0, [t, 'lam.mu_l']) 
211      t_ok(len(lam['mu_u']) == 0, [t, 'lam.mu_u']) 
212      t_is(lam['lower'], zeros(x.shape), 13, [t, 'lam[\'lower\']']) 
213      t_is(lam['upper'], zeros(x.shape), 13, [t, 'lam[\'upper\']']) 
214   
215      t = 'unconstrained 3-d quadratic : ' 
216      ## from http://www.akiti.ca/QuadProgEx0Constr.html 
217      f_fcn = f3 
218      x0 = array([0, 0, 0], float) 
219      # solution = pips(f_fcn, x0, opt={'verbose': 2}) 
220      solution = pips(f_fcn, x0) 
221      x, f, s, lam, out = solution["x"], solution["f"], solution["eflag"], \ 
222              solution["lmbda"], solution["output"] 
223      t_is(s, 1, 13, [t, 'success']) 
224      t_is(x, [3, 5, 7], 13, [t, 'x']) 
225      t_is(f, -244, 13, [t, 'f']) 
226      t_is(out['hist'][-1]['compcond'], 0, 6, [t, 'compcond']) 
227      t_ok(len(lam['mu_l']) == 0, [t, 'lam.mu_l']) 
228      t_ok(len(lam['mu_u']) == 0, [t, 'lam.mu_u']) 
229      t_is(lam['lower'], zeros(x.shape), 13, [t, 'lam[\'lower\']']) 
230      t_is(lam['upper'], zeros(x.shape), 13, [t, 'lam[\'upper\']']) 
231   
232      t = 'constrained 4-d QP : ' 
233      ## from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm 
234      f_fcn = f4 
235      x0 = array([1.0, 0.0, 0.0, 1.0]) 
236      A = array([ 
237          [1.0,  1.0,  1.0,  1.0 ], 
238          [0.17, 0.11, 0.10, 0.18] 
239      ]) 
240      l = array([1,  0.10]) 
241      u = array([1.0, Inf]) 
242      xmin = zeros(4) 
243      # solution = pips(f_fcn, x0, A, l, u, xmin, opt={'verbose': 2}) 
244      solution = pips(f_fcn, x0, A, l, u, xmin) 
245      x, f, s, lam, out = solution["x"], solution["f"], solution["eflag"], \ 
246              solution["lmbda"], solution["output"] 
247      t_is(s, 1, 13, [t, 'success']) 
248      t_is(x, array([0, 2.8, 0.2, 0]) / 3, 6, [t, 'x']) 
249      t_is(f, 3.29 / 3, 6, [t, 'f']) 
250      t_is(out['hist'][-1]['compcond'], 0, 6, [t, 'compcond']) 
251      t_is(lam['mu_l'], array([6.58, 0]) / 3, 6, [t, 'lam.mu_l']) 
252      t_is(lam['mu_u'], array([0, 0]), 13, [t, 'lam.mu_u']) 
253      t_is(lam['lower'], array([2.24, 0, 0, 1.7667]), 4, [t, 'lam[\'lower\']']) 
254      t_is(lam['upper'], zeros(x.shape), 13, [t, 'lam[\'upper\']']) 
255   
256      # H = array([ 
257      #     [1003.1, 4.3, 6.3,  5.9], 
258      #     [   4.3, 2.2, 2.1,  3.9], 
259      #     [   6.3, 2.1, 3.5,  4.8], 
260      #     [   5.9, 3.9, 4.8, 10.0] 
261      # ]) 
262      # c = zeros(4) 
263      # ## check with quadprog (for dev testing only) 
264      # x, f, s, out, lam = quadprog(H,c,-A(2,:), -0.10, A(1,:), 1, xmin) 
265      # t_is(s, 1, 13, [t, 'success']) 
266      # t_is(x, [0 2.8 0.2 0]/3, 6, [t, 'x']) 
267      # t_is(f, 3.29/3, 6, [t, 'f']) 
268      # t_is(lam['eqlin'], -6.58/3, 6, [t, 'lam.eqlin']) 
269      # t_is(lam.['ineqlin'], 0, 13, [t, 'lam.ineqlin']) 
270      # t_is(lam['lower'], [2.24001.7667], 4, [t, 'lam[\'lower\']']) 
271      # t_is(lam['upper'], [0000], 13, [t, 'lam[\'upper\']']) 
272   
273      t = 'constrained 2-d nonlinear : ' 
274      ## from http://en.wikipedia.org/wiki/Nonlinear_programming#2-dimensional_example 
275      f_fcn = f5 
276      gh_fcn = gh5 
277      hess_fcn = hess5 
278      x0 = array([1.1, 0.0]) 
279      xmin = zeros(2) 
280      # xmax = 3 * ones(2, 1) 
281      # solution = pips(f_fcn, x0, xmin=xmin, gh_fcn=gh_fcn, hess_fcn=hess_fcn, opt={'verbose': 2}) 
282      solution = pips(f_fcn, x0, xmin=xmin, gh_fcn=gh_fcn, hess_fcn=hess_fcn) 
283      x, f, s, lam, out = solution["x"], solution["f"], solution["eflag"], \ 
284              solution["lmbda"], solution["output"] 
285      t_is(s, 1, 13, [t, 'success']) 
286      t_is(x, [1, 1], 6, [t, 'x']) 
287      t_is(f, -2, 6, [t, 'f']) 
288      t_is(out['hist'][-1]['compcond'], 0, 6, [t, 'compcond']) 
289      t_is(lam['ineqnonlin'], array([0, 0.5]), 6, [t, 'lam.ineqnonlin']) 
290      t_ok(len(lam['mu_l']) == 0, [t, 'lam.mu_l']) 
291      t_ok(len(lam['mu_u']) == 0, [t, 'lam.mu_u']) 
292      t_is(lam['lower'], zeros(x.shape), 13, [t, 'lam[\'lower\']']) 
293      t_is(lam['upper'], zeros(x.shape), 13, [t, 'lam[\'upper\']']) 
294      # ## check with fmincon (for dev testing only) 
295      # # fmoptions = optimset('Algorithm', 'interior-point') 
296      # # [x, f, s, out, lam] = fmincon(f_fcn, x0, [], [], [], [], xmin, [], gh_fcn, fmoptions) 
297      # [x, f, s, out, lam] = fmincon(f_fcn, x0, [], [], [], [], [], [], gh_fcn) 
298      # t_is(s, 1, 13, [t, 'success']) 
299      # t_is(x, [1 1], 4, [t, 'x']) 
300      # t_is(f, -2, 6, [t, 'f']) 
301      # t_is(lam.ineqnonlin, [00.5], 6, [t, 'lam.ineqnonlin']) 
302   
303      t = 'constrained 3-d nonlinear : ' 
304      ## from http://en.wikipedia.org/wiki/Nonlinear_programming#3-dimensional_example 
305      f_fcn = f6 
306      gh_fcn = gh6 
307      hess_fcn = hess6 
308      x0 = array([1.0, 1.0, 0.0]) 
309      # solution = pips(f_fcn, x0, gh_fcn=gh_fcn, hess_fcn=hess_fcn, opt={'verbose': 2, 'comptol': 1e-9}) 
310      solution = pips(f_fcn, x0, gh_fcn=gh_fcn, hess_fcn=hess_fcn) 
311      x, f, s, lam, out = solution["x"], solution["f"], solution["eflag"], \ 
312              solution["lmbda"], solution["output"] 
313      t_is(s, 1, 13, [t, 'success']) 
314      t_is(x, [1.58113883, 2.23606798, 1.58113883], 6, [t, 'x']) 
315      t_is(f, -5 * sqrt(2), 6, [t, 'f']) 
316      t_is(out['hist'][-1]['compcond'], 0, 6, [t, 'compcond']) 
317      t_is(lam['ineqnonlin'], array([0, sqrt(2) / 2]), 7, [t, 'lam.ineqnonlin']) 
318      t_ok(len(lam['mu_l']) == 0, [t, 'lam.mu_l']) 
319      t_ok(len(lam['mu_u']) == 0, [t, 'lam.mu_u']) 
320      t_is(lam['lower'], zeros(x.shape), 13, [t, 'lam[\'lower\']']) 
321      t_is(lam['upper'], zeros(x.shape), 13, [t, 'lam[\'upper\']']) 
322      # ## check with fmincon (for dev testing only) 
323      # # fmoptions = optimset('Algorithm', 'interior-point') 
324      # # [x, f, s, out, lam] = fmincon(f_fcn, x0, [], [], [], [], xmin, [], gh_fcn, fmoptions) 
325      # [x, f, s, out, lam] = fmincon(f_fcn, x0, [], [], [], [], [], [], gh_fcn) 
326      # t_is(s, 1, 13, [t, 'success']) 
327      # t_is(x, [1.58113883 2.23606798 1.58113883], 4, [t, 'x']) 
328      # t_is(f, -5*sqrt(2), 8, [t, 'f']) 
329      # t_is(lam.ineqnonlin, [0sqrt(2)/2], 8, [t, 'lam.ineqnonlin']) 
330   
331      t = 'constrained 3-d nonlinear (dict) : ' 
332      p = {'f_fcn': f_fcn, 'x0': x0, 'gh_fcn': gh_fcn, 'hess_fcn': hess_fcn} 
333      solution = pips(p) 
334      x, f, s, lam, out = solution["x"], solution["f"], solution["eflag"], \ 
335              solution["lmbda"], solution["output"] 
336      t_is(s, 1, 13, [t, 'success']) 
337      t_is(x, [1.58113883, 2.23606798, 1.58113883], 6, [t, 'x']) 
338      t_is(f, -5 * sqrt(2), 6, [t, 'f']) 
339      t_is(out['hist'][-1]['compcond'], 0, 6, [t, 'compcond']) 
340      t_is(lam['ineqnonlin'], [0, sqrt(2) / 2], 7, [t, 'lam.ineqnonlin']) 
341      t_ok(len(lam['mu_l']) == 0, [t, 'lam.mu_l']) 
342      t_ok(len(lam['mu_u']) == 0, [t, 'lam.mu_u']) 
343      t_is(lam['lower'], zeros(x.shape), 13, [t, 'lam[\'lower\']']) 
344      t_is(lam['upper'], zeros(x.shape), 13, [t, 'lam[\'upper\']']) 
345   
346      t = 'constrained 4-d nonlinear : ' 
347      ## Hock & Schittkowski test problem #71 
348      f_fcn = f7 
349      gh_fcn = gh7 
350      hess_fcn = hess7 
351      x0 = array([1.0, 5.0, 5.0, 1.0]) 
352      xmin = ones(4) 
353      xmax = 5 * xmin 
354      # solution = pips(f_fcn, x0, xmin=xmin, xmax=xmax, gh_fcn=gh_fcn, hess_fcn=hess_fcn, opt={'verbose': 2, 'comptol': 1e-9}) 
355      solution = pips(f_fcn, x0, xmin=xmin, xmax=xmax, gh_fcn=gh_fcn, hess_fcn=hess_fcn) 
356      x, f, s, lam, _ = solution["x"], solution["f"], solution["eflag"], \ 
357              solution["lmbda"], solution["output"] 
358      t_is(s, 1, 13, [t, 'success']) 
359      t_is(x, [1, 4.7429994, 3.8211503, 1.3794082], 6, [t, 'x']) 
360      t_is(f, 17.0140173, 6, [t, 'f']) 
361      t_is(lam['eqnonlin'], 0.1614686, 5, [t, 'lam.eqnonlin']) 
362      t_is(lam['ineqnonlin'], 0.55229366, 5, [t, 'lam.ineqnonlin']) 
363      t_ok(len(lam['mu_l']) == 0, [t, 'lam.mu_l']) 
364      t_ok(len(lam['mu_u']) == 0, [t, 'lam.mu_u']) 
365      t_is(lam['lower'], [1.08787121024, 0, 0, 0], 5, [t, 'lam[\'lower\']']) 
366      t_is(lam['upper'], zeros(x.shape), 7, [t, 'lam[\'upper\']']) 
367   
368      t_end() 
369   
370   
371      # ##-----  eg99 : linearly constrained fmincon example, pips can't solve  ----- 
372      # function [f, df, d2f] = eg99(x) 
373      # f = -x(1)*x(2)*x(3) 
374      # df = -[ x(2)*x(3) 
375      #         x(1)*x(3) 
376      #         x(1)*x(2)   ] 
377      # d2f = -[    0       x(3)    x(2) 
378      #             x(3)    0       x(1) 
379      #             x(2)    x(1)    0   ] 
380      # end 
381      # 
382      # x0 = [101010] 
383      # A = [1 2 2] 
384      # l = 0 
385      # u = 72 
386      # fmoptions = optimset('Display', 'testing') 
387      # fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point') 
388      # [x, f, s, out, lam] = fmincon(f_fcn, x0, [-A A], [-l u], [], [], [], [], [], fmoptions) 
389      # t_is(x, [24 12 12], 13, t) 
390      # t_is(f, -3456, 13, t) 
391   
392   
393  if __name__ == '__main__': 
394      t_pips(quiet=False) 
395