%------------------------------------------------------------------ % embedding branch and bound into Prolog: 0-1 knapsack problem. %------------------------------------------------------------------ % knapsack0.pl % By Kenryo Indo (Kanto Gakuen University) % 23,25 Aug 2003. % modified: 26-31 Aug 2003. %--------------------------- % Reference: % [1] M. Sakawa (2000). Risan Sistemu no Saitekika. % (Optimization of Discrete Systems.) % Morikita Shuppan.(Japanese) % main: solve_knapsack/4 % (with auto save data to file by tell_data/0.) % output data: solved_data/4, fathomed_data/3. % an integer (not 0-1) programming problem. ilp1(max,_Z=2*X+5*Y,[X,Y], [ 2*X+6*Y =< 27, 8*X+6*Y =< 45, 3*X+1*Y =< 15 ] ). % -- a naive enumerating method to solve ILP problem. solve_ilp(max, Z=C1*X+C2*Y,[X,Y]):- ilp1(max, Z=C1*X+C2*Y,[X,Y],CON), boundary_of_feasible_solution(CON,[BX,BY]), max(Z,feasible_solution([X,Y],CON,[BX,BY],Z=C1*X+C2*Y)). boundary_of_feasible_solution(CON,[BX,BY]):- findall((BX,BY), ( member((A1*_X+A2*_Y =< B),CON), BX is integer(B/A1), BY is integer(B/A2) ), BS), member((BX,_),BS),\+ (member((BX1,_),BS),BX1 < BX), member((_,BY),BS),\+ (member((_,BY1),BS),BY1 < BY). feasible_solution_1([X,Y],CON,[BX,BY]):- length(LX,BX), length(LY,BY), nth0(X,[BX|LX],_), nth0(Y,[BY|LY],_), forall(member(C,CON),C). feasible_solution([X,Y],CON,[BX,BY],Z=C1*X+C2*Y):- ilp1(max, Z=C1*X+C2*Y,[X,Y],CON), feasible_solution_1([X,Y],CON,[BX,BY]), Z is C1*X+C2*Y. % utilities (mathematics) %-------------------------------------------------------- min(X,Goal):- max(-X,Goal). max(X,Goal):- % X: the objective variable, % Goal: the objective function and constraints, setof((X,Goal),Goal,Z), member((X,Goal),Z), \+ ( member((Y,_),Z), Y > X ). product_sum([],[],[],0). product_sum([P|Q],[A|B],[E|F],V):- length(Q,N), length(B,N), product_sum(Q,B,F,V1), E is P * A, V is V1 + E. %-------------------------------------------------------- % 4 variables 0-1 knapsack problem %-------------------------------------------------------- % -- a naive enumerating method to solve 0-1 knapsack problem. solve_knapsack( min(Z = -4 * X1 -9 * X2 -11 * X3 -3 * X4), st(3 * X1 + 8 * X2 + 5 * X3 + 2 * X4 =< 9), [X1,X2,X3,X4] ):- min(Z, knapsack_0( min(Z = -4 * X1 -9 * X2 -11 * X3 -3 * X4), st(3 * X1 + 8 * X2 + 5 * X3 + 2 * X4 =< 9), [X1,X2,X3,X4] ) ). knapsack_0(min(Z = OBJECTIVE), st(CONSTRAINT), [X1,X2,X3,X4]):- member(X1,[0,1]), member(X2,[0,1]), member(X3,[0,1]), member(X4,[0,1]), CONSTRAINT = (3 * X1 + 8 * X2 + 5 * X3 + 2 * X4 =< 9), CONSTRAINT, OBJECTIVE = -4 * X1 -9 * X2 -11 * X3 -3 * X4, Z is OBJECTIVE. knapsack_1( min(_Z = [-4, -9, -11, -3] * [X1,X2,X3,X4] ), st( [3, 8, 5, 2] * [X1,X2,X3,X4] =< 9 ) ). %-------------------------------------------------------- % solver for the continuous relaxed 0-1 knapsack problems %-------------------------------------------------------- solve_continuous_relaxed(P,min(Z=C*X),st(A*X= F = feasible ; F = infeasible ). % optimality for sorted continuos relaxed 0-1 knapsack problem. %-------------------------------------------------------- % See Theorem 3.1 of Ref.1, p.27. solve_continuous_relaxed_knapsack_1(0,[],_B,[],0). solve_continuous_relaxed_knapsack_1(P,[AJ|A],B,[XJ|X],S):- solve_continuous_relaxed_knapsack_1(P0,A,B,X,S0), S is S0 + AJ, pattern_of_solve_relaxed_knapsack(B,AJ,S0,S,P0,P,XJ). pattern_of_solve_relaxed_knapsack(B,_,S0,_S,P0,P,0):- B < S0, P = P0. pattern_of_solve_relaxed_knapsack(B,_,_S0,S,P0,P,1):- B > S, P is P0 + 1. pattern_of_solve_relaxed_knapsack(B,AJ,S0,S,P0,P,XJ):- B >= S0, B < S, XJ=(B-S0)/AJ, P is P0 + 1. % sorted problem by the efficiency. %-------------------------------------------------------- efficiency((PK,K,CK,AK,XK),(C0,A0,X0)):- nth1(K,C0,CK), nth1(K,A0,AK), nth1(K,X0,XK), PK is CK/AK. knapsack(min(Z=C*X),st(A*X= P is P0; P is P0 + 1). %-------------------------------------------------------- % branch-and-bound method for 0-1 knapsack problem %-------------------------------------------------------- % cf. serially satisfying plan of 0-1 constraints. %----------------------------------------------- % It may be seen as the abbreviated search or prior planning % which will be extended to the branch and bound algorithm % next in the program. :- dynamic plan_data/3. plan_0(N,X,PLAN):- abolish(plan_data/3), assert(plan_data(current(0),0,[])), length(X,N), plan_constraint_satisfaction_0(N,X,PLAN). plan_constraint_satisfaction_0(0,_,[]):- update_plan_0(0,[]). plan_constraint_satisfaction_0(K,X,[(P,XP)|PLAN]):- length(X,N), length(PLAN,K0), (K0 >= N -> !, fail;true), length([(P,XP)|PLAN],K), plan_constraint_satisfaction_0(K0,X,PLAN), K is K0 + 1, set_nonrecurrent_value_0(P,X,XP,PLAN), update_plan_0(K,[(P,XP)|PLAN]). set_nonrecurrent_value_0(P,X,XP,PLAN):- member(XP,[0,1]), nth1(P,X,XP), \+ member((P,_),PLAN). update_plan_0(K,PLAN):- retract(plan_data(current(L0),D1,D2)), assert(plan_data(log(L0),D1,D2)), L is L0 + 1, assert(plan_data(current(L),K,PLAN)). % solver of constrained continuous relaxed subproblems %----------------------------------------------- % a branch and bound algorithm (Ref.1, pp.26-29) solve_knapsack(K,PROBLEM,SOLVED,CONSTRAINTS):- PROBLEM=knapsack(min(_Z=_C*X),st(_A*X=<_B)), PROBLEM, (var(K)->length(X,K);integer(K)), initialize_experimental_data(PROBLEM,X), plan_constraint_satisfaction(K,PROBLEM,[SOLVED|_],CONSTRAINTS). solve_knapsack(_,_,_,_):- tell_data. initialize_experimental_data(PROBLEM,X):- abolish(solved_data/4), abolish(fathomed_data/3), assert(solved_data(current(0),PROBLEM,init,[])), assert(fathomed_data(current(0),[_,infinite,X],_)). % PLAN = Constraints. plan_constraint_satisfaction(0,PROBLEM,[SOLVED],[]):- %PROBLEM=knapsack(min(Z=C*X),st(A*X== N -> !, fail;true), length([_|PLAN],K), plan_constraint_satisfaction(K0,PROBLEM,SEARCHED,PLAN), %K is K0 + 1, non_integer_element(P,SEARCHED), generate_a_new_constraint((P,XP),X1,PLAN), solve_and_update(Q,PROBLEM1,_SUBPROBLEM,F,[(P,XP)|PLAN],_CASE), SOLVED=(PROBLEM1,F,Q). copy_frame_of_problem(PROBLEM,PROBLEM1,X,X1,N):- PROBLEM=knapsack(min(_Z=C*X),st(A*X=fail;true). :- dynamic solved_data/4. generate_a_new_constraint((P,XP),X,PLAN):- set_nonrecurrent_01_value_at(P,X,XP,PLAN). set_nonrecurrent_01_value_at(P,X,XP,PLAN):- member(XP,[0,1]), nth1(P,X,XP), \+ member((P,_),PLAN). transform_constraints_into_projection(X,PLAN,PROJECT,USE):- findall(XP,(nth1(P,X,XP),(member((P,XP),PLAN)->true;true)),X), findall(Y,(nth1(P,X,_),(member((P,_),PLAN)->Y=1;Y=0)),PROJECT), findall(Y,(nth1(P,X,_),((member((P,XP),PLAN),XP==1)->Y=1;Y=0)),USE). project_subproblem(PROBLEM,PROBLEM,[]):-!. project_subproblem(PROBLEM,SUBPROBLEM,PLAN,[Used,Costed]):- PROBLEM=knapsack(min(_Z=C*X),st(A*X= (nl,write(infeasible(PROJECT,B-Used)));true). update_solved_data(PROBLEM,FB,PLAN):- retract(solved_data(current(K0),D1,D2,D3)), K is K0 + 1, assert(solved_data(log(K0),D1,D2,D3)), assert(solved_data(current(K),PROBLEM,FB,PLAN)). % fathoming (with bounding) %---------------------------- % ``Zantei Kai" means ``incumbent solution" and % ``Soku-shin"(Measure of the depth) means ``fathoming" % respectively. These use of words appear in Ref. 1. :- dynamic fathomed_data/3. fathomed_case(infeasible,[infeasible,_,_,_,_]):-!. fathomed_case(dominated,[_,Z,ZANTEI,_,_]):- ZANTEI \= infinite, Z > ZANTEI,!. fathomed_case(update,[_,_,_,X,_]):- forall((nth1(_P,X,XP), XP1 is XP), integer(XP1)), !. fathomed_case(branch,[_,_,_,X,P]):- nth1(P,X,XP), XP1 is XP, \+ integer(XP1),!. fathomed_case(anormal,[_,_,_,_,_]). %bounding(CASE,[ZANTEI,Y],[ZANTEI0,Y0],[Z,X]). bounding(infeasible,[Z0,Y0],[Z0,Y0],_). bounding(dominated,[Z0,Y0],[Z0,Y0],_). bounding(update,[Z,X],_,[Z,X]). bounding(branch,[Z0,Y0],[Z0,Y0],_). bounding(anormal,[Z0,Y0],[Z0,Y0],_). fathoming([CASE,ZANTEI,Y],[F,Z,ZANTEI0,Y0,X,P]):- fathomed_data(current(_),D1,_), D1 = [_,ZANTEI0,Y0], fathomed_case(CASE,[F,Z,ZANTEI0,X,P]), bounding(CASE,[ZANTEI,Y],[ZANTEI0,Y0],[Z,X]). update_fathomed_data([CASE,ZANTEI,Y],[FB,ZB,ZANTEI0,Y0,X1,P]):- retract(fathomed_data(current(K0),D1,D2)), K is K0 + 1, assert(fathomed_data(log(K0),D1,D2)), assert(fathomed_data(current(K),[CASE,ZANTEI,Y],[FB,ZB,ZANTEI0,Y0,X1,P])). % review and save the experimental data %-------------------------------------------------------- output_data(A,B,C,D,E,F,G):- solved_data(A,B,C,D), fathomed_data(A,[G|E],[C|F]). tell_data:- tell_goal('knapsack_exp.txt', forall_such_that, (problem(A,C,B),J,G,X,Z,F,I,Plan), output_data(J, knapsack(min(Z=C*X),st(A*X=close(S0);true), open(File,write,S), tell(File), nl, time_stamp('% file output start time ',_), nl, write('%---------- start from here ------------%'), nl, G, nl, write('%---------- end of data ------------%'), nl, time_stamp('% file output end time ',_), tell(user), close(S), % The following is to cope with the duplicated stream problem. (current_stream(File,write,S1)->close(S1);true). % 成功するゴールをすべて保存 %-------------------------------- tell_goal(File,forall,G):- G0 = (nl,write(G),write('.')), G1 = forall(G,G0), tell_goal(File,G1). tell_goal(File,forall_such_that,G,Condition):- % G should be occurred in the Condition. WRITE = (nl,write(G),write('.')), G1 = forall(Condition,WRITE), tell_goal(File,G1). % time stamp %-------------------------------- time_stamp(no,T):- get_time(U), convert_time(U,A,B,C,D,E,F,_G), T = [date(A/B/C), time(D:E:F)]. time_stamp(Word,T):- \+ var(Word), Word \= no, get_time(U), convert_time(U,A,B,C,D,E,F,_G), T = [date(A/B/C), time(D:E:F)], % format('~`.t~t~a~30|~`.t~t~70|~n',[Word]), write((Word,T)), nl. /************************************************* sample of the output **************************************************/ /* ?- solve_knapsack(K,PROBLEM,SOLVED,CONSTRAINTS). qqqqqqpqqpq No ?- tell_data. Yes % file: knapsack_exp.txt % file output start time , [date(2003/9/1), time(21:22:0)] %---------- start from here ------------% problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(0), infinite, [_G378, _G381, _G384, _G387], _G281, init, _G287, []. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(1), infinite, [1, 1, (9-7)/3, 0], -16.6667, feasible, branch, []. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(2), infinite, [1, 1, 0, (9-7)/8], -16.25, feasible, branch, [ (3, 0)]. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(3), -14, [1, 1, 0, 0], -14, feasible, update, [ (4, 0), (3, 0)]. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(4), -14, [ (1-0)/5, 0, 0, 1], -11.2, feasible, dominated, [ (4, 1), (3, 0)]. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(5), -14, [1, (6-5)/2, 1, 0], -16.5, feasible, branch, [ (3, 1)]. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(6), -14, [1, 0, 1, (6-5)/8], -16.125, feasible, branch, [ (2, 0), (3, 1)]. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(7), -15, [1, 0, 1, 0], -15, feasible, update, [ (4, 0), (2, 0), (3, 1)]. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(8), -15, [0, 0, 1, 1], -13, infeasible, infeasible, [ (4, 1), (2, 0), (3, 1)]. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(9), -15, [ (4-0)/5, 1, 1, 0], -15.8, feasible, branch, [ (2, 1), (3, 1)]. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), log(10), -15, [0, 1, 1, (4-0)/8], -11.5, feasible, dominated, [ (1, 0), (2, 1), (3, 1)]. problem([5, 2, 3, 8], [-11, -3, -4, -9], 9), current(11), -15, [1, 1, 1, 0], -18, infeasible, infeasible, [ (1, 1), (2, 1), (3, 1)]. %---------- end of data ------------% % file output end time , [date(2003/9/1), time(21:22:0)] */ %---- end