Optimization: principles and algorithms, by Michel Bierlaire
predictorCorrector.m
Go to the documentation of this file.
1 %> \file
2 %> Algorithm 18.3: Predictor-corrector interior point algorithm. Implementation of algorithm 18.3 of \cite Bier15-book
3 %>
4 %> @note Tested with \ref run1807predictorCorrector.m
5 %>
6 %> @author Michel Bierlaire
7 %> @date Mon Mar 23 12:10:29 2015
8 %> @ingroup Algorithms
9 %> @ingroup chap18
10
11 %> Applies the predictor-corrector interior point algorithm to solve \f[\min_x c^Tx \f] subject to \f[Ax = b\f] and \f[ x \geq 0 \f].
12 %> @param A the constraint matrix
13 %> @param b the constraint right hand side
14 %> @param c the cost vector for the objective function
15 %> @param x0 starting primal point (nx1)
16 %> @param lambda0 starting dual point for equality constraints (mx1)
17 %> @param mu0 starting dual point for inequality constraints (nx1)
18 %> @param eps algorithm stops if \f$\|d_k\| \leq \varepsilon \f$.
19 %> @param thetapred parameter such that 0 <= theta <= 1 (default: 0.5)
20 %> @param thetacorr parameter such that 0 <= theta <= 1 (default: 0.25)
21 %> @param maxiter maximum number of iterations (default: 100)
22 %> @return [x,lambda,mu]
23 %> @return x: primal solution
24 %> @return lambda: dual solution for equality constraints
25 %> @return mu: dual solution for inequality constraints
26 function [x,lambda,mu] = predictorCorrector(A,b,c,x0,lambda0,mu0,eps,thetapred=0.5,thetacorr=0.25,maxiter=100)
27
28  [m,n] = size(A) ;
29  if (m != size(b))
30  error("The number of rows in A and b do not match") ;
31  endif
32  if (n != size(c))
33  error("The number of columns in A and do not match the size of c") ;
34  endif
35  if (n != size(x0))
36  error("The number of columns in A and do not match the size of x0") ;
37  endif
38  if (m != size(lambda0))
39  error("The number of rows in A and do not match the size of lambda0") ;
40  endif
41  if (n != size(mu0))
42  error("The number of columns in A and do not match the size of mu0") ;
43  endif
44
45  primalcons = norm(A*x0 - b) ;
46  if (primalcons > eps)
47  error("Initial point is not primal feasible") ;
48  endif
49  dualcons = norm(A' * lambda0 + mu0 - c) ;
50  if (dualcons > eps)
51  error("Initial point is not dual feasible") ;
52  endif
53
54  if (min(x0) <= 0)
55  error("x0 must be strictly positive")
56  endif
57  if (min(mu0) <= 0)
58  error("mu0 must be strictly positive")
59  endif
60
61
62
63  x = x0 ;
64  lambda = lambda0 ;
65  mu = mu0 ;
66  e = ones(n,1) ;
67  sigma = 0 ;
68
69  k=0 ;
70  do
71  nu = (mu' * x ) / n ;
72  S = diag(mu) ;
73  X = diag(x) ;
74  voisinage = norm(X*S*e-nu*e) / nu ;
75  if (voisinage > thetacorr)
76  error("Iterate outside the neighborhood: %e > %f",voisinage,thetacorr)
77  endif
78  if (k==0)
79  printf("%d\t%e\t\t%e\n",k,voisinage,nu) ;
80  else
81  if (sigma == 0)
82  printf("%d\t%e\t%e\t%e\t%e\tCorr\n",k,voisinage,norm(alpha*d),nu,alpha);
83  else
84  printf("%d\t%e\t%e\t%e\t%e\tPred\n",k,voisinage,norm(alpha*d),nu,alpha);
85  endif
86  endif
87  Q = [ A zeros(m,m) zeros(m,n) ; zeros(n,n) A' eye(n,n) ; S zeros(n,m) X] ;
88  lefthand = [ zeros(m+n,1) ; -X*S*e + nu * sigma * e] ;
89  d = Q \ lefthand ;
90  alpha = 1.0 ;
91  if (sigma == 0)
92  # Prediction step
93  ok = 0 ;
94  while (ok != 1)
95  # Candidate
96  cx = x + alpha * d(1:n) ;
97  clambda = lambda + alpha * d(n+1:n+m) ;
98  cmu = mu + alpha * d(n+m+1:n+m+n) ;
99  if (min(cx) <= 0 || min(cmu) <= 0)
100  ok = 0 ;
101  alpha = alpha / 2.0 ;
102  else
103  cnu = (cmu' * cx ) / n ;
104  if (cnu <= eps)
105  # Solution found
106  x = cx ;
107  lambda = clambda ;
108  mu = cmu ;
109  return ;
110  endif
111  cS = diag(cmu) ;
112  cX = diag(cx) ;
113  cvoisinage = norm(cX*cS*e-cnu*e) / cnu ;
114  if (cvoisinage <= thetapred)
115  x = cx ;
116  lambda = clambda ;
117  mu = cmu ;
118  ok = 1 ;
119  else
120  ok = 0 ;
121  alpha = alpha / 2 ;
122  endif
123  endif
124  endwhile
125  else
126  # Correction step
127  x = x + d(1:n) ;
128  lambda = lambda + d(n+1:n+m) ;
129  mu = mu + d(n+m+1:n+m+n) ;
130  endif
131  k=k+1 ;
132  sigma = 1 - sigma ;
133  until (nu <= eps || k >= maxiter)
134 endfunction
135
function predictorCorrector(in A, in b, in c, in x0, in lambda0, in mu0, in eps, in thetapred)
Applies the predictor-corrector interior point algorithm to solve subject to and ...