Different results when using library(simplex) or R-package lpSolve / lpSolveAPI

Hi,

One of the exercises from the course material for the Linear Programming course, involved baking cakes. For chocolate, vanilla, and lemon cakes it is known how many eggs and how much flour are needed and what is the baking time. The question is to maximize the profit.

This is the model in lpSolveAPI:

library (lpSolveAPI)

result <- make.lp(3,3)
set.column(result, 1, c(4, 4, 20))
set.column(result, 2, c(2, 5, 40))
set.column(result, 3, c(6, 6, 50))
set.objfn(result, c(3, 5, 7))
set.constr.type(result, rep("<=", 3))
set.rhs(result, c(32, 40, 260))
lp.control(result, sense = "max")
set.type(result, 1, "integer")
set.type(result, 2, "integer")
set.type(result, 3, "integer")
dimnames(result) <- list(c("Ei", "Bloem", "Tijd"), c("Chocola", "Vanille", "Citroen"))

4 Chocolate, 2 vanilla and 2 lemmon cakes give a profit of 36

I have defined the same problem with library(simplex):

:- use_module(library(simplex)).

cakes_3(S) :-
        gen_state(S0),
        post_constraints_v3(S0, S1),
        maximize([3*chocola, 5*vanille, 7*citroen], S1, S).

post_constraints_v3 -->
        constraint([4*chocola, 2*vanille, 6*citroen] =< 32),
        constraint([4*chocola, 5*vanille, 6*citroen] =< 40),
        constraint([20*chocola, 40*vanille, 50*citroen] =< 260),
        constraint(integral(chocola)),
        constraint(integral(vanille)),
        constraint(integral(citroen)),
        constraint([chocola] >= 0),
        constraint([vanille] >= 0),
        constraint([citroen] >= 0).


test_cakes_v3 :- cakes_3(S),
        variable_value(S, chocola, Val1),
        variable_value(S, vanille, Val2),
        variable_value(S, citroen, Val3),
        writeln(Val1),
        writeln(Val2),
        writeln(Val3), nl,
        result(Val1, Val2, Val3), nl.

result(A, B, C) :-
        Ei    is A*4 +  B*2 +  C*6, writeln(Ei),
        Bloem is A*4 +  B*5 +  C*6, writeln(Bloem),
        Tijd  is A*20 + B*40 + C*50, writeln(Tijd),
        Winst is A*3 +  B*5 +  C*7, writeln(Winst).

Now 1 Chocolate, 1 vanilla and 4 lemmon cakes give a profit of 36

Both solutions are correct.

My question is why both environments produce diffent results? And since both solutions are correct, is it possible to have simplex return all the solutions?

Ben

They get different results because there are different methods to find the optimum. When an optimum is found, the algorithm stops.

In most cases, there is only one optimum, but in some cases there can be many, or when the integral constraint is gone, infinite numbers. For 3 variable problems, it’s reasonable to enumerate all the solutions, but for large integer problems, each possible solution could be very slow to find.

FWIW, when using pack(clpBNR) this is done in two steps for this reason. First find the global maximum, then the set of “maximizers” corresponding to that maximum. In general, the second step is non-deterministic.

?- [Choc,Van,Cit]::integer(0,_),
{4*Choc+2*Van+6*Cit=<32, 4*Choc+5*Van+6*Cit=<40, 20*Choc+40*Van+50*Cit=<260},
global_maximum(3*Choc+5*Van+7*Cit, Profit).
Profit = 36,
Choc::integer(0, 8),
Van::integer(0, 6),
Cit::integer(0, 5).

?- [Choc,Van,Cit]::integer(0,_),
{4*Choc+2*Van+6*Cit=<32, 4*Choc+5*Van+6*Cit=<40, 20*Choc+40*Van+50*Cit=<260},
global_maximum(3*Choc+5*Van+7*Cit, Profit),
enumerate([Choc,Van,Cit]).
Choc = Van, Van = 1,
Cit = 4,
Profit = 36 ;
Choc = 4,
Van = Cit, Cit = 2,
Profit = 36 ;
false.

Finding the global maximum partially narrows the domains, but a subsequent labelling step is required to separate the solutions.

2 Likes

Thaks for the hint. I have installed the package and start experimenting with it.

If the focus of the course is linear programming, just note that linear programming techniques are not used in in this clpBNR solution.

IMO, it’s a bit unfortunate that the simplex library chooses to bind (narrow) the “maximizer” values when determining the maximum value for the objective function. If it just treated it as another constraint to be applied, standard labelling could generate all the solutions.

Some further thoughts:

The clpBNR approach is general enough to handle global optimization of linear and non-linear problems, but it does so via a branch-and-bound search. This is never going to be as efficient as an algorithmic approach that generates a solution, such as that implemented by library(simplex). For many problems it doesn’t matter, but often the search space is large enough that efficiency becomes an issue.

For that reason, it’s advantageous to use known algorithms to generate new constraints for your favourite CLP which then can help reduce the search space for the labelling step. Looking at library(simplex) closer, I note that it really doesn’t bind any (constrained) variables until variable_value/3 is called. So an interesting approach is to use library(simplex) to generate a maximum/minimum value which can be used to define a new constraint for the CLP, and then use labelling to generate all solutions. (If you only want one solution, just use variable_value.)

This would be fairly straight forward if the CLP and library(simplex) shared a common representation of constraints and the objective expression to be optimized. But it doesn’t take a lot of Prolog code to implement wrapper functions to provide a “friendlier” API. I’ve done a small prototype for clpBNR which implements lin_minimum(ObjF,Constraints,MinValue) and lin_maximum(ObjF,Constraints,MaxValue) where ObjF and Constraints are an expression and list of constraints as defined by the particular CLP syntax and Min/MaxValue is unified with the requested optimal value. So the cakes example (clpBNR version) would look something like:

?- [Choc,Van,Cit]::integer(0,_),
Profit = 3*Choc+5*Van+7*Cit,
Constraints = [4*Choc+2*Van+6*Cit=<32, 4*Choc+5*Van+6*Cit=<40, 20*Choc+40*Van+50*Cit=<260],
lin_maximum(Profit,Constraints,Max),
{Profit==Max,Constraints}, enumerate([Choc,Van,Cit]),nl.

Choc = Van, Van = 1,
Cit = 4,
Profit = 3*1+5*1+7*4,
Constraints = [4*1+2*1+6*4=<32, 4*1+5*1+6*4=<40, 20*1+40*1+50*4=<260],
Max = 36 ;

Choc = 4,
Van = Cit, Cit = 2,
Profit = 3*4+5*2+7*2,
Constraints = [4*4+2*2+6*2=<32, 4*4+5*2+6*2=<40, 20*4+40*2+50*2=<260],
Max = 36 ;
false.

I don’t think it would take much to change the code to handle clpfd syntax; mainly an issue of using different arithmetic operators. Here’s the unpolished source - sorry for the lack of comments:

Linear optimization code
:- use_module(library(simplex)).
:- use_module(library(clpBNR),[domain/2]).

lin_minimum(ObjF,Constraints,MinValue) :-
	init_simplex(ObjF,Constraints,Objective,S0),
	minimize(Objective,S0,S),
	objective(S,MinValue).
	
lin_maximum(ObjF,Constraints,MaxValue) :-
	init_simplex(ObjF,Constraints,Objective,S0),
	maximize(Objective,S0,S),
	objective(S,MaxValue).

init_simplex(ObjF,Constraints,Objective,S) :-
	gen_state(S0),
	sim_constraints(Constraints,S0,S1),
	term_variables((ObjF,Constraints),Vs),
	constrain_ints(Vs,S1,S),
	map_simplex(ObjF,T/T,Objective/[]).
	
constrain_ints([],S,S).
constrain_ints([V|Vs],Sin,Sout) :-
	(domain(V,integer(_,_))
	 -> simplex_var(V,SV), constraint(integral(SV),Sin,Snxt)
	 ;  Snxt = Sin
	),
	constrain_ints(Vs,Snxt,Sout).

sim_constraints([],S,S).
sim_constraints([C|Cs],Sin,Sout) :-
	sim_constraint(C,SC),
	constraint(SC,Sin,Snxt),
	sim_constraints(Cs,Snxt,Sout).

sim_constraint(C,SC) :-
	C=..[Op,LHS,RHS],
	constraint_op(Op,COp),
	map_simplex(LHS,T/T,Sim/[]),
	SC=..[COp,Sim,RHS].

map_simplex(T,CT/[S|Tail],CT/Tail) :-
	map_simplex_term(T,S),
	!.
map_simplex(A+B,Cin,Cout) :- !,
	map_simplex(A,Cin,Cnxt),
	map_simplex(B,Cnxt,Cout).
map_simplex(A-B,Cin,Cout) :- !,
	map_simplex(A,Cin,Cnxt),
	map_simplex(-B,Cnxt,Cout).
			
map_simplex_term(V,1*SV) :- var(V), simplex_var(V,SV), !.
map_simplex_term(-T,NN*V) :- !,
	map_simplex_term(T,N*V),
	NN is -N.
map_simplex_term(N*V,N*SV) :- number(N), var(V), simplex_var(V,SV), !.
map_simplex_term(V*N,N*SV) :- number(N), var(V), simplex_var(V,SV).

constraint_op(==,=).
constraint_op(=<,=<).
constraint_op(>=,>=).

simplex_var(V,var(Vstr)) :- term_string(V,Vstr).

So this hybrid approach seems to nicely combine the use of a known algorithm to produce the optimal value with the ability of a CLP to generate all solutions corresponding to that value.

2 Likes

FWIW, an up-to-date version of clpBNR now includes a “toolkit” module which exports CLP friendly wrapper predicates for library(simplex) including lin_maximum and lin_maximize (the latter deterministically provides a single set of maximizing values):

?- use_module(library(clpBNR_toolkit)).

?- [Choc,Van,Cit]::integer(0,_),
Profit = 3*Choc+5*Van+7*Cit,
Constraints = {4*Choc+2*Van+6*Cit=<32, 4*Choc+5*Van+6*Cit=<40, 20*Choc+40*Van+50*Cit=<260},
lin_maximum(Profit,Constraints,Max),enumerate([Choc,Van,Cit]),nl.

Choc = Van, Van = 1,
Cit = 4,
Profit = 3*1+5*1+7*4,
Constraints = {4*1+2*1+6*4=<32, 4*1+5*1+6*4=<40, 20*1+40*1+50*4=<260},
Max = 36 ;

Choc = 4,
Van = Cit, Cit = 2,
Profit = 3*4+5*2+7*2,
Constraints = {4*4+2*2+6*2=<32, 4*4+5*2+6*2=<40, 20*4+40*2+50*2=<260},
Max = 36 ;
false.

?- [Choc,Van,Cit]::integer(0,_),
Profit = 3*Choc+5*Van+7*Cit,                                                                                             Constraints = {4*Choc+2*Van+6*Cit=<32, 4*Choc+5*Van+6*Cit=<40, 20*Choc+40*Van+50*Cit=<260}, 
lin_maximize(Profit,Constraints,Max).
Choc = Van, Van = 1,
Cit = 4,
Profit = 3*1+5*1+7*4,
Constraints = {4*1+2*1+6*4=<32, 4*1+5*1+6*4=<40, 20*1+40*1+50*4=<260},
Max = 36.

lin_minimize/3,and lin_maximize/3 can be quite helpful for linear programming problems when the general clpBNR global optimization of non-linear systems is infeasible due to the extensive search space. (E.g., 77 real variables for The Stigler Diet Problem  |  OR-Tools  |  Google Developers. Discovering all solutions for such a large search space is still an issue.)

2 Likes