Starter Guide of PRINTEMPS

Introduction

PRINTEMPS can compute approximate solutions for integer linear programming problems using Weighted Tabu Search [1]. To solve the problem, the solver replaces constraints with penalty functions that return violations of the corresponding constraints, and these penalty functions, scaled by positive penalty coefficients, are added to the objective function. For more information on the algorithm, please refer to the paper [2]. The figure below shows the flow-chart of the algorithm of PRINTEMPS.

The flow-chart of the optimization algorithm of PRINTEMPS.

Let us consider a simple linear integer optimization problem [3]:

(P):   minimize       x_1 + 10 x_2
          x
      subject to   66 x_1 + 14 x_2 >= 1430,
                  -82 x_1 + 28 x_2 >= 1306,
                      x_1 and x_2 are integer.

The following code shows an implementation to solve the problem (P) using PRINTEMPS.

#include <printemps.h>
int main(void) {
    // (1) Modeling
    printemps::model::IPModel model;

    auto& x = model.create_variables("x", 2);
    auto& g = model.create_constraints("g", 2);

    g(0) =  66 * x(0) + 14 * x(1) >= 1430;
    g(1) = -82 * x(0) + 28 * x(1) >= 1306;
    model.minimize(x(0) + 10 * x(1));

    // (2) Running Solver
    auto result = printemps::solver::solve(&model);

    // (3) Accessing the Result
    std::cout << "objective = " << result.solution.objective() << std::endl;
    std::cout << "x(0) = "      << result.solution.variables("x").values(0) << std::endl;
    std::cout << "x(1) = "      << result.solution.variables("x").values(1) << std::endl;

    return 0;
}

By compiling the code above and then running the generated executable, we will obtain the following result:

objective = 707
x(0) = 7
x(1) = 70

An optimization using PRINTEMPS consists of three steps: (1) modeling, (2) optimization, and (3) accessing the result. Each step is detailed in the following sections of this document.

Modeling

Basics

Let us consider the following optimization problem:

(P2):  minimize  f_1 + f_2
         x,y
      subject to x_1 + x_2 + x_3              >= 2,
                                   y_1 + y_2   = 1,
                 x_1             + y_1         = 1,
                 x_1, x_2, x_3, y_1, y_2 \in {0, 1},

        where    f_1 = 2 x_1 + 7 x_2 + 9 x_3 + 1,
                 f_2 = 5 y_1 + 6 y_2.

The following code shows an implementation of the mathematical model of the problem (P2) using PRINTEMPS.

printemps::model::IPModel model;

// Create binary decision variables.
auto& x_1 = model.create_variable("x_1", 0, 1);
auto& x_2 = model.create_variable("x_2", 0, 1);
auto& x_3 = model.create_variable("x_3", 0, 1);

auto& y_1 = model.create_variable("y_1", 0, 1);
auto& y_2 = model.create_variable("y_2", 0, 1);

// Create expressions.
auto& f_1 = model.create_expression("f_1");
auto& f_2 = model.create_expression("f_2");

// Create constraints.
auto& g_1 = model.create_constraint("g_1");
auto& g_2 = model.create_constraint("g_2");
auto& g_3 = model.create_constraint("g_3");

// Define expressions.
f_1 = 2 * x_1 + 7 * x_2 + 9 * x_3 + 1;
f_2 = 5 * y_1 + 6 * y_2;

// Define constraints.
g_1 = x_1 + x_2 + x_3 >= 2;
g_2 = y_1 + y_2 == 1;
g_3 = x_1 + y_1 == 1;

// Define the objective function to be minimized.
model.minimize(f_1 + f_2);

An optimization problem can be described using arithmetic operations on decision variables, expressions, constraints, and related objects. An expression can represent a linear function of decision variables, such as 2 x_1 + 7 x_2 + 9 x_3 + 1, and can be used to simplify the model. Of course, the statement model.minimize(f_1 + f_2); in the code above can be replaced with model.minimize(2 * x_1 + 7 * x_2 + 9 * x_3 + 5 * y_1 + 6 * y_2 + 1);. Please note that instances of decision variables, expressions, and constraints created by the methods create_variable(), create_expression(), and create_constraint() are managed by the IPModel object. These methods return references to the created instances, which should be assigned to references (aliases) representing the symbols of the problem.

PRINTEMPS supports one- or multi-dimensional indexed variables such as x(0), x(1), x(0,0), x(0,1), ... for intuitive modeling. The problem (P2) can be equivalently modeled as follows:

printemps::model::IPModel model;

// Create indexed binary decision variables x(0), x(1), x(2), y(0), y(1).
auto& x = model.create_variables("x", 3, 0, 1);
auto& y = model.create_variables("y", 2);

// Create indexed expressions f(0), f(1).
auto& f = model.create_expressions("f", 2);

// Create indexed constraints g(0), g(1), g(2).
auto& g = model.create_constraints("g", 3);

// Define expressions.
f(0) = 2 * x(0) + 7 * x(1) + 9 * x(2) + 1;
f(1) = 5 * y(0) + 6 * y(1);

// Define constraints.
g(0) = x(0) + x(1) + x(2) >= 2;
g(1) = y(0) + y(1) == 1;
g(2) = x(0) + y(0) == 1;

// Define the objective function to be minimized.
model.minimize(f(0) + f(1));

In addition, the following special methods are provided to simplify problem formulation:

With these methods, the problem (P2) can be simplified as follows:

printemps::model::IPModel model;

// Create indexed binary decision variables x(0), x(1), x(2), y(0), y(1).
auto& x = model.create_variables("x", 3, 0, 1);
auto& y = model.create_variables("y", 2);

// Create indexed expressions f(0), f(1).
auto& f = model.create_expressions("f", 2);

// Create indexed constraints g(0), g(1), g(2).
auto& g = model.create_constraints("g", 3);

// Create coefficient vectors.
std::vector<int> c = {2, 7, 9};
std::vector<int> d = {5, 6};

// Define expressions.
f(0) = x.dot(c) + 1;
f(1) = y.dot(d);

// Define constraints.
g(0) = x.sum() >= 2;
g(1) = y.selection();
g(2) = x(0) + y(0) == 1;

// Define the objective function to be minimized.
model.minimize(f.sum());

Decision Variable

Creating Decision Variables

Decision variables to describe the problem can be created as follows:

printemps::model::IPModel model;

// Create scalar (non-indexed) variable x.
auto& x = model.create_variable("x");

// Create scalar (non-indexed) variable y with the lower bound 0 and the upper
// bound 1.
auto& y = model.create_variable("y", 0, 1);

// Create one-dimensional indexed variables z(0), ..., z(9).
auto& z = model.create_variables("z", 10);

// Create two-dimensional indexed variables w(0, 0), ..., w(19, 29).
auto& w = model.create_variables("w", {20, 30});

The methods create_variable(name) and create_variable(name, lb, ub) create a scalar (non-indexed) decision variable named name. Optionally, lower and upper bounds can be set through the arguments lb and ub, respectively.

The methods create_variables(name, n) and create_variables(name, n, lb, ub) create n decision variables named name with one-dimensional index. Each decision variable can be accessed as z(0). Optionally, lower and upper bounds can be set through the arguments lb and ub, respectively, and they will be commonly set for all variables.

The methods create_variables("name", {n_1, ..., n_m}) and create_variables("name", {n_1, ..., n_m}, lb, ub) create n_1 * ... * n_m decision variables named name with m-dimensional index. Each variable can be accessed as w(0, 0). Optionally, lower and upper bounds can be set through the arguments lb and ub, respectively, and they will be commonly set for all variables.

Setting Lower and Upper Bounds for Decision Variables

The lower and upper bounds of a decision variable can be specified after creation of the decision variable by the method set_bound(lb, ub). In the following code, the lower and upper bounds of x(1) to x(9) are commonly set as -5 and 5, respectively, while those of x(0) are set as 0 and 20, respectively.

printemps::model::IPModel model;

// Create decision variables x(0), ..., x(9) with the lower bound -5 and the
// upper bound 5.
auto& x = model.create_variables("x", 10, -5, 5);

// The lower and upper bounds of x(0) are changed to 0 and 20, respectively.
x(0).set_bound(0, 20);

Setting Initial Values for Decision Variables

The default initial value of a decision variable is 0. Users can specify initial values using the = operator, and the current value of the decision variable can be accessed using the value() method. In the following code, the initial value for x(0) is set to 1, while those of x(1) to x(9) are still 0.

printemps::model::IPModel model;

// Create decision variables x(0), ..., x(9). 
// The initial values are set to the default value 0.
auto& x = model.create_variables("x", 10);

// The initial value of x(0) is changed to 1.
x(0) = 1;

// x_value will contain 1.
auto x_value = x(0).value();

Fixing Decision Variable Values

Users can also "fix" the values of decision variables by the method fix_by(value). The value of a fixed variable will not be changed during optimization. In the following code, the value for x(0) is fixed to 1, while those of x(1) to x(9) can vary.

printemps::model::IPModel model;

// Create decision variables x(0), ..., x(9). 
auto& x = model.create_variables("x", 10);

// The value of x(0) is fixed by 1.
x(0).fix_by(1);

There are some other methods related to fixing. The method unfix() cancels the fixing, and the method set_value_if_not_fixed(value) overrides the value of the decision variable if it is not fixed. If the decision variable is fixed, assignment using the = operator discussed above will throw an error. Users can check whether the decision variable is fixed or not by the method is_fixed(), which returns true if and only if the decision variable is fixed.

Expression

Creating Expression

Expressions to describe the problem can be created as follows:

printemps::model::IPModel model;

auto& x = model.create_variable("x");
auto& y = model.create_variable("y");

// Create scalar (non-indexed) expression p.
auto& p = model.create_expression("p");

// Create scalar (non-indexed) expression q by given expression 2 * x + 3 * y.
auto& q = model.create_expression("q", 2 * x + 3 * y);

// Create one-dimensional expressions r(0), ..., r(9).
auto& r = model.create_expressions("r", 10);

// Create two-dimensional expressions r(0,0), ..., r(19, 29).
auto& s = model.create_expressions("s", {20, 30});

The methods create_expression(name) and create_expression(name, expression) create a scalar (non-indexed) expression named name. Optionally, the created expression can be directly defined by an expression.

The methods create_expressions(name, n) creates n expressions named name with one-dimensional index. Each expression can be accessed as r(0).

The create_expressions("name", {n_1, ..., n_m}) method creates n_1 * ... * n_m expressions named name with m-dimensional index. Each expression can be accessed as s(0, 0).

Defining Expression

An expression can be defined as a linear function of decision variables and other expressions. The following code shows examples of expression definitions.

printemps::model::IPModel model;

// Create decision variables.
auto &x = model.create_variable("x");
auto &y = model.create_variable("y");
auto &z = model.create_variables("z", 2);

// Create and define expressions.

// OK. A first-order function on decision variables can be accepted as an
// expression.
auto &p = model.create_expression("p");
p       = 2 * x + 3 * y + 4;

// OK. A decision variable itself can be accepted as an expression.
auto &q = model.create_expression("q");
q       = x;

// OK. A numerical constant can be accepted as an expression.
auto &r = model.create_expression("r");
r       = 4;

// OK. An expression can include a term of a decision variable divided by a
// constant. An expression can be defined in its creation.
auto &s = model.create_expression("s", 2 * x + y / 3);

// OK. Addition and subtraction of the predefined expression are also possible.
s += 4 * x + 1; // s = 6 * x + y / 3 + 1 
s -= 5 * x;     // s = x + y / 3 + 1

auto &t = model.create_expressions("t", 2);
// OK. Parentheses can be used in the definition of an expression.
t(0) = x + 3 * (y * z(0));
// OK. A first-order function on expressions can be accepted as an expression.
t(1) = 5 * s(0) + 6 * s(1) + 6;

auto &u = model.create_expression("u");
// Invalid. Multiplication and division that produce nonlinear terms are not
// supported. 
// u = x * x + x / y;

// Invalid. Nonlinear functions are not supported.
// u = sin(x);

Constraint

Creating Constraint

Constraints to describe the problem can be created as follows:

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variable("x");
auto& y = model.create_variable("y");

// Create scalar (no-indexed) constraint g.
auto& g = model.create_constraint("g");

// Create scalar (no-indexed) constraint h by given constraint x + y = 1.
auto& h = model.create_constraint("h", x + y == 1);

// Create one-dimensional constraints u(0), ..., u(9).
auto& u = model.create_constraints("u", 10);

// Create two-dimensional constraints v(0,0), ..., v(19, 29).
auto& v = model.create_constraints("v", {20, 30});

The methods create_constraint(name) and create_constraint(name, constraint) create scalar (no-indexed) constraint named name. As optional, the created constraint can be directly defined by constraint.

The methods create_constraints(name, n) creates n constraints named name with one-dimensional index. Each constraint can be accessed by like as u(0).

The create_constraints("name", {n_1, ..., n_m}) method creates n_1 * ... * n_m constraints named name with m-dimensional index. Each constraint can be accessed by like as v(0, 0).

Defining Constraint

A constraint can be defined as an equal or inequality between decision variables, expressions, and constants (numerical values of primitive types such as int, double) with operators ==, <=, and >=.
Note that a constraint whose sides are both constants is not acceptable, meaning that at least one side must contain a decision variable or an expression. The following code shows examples of constraint definitions.

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variable("x");
auto& y = model.create_variable("y");
auto& z = model.create_variables("z", {10, 10});

// Create and define an expression.
auto& p = model.create_expression("p", x + y);

// Create and define constraints.
auto& g = model.create_constraint("g");
g       = (2 * x + 3 * y + 5 <= 10); // OK

auto& h = model.create_constraint("h", p >= 7 + 2); // OK

auto& u = model.create_constraint("u");
// Invalid. Both sides are constants.
// u = 4 <= 5;

auto& v = model.create_constraints("v", 10);
for (auto i = 0; i < 10; i++) {
    // OK. This code is identical to the following code:
    // u(i) = z.sum({i, printemps::model_component::All}) == 1;
    v(i) = z.selection({i, printemps::model_component::All});
}

Objective Function

The objective function and the optimization sense (minimization or maximization) can be specified by the methods minimize() and maximize() of IPModel object. An expression or a decision variable can be objective function by specifying them as the argument of the methods minimize() or maximize(). If there is no objective function definition, PRINTEMPS solves the problem as a constraint satisfaction problem, which aims to find a feasible solution. If the minimize() or maximize() methods are called more than once, only the last definition will take effect. The following code shows an example of objective function definitions.

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variable("x");
auto& y = model.create_variable("y");

// Create and define an expression.
auto& p = model.create_expression("p", 2 * x + 3 * y + 4);

// The part of constraint creation and definition is omitted.

// OK. A decision variable can be the objective function.
model.minimize(x);
// OK. An expression can be the objective function.
model.maximize(p);
// OK. A formula to define an expression can also be the objective function.
model.minimize(5 * x - 6 * y + 7);

// Invalid. The arguments of minimize() and maximize() must not be a constant.
// model.maximize(5);

Special Methods for Practical Modeling

Sum of Indexed Components

The method sum() returns an expression defined as the sum of all components of the indexed decision variables or expressions. The code

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variables("x", {10, 20});

// Create and define an expression.
auto& p = model.create_expression("p", x.sum());

is identical to the following code:

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variables("x", {10, 20});

// Create and define an expression.
auto& p = model.create_expression("p");

for (auto i = 0; i < 10; i++) {
    for (auto j = 0; j < 20; j++) {
        p += x(i, j);
    }
}

Partial sums of components of indexed decision variables or expressions can be obtained by specifying which dimensions are fixed and which are summed over. Summation is calculated along the dimension where the index is specified by the constant printemps::model_component::All. The code

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variables("x", {10, 20});

// Create and define expressions.
auto& p = model.create_expressions("p", 10);
for (auto i = 0; i < 10; i++) {
    p(i) += x.sum({i, printemps::model_component::All});
}

is identical to the following code:

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variables("x", {10, 20});

// Create and define expressions.
auto& p = model.create_expressions("p", 10);
for (auto i = 0; i < 10; i++) {
    for (auto j = 0; j < 20; j++) {
        p(i) += x(i, j);
    }
}

Inner Product

The method dot() returns an expression representing the inner product of the indexed variables or expressions and a numerical array with the same number of components. The code

std::vector<int> a = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variables("x", 10);

// Create and define expressions.
auto& p = model.create_expression("p", x.dot(a));

is identical to the following code:

std::vector<int> a = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variables("x", 10);

// Create and define expressions.
auto& p = model.create_expression("p");
for (auto i = 0; i < 10; i++) {
    p += a[i] * x(i);
}

The method dot() is also applicable to two- or higher-dimensional indexed decision variables and expressions by specifying fixed indices and the index to be summed over. Note that only one index can be summed over (using Range::All) at a time. The following code shows example of allowable and non-allowable usages.

std::vector<int> a = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variables("x", {10, 10});

// Create and define expressions.

auto& p = model.create_expressions("p", 10);
// OK. The following part is identical to
// for (auto i=0; i<10; i++){
//    for (auto j=0; j<10; j++){
//        p(i) += a[j] * x(i, j);
//    }
// }
for (auto i = 0; i < 10; i++) {
    p(i) = x.dot(a, {i, printemps::model_component::All});
}

auto& q = model.create_expressions("q", 10);
// OK. The following part is identical to
// for (auto i=0; i<10; i++){
//    for (auto j=0; j<10; j++){
//        q(i) += a[j] * x(j, i);
//    }
// }
for (auto i = 0; i < 10; i++) {
    q(i) = x.dot(a, {printemps::model_component::All, i});
}

auto& r = model.create_expressions("r", 10);
// Invalid. The number of indices for summation must be one.
// for(auto i=0; i<10; i++){
//     r[i] = x.dot(a, {printemps::model_component::All, printemps::model_component::All});
// }

Selection Constraint

The method selection() returns a constraint forcing exactly one of the indexed binary decision variables to take the value 1. Such constraints are frequently used in practice. The code

printemps::model::IPModel model;

// Create binary decision variables.
auto& x = model.create_variables("x", {10, 20}, 0, 1);

// Create and define a constraint.
auto& g = model.create_constraint("g", x.selection());

is identical to the following code:

printemps::model::IPModel model;

// Create binary decision variables.
auto& x = model.create_variables("x", {10, 20}, 0, 1);

// Create and define an expression for temporary.
auto& p = model.create_expression("p");
for (auto i = 0; i < 10; i++) {
    for (auto j = 0; j < 20; j++) {
        p += x(i, j);
    }
}
// Create and define a constraint.
auto& g = model.create_constraint("g");
g = (p == 1);

Of course, a selection constraint for partial sums of indexed binary decision variables can be obtained by specifying fixed indices and the index to be summed over. The code

printemps::model::IPModel model;

auto& x = model.create_variables("x", {10, 20}, 0, 1);
auto& g = model.create_constraints("g", 10);

for (auto i = 0; i < 10; i++) {
    g(i) = x.selection({i, printemps::model_component::All});
}

is identical to the following code:

printemps::model::IPModel model;

auto& x = model.create_variables("x", {10, 20}, 0, 1);
auto& g = model.create_constraints("g", 10);

// Create and define an expression for temporary.
auto& p = model.create_expressions("p", 10);
for (auto i = 0; i < 10; i++) {
    for (auto j = 0; j < 20; j++) {
        p(i) += x(i, j);
    }
}

// Create and define constraints.
for (auto i = 0; i < 10; i++) {
    g(i) = (p(i) == 1);
}

Optimization

The solver can be run by the printemps::solver::solve() function with the pointer to the IPModel object as the argument:

printemps::model::IPModel model;

// The modeling part is omitted.

auto result = printemps::solver::solve(&model);

Users can configure the behavior of the solver by setting options. The options include those related to optimization parameters, log verbosity, multi-thread computation enabling, etc. Solver Option Guide gives a detailed description of the options and their default values.

Accessing the Result

The result of the optimization is stored in the solution member of the result object, which is returned by the solve() function. If one or more feasible solutions were found during the optimization, the best solution among them is selected as the final best solution. Otherwise, the solution with the minimum augmented objective function value is selected. The augmented objective function is the sum of the original objective function value and the penalty values of constraint violations. The following code shows an example of accessing the optimization result.

printemps::model::IPModel model;

// Create decision variables.
auto& x = model.create_variable("x");
auto& y = model.create_variables("y", 1000);

// Create expressions.
auto& p = model.create_expression("p");
auto& q = model.create_expressions("q", {5, 5});

// Create constraints.
auto& g = model.create_constraint("g");
auto& h = model.create_constraints("h", {10, 10, 10});

// The formulation part is omitted.

auto result = printemps::solver::solve(&model);

// Print whether a feasible solution was obtained or not and the objective
// function value of the final best solution.
std::cout << result.solution.is_feasible() << std::endl;
std::cout << result.solution.objective() << std::endl;

// Print the value of the decision variables "x" and "y(0)" in the final
// best solution.
std::cout << result.solution.variables("x").value() << std::endl;
std::cout << result.solution.variables("y").values(0) << std::endl;

// Print the value of the expressions "p" and "q(0, 0) "in final best solution.
std::cout << result.solution.expressions("p").value() << std::endl;
std::cout << result.solution.expressions("q").values(0, 0) << std::endl;

// Print the value of the constraints "g" and "h(0, 0, 0)" and the violations
// for them in final best solution.
std::cout << result.solution.constraints("g").value() << std::endl;
std::cout << result.solution.constraints("h").values(0, 0, 0) << std::endl;
std::cout << result.solution.violations("g").value() << std::endl;
std::cout << result.solution.violations("h").values(0, 0, 0) << std::endl;

For accessing the result, the following methods are provided by the member solution in the Result object.

To make it easier to inspect the values, the utility functions print_variable_values(), print_expression_values(), print_constraint_values(), and print_violation_values() are provided. They respectively print the values of all variables, expressions, constraints, and violations to the standard output. Please refer to the example.

References

Python Bindings

Official Python bindings built with pybind11 are supported from v2.9.0.

Installation

You can install the package from the repository root:

pip install ./python

Quick Start

The Python API mirrors the C++ API closely. Here is a simple example solving the same problem (P):

import printemps

# (1) Modeling
model = printemps.IPModel()
x = model.create_variables("x", 2, -100, 100)
g = model.create_constraints("g", 2)

g[0] = 66 * x(0) + 14 * x(1) >= 1430
g[1] = -82 * x(0) + 28 * x(1) >= 1306
model.minimize(x(0) + 10 * x(1))

# (2) Running Solver
result = printemps.solve(model)

# (3) Accessing the Result
print("objective =", result.solution.objective())
print("x(0) =", result.solution.variables("x")(0))
print("x(1) =", result.solution.variables("x")(1))

Key Differences from C++

The main syntax differences are summarised as follows:

C++ idiom Python idiom
auto& x = model.create_variables(...) x = model.create_variables(...)
g(0) = expr >= 1430 g[0] = expr >= 1430
f(0) = expr f[0] = expr
objective += term objective += term
expression_block(i, j, k) += x(n, m, k) expression_block[i, j, k] += x(n, m, k)
x.selection({n, Range::All, k}) x.selection([n, printemps.ALL, k])
option.output.verbose = verbose::Full option.output.verbose = printemps.Verbose.Full
result.solution.variables("x").values(n,m) result.solution.variables("x")(n, m)

For details, please refer to the Python binding README.md in the repository.