Here's my solution for the water jugs problem
:- use_module(library(clpfd)).
initial(state(8, 0, 0)).
final(state(4, 4, 0)).
constraint(A0, A1, B0, B1, C0, C1, V) :-
A0 #< A1,
( B0 #> 0, T #= min(V - A0, B0), A1 #= A0 + T, B1 #= B0 - T, C1 #= C0
; C0 #> 0, T #= min(V - A0, C0), A1 #= A0 + T, C1 #= C0 - T, B1 #= B0
).
transition(state(A0, B0, C0), state(A1, B1, C1)) :-
( constraint(A0, A1, B0, B1, C0, C1, 8)
; constraint(B0, B1, A0, A1, C0, C1, 5)
; constraint(C0, C1, A0, A1, B0, B1, 3)
).
solve(A, A, _, [A]).
solve(A, B, P, [A|Q]) :-
transition(A, A1),
\+ member(A1, P),
solve(A1, B, [A|P], Q).
path(P) :-
initial(S0),
final(S),
solve(S0, S, [], P).
Is there a way to find the P of minimal length without traversing all options?
Here is a solution that makes more use of the power of clpfd: First state the problem, then try to solve it (using labeling/2 or similar). Given that we do not know the length of the (shortest) path, this will generate larger and larger problems until a solution is found. In my code, I do not prevent visiting the same state twice (but this could be added in the same way as in the MiniZinc model written by @DavidTonhofer, or as some post-processing). However, in order to ensure a finite search space, I've added code to stop the problem generation if the length of the path is longer than (5+1)*(3+1), as this is an upper bound on the number of different states (assuming we have do not add or remove water outside of the 3 jugs).
:- use_module(library(clpfd)).
initial(state(8, 0, 0)).
final(state(4, 4, 0)).
constraint(A0,A1,B0,B1,C0,C1,R,Max):-
T#=min(Max-B0,A0),
R in 0..1,
R#==>T#>0,
R#==>A1#=A0-T,
R#==>B1#=B0+T,
R#==>C1#=C0.
transition(state(A0, B0, C0), state(A1, B1, C1)) :-
A0+B0+C0#=A1+B1+C1,
A0 in 0..8,
B0 in 0..5,
C0 in 0..3,
A1 in 0..8,
B1 in 0..5,
C1 in 0..3,
constraint(A0,A1,B0,B1,C0,C1,RAB,5),
constraint(B0,B1,A0,A1,C0,C1,RBA,8),
constraint(A0,A1,C0,C1,B0,B1,RAC,3),
constraint(C0,C1,A0,A1,B0,B1,RCA,8),
constraint(C0,C1,B0,B1,A0,A1,RCB,5),
constraint(B0,B1,C0,C1,A0,A1,RBC,3),
RAB+RBA+RAC+RCA+RCB+RBC#=1.
solve(A, A, Xs, [A]):-
labeling([],Xs).
solve(A, B, Xs, [A|Q]) :-
length(Xs, L),
L < 24*3,
transition(A, A1),
A=state(X1,X2,X3),
solve(A1, B, [X1,X2,X3|Xs], Q).
path(P) :-
initial(S0),
final(S),
solve(S0, S, [], P).
I tried to keep the code relatively close to the one in the question. The main difference is that all the prolog-level disjunctions in transition/2 and constraint/7 have been removed and replaced by reification. In particular, I added the parameter R to constraint/8 which is equal to 1 if that specific transition is taken. Then I state in transition/2 that exactly one of the transitions must take place.
I must add that this formulation is not particularly efficient and I would not be surprised to find out that one can solve the problem more efficiently with either a different clpfd formulation or without using clpfd at all.