2 %> Algorithm 18.4:
long steps interior point algorithm. Implementation of algorithm 18.4 of \cite Bier15-book
4 %> @note Tested with \ref run1807longSteps.m
6 %> @author Michel Bierlaire
7 %> @date Mon Mar 23 14:52:36 2015
11 %> Applies the
long step 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 gamma parameter such that gamma > 0 (
default: 1.0e-3)
20 %> @param sigma parameter such that 0 < sigma < 1 (
default: 0.1)
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] =
longSteps(A,b,c,x0,lambda0,mu0,eps,gamma=1.0e-3,sigma=0.1,maxiter=100)
30 error(
"The number of rows in A and b do not match") ;
33 error(
"The number of columns in A and do not match the size of c") ;
36 error(
"The number of columns in A and do not match the size of x0") ;
38 if (m != size(lambda0))
39 error(
"The number of rows in A and do not match the size of lambda0") ;
42 error(
"The number of columns in A and do not match the size of mu0") ;
45 primalcons = norm(A*x0 - b) ;
47 error(
"Initial point is not primal feasible") ;
49 dualcons = norm(A
' * lambda0 + mu0 - c) ; 51 error("Initial point is not dual feasible") ; 55 error("x0 must be strictly positive") 58 error("mu0 must be strictly positive") 70 voisinage = max(x .* mu) / nu;
71 if (voisinage < gamma)
72 error(
"Iterate not in the vicinity of the central path: %f",voisinage)
75 printf(
"%d\t%e\t\t%e\n",k,voisinage,nu) ;
77 printf(
"%d\t%e\t%e\t%e\t%e\n",k,voisinage,norm(alpha*d),nu,alpha);
79 Q = [ A zeros(m,m) zeros(m,n) ; zeros(n,n) A
' eye(n,n) ; S zeros(n,m) X] ; 80 lefthand = [ zeros(m+n,1) ; -X*S*e + nu * sigma * e] ; 86 cx = x + alpha * d(1:n) ; 87 clambda = lambda + alpha * d(n+1:n+m) ; 88 cmu = mu + alpha * d(n+m+1:n+m+n) ; 89 if (min(cx) <= 0 || min(cmu) <= 0) 93 cnu = (cmu' * cx ) / n ;
101 cvoisinage = max(cx .* cmu) / cnu;
102 if (cvoisinage >= gamma)
109 alpha = alpha / 2.0 ;
114 until (nu <= eps || k >= maxiter)
function longSteps(in A, in b, in c, in x0, in lambda0, in mu0, in eps, in gamma)
Applies the long step interior point algorithm to solve subject to and .