NOPT042 Constraint programming: Tutorial 5 – Advanced constraint modelling¶
The element constraint¶
Picat doesn't support indexation in constraints, e.g. X #= L[I]. Instead, it implements the element constraint:
element(I, L, X)
("X is on the Ith position in the list L"). But this constraint can also be used for "reverse indexing": when we have X and want to find its position in L. The constraint doesn't care about the direction; this is called bidirectionality.
For matrices (2D arrays), Matrix[I,J]#=X is expressed using the following constraint:
matrix_element(Matrix,I,J,X)
%load_ext ipicat
Picat version 3.7
%%picat
import cp.
main =>
L = new_array(10),
L :: 0..9,
all_different(L),
% X #= L[I], % this does not work, use element instead
element(I, L, X),
X #= 5,
solve($[max(I)],L),
println(L).
{0,1,2,3,4,6,7,8,9,5}
Example: Langford's number problem¶
Consider the following problem (formulation slightly modified from the book):
Consider two sets of the numbers from 1 to N. The problem is to arrange the 2N numbers in the two sets into a single sequence in which the two 1’s appear one number apart, the two 2’s appear two numbers apart, the two 3’s appear three numbers apart, etc.
Try to formulate a model for this problem.
!picat langford/langford.pi 8
CPU time 0.117 seconds. Backtracks: 9160
{2,4,7,2,8,6,4,1,5,1,7,3,6,8,5,3}
!picat langford/langford2.pi 12
CPU time 0.034 seconds. Backtracks: 1346 solution = [1,2,1,3,2,8,9,3,10,11,12,5,7,4,8,6,9,5,4,10,7,11,6,12] position = [1,2,4,14,12,16,13,6,7,9,10,11,3,5,8,19,18,23,21,15,17,20,22,24]
!picat langford/langford-primal.pi 8
!picat langford/langford-dual.pi 12
!picat langford/langford-channeling.pi 8
CPU time 2.06 seconds. Backtracks: 250746
CPU time 0.002 seconds. Backtracks: 970
{1,2,4,10,11,13,16,14,12,7,9,6,3,5,8,15,17,20,24,23,22,18,21,19}
CPU time 3.12 seconds. Backtracks: 250746
{1,3,1,6,7,3,8,5,2,4,6,2,7,5,4,8}
{1,9,2,10,8,4,5,7,3,12,6,15,14,11,13,16}
The assignment problem¶
From Wikipedia:
The assignment problem is a fundamental combinatorial optimization problem. In its most general form, the problem is as follows:
The problem instance has a number of agents and a number of tasks. Any agent can be assigned to perform any task, incurring some cost that may vary depending on the agent-task assignment. It is required to perform as many tasks as possible by assigning at most one agent to each task and at most one task to each agent, in such a way that the total cost of the assignment is minimized.
Alternatively, describing the problem using graph theory:
The assignment problem consists of finding, in a weighted bipartite graph, a matching of a given size, in which the sum of weights of the edges is minimum.
If the numbers of agents and tasks are equal, then the problem is called balanced assignment.
Example: Swimmers¶
(From: W. Winston, Operations Research: Applications & Algorithms.)
In medley swimming relay, a team of four swimmers must swim 4x100m, each swimmer using a different style: breaststroke, backstroke, butterfly, or freestyle. The table below gives their average times for 100m in each style. Which swimmer should swim which stroke to minimize total time?
| Swimmer | Free | Breast | Fly | Back |
|---|---|---|---|---|
| A | 54 | 54 | 51 | 53 |
| B | 51 | 57 | 52 | 52 |
| C | 50 | 53 | 54 | 56 |
| D | 56 | 54 | 55 | 53 |
Write a general model, generate larger instances, and try to make your model as efficient as possible.
!cd swimmers && picat instances.pi
Sample instance:
{A,B,C,D}
{Free,Breast,Fly,Back}
{{54,54,51,53},{51,57,52,52},{50,53,54,56},{56,54,55,53}}
Random instance:
{Swimmer1,Swimmer2,Swimmer3,Swimmer4,Swimmer5,Swimmer6,Swimmer7,Swimmer8}
{Style1,Style2,Style3,Style4,Style5,Style6,Style7,Style8}
{{46,49,54,64,53,55,55,54},{60,55,47,64,65,49,65,52},{48,60,61,61,62,59,57,54},{47,50,50,58,46,64,50,45},{48,57,62,54,46,52,61,61},{60,63,57,59,65,55,65,47},{47,60,62,64,52,53,53,54},{56,56,46,55,54,51,56,57}}
!cd swimmers && picat swimmers.pi
CPU time 57.455 seconds. Backtracks: 417964
!cd swimmers && cat swimmers.pi
import cp.
main =>
cl(sample_instance),
sample_instance(SwimmerNames, StyleNames, Times),
% go(SwimmerNames, StyleNames, Times),
M = 1000,
time2(test(M)).
test(M) =>
cl(random_instance),
MinTime = 45,
MaxTime = 65,
N = 12,
K = 12,
foreach(_ in 1..M)
random_instance(N, K, MinTime, MaxTime, SwimmerNames, StyleNames, Times),
primal_model(StyleOfSwimmer, Times, TotalTime),
solve([$min(TotalTime), StyleOfSwimmer]),
% dual_model(SwimmerOfStyle, Times, TotalTime),
% solve([$min(TotalTime), SwimmerOfStyle]),
% channeling_model(StyleOfSwimmer, SwimmerOfStyle, Times, TotalTime),
% solve([$min(TotalTime), StyleOfSwimmer ++ SwimmerOfStyle])
end.
go(SwimmerNames, StyleNames, Times) =>
primal_model(StyleOfSwimmer, Times, TotalTime),
time2(solve([$min(TotalTime), StyleOfSwimmer])),
println("\nPrimal model:"),
foreach(I in 1..SwimmerNames.length)
printf("Swimmer %w is swims %w\n", SwimmerNames[I], StyleNames[StyleOfSwimmer[I]])
end,
dual_model(SwimmerOfStyle, Times, TotalTime),
time2(solve([$min(TotalTime), SwimmerOfStyle])),
println("\nDual model:"),
foreach(J in 1..StyleNames.length)
printf("Style %w is swum by %w\n", StyleNames[J], SwimmerNames[SwimmerOfStyle[J]])
end,
channeling_model(StyleOfSwimmer, SwimmerOfStyle, Times, TotalTime),
time2(solve([$min(TotalTime), StyleOfSwimmer ++ SwimmerOfStyle])),
println("\nChanneling model:"),
foreach(I in 1..SwimmerNames.length)
printf("Swimmer %w is swims %w\n", SwimmerNames[I], StyleNames[StyleOfSwimmer[I]])
end,
println("or in the dual view"),
foreach(J in 1..StyleNames.length)
printf("Style %w is swum by %w\n", StyleNames[J], SwimmerNames[SwimmerOfStyle[J]])
end.
primal_model(StyleOfSwimmer, Times, TotalTime) =>
N = Times.length,
%K = Times[1].length,
K = N,
StyleOfSwimmer = new_array(N),
StyleOfSwimmer :: 1..K,
all_different(StyleOfSwimmer),
TimeOfSwimmer = new_array(N),
foreach(I in 1..N)
matrix_element(Times, I, StyleOfSwimmer[I], TimeOfSwimmer[I])
end,
TotalTime #= sum(TimeOfSwimmer).
dual_model(SwimmerOfStyle, Times, TotalTime) =>
N = Times.length,
K = N,
SwimmerOfStyle = new_array(K),
SwimmerOfStyle :: 1..N,
all_different(SwimmerOfStyle),
TimeOfStyle = new_array(K),
foreach(J in 1..K)
matrix_element(Times, J, SwimmerOfStyle[J], TimeOfStyle[J])
end,
TotalTime #= sum(TimeOfStyle).
channeling_model(StyleOfSwimmer, SwimmerOfStyle, Times, TotalTime) =>
primal_model(StyleOfSwimmer, Times, TotalTime),
dual_model(SwimmerOfStyle, Times, TotalTime),
assignment(StyleOfSwimmer, SwimmerOfStyle). %chanelling constraint
Modelling functions¶
In general, how to model a function (mapping) f:A→B? Let's say A={1,…,n} and B={1,…,k}.
- as an array:
F = new_array(N},
F :: 1..K.
- injective:
all_different(F) - surjective: a partition of A into classes labelled by B, to each element of B map a set of elements of A. In Picat we can model set as their characteristic vectors. More on modelling with sets later.
- partial function: a dummy value for undefined inputs
- dual model: switch the role of variables and values (not a function unless F injective, see above)
- channelling: combine the primal and dual models, if it is a bijection, then use
assignment(F, FInv)
!picat functions.pi 4 4
{1,2,3,4}
{1,2,3,4}
!cat functions.pi
import cp.
main([N, K]) =>
N := N.to_int,
K := K.to_int,
% function
F = new_array(N),
F :: 1..K,
% injective
all_different(F),
% dual model if it is a bijection (K=N and injective)
FInv = new_array(K),
FInv :: 1..N,
% channeling if it is a bijection (K=N and injective)
assignment(F, FInv),
% % dual model in general
% FInv = new_array(K, N),
% FInv :: 0..1,
% % surjective in general
% foreach(J in 1..N)
% sum([FInv[I, J]: I in 1..K]) #>= 1
% end,
% % channeling in general
% foreach(I in 1..N)
% foreach(J in 1..N)
% (FInv[J, I] #= 1) #<=> (F[I] #= J)
% end
% end,
solve(F ++ FInv),
println(F),
println(FInv).