# 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