`constrained_optimization`

This library complements local_optimization (which only handles box constraints via position_bounds/1) with general equality and inequality constraints.

API documentation

Open the ../../apis/library_index.html#constrained-optimization link in a web browser.

Loading

To load all entities currently in this library, load the loader.lgt file:

| ?- logtalk_load(constrained_optimization(loader)).

Testing

To test this library’s predicates, load the tester.lgt file:

| ?- logtalk_load(constrained_optimization(tester)).

Solvers

`qp_active_set`

qp_active_set is a dense primal active-set solver for convex quadratic programs, following Nocedal and Wright, section 16.5. It repeatedly solves the KKT system for the current working set. A nonzero step is taken either fully or only as far as the nearest newly binding inactive inequality. At a zero step, the first active inequality with a negative multiplier in row order is dropped, following Bland’s rule. Equalities are always active and are never dropped.

The initial feasible point is found heuristically. The solver starts from the minimum-norm point satisfying the equalities, or the origin when there are no equalities, then repeatedly adds violated inequalities to that system. It verifies feasibility both after this phase and before returning. If no feasible point is found, the active-set iteration limit is reached, or a working-set KKT system is singular, solve/8 fails rather than returning an infeasible or non-optimal point.

Linear systems use linear_algebra::solve_linear_system/3; the phase-1 minimum-norm point uses linear_algebra::pseudo_inverse/2. The solver is intended for small dense QPs and is also used internally by sqp_active_set(_).

Current limitation: phase 1 is not a complete two-phase or elastic-mode method and can fail to find a feasible point even when one exists.

`sqp_active_set(Problem)`

This Sequential Quadratic Programming solver builds and solves a qp_active_set subproblem at each iteration. It maintains a damped-BFGS approximation to the Hessian of the Lagrangian. The secant pair uses Lagrangian gradients and Powell damping; multipliers for box constraints are omitted because those constraints are linear.

Globalization uses the L1 exact-penalty merit function f(x) + mu*(||g(x)||_1 + ||max(0,h(x))||_1) with Armijo backtracking. The penalty weight is increased as required by the current QP multipliers and is never decreased. A failed QP solve or exhausted line search causes run/4 to fail.

Bounds from position_bounds/1 are converted into exact linear inequality rows on the QP step. Convergence requires both a sufficiently small QP step and acceptable constraint violation. Statistics include final_step_norm/1, final_constraint_violation/1, and termination_reason/1; possible reasons are converged, target_reached, stop_condition, and max_iterations.

This solver is generally the best choice for problems close to a quadratic objective with linear constraints. Current limitations: there is no feasibility-restoration phase when a nonlinear constraint linearization produces an infeasible QP, and the QP working set is not warm-started between iterations or runs.

`augmented_lagrangian(Problem, InnerSolver)`

This method of multipliers delegates each outer subproblem to an existing local_optimization solver such as bfgs or lbfgs. Equality constraints use the Hestenes-Powell term lambda.g(x) + (rho/2)*||g(x)||^2. Inequalities use the smooth Rockafellar-Bertsekas term based on max(0, lambda_j + rho*h_j(x)). After each inner solve, equality multipliers are updated by lambda + rho*g(x) and inequality multipliers by max(0, lambda + rho*h(x)).

The penalty parameter grows only when the new constraint violation fails to decrease by the configured factor. Bounds are forwarded unchanged to the inner solver instead of being penalized. The inner solver always minimizes because the subproblem already applies the sign required for the requested objective direction.

max_iterations, tol_x, tol_f, and tol_g apply to each inner solve; max_outer_iterations and outer_tolerance control the outer loop. updates(N) controls outer progress/5 callbacks, while inner_updates(N) controls delegated inner_progress/6 callbacks. Inner callbacks use the stage outer(K) and report the transformed subproblem objective and convergence measure. Both options default to zero, disabling their respective callbacks. Statistics include final_violation/1 and termination_reason/1, with converged, target_reached, or max_iterations as possible reasons.

This is the usual default among the delegated outer-loop methods. Current limitation: inner-solver state such as a BFGS approximation is rebuilt for every outer iteration.

`quadratic_penalty(Problem, InnerSolver)`

This solver delegates minimization of f(x) + (rho/2)*(||g(x)||^2 + ||max(0,h(x))||^2) to an existing local_optimization solver. It has the same outer-loop and bound handling as augmented_lagrangian(_, _), but maintains no multiplier estimates. The penalty parameter grows when constraint violation does not decrease by the configured factor.

The inner and outer option scopes and termination reasons are the same as for augmented_lagrangian(_, _), including the independent outer and inner progress controls. Statistics include final_violation/1 and termination_reason/1.

The method is simple and can suit small, well-conditioned problems, but augmented_lagrangian(_, _) is normally preferable. Current limitation: exact feasibility is reached only as rho tends to infinity. Large penalties increasingly ill-condition the inner problem, so an insufficient outer budget or unsuitable penalty schedule can leave a material residual.

`log_barrier(Problem, InnerSolver)`

This primal interior-point method handles inequalities using -mu*sum(log(-h_j(x))) and handles equalities using a quadratic penalty. Each outer iteration delegates the resulting subproblem to an existing local_optimization solver, then unconditionally shrinks mu and grows the equality penalty only when equality violation is not decreasing fast enough. Bounds are forwarded to the inner solver.

If the initial point is not strictly inequality-feasible, the solver first minimizes a squared-violation phase-1 subproblem. If that heuristic does not produce a strictly feasible point, run/4 raises domain_error(strictly_feasible_initial_point, Point). During the main solve, trial points outside the barrier domain use a smooth quadratic fallback so an inner line search cannot trigger an arithmetic exception. Outer and inner progress use the same independent controls as the other delegated solvers. Main inner solves use the stage outer(K); phase-1 callbacks use phase1.

This solver is useful when the objective is defined only in the interior of the feasible region. Statistics report final_mu/1 and final_equality_violation/1 separately, plus termination_reason/1 with converged, target_reached, or max_iterations.

Current limitation: phase 1 is heuristic and may fail on a feasible nonlinear problem. The inner solver is also not domain-aware; the smooth fallback protects trial evaluations but is not a replacement for a feasibility-preserving line search.

`primal_dual_interior_point(Problem)`

This infeasible-start method implements Mehrotra’s predictor-corrector (Nocedal and Wright, chapter 19, algorithm 19.3). At each iteration it solves an affine-scaling predictor system, derives the cubic centering parameter, then solves a corrector system with the centering and second-order correction terms. Slack and inequality multiplier steps use the fraction-to-boundary rule. A damped-BFGS approximation represents the Hessian of the Lagrangian, and an L1 merit function globalizes the step.

Slack variables allow the problem’s initial point to violate equality, inequality, and box constraints. Bounds are folded into the same reduced KKT system as constant-Jacobian inequalities. The centering quantity sigma*mu is distinct from the merit-function penalty and is recomputed from the predictor step every iteration.

Statistics report final_stationarity_residual/1, final_primal_infeasibility/1, final_complementarity_gap/1, and termination_reason/1. Reasons are converged, target_reached, stop_condition, max_iterations, or kkt_singular.

Current limitations: KKT solves are dense and have no scaling or iterative refinement; an accurate point can therefore be returned with termination_reason(kkt_singular) before every requested residual tolerance is met. The merit line search retains a safeguarded fallback because strict rejection of all non-descent or non-Armijo trials is not robust on the supported nonlinear and maximization cases. Like the other nonlinear solvers, this is a local method and may converge to any reachable KKT point rather than a global optimum.

Architecture

constrained_optimization_solver - a category providing shared numeric, objective-direction, list, inner-solver construction, penalty-update, and validation auxiliary predicates used by the solvers and their internal problem wrappers. Imported alongside local_optimization_solver(_Problem_) wherever both are needed, via imports([constrained_optimization_solver, local_optimization_solver(_Problem_)]); called via ^^, like any other category predicate.
constrained_optimization_problem_protocol - extends local_optimization_problem_protocol with equality_constraints/2, equality_jacobian/2, inequality_constraints/2, and inequality_jacobian/2. Box constraints continue to be expressed via the inherited position_bounds/1. Implemented by problems passed to any of the five problem-object solvers.
qp_solver_protocol - protocol for QP subroutines operating on plain matrices/vectors (minimize 0.5*x^T*H*x + c^T*x s.t. Aeq*x = beq, Aineq*x =< bineq), as opposed to the problem-object-based protocols above. Implemented by qp_active_set.
The three solvers that delegate to an inner local solver use internal problem wrappers: sub_problem/6 for augmented_lagrangian(_, _), penalty_sub_problem/4 for quadratic_penalty(_, _), and barrier_sub_problem/5 plus phase1_sub_problem/2 for log_barrier(_, _). The wrappers expose the transformed objective and constraints expected by the selected inner solver.

Defining a QP directly

qp_active_set does not use a problem object; solve/8 takes the problem in standard form directly:

minimize    0.5 * x^T * H * x + c^T * x
subject to  Aeq * x = beq
            Aineq * x =< bineq

H - dense N x N matrix (list of rows), symmetric positive semi-definite over the feasible set.
C - length-N vector.
Aeq, Beq - equality constraint matrix (Meq x N, possibly []) and right-hand side (length Meq).
Aineq, Bineq - inequality constraint matrix (Mineq x N, possibly [], rows in the form row.x =< b) and right-hand side (length Mineq).

solve/8 fails (rather than raising an error) when no feasible point is found, the problem is unbounded below, or the current working set becomes numerically singular; it raises a domain_error/2 when the input matrices/vectors have inconsistent dimensions.

| ?- qp_active_set::solve(
         [[1.0,0.0],[0.0,1.0]], [0.0,0.0],
         [], [],
         [[-1.0,-1.0]], [-2.0],
         X, Lambda
     ).
X = [1.0, 1.0],
Lambda = [1.0].

% minimize 0.5*(x1^2+x2^2) subject to x1+x2 >= 2: the closest point
% on that line to the origin.

Solving a nonlinear constrained problem

All five nonlinear solvers take a problem object implementing constrained_optimization_problem_protocol, in the same way bfgs(_)/lbfgs(_) take a local_optimization_problem_protocol object:

| ?- sqp_active_set(my_problem)::run(BestPoint, BestValue).

| ?- augmented_lagrangian(my_problem, bfgs)::run(BestPoint, BestValue).

sqp_active_set(_) and primal_dual_interior_point(_) require the corresponding Jacobian whenever a constraint predicate is defined. The three delegated solvers require Jacobians when their selected inner solver uses gradients. A missing required Jacobian raises existence_error/2 rather than silently falling back to finite differences.

Library-wide limitations

All solvers use dense matrices. Problem sizes are therefore expected to remain within the range where dense KKT systems and dense inner-solver Hessian approximations are practical. An initial_point(Point) option can warm-start the point, but no solver preserves internal state such as a QP working set, dual and slack variables, or an inner solver’s Hessian approximation between runs.

Usage

Solving a QP directly

| ?- qp_active_set::solve(
         [[1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]], [0.0,0.0,-1.0],
         [[1.0,1.0,1.0]], [3.0],
         [[-1.0,0.0,0.0],[0.0,-1.0,0.0],[0.0,0.0,1.0]], [0.0,0.0,1.0],
         X, _Lambda
     ).
X = [1.0, 1.0, 1.0].

% minimize 0.5*(x1^2+x2^2+x3^2) - x3
% subject to x1+x2+x3 = 3, x1 >= 0, x2 >= 0, x3 =< 1

Solving the same problem with SQP

:- object(my_problem,
    implements(constrained_optimization_problem_protocol)).

    initial_point([0.0, 0.0, 0.0]).
    objective(X, V) :- X = [X1,X2,X3], V is 0.5*(X1*X1+X2*X2+X3*X3) - X3.
    gradient(X, [X1,X2,G3]) :- X = [_,_,X3], G3 is X3 - 1.0.
    equality_constraints(X, [G]) :- X = [X1,X2,X3], G is X1+X2+X3-3.0.
    equality_jacobian(_, [[1.0,1.0,1.0]]).
    inequality_constraints(X, [H1,H2,H3]) :- X = [X1,X2,X3],
        H1 is -X1, H2 is -X2, H3 is X3-1.0.
    inequality_jacobian(_, [[-1.0,0.0,0.0],[0.0,-1.0,0.0],[0.0,0.0,1.0]]).

:- end_object.

| ?- sqp_active_set(my_problem)::run(X, V).
X = [1.0, 1.0, 1.0].

Solving it with augmented_lagrangian instead

| ?- augmented_lagrangian(my_problem, bfgs)::run(X, V).
X = [1.0, 1.0, 1.0].

Solving it with quadratic_penalty instead

| ?- quadratic_penalty(my_problem, bfgs)::run(X, V).
X = [1.0, 1.0, 1.0].

Solving it with log_barrier instead

initial_point/1 need not be strictly feasible - a heuristic phase 1 runs first if it is not, but supplying one directly, when convenient, skips that extra work:

| ?- log_barrier(my_problem, bfgs)::run(X, V).
X = [1.0, 1.0, 1.0].

Solving it with primal_dual_interior_point instead

No inner solver to name, and initial_point/1 need not be feasible at all - [0.0, 0.0, 0.0] (which violates x1 >= 0 exactly, on the boundary, and doesn’t satisfy the equality either) works fine:

| ?- primal_dual_interior_point(my_problem)::run(X, V).
X = [1.0, 1.0, 1.0].

Warm-starting from a previous solution

Every solver in this library accepts an initial_point(Point) option (from local_optimization_solver(_Problem_)) that overrides Problem’s own initial_point/1 - useful when re-solving a problem that changed only slightly from one already solved:

| ?- sqp_active_set(my_problem)::run(X, V, _Statistics, [initial_point([1.0, 1.0, 1.0])]).
X = [1.0, 1.0, 1.0].

This overrides only the starting point, not any solver-internal state (a QP working set, dual variables, or an inner solver’s Hessian approximation).

constrained_optimization