Optimization: principles and algorithms, by Michel Bierlaire
restrictedSteps.m
Go to the documentation of this file.
1 %> \file
2 %> Algorithm 18.2: Interior point algorithm with restricted steps. Implementation of algorithm 18.2 of \cite Bier15-book
3 %>
4 %> @note Tested with \ref run1807restrictedSteps.m
5 %>
6 %> @author Michel Bierlaire
7 %> @date Mon Mar 23 10:59:13 2015
8 %> @ingroup Algorithms
9 %> @ingroup chap18
10
11 %> Applies the interior point algorithm with restricted steps 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 theta parameter such that 0 <= theta <= 1 (default: 0.4)
20 %> @param maxiter maximum number of iterations (default: 100)
21 %> @return [x,lambda,mu]
22 %> @return x: primal solution
23 %> @return lambda: dual solution for equality constraints
24 %> @return mu: dual solution for inequality constraints
25 function [x,lambda,mu] = restrictedSteps(A,b,c,x0,lambda0,mu0,eps, theta=0.4,maxiter=100)
26
27  [m,n] = size(A) ;
28  if (m != size(b))
29  error("The number of rows in A and b do not match") ;
30  endif
31  if (n != size(c))
32  error("The number of columns in A and do not match the size of c") ;
33  endif
34  if (n != size(x0))
35  error("The number of columns in A and do not match the size of x0") ;
36  endif
37  if (m != size(lambda0))
38  error("The number of rows in A and do not match the size of lambda0") ;
39  endif
40  if (n != size(mu0))
41  error("The number of columns in A and do not match the size of mu0") ;
42  endif
43
44  primalcons = norm(A*x0 - b)
45  if (primalcons > eps)
46  error("Initial point is not primal feasible") ;
47  endif
48  dualcons = norm(A' * lambda0 + mu0 - c) ;
49  if (dualcons > eps)
50  error("Initial point is not dual feasible") ;
51  endif
52
53  sigma = 1.0 - theta / sqrt(n) ;
54  x = x0 ;
55  lambda = lambda0 ;
56  mu = mu0 ;
57
58  e = ones(n,1) ;
59
60  k=0 ;
61  do
62  nu = (mu' * x ) / n ;
63  S = diag(mu) ;
64  X = diag(x) ;
65  voisinage = norm(X*S*e-nu*e) / nu ;
66  if (voisinage > theta)
67  error("Iterate not in the vicinity of the central path")
68  endif
69  if (k==0)
70  printf("%d\t%e\t\t%e\n",k,voisinage,nu) ;
71  else
72  printf("%d\t%e\t%e\t%e\n",k,voisinage,norm(d),nu) ;
73  endif
74  Q = [ A zeros(m,m) zeros(m,n) ; zeros(n,n) A' eye(n,n) ; S zeros(n,m) X] ;
75  lefthand = [ zeros(m+n,1) ; -X*S*e + nu * sigma * e] ;
76  d = Q \ lefthand ;
77  x = x + d(1:n) ;
78  lambda = lambda + d(n+1:n+m) ;
79  mu = mu + d(n+m+1:n+m+n) ;
80  k=k+1 ;
81  until (nu <= eps || k >= maxiter)
82 endfunction
function restrictedSteps(in A, in b, in c, in x0, in lambda0, in mu0, in eps, in theta)
Applies the interior point algorithm with restricted steps to solve subject to and ...