Linear Programming
Mathwrist C++ API Series

Copyright ©Mathwrist LLC 2023
(January 20, 2023)
Abstract

This document details the C++ programming interface of linear programming (LP) features in Mathwrist’s Numerical Programming Library (NPL). Please refer to our linear programming white paper for technical overview.

1 Introduction

Linear programming (LP) solves the following optimization problem in a general form.

min𝐱∈ℝn⁡𝐜T⁢𝐱, s.t. ⁢𝐛l≤𝐀𝐱≤𝐛u⁢ and ⁢𝐥≤𝐱≤𝐮 (1)

, for some cost vector 𝐜 to minimize. 𝐛l≤𝐀𝐱≤𝐛u is referred as general linear constraints. 𝐥≤𝐱≤𝐮 is referred as simple bound constraints. This formulation is represented by the LinearProg base class in header file linear_prog.h. NPL provides two LP solvers based on the classic Simplex method and active-set method respectively. Both solvers are derived from the LinearProg base class.

1.1 Constraint Specification

The constraint specification in the general form (1) is represented by a nested structure LinearProg::Constraints, which has the following data members. Let m denote the total number of general constraints. Let n denote the total number of unknowns.

1.1.1 Data Members

1 Matrix A; // Constraint coefficient matrix
2 ftl::Buffer<Bound> simple_bounds;
3 ftl::Buffer<Bound> general_bounds;

Description:

  1. 1.

    A is the m×n coefficient matrix 𝐀 of the general linear constraints.

  2. 2.

    simple_bounds holds n Bound objects representing (𝐥,𝐮). REQUIRED: Each simple bound cannot have its lower bound equal to the upper bound.

  3. 3.

    general_bounds holds m Bound objects representing (𝐛l,𝐛u). A general bound could have its lower bound equal to the upper bound. In which case it indicates the corresponding constraint is an equality constraint, 𝐀i⁢𝐱=𝐛i.

Either one of the simple_bounds buffer or general_bounds buffer could be empty to indicate unknowns are unbounded or there are no general linear constraints.

1.1.2 Constructors

1 Constraints();
2 Constraints(size_t m, size_t n);
3 explicit Constraints(const Matrix& _A);
4 Constraints(const Matrix& _A, const ftl::Buffer<Bound>& _b,
5             bool non_neg_x = false);

Description:

  1. 1.

    Default constructor, constraints are empty.

  2. 2.

    Create a constraint specification that has m×n general linear constraints and n simple bounds. Users need set the coefficient matrix and bounds at a later time.

  3. 3.

    Initialize with a given m×n coefficient matrix 𝐀 passed from parameter _A. Users need set the bounds at a later time.

  4. 4.

    Initialize with a given a m×n coefficient matrix 𝐀 passed from parameter _A and the general constraint bounds passed from parameter _b. Parameter non_neg_x indicates if the unknowns are subject to non-negative constraints 𝐱≥0. REQUIRED: Parameter _b needs have buffer size m.

1.1.3 Resize and Validate

1 void resize(size_t m, size_t n);
2 void validate() const;

Description:

  1. 1.

    Resize the constraint specification to take m×n general linear constraints and n simple bounds.

  2. 2.

    Validate the constraint specification. This includes: 1) if both general constraints and simple bounds are provided, their sizes are consistent; 2) if simple bounds are provided, they cannot be equality constraints. An Exception will be thrown if the validation fails.

1.1.4 Feasibility Functions

1 void feasible_point(Matrix& x, size_t max_iter = 0) const;
2 bool is_feasible(
3    const Matrix& x,
4    const Matrix& x_tol,
5    const Matrix& gc_tol) const;

Description:

  1. 1.

    Use phase-I Simplex method to find the initial feasible point. Parameter x holds the initial guess on call and feasible point solution on return. Parameter max_iter is the max number of iterations allowed. If max_iter is not provided, it is taken from the IterativeMethod’s default value.

  2. 2.

    Returns if a given point, parameter x, is a feasible point with respect to the current constraints. Parameter x_tol holds n elements of tolerance to test the feasibility of simple bounds. Parameter gc_tol holds m elements of tolerance to test the feasibility of general linear constraints.

1.2 LinearProg Class Member Functions

1 protected:
2    LinearProg();
3    LinearProg(const Constraints& constraint);
4 public:
5    virtual ~LinearProg();
6    const Constraints& constraint() const;
7    virtual double operator() (const Matrix& c, Matrix& x) = 0;

Description:

LinearProg class is a pure virtual base class for LP solvers. Its constructors are all protected.

  1. 1.

    Default constructor.

  2. 2.

    Initialize the object with a given constraint specification.

  3. 3.

    Destructor.

  4. 4.

    Return a reference to the constraint specification held inside the LP solver.

  5. 5.

    Interface function to solve the LP problem with a given cost vector c and initial guess x. It returns the objective function value and saves the solution to parameter x.

2 Simplex Method

The Simplex solver is defined in header file simplex.h. It is derived from the LinearProg and IterativeMethod base classes. Please refer to LinearProg for linear constraint specification. Please refer to IterativeMethod for iteration and convergency control.

1 class Simplex : public LinearProg, public IterativeMethod
2 {
3   ...
4 };

2.1 Constructor and Destructor

1 explicit Simplex(const Constraints& constraint);
2 ~Simplex();

Description:

  1. 1.

    Initialize the Simplex solver object with a given constraint specification.

  2. 2.

    Free all internal resources.

2.2 Phase-I Feasibility

1 bool phaseI(Matrix& x);

Description:

Given an initial guess passed from parameter x, solve the phase-I feasibility problem with the solution being saved to parameter x on return. This function returns true if the Simplex algorithm finds a basic feasible point before the max number of iterations is reached.

2.3 Solver Operator

1 double operator() (const Matrix& c, Matrix& x);

Description:

Please refer to LinearProg::operator() for function details.

3 Active Set Method

The ActiveSet base class in header file active_set.h is derived from IterativeMethod class. It provides a framework implementation that is shared by various optimization solvers. Some public functions in this class are for internal use only, hence not documented in this section.

3.1 Crash Start Radius

1 void set_crash_start_radius(double r);
2 double crash_start_radius() const;

Description:

Please refer to our LP white paper for the details of active set “crash start” procedure. This function sets the “crash start” radius r that is used to identify the initial working set of constraints. Parameter r is a positive number. The larger the radius, the more constraints are potentially included in the initial working set.

  1. 1.

    Set the crash start radius.

  2. 2.

    Return the crash start radius used by the active-set method.

3.2 Stationary Tolerance

1 void set_stationary_tolerance(double eps);
2 double stationary_tolerance() const;

Description:

  1. 1.

    Set the stationary tolerance ϵ that is used to test if a search direction 𝐩 is a zero vector, |𝐩|∞<ϵ.

  2. 2.

    Return the stationary tolerance used by the active-set method.

3.3 Feasibility Tolerance

1 void set_feasibility_tolerance(double eps);
2 double feasibility_tolerance() const;

Description:

  1. 1.

    Set the tolerance control level that is used to test if the problem has a feasibile solution.

  2. 2.

    Return the current feasibility control tolerance.

3.4 Elastic Programming Control

1 double elastic_penalty() const;
2 bool in_ep_mode() const;
3 void set_elastic_penalty(double p);

Description:

This set of functions are only used by Sequential Quadratic Programming (SQP) solver at the moment. When the active set method is in elastic programming mode, it penalizes the violated elastic constraints. The penalty amount is a part of the minimization objective.

  1. 1.

    Return the factor to penalize elastic constraints that are violated.

  2. 2.

    Return if the active-set method is currently in elastic programming mode.

  3. 3.

    Set the penalty factor for elastic programming control.

3.5 Abnormal Lagrange Threshold

1 void set_abnormal_lagrange_threshold(double l);
2 double abnormal_lagrange_threshold() const;

Description:

At a stationary point, Lagrange multipliers are tested. If the magnitude of Lagrange multiplers is extremely large, it is very likely that the algorithm has accumulated enough numerical error. The abnormal Lagrange threshold is used to identify this situation.

  1. 1.

    Set the abnormal Lagrange threshold. Parameter l is the threshold provided by users. The actual threshold being used is max⁡(l,1.0⁢e8).

  2. 2.

    Return the abnormal Lagrange threshold used by the algorithm.

4 LP Active Set

The active set framework provides a solver LP_ActiveSet class in lp_active_set.h for LP problems. It is derived from LinearProg and ActiveSetWithLagrange classes.

1 class LP_ActiveSet : public LinearProg,
2                      public ActiveSetWithLagrange
3 {
4   ...
5 };

ActiveSetWithLagrange is an intermediate derived class from the ActiveSet base class. It serves for internal implementation and doesn’t have public functions, hence not documented.

4.1 Constructor

1 explicit LP_ActiveSet(const LinearProg::Constraints& constr);

Description:

Initialize the active set LP solver with a given constraint specification.

4.2 Solver Operator

1 double operator() (const Matrix& c, Matrix& x);

Description:

Please refer to LinearProg::operator() for function details.