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

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