Penalty- and Barrier - Method

Directly to the input form: predefined testcases

Directly to the input form: userdefined testcase

• Here you find a code for solving a constrained optimization problem using either the exterior penalty (quadratic loss) function or the logarithmic barrier method (mixed form for equality constraints). For basic facts on optimization problems of this kind: f(x) = min, h(x)=0, g(x)>=0
  look here.
• These two are the oldest methods in use for nonlinear optimization problems and ''state of the art'' up to the end of the 60's of the last century.
  Both approaches build a function whose local unconstrained minima approximate the solution of the original problem and obtain (at least theoretically) arbitrary precision as the so called penalty parameter (we name it ρ here) tends to infinity.
  If the approximated solution is a point satisfying the strong second order sufficiency conditions, then the final precision in x and the multipliers will be O(1/ρ). In other cases the final precision might be much worse. However, unfortunately the condition number of the Hessian of these unconstrained functions also is O(1/ρ) (with the exception that the solution x^* is a vertex), hence the unconstrained minimizer will run into more and more trouble and finally will probably fail. This limits the obtainable precision drastically and also degrades the reliability of these solvers.
• For the quadratic loss function one knows that the objective function increases during the the ρ-iteration, whereas the penalty term decreases, that means that intermediate values are infeasible and feasible only in the limit.
• We present two estimates for the multipliers: the least squares multipliers, computed using the final x^*(ρ), and the ''asymptotic multipliers'', obtained from the convergence theorem. These are
  λ_i(ρ) = - (2/ρ)min{0,g_i(x^*(ρ))}
  μ_j(ρ) = -(2/ρ)h_j(x^*(ρ)).
  Typically these asymptotic estimates are worse than the least squares multipliers.
• The quadratic loss function reads
  φ(x;ρ) = f(x)+ρ*(||h(x)||²+||g^(-)(x)||²)
  whereas the logarithmic barrier function is
  φ(x;ρ) = f(x)+ρ*(||h(x)||²)-(1/ρ)*∑_i=1,..,m log(g_i(x))
  
• Here ρ is held fixed for one unconstrained minimization and (theoretically) takes values from a sequence tending to infinity. Practically ρ is limited by the user defined parameter rhomax.
• The unconstrained minimization of φ uses the BFGS-minimizer, for its details see there.
• The parameter ρ takes values from rhomin up to rhomax, increased by the userdefined parameter rhofac.
• Termination occurs if either the desired precision is reached within the sequence of minimization steps or if the unconstrained minimizer stops unsuccessfully. The termination messages you might obtain are
  1. ''function evaluation signals error'' (you can set ''err=.true.'' inside your function programs, if for example the optimizer provides you with x(1)=-1, but you need log(x(1)))
  2. ''too many steps for rho fixed (slow convergence)'' (the unconstrained minimizer did not reach its normal termination within 4000 steps)
  3. ''backtracking failed (possibly rho too large?)'' (the stepsize scheme of the BFGS minimizer failed. Since we compute gradients analytically this occurs typically if the conditioning of the problem has become so bad that the ''valleys'' become too narrow to identify an unidimensional minimizer)
  4. ''inner iteration : sufficient precision'' (successful termination w.r.t. x),
  5. ''inner iteration : successful termination'' (successful termination w.r.t. ∇ φ)
  6. ''inner iteration : directional derivative small'' (the directional derivative has been reduced to the roundoff level of φ)
  7. ''inner iteration : no decrease in penalty obtained'' (BFGS step failed, surely ρ already too large)
  8. ''too many steps in total '' (since we store some information of all steps for providing the graphical output, the total number of steps is limited to 99999 )
• you can play with predefined testcases, see testcases for this, or define a problem of your own.
  For the log-barrier function the initial point must be strictly feasible for all constraints, which is not the case for many of the predefined testcases. Hence you must use the feature of modifying the initial point in these cases.

There are two different input pages: Predefined und Userdefined

Input - predefined variant

• Select a problem number from testcases The problem number always looks like hs followed by a number and then sometimes followed by scal (=scaled). You must write the problem number exactly as seen there. The scaled versions allow you to study the effects of bad or good scaling of a problem, which in real life applications is of outmost importance.
• you choose between the penalty and the log barrier function.
• You have the choice to see the intermediate results displayed or only the final result.
• There is a standard initialization for the algorithmic parameters but you can use or own choice:
  1. epsx = ε: the inner minimization terminates successfully if
     ||∇ φ(x)|| <= ε*(1+ε*|φ|)
     or
     ||x-x_previous|| <= ε*(1+||x||).
     Important: epsx > 0 !
  2. rhomin, rhomax: the minimal and the maximal value of the parameter ρ.
  3. rhofac: factor by which ρ is increased. Important : rhofac > 1 !
  4. you may change the initial point. If so, then of course you must type n values here.
  The outer iteration w.r.t. ρ is terminated successfully, if the primal and dual infeasibility are less than ε and the norm of the gradient of the Lagrange function is <= max(ε,epsmach*100*cond), but at least if rhomax exceeded, which is interpreted as failure. cond is an estimate of the condition number of the Hessian of the penalty- resp. barrier function, obtained from the BFGS update. This condition number grows with ρ growing to infinity itself to infinity. The Lagrange multiplier used in the termination test are the least squares multipliers.
  ||∇ f(x) - N (μ ,λ)||₂ = min_{μ , λ} .
  Here N is the matrix of the gradients of binding constraints at x.
  Inequalites with values <= ε are considered as binding and of course all equality constraints.
• You might require a contour plot of the penalty- resp. barrier function w.r.t. two components of x with the other components held at their optimal values. In this case you specify these two components and the width of the enclosing box you want to see. you might also require that bounds on the variables will not be violated in the plot. For the barrier function we use outside the feasible domain a very large value in order to indicate ''undefined''.

Input: user defined variant

• You first specify the dimensions of your problem: n the number of variables and
• nh, the number of equality constraints (can be zero) and ng the number of inequality constraints, can also be zero but (nh+ng > 0 )!
• Then you specify the objective function fx and the functions decribing the equality constraints hxi and the inequality constraints gxi as described here.
• You choose between the penalty or the log barrier function.
• you must specify your initial guess : of course you must type here n numbers.
• you choose the amount of displayed output.
• There is a standard initialization for the algorithmic parameters but you can use or own choice:
  1. epsx = ε: the inner minimization terminates successfully if
     ||∇ φ(x)|| <= ε*(1+ε*|φ|)
     or
     ||x-x_previous|| <= ε*(1+||x||).
     Important: epsx > 0 !
  2. rhomin, rhomax: the minimal and the maximal value of the parameter ρ.
  3. rhofac: factor by which ρ is increased. Important : rhofac > 1 !
• You might require a contour plot of the penalty- resp. barrier function w.r.t. two components of x with the other components held at their optimal values. In this case you specify these two components and the width of the enclosing box you want to see. you might also require that bounds on the variables will not be violated in the plot. For the barrier function we use outside the feasible domain a very large value in order to indicate ''undefined''.

Output

• You get the final result x^*, the two estimates of the Lagrange multipliers and the optimal function value as well as the gradient of the Lagrange function.
• If the iteration protocol had been required then you get for the inner iteration, that is for ρ fixed, a table with the values of the penalty- resp. barrier function, the norm of its gradient and the value of the infeasibility ||(h(x),g^(-)(x)||. Since ρ is growing in the outer iteration, the gradient doesn't become really ''small''.
• A plot showing the decadic logarithm of the penalty/barrier function (minus its minimal value) and of the penalty / barrier term (ψ). The jumps in these graphs appear at the steps where ρ is increased.
• If required, the contour plot.

To the input form: predefined testcases

To the input form: userdefined testcase