A.31 library(simplex): Solve linear programming problems

Author: Markus Triska

A linear programming problem consists of a set of (linear) constraints, a number of variables and a linear objective function. The goal is to assign values to the variables so as to maximize (or minimize) the value of the objective function while satisfying all constraints.

Many optimization problems can be modeled in this way. Consider having a knapsack with fixed capacity C, and a number of items with sizes s(i) and values v(i). The goal is to put as many items as possible in the knapsack (not exceeding its capacity) while maximizing the sum of their values.

As another example, suppose you are given a set of coins with certain values, and you are to find the minimum number of coins such that their values sum up to a fixed amount. Instances of these problems are solved below.

The library(simplex) module provides the following predicates:

assignment(+Cost, -Assignment)

Solves an assignment problem. Cost is a list of lists representing the quadratic cost matrix, where element (i,j) denotes the integer cost of assigning entity i to entity j. An assignment with minimal cost is computed and unified with Assignment as a list of lists, representing an adjacency matrix.

constraint(+Constraint, +S0, -S)

Adds a linear or integrality constraint to the linear program corresponding to state S0. A linear constraint is of the form "Left Op C", where "Left" is a list of Coefficient*Variable terms (variables in the context of linear programs can be atoms or compound terms) and C is a non-negative numeric constant. The list represents the sum of its elements. Op can be =, =< or >=. The coefficient "1" can be omitted. An integrality constraint is of the form integral(Variable) and constrains Variable to an integral value.

constraint(+Name, +Constraint, +S0, -S)

Like constraint/3, and attaches the name Name (an atom or compound term) to the new constraint.

constraint_add(+Name, +Left, +S0, -S)

Left is a list of Coefficient*Variable terms. The terms are added to the left-hand side of the constraint named Name. S is unified with the resulting state.

gen_state(-State)

Generates an initial state corresponding to an empty linear program.

maximize(+Objective, +S0, -S)

Maximizes the objective function, stated as a list of "Coefficient*Variable" terms that represents the sum of its elements, with respect to the linear program corresponding to state S0. S is unified with an internal representation of the solved instance.

minimize(+Objective, +S0, -S)

Analogous to maximize/3.

objective(+State, -Objective)

Unifies Objective with the result of the objective function at the obtained extremum. State must correspond to a solved instance.

shadow_price(+State, +Name, -Value)

Unifies Value with the shadow price corresponding to the linear constraint whose name is Name. State must correspond to a solved instance.

transportation(+Supplies, +Demands, +Costs, -Transport)

Solves a transportation problem. Supplies and Demands must be lists of non-negative integers. Their respective sums must be equal. Costs is a list of lists representing the cost matrix, where an entry (i,j) denotes the integer cost of transporting one unit from i to j. A transportation plan having minimum cost is computed and unified with Transport in the form of a list of lists that represents the transportation matrix, where element (i,j) denotes how many units to ship from i to j.

variable_value(+State, +Variable, -Value)

Value is unified with the value obtained for Variable. State must correspond to a solved instance.

A.31.1 Example 1

This is the "radiation therapy" example, taken from "Introduction to Operations Research" by Hillier and Lieberman. DCG notation is used to implicitly thread the state through posting the constraints:

:- use_module(library(simplex)).

post_constraints -->
        constraint([0.3*x1, 0.1*x2] =< 2.7),
        constraint([0.5*x1, 0.5*x2] = 6),
        constraint([0.6*x1, 0.4*x2] >= 6),
        constraint([x1] >= 0),
        constraint([x2] >= 0).

radiation(S) :-
        gen_state(S0),
        post_constraints(S0, S1),
        minimize([0.4*x1, 0.5*x2], S1, S).

An example query:

?- radiation(S), variable_value(S, x1, Val1),
                 variable_value(S, x2, Val2).

Val1 = 15 rdiv 2
Val2 = 9 rdiv 2 ;

A.31.2 Example 2

Here is an instance of the knapsack problem described above, where C = 8, and we have two types of items: One item with value 7 and size 6, and 2 items each having size 4 and value 4. We introduce two variables, x(1) and x(2) that denote how many items to take of each type.

knapsack_constrain(S) :-
        gen_state(S0),
        constraint([6*x(1), 4*x(2)] =< 8, S0, S1),
        constraint([x(1)] =< 1, S1, S2),
        constraint([x(2)] =< 2, S2, S).

knapsack(S) :-
        knapsack_constrain(S0),
        maximize([7*x(1), 4*x(2)], S0, S).

An example query yields:

?- knapsack(S), variable_value(S, x(1), X1),
                variable_value(S, x(2), X2).

X1 = 1
X2 = 1 rdiv 2 ;

That is, we are to take the one item of the first type, and half of one of the items of the other type to maximize the total value of items in the knapsack.

If items can not be split, integrality constraints have to be imposed:

knapsack_integral(S) :-
        knapsack_constrain(S0),
        constraint(integral(x(1)), S0, S1),
        constraint(integral(x(2)), S1, S2),
        maximize([7*x(1), 4*x(2)], S2, S).

Now the result is different:

?- knapsack_integral(S), variable_value(S, x(1), X1),
                         variable_value(S, x(2), X2).

X1 = 0
X2 = 2

That is, we are to take only the two items of the second type. Notice in particular that always choosing the remaining item with best performance (ratio of value to size) that still fits in the knapsack does not necessarily yield an optimal solution in the presence of integrality constraints.

A.31.3 Example 3

We are given 3 coins each worth 1, 20 coins each worth 5, and 10 coins each worth 20 units of money. The task is to find a minimal number of these coins that amount to 111 units of money. We introduce variables c(1), c(5) and c(20) denoting how many coins to take of the respective type:

coins -->
        constraint([c(1), 5*c(5), 20*c(20)] = 111),
        constraint([c(1)] =< 3),
        constraint([c(5)] =< 20),
        constraint([c(20)] =< 10),
        constraint([c(1)] >= 0),
        constraint([c(5)] >= 0),
        constraint([c(20)] >= 0),
        constraint(integral(c(1))),
        constraint(integral(c(5))),
        constraint(integral(c(20))),
        minimize([c(1), c(5), c(20)]).

coins(S) :-
        gen_state(S0),
        coins(S0, S).

An example query:

?- coins(S), variable_value(S, c(1), C1),
             variable_value(S, c(5), C5),
             variable_value(S, c(20), C20).

C1 = 1
C5 = 2
C20 = 5