Optimization: principles and algorithms, by Michel Bierlaire
globalSqp.m
Go to the documentation of this file.
1 %> \file
2 %> Algorithm 20.2: global SQP algorithm. Implementation of algorithm 20.2 of \cite Bier15-book
3 %>
4 %> @note Tested with \ref run2007.m
5 %>
6 %> @author Michel Bierlaire
7 %> @date Fri Mar 27 15:54:09 2015
8 %> @ingroup Algorithms
9 %> @ingroup chap20
10
11 %> Applies the global SQP method to solve \f[\min_x f(x) \f] subject to \f[h(x)=0,\f] where \f$f:\mathbb{R}^n \to \mathbb{R}\f$ and \f$h:\mathbb{R}^n \to \mathbb{R}^m \f$.
12 %> @param problem the name of the Octave function defining f(x), h(x) and their derivatives. The funtion has two arguments: x and index. If index=0, the objective function \f$f\f$ and its derivatives are evaluated. If index=\f$i\f$, the constraint \f$h_i\f$ and its derivtives are evaluated.
13 %> @param x0 starting primal point (nx1)
14 %> @param lambda0 starting dual point (mx1)
15 %> @param eps algorithm stops if \f$\|\nabla L(x_k,\lambda_k\| \leq \varepsilon \f$ and \f$\|h(x_k)\|^2\f$.
16 %> @param maxiter maximum number of iterations (default: 100)
17 %> @return [solution,lambda]
18 %> @return x: primal solution
19 %> @return lambda: dual solution
20 function [solution,lambda] = globalSqp(problem,x0,lambda0,eps,maxiter)
22  n = size(x0,1) ;
23  m = size(lambda0,1) ;
24  xk = x0 ;
25  beta1 = 0.3 ;
26  cconst = 0.1 ;
27  lambdak = lambda0 ;
28  c = norm(lambda0,Inf) + cconst ;
29  dx = [ 0 ; 0] ;
30 # printf("x(1)\tx(2)\tlambda\tdx1\tdx2\tnorm\n") ;
31  k = 0 ;
32  do
33  printf("%d\t%12.5e\t%12.5e\t%12.5e",k,xk(1),xk(2),lambdak) ;
34  [f,g,H] = feval(problem,xk,0) ;
35  contr = [] ;
37  hesslagrangien = H ;
38  for i=1:m
39  [fi,gi,Hi] = feval(problem,xk,i) ;
40  contr = [contr ; fi] ;
42  hesslagrangien += lambdak(i) * Hi ;
43  endfor
44
45  [L,tau] = modifiedCholesky(hesslagrangien) ;
46
47 # printf(" tau = %f\n", tau) ;
48  hesslagrangien = L * L' ;
49
52
53  %Update the penalty parameter
54  target = norm(dl,Inf) + cconst ;
55  if (c >= 1.1 * target)
56  % reduce c
57  c = (c + target) / 2.0 ;
58  else
59  if ( c < target)
60  c = max(1.5 * c, target) ;
61  endif
62  endif
63
64  normConstraint = norm(contr,1) ;
65  derivdirect = g' * dx - c * normConstraint ;
66
67  if (derivdirect < 0)
68  phi = f + c * normConstraint ;
69  alpha = 2.0 ;
70  do
71  alpha /= 2.0 ;
72  cand = xk + alpha * dx ;
73  fcand = feval(problem,cand,0) ;
74  contrcand = [] ;
75  for i=1:m
76  fi = feval(problem,cand,i) ;
77  contrcand = [contrcand ; fi] ;
78  endfor
79  phicand = fcand + c * norm(contrcand,1) ;
80  until (phicand <= phi + alpha * beta1 * derivdirect )
81  xk = xk + alpha * dx ;
82  lambdak = dl ;
83  k = k + 1 ;
84  endif