%% bouncingBall.pl %% %% Using Constraint Logic Programming to solve a planning problem for bouncing balls. %% This example adapts a DFS planner without heuristics and reduces the planning %% problem to solving a Linear Programming problem with an external EPLEX library. %% This library is designed to interface with Linear Programming solvers only. %% For this reason, the trajectory of the bouncing ball is approximated as a %% linear function of time, since otherwise a different interface to an external %% non-linear programming solver would be required. %% %% Author: Mikhail Soutchanski. %% Contact Email: soutchanski (at) protonmail (dot) com %% Copyright (c) 2020, 2023 Mikhail Soutchanski %% %% This software is available with the MIT Licence https://opensource.org/license/MIT %% See details about physics at %% https://www.physicsclassroom.com/class/vectors/Lesson-2/Characteristics-of-a-Projectile-s-Trajectory %% This example works with ECLiPSe Prolog 6.1 and more recent releases. /* Eplex: Harnessing Mathematical Programming Solvers for Constraint Logic Programming by Kish Shen and Joachim Schimpf. In: Proceedings of the 11th International Conference on Principles and Practice of Constraint Programming, pp 622–636. CP-2005, LNCS 3709, Springer Verlag, 2005. https://link.springer.com/chapter/10.1007/11564751_46 http://eclipseclp.org/reports/eplex_cp2005.pdf */ :- lib(eplex). /* EPLEX library with an external linear constraints solver */ /* The following simple iterative deepening DFS planner is taken from the Hector Levesque's course notes prepared for the 1st year SCI-199 course. A similar program is in Chapter 9 of his textbook "Thinking as Computation", MIT Press, 2017. It turns out that his program can be adapted to solve planning problems for the mixed discrete-continuous (hybrid) systems. */ solve_problem(N,Plan) :- C0 is cputime, max_length(L,N), reachable(S,L), goal_state(S), myReverse(L,Plan), Cf is cputime, Duration is Cf - C0, nl, write('Elapsed time (sec): '), write(Duration), nl. max_length([],N) :- N >= 0. max_length([_|L],N1) :- N1 > 0, N is N1 -1, max_length(L,N). myReverse(L1,L2) :- myRev(L1, [], L2). myRev([], L1, L1). myRev([H|L1], L2, L3) :- myRev(L1, [H|L2], L3). reachable(S, []) :- initial_state(S). reachable(S2, [M|ListOfActions]) :- reachable(S1, ListOfActions), legal_move(S2, M, S1). %using legal_move predicate with natural actions legal_move([A|S], A, S) :- natural(Nat), ( poss(Nat,S) -> A=Nat ; (agent(A), poss(A,S) ) ). /* % The above rule can be simplified as follows. legal_move([A|S], A, S) :- natural(A), poss(A,S). legal_move([A|S], A, S) :- agent(A), poss(A,S). */ ball(b1). ball(b2). epsilon(E) :- E is (1/10). %approximation natural(bounce(Ball,Time)). natural(atPeak(Ball,Time)). agent(drop(Ball,Time)). agent(catch(Ball, Time)). %--------- Precondition axioms ---------------% poss(drop(Ball,Time),S) :- ball(Ball), not falling(Ball,S), not flying(Ball,S), start(S, T), T $=< Time. poss(catch(Ball, Time), S) :- ball(Ball), ( falling(Ball, S); flying(Ball,S) ), start(S,T), T $=< Time. poss(bounce(Ball, Time), S) :- ball(Ball), falling(Ball, S), epsilon(E), dist(Ball, Distance, Time, S), Distance $= 0, vel(Ball, Velocity, Time, S), (Velocity $>= E ), start(S,T), T $=< Time. poss(atPeak(Ball, Time), S) :- ball(Ball), flying(Ball, S), epsilon(E), dist(Ball, Distance, Time, S), Distance $>= E, vel(Ball, Velocity, Time, S), Velocity $=0, start(S,T), T $=< Time. /* falling(B,S): ball B is falling down and accelerating under the Earth gravity. flying(B,S): ball B bounced and it is flying up and decelerating due to gravity. Assume the vertical axis is oriented downwards, i.e., if a ball is falling down, then its speed is positive and increases. But when it bounces, its speed is negative and decreases. For simplicity: elastic ball, no damping. */ %-----Successor state axioms ------------% falling(Ball, [A|S]) :- A = drop(Ball, Time). falling(Ball, [A|S]) :- A = atPeak(Ball, Time). falling(Ball, [A|S]) :- not A = catch(Ball,Time), not A=bounce(Ball, Time), falling(Ball, S). flying(Ball, [A|S]) :- A = bounce(Ball, Time). flying(Ball, [A|S]) :- not A = catch(Ball,Time), not A=atPeak(Ball, Time), flying(Ball, S). numberFlights(0,B, []). numberFlights(N,B, [A | S]) :- A = atPeak(B,_), numberFlights(M,B,S), N is M+1. numberFlights(N,B, [A | S]) :- not A = atPeak(B,_), numberFlights(N,B,S). /* The height of the ball from the floor dist(B,Dist,T,S) changes continuously, no matter what actions or events happen. */ init_dist(Ball,Dist,[Action|S]) :- time(Action, Time), dist(Ball,Dist,Time,S). init_vel(b1,0,[]). init_vel(b2,0,[]). init_vel(Ball,Vel, [Action|S]) :- Action = catch(Ball, Time), Vel $= 0. init_vel(Ball,Vel, [Action|S]) :- Action = bounce(Ball, Time), time(Action, Time), vel(Ball, Old_Vel, Time, S), Vel $= -1*Old_Vel. init_vel(Ball,Vel,[Action|S]) :- time(Action,T), not Action = catch(Ball, _), not Action = bounce(Ball, _), vel(Ball, Vel, T, S). /* ------------ State evolution axioms -------------- See details at https://arxiv.org/abs/1807.04861 Hybrid Temporal Situation Calculus Vitaliy Batusov, Giuseppe De Giacomo, Mikhail Soutchanski */ vel(b1, 0, 0, []). vel(b2, 0, 0, []). vel(Ball, Vel, Time, S) :- not falling(Ball, S), not flying(Ball, S), Vel $= 0. vel(Ball, Vel, Time, S) :- falling(Ball, S), init_vel(Ball, Old_Vel, S), start(S, T1), Time $>= T1, Vel $= Old_Vel + 9.81*(Time - T1). /* since Old_vel >0, Ball accelerates */ vel(Ball, Vel, Time, S) :- flying(Ball, S), init_vel(Ball, Old_Vel, S), start(S, T1), Time $>= T1, Vel $= Old_Vel + 9.81*(Time - T1). /* since Old_vel < 0, Ball decelerates */ dist(b1, 100, 0, []). dist(b2, 150, 0, []). dist(Ball, Dist, Time, S) :- not falling(Ball, S), not flying(Ball, S), init_dist(Ball, Dist, S). /* Since eplex solves linear constraints only, approximate the law for distance with a linear function of time, instead of a quadratic function. */ dist(Ball, Dist, Time, S) :- falling(Ball, S), init_dist(Ball, Old_Dist, S), start(S,T1), Time $>= T1, Dist $= Old_Dist - 0.5*9.81*(Time - T1). /* Dist decreases if Ball falls down */ dist(Ball, Dist, Time, S) :- flying(Ball, S), init_dist(Ball, Old_Dist, S), start(S,T1), Time $>= T1, Dist $= Old_Dist + 0.5*9.81*(Time - T1). /* Dist increases if Ball moves up */ %--------------Goal conditions and time minimization-------% firstAction([A], [A]). /* The following rule makes sure there is no other solution upon backtracking */ firstAction([A1, A2 |S], Sit) :- firstAction([A2 | S], Sit). duration([A|S], Duration) :- firstAction([A|S], [FirstAct]), time(FirstAct, Tinit), time(A, Tlast), Duration $= (Tlast - Tinit), Duration $>= 0. /* Linear Programming: schedule actions to execute as soon as possible */ minimizeDuration(TotalDur, Cost) :- eplex_solver_setup(min(TotalDur)), eplex_solve(Cost). goal_state(S) :- satisfy(S,TotalTime), duration(S,Time), /* Let Delta be extra time that passed since the start of the last action*/ Delta $>= 0, TotalTime $= Time + Delta, minimizeDuration(TotalTime, Cost), chooseTimes(S), eplex_get(cost,EndTime), nl, write("EndTime= "), writeln(EndTime). satisfy(S,Time) :- start(S,Tinit), Time $>= Tinit, falling(b1,S), falling(b2,S), numberFlights(N1,b1,S), N1 >= 1, numberFlights(N2,b2,S), N2 >= 1, vel(b1, V1, Time, S), vel(b2, V2, Time, S), V1 - V2 $= 0, dist(b1, D1, Time, S), dist(b2, D2, Time, S), D1 - D2 $= 0 . %--------------Auxiliary axioms ---------------------------% time( drop(Ball,Time), Time). time( catch(Ball,Time), Time). time( bounce(Ball, Time), Time). time( atPeak(Ball, Time), Time). initial_state([]). %Reiter domain independent axioms about start time of situations. start([], 0). start([A|S], T) :- time(A,T). chooseTimes([]). chooseTimes([A|S]) :- chooseTimes(S), time(A,T), eplex_var_get(T, typed_solution, T). /* ECLiPSe Constraint Logic Programming System [kernel threads] Kernel and basic libraries copyright Cisco Systems, Inc. and subject to the Cisco-style Mozilla Public Licence 1.1 (see legal/cmpl.txt or http://eclipseclp.org/licence) Source available at www.sourceforge.org/projects/eclipse-clp GMP library copyright Free Software Foundation, see legal/lgpl.txt For other libraries see their individual copyright notices Version 7.0 #63 (x86_64_linux), Sun Apr 24 18:46 2022 [eclipse 1]: solve_problem(1, Plan). No (0.00s cpu) [eclipse 2]: solve_problem(7, Plan). No (0.08s cpu) [eclipse 3]: solve_problem(8, Plan). EndTime= 91.743119266055 Elapsed time (sec): 0.180136583 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), drop(b1, 50.9683995922528), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 71.35575942915392), drop(b2, 91.743119266055047), atPeak(b1, 91.743119266055047)] Yes (0.18s cpu, solution 1, maybe more) ? [eclipse 4]: solve_problem(8, Plan). EndTime= 91.743119266055 Elapsed time (sec): 0.181632581 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), drop(b1, 50.9683995922528), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 71.35575942915392), drop(b2, 91.743119266055047), atPeak(b1, 91.743119266055047)] Yes (0.18s cpu, solution 1, maybe more) ? ; EndTime= 91.743119266055 Elapsed time (sec): 0.185122109 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), drop(b1, 50.968399592252787), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 71.3557594291539), atPeak(b1, 91.743119266055032), drop(b2, 91.743119266055032)] Yes (0.19s cpu, solution 2, maybe more) ? ; EndTime= 101.936799184506 Elapsed time (sec): 0.213021238 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), drop(b1, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 81.549439347604476), drop(b2, 101.9367991845056), atPeak(b1, 101.9367991845056)] Yes (0.21s cpu, solution 3, maybe more) ? ; EndTime= 101.936799184506 Elapsed time (sec): 0.216116812 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), drop(b1, 61.162079510703364), catch(b2, 71.355759429153949), bounce(b1, 81.549439347604476), atPeak(b1, 101.9367991845056), drop(b2, 101.9367991845056)] Yes (0.22s cpu, solution 4, maybe more) ? ; EndTime= 112.130479102956 Elapsed time (sec): 0.222204706 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), drop(b2, 112.13047910295616), atPeak(b1, 112.13047910295616)] Yes (0.22s cpu, solution 5, maybe more) ? ; EndTime= 112.130479102956 Elapsed time (sec): 0.22502153 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), atPeak(b1, 112.13047910295616), drop(b2, 112.13047910295616)] Yes (0.23s cpu, solution 6, maybe more) ? ; EndTime= 112.130479102956 Elapsed time (sec): 0.231816381 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), drop(b2, 112.13047910295616), atPeak(b1, 112.13047910295616)] Yes (0.23s cpu, solution 7, maybe more) ? ; EndTime= 112.130479102956 Elapsed time (sec): 0.234644073 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), atPeak(b1, 112.13047910295616), drop(b2, 112.13047910295616)] Yes (0.23s cpu, solution 8, maybe more) ? ; EndTime= 50.9683995922528 Elapsed time (sec): 0.271247942 Plan = [drop(b2, 0), drop(b1, 10.19367991845056), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)] Yes (0.27s cpu, solution 9, maybe more) ? ; EndTime= 50.9683995922528 Elapsed time (sec): 0.296041525 Plan = [drop(b2, 0), catch(b2, 10.193679918450554), drop(b1, 10.193679918450554), drop(b2, 10.193679918450554), bounce(b1, 30.581039755351679), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.968399592252808)] Yes (0.30s cpu, solution 10, maybe more) ? ; EndTime= 50.9683995922528 Elapsed time (sec): 0.322303164 Plan = [drop(b2, 0), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), drop(b1, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)] Yes (0.32s cpu, solution 11, maybe more) ? ; EndTime= 50.9683995922528 Elapsed time (sec): 0.331785495 Plan = [drop(b2, 0), catch(b2, 10.193679918450554), drop(b1, 10.193679918450554), drop(b2, 10.193679918450554), bounce(b1, 30.581039755351679), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.968399592252808)] Yes (0.33s cpu, solution 12, maybe more) ? ; EndTime= 50.9683995922528 Elapsed time (sec): 0.358241188 Plan = [drop(b2, 0), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), drop(b1, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)] Yes (0.36s cpu, solution 13, maybe more) ? ; EndTime= 91.743119266055 Elapsed time (sec): 0.485687106 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), drop(b1, 50.9683995922528), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 71.35575942915392), drop(b2, 91.743119266055047), atPeak(b1, 91.743119266055047)] Yes (0.49s cpu, solution 14, maybe more) ? ; EndTime= 91.743119266055 Elapsed time (sec): 0.489053124 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), drop(b1, 50.968399592252787), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 71.3557594291539), atPeak(b1, 91.743119266055032), drop(b2, 91.743119266055032)] Yes (0.49s cpu, solution 15, maybe more) ? ; EndTime= 101.936799184506 Elapsed time (sec): 0.517981219 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), drop(b1, 61.162079510703364), catch(b2, 71.35575942915392), bounce(b1, 81.549439347604476), drop(b2, 101.9367991845056), atPeak(b1, 101.9367991845056)] Yes (0.52s cpu, solution 16, maybe more) ? ; EndTime= 101.936799184506 Elapsed time (sec): 0.521125555 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), drop(b1, 61.162079510703364), catch(b2, 71.355759429153949), bounce(b1, 81.549439347604476), atPeak(b1, 101.9367991845056), drop(b2, 101.9367991845056)] Yes (0.52s cpu, solution 17, maybe more) ? ; EndTime= 112.130479102956 Elapsed time (sec): 0.52731797 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), drop(b2, 112.13047910295616), atPeak(b1, 112.13047910295616)] Yes (0.53s cpu, solution 18, maybe more) ? ; EndTime= 112.130479102956 Elapsed time (sec): 0.53014014 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), atPeak(b1, 112.13047910295616), drop(b2, 112.13047910295616)] Yes (0.53s cpu, solution 19, maybe more) ? ; EndTime= 112.130479102956 Elapsed time (sec): 0.536990292 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), drop(b2, 112.13047910295616), atPeak(b1, 112.13047910295616)] Yes (0.54s cpu, solution 20, maybe more) ? ; EndTime= 112.130479102956 Elapsed time (sec): 0.539800036 Plan = [drop(b2, 0), bounce(b2, 30.581039755351682), atPeak(b2, 61.162079510703364), catch(b2, 71.35575942915392), drop(b1, 71.35575942915392), bounce(b1, 91.743119266055047), atPeak(b1, 112.13047910295616), drop(b2, 112.13047910295616)] Yes (0.54s cpu, solution 21, maybe more) ? ; EndTime= 50.9683995922528 Elapsed time (sec): 0.576289241 Plan = [drop(b2, 0), drop(b1, 10.19367991845056), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)] Yes (0.58s cpu, solution 22, maybe more) ? ; EndTime= 50.9683995922528 Elapsed time (sec): 0.586033084 Plan = [drop(b2, 0), catch(b2, 10.193679918450554), drop(b1, 10.193679918450554), drop(b2, 10.193679918450554), bounce(b1, 30.581039755351679), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.968399592252808)] Yes (0.59s cpu, solution 23, maybe more) ? ; EndTime= 50.9683995922528 Elapsed time (sec): 0.596641336 Plan = [drop(b2, 0), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), drop(b1, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)] Yes (0.60s cpu, solution 24, maybe more) ? ; EndTime= 50.9683995922528 Elapsed time (sec): 0.606203237 Plan = [drop(b2, 0), catch(b2, 10.193679918450554), drop(b1, 10.193679918450554), drop(b2, 10.193679918450554), bounce(b1, 30.581039755351679), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.968399592252808)] Yes (0.61s cpu, solution 25, maybe more) ? ; EndTime= 50.9683995922528 Elapsed time (sec): 0.616902587 Plan = [drop(b2, 0), catch(b2, 10.19367991845056), drop(b2, 10.19367991845056), drop(b1, 10.19367991845056), bounce(b1, 30.581039755351682), bounce(b2, 30.581039755351682), atPeak(b1, 50.9683995922528), atPeak(b2, 50.9683995922528)] Yes (0.62s cpu, solution 26, maybe more) ? ; No (0.62s cpu) [eclipse 5]: halt. */