(!****************************************************************
CP example problems
===================
file runfapp.mos
````````````````
Solution algorithm for solving the
Frequency Assignment Problem with Polarity
- Master model executing the submodel fappmod.mos -
This problem, taken from the ROADEF'2001 challenge, is a Frequency
Assignment Problem with Polarization constraints. It originates
from the CALMA European project (Combinatorial Algorithms
for Military Applications), from 1993 to 1995. The formulation
has been extended to take into account the polarisation
constraints and the undercontrolled relaxation of
electromagnetic compatibility constraints.
You can find more information about the challenge on
the following URL:
http://challenge.roadef.org/2001/en/
A complete description of the problem can be dowloaded here:
http://challenge.roadef.org/2001/en/sujet.php
The algorithm implemented is a two phase algorithm:
- First find a good lower bound on the relaxation level 'k' of
the CEM (ElectroMagnetic Compatibility constraints). This is
done by successively propagation CEM constraints by decreasing
relaxation level. When the problem becomes infeasible we obtain
a lower bound on 'k'
- We then apply a simple solution algorithm to find (on some instances)
an optimal solution relative to the relaxation level 'k'.
The lower bound algorithm enables the finding of an exact lower bound
just using constraint propagation for 37 of the 44 instances published
by the challenge.
(c) 2008 Artelys S.A. and Fair Isaac Corporation
Creation: 2005
*****************************************************************!)
model "Run FAPP"
uses "mmsystem", "mmjobs"
parameters
FILENAME = "Data/fapp_instance_01.in"
! You can alternatively choose another source data set with the expected results below.
!FILENAME = "Data/fapp_instance_02.in"
!FILENAME = "Data/fapp_instance_03.in"
!FILENAME = "Data/fapp_instance_04.in"
end-parameters
(!
===========================
Instance Kopt Bound
===========================
fapp_instance_01 4 4
fapp_instance_02 8 8
fapp_instance_03 10 10
fapp_instance_04 4 4
===========================
!)
! Declaration of procedures
forward procedure read_data
forward procedure clean_data
! Data model declaration
declarations
nbTravel: integer ! Number of travels
nbFreqDom: integer ! Number of frequency domains
nbCI: integer ! Number of hard constraints
nbCEME: integer ! Number of CEME soft constraints
nbCEMD: integer ! Number of CEMD soft constraints
FREQDOMS: array(range,range) of integer ! Frequency domains
SIZEDOMS: array(range) of integer ! Frequency domain sizes
VKEV: array(range,range) of integer ! Soft constraint indices
VKE : array(range,range) of integer ! CEME coefficients
VKD : array(range,range) of integer ! CEMD coefficients
TR: array(range, {"F","P"}) of integer ! Travel data
CI: dynamic array(set of integer, set of integer, {"F","P"}, {"E","I"}) of integer
! Hard constraints data
fappmodel: Model ! Reference to the CP model
lb: integer ! Lower bound value from Phase 1
end-declarations
! **** Compile and load the CP model ****
res:=compile("fappmod.mos") ! Compile the CP model
load(fappmodel, "fappmod.bim") ! Load the CP model
read_data ! Read problem data from file
! **** Solution algorithm ****
starttime:= gettime
run(fappmodel, "PHASE=1") ! Solving Phase 1: lower bound
wait ! Wait for model to finish
dropnextevent ! Ignore termination event
initializations from "raw:" ! Retrieve lower bound value
lb as "shmem:lb"
end-initializations
writeln("Lower bound = " + lb);
run(fappmodel, "PHASE=2,BOUND=" + lb) ! Solving Phase 2: solve the prob.
wait ! Wait for model to finish
dropnextevent ! Ignore termination event
writeln("Total time: ", gettime-starttime, "sec")
clean_data ! Delete data from shared memory
!*************************************************************************
! Read data from file and pass them to the CP model (via shared memory)
!*************************************************************************
procedure read_data
declarations
SC: string ! Marker of data table names
val,freq,v0,v1,v2:integer
s0,s1:string
end-declarations
writeln("***************************************************************")
fopen(FILENAME, F_INPUT) ! Open the data file for reading
writeln("Reading file `", FILENAME, "'...")
nbFreqDom := 0 ! number of frequency domains
nbCI := 0 ! number of imperative constraints
nbCEME := 0 ! number of CEM constraints
nbCEMD := 0
nbTravel := 0 ! number of travels involved
while (not iseof) do
read(SC) ! Read the name of the table on a line
case SC of
"DM": do ! Reading a domain for a travel
readln(val,freq) ! reading the domain number and the frequency to add to this domain
FREQDOMS(val,SIZEDOMS(val)) := freq
SIZEDOMS(val) := SIZEDOMS(val) + 1
if (val+1 > nbFreqDom) then
nbFreqDom := val+1 ! increase the counter of frequency domains
end-if
end-do
"TR": do ! Reading a travel definition
readln(v0,v1,v2)
TR(nbTravel, "F"):= v1 ! Travel number 'nbTravel' takes values in the frequency domain number 'v1'
TR(nbTravel, "P"):= v2 ! Travel number 'nbTravel' takes values in the polarity domain number 'v2'
nbTravel := nbTravel + 1 ! Increase the counter of travels
end-do
"CI": do ! Reading an imperative constraint
readln(v0,v1,s0,s1,v2) !
CI(v0,v1,s0,s1):= v2 !
nbCI := nbCI + 1
end-do
"CE": do ! Reading the relaxation table for CEM constraints when polarities are equal
read(VKEV(nbCEME,0),VKEV(nbCEME,1)) ! get the frequency involved in the CEM constraint
forall (kk in 0..9) read(VKE(nbCEME,kk)) ! read the relaxation table
readln(VKE(nbCEME,10))
nbCEME := nbCEME + 1 ! increase the counter of CEM constraints
end-do
"CD": do ! Reading the relaxation table for CEM constraints when polarities are different
read(v0,v1) ! get the frequency involved in the CEM constraint
forall (kk in 0..9) read(VKD(nbCEMD,kk)) ! read the relaxation table
readln(VKD(nbCEMD,10))
nbCEMD := nbCEMD + 1 ! increase the counter of CEM constraints
end-do
end-case
end-do
fclose(F_INPUT)
! Write some problem statistics
writeln("Problem Stats :")
writeln(nbTravel ," Travels")
writeln(nbFreqDom," Frequency domains")
writeln(nbCI ," Hard Constraints")
writeln(nbCEME ," CEME Soft Constraints")
writeln(nbCEMD ," CEMD Soft Constraints")
! Write data to memory
initializations to "raw:"
nbTravel as "shmem:nbTravel"
nbFreqDom as "shmem:nbFreqDom"
nbCI as "shmem:nbCI"
nbCEME as "shmem:nbCEME"
nbCEMD as "shmem:nbCEMD"
FREQDOMS as "shmem:FREQDOMS"
SIZEDOMS as "shmem:SIZEDOMS"
TR as "shmem:TR"
CI as "shmem:CI"
VKEV as "shmem:VKEV"
VKE as "shmem:VKE"
VKD as "shmem:VKD"
end-initializations
writeln("***************************************************************")
end-procedure
!*************************************************************************
! Delete data stored in shared memory
!*************************************************************************
procedure clean_data
fdelete("shmem:nbTravel")
fdelete("shmem:nbFreqDom")
fdelete("shmem:nbCI")
fdelete("shmem:nbCEME")
fdelete("shmem:nbCEMD")
fdelete("shmem:FREQDOMS")
fdelete("shmem:SIZEDOMS")
fdelete("shmem:TR")
fdelete("shmem:CI")
fdelete("shmem:VKEV")
fdelete("shmem:VKE")
fdelete("shmem:VKD")
fdelete("shmem:lb")
end-procedure
end-model
|
(!****************************************************************
CP example problems
===================
file fappmod.mos
````````````````
Frequency Assignment Problem with Polarity
- Submodel executed by runfapp.mos -
(c) 2008 Artelys S.A. and Fair Isaac Corporation
Creation: 2005, rev. Mar. 2013
*****************************************************************!)
model FAPP
uses "kalis", "mmjobs"
parameters
PHASE=1 ! 1: lower bound, 2: solve problem
BOUND=0 ! lower bound value (phase 2 only)
end-parameters
! Evaluate the relaxation level of the current solution
forward function evaluate_solution(k:integer): integer
! Post the electromagnetic compatibility constraints for a certain level 'k'
forward function post_CEMX(k:integer) : boolean
! Find a solution minimizing 'k' knowing that kmin is a lower bound on 'k'
forward procedure solve_FAPP(kmin:integer)
! Find a lower bound on 'k'
forward function find_lb_FAPP : integer
declarations
nbTravel:integer ! Number of travels
nbFreqDom:integer ! Number of frequency domains
nbCI:integer ! Number of hard constraints
nbCEME:integer ! Number of CEME soft constraints
nbCEMD:integer ! Number of CEMD soft constraints
Travels: range ! Set of Travels
Freqs: range ! Set of frequency domains
FREQDOMS: array(Freqs,range) of integer ! Frequency domains
SIZEDOMS: array(Freqs) of integer ! Frequency domain sizes
VKEV: array(range,range) of integer ! Soft constraint indices
VKE : array(range,range) of integer ! CEME coefficients
VKD : array(range,range) of integer ! CEMD coefficients
Types={"F","P"} ! Type of constraint ('F' stands for frequency and 'P' for polarity)
Ops={"E","I"} ! Type of constraint ('E' stands for equality and 'I' for inequality)
TR: array(Travels, Types) of integer ! Travel data
CI:dynamic array(Travels, Travels, Types, Ops) of integer ! Hard constraints data
end-declarations
! ***** Read data from memory *****
initializations from "raw:"
nbTravel as "shmem:nbTravel" nbFreqDom as "shmem:nbFreqDom"
nbCI as "shmem:nbCI" nbCEME as "shmem:nbCEME" nbCEMD as "shmem:nbCEMD"
FREQDOMS as "shmem:FREQDOMS" SIZEDOMS as "shmem:SIZEDOMS"
TR as "shmem:TR" CI as "shmem:CI"
VKEV as "shmem:VKEV" VKE as "shmem:VKE" VKD as "shmem:VKD"
end-initializations
finalize(Travels)
! ***** Define the problem *****
declarations
freq: array(Travels) of cpvar ! Frequency variables
pols: array(Travels) of cpvar ! Polarity variables
strategy: array(range) of cpbranching ! Branching strategy
BACKTRACKSLIMIT: integer ! Limit on no. of backtracks
end-declarations
! Domains of decision variables
forall(t in Travels) do ! for all travels:
setdomain(freq(t), union(f in 0..SIZEDOMS(TR(t,"F"))-1) {FREQDOMS(TR(t,"F"),f)}) ! set the corresponding frequency domain
! set the corresponding polarity domain
case TR(t,"P") of
-1: pols(t)=-1
0: setdomain(pols(t), {-1,1})
1: pols(t)=1
end-case
end-do
! Building Hard constraints
forall(i,j in Travels, s in Types, o in Ops | exists(CI(i,j,s,o))) do
case s of
"F": if (CI(i,j,s,o) = 0) then ! a constraint between frequency of travels
if (o= "E") then ! an equality constraint
freq(i) = freq(j) ! posting the constraint Fi = Fj
elif (o= "I") then ! a difference constraint
freq(i) <> freq(j) ! posting the constraint Fi <> Fj
end-if
elif (CI(i,j,s,o) <> -1 ) then
if (o= "E") then ! an equality constraint
CI(i,j,s,o) = distance(freq(i),freq(j)) ! Hard Constraint: |Fi - Fj| == CI(i,j,s,o)
elif (o = "I") then ! a difference constraint
CI(i,j,s,o) <> distance(freq(i),freq(j)) ! Hard Constraint: |Fi - Fj| <> CI(i,j,s,o)
end-if
end-if
"P": if (CI(i,j,s,o) = 0) then ! a constraint between polarity of travels
if (o= "E") then ! an equality constraint
pols(i) = pols(j) ! Hard Constraint: Pi == Pj
elif (o= "I") then ! a difference constraint
pols(i) <> pols(j) ! Hard Constraint: Pi <> Pj
end-if
end-if
end-case
end-do
! ***** Start solution algorithm *****
if PHASE=1 then
lb := find_lb_FAPP ! Determine a lower bound
initializations to "raw:" ! Send solution to master problem
lb as "shmem:lb"
end-initializations
else
solve_FAPP(lb) ! Solve the problem
end-if
!============================================================================
! Calculating a lower bound by posting soft constraints
!============================================================================
function find_lb_FAPP: integer
k := 11
while ((k >= 0) and post_CEMX(k)) k := k - 1
writeln("==================================================================")
writeln("LOWER BOUND : K = " + (k+1))
writeln("==================================================================")
returned := k+1
end-function
!============================================================================
! Solving algorithm
!============================================================================
procedure solve_FAPP(kmin:integer)
currentStrategy := 0 ! current strategy
strategy(0) := assign_var(KALIS_SMALLEST_DOMAIN,KALIS_RANDOM_VALUE) ! first strategy
strategy(1) := assign_var(KALIS_SDOMDEG_RATIO,KALIS_NEAREST_VALUE) ! second strategy
strategy(2) := assign_var(KALIS_RANDOM_VARIABLE,KALIS_NEAREST_VALUE) ! third strategy
k := 11 ! initial relaxation level
BACKTRACKSLIMIT := integer(nbTravel) ! initial backtracks limit
while ((k >= kmin)) do ! while the relaxation level of the current solution is greater than our lower bound.
cp_set_branching(strategy(currentStrategy)) ! set the search strategy
setparam("KALIS_MAX_COMPUTATION_TIME",2) ! limit the computation time to two seconds
writeln("### Now searching for a solution of level ", k, " / ", kmin, " (", BACKTRACKSLIMIT, ")")
if (not cp_find_next_sol) then ! if no solution was found in the given time
currentStrategy := (currentStrategy + 1) mod 3 ! switch the search strategy for diversification of search
if (currentStrategy = 0) then ! all the strategy have been used
BACKTRACKSLIMIT := BACKTRACKSLIMIT + 10 ! increase backtracks limit
end-if
writeln("No solution has been found ... diversification of the search: ", currentStrategy)
else ! a solution has been found
BACKTRACKSLIMIT := integer(nbTravel) ! reset the backtracks limit
currentStrategy := 0 ! reset the current strategy
k := evaluate_solution(11) ! evaluate solution
writeln( "==================================================================")
writeln("a solution of level k = ", k, " has been found with strategy ", currentStrategy)
writeln( "==================================================================")
cp_reset_search ! resetting the search tree
if not post_CEMX(k-1) then ! posting the soft CEM constraints at level k-1
writeln("optimal") ! a contradiction has occured so the minimal k level i
break
end-if
k := k - 1 ! solution current relaxation level is k-1
end-if
cp_reset_search ! reset the search tree
end-do
writeln("==================================================================")
writeln("OPTIMAL SOLUTION FOUND : K = " + (k))
writeln("==================================================================")
end-procedure
!============================================================================
! Posting soft constraints
!============================================================================
function post_CEMX(k:integer): boolean
writeln("### Now Posting " + (nbCEME) + " CEMX Constraints of level " + k)
returned := true
if (k <>11) then
forall(c in 0..nbCEME-1) do ! for all CEM soft constraints:
f0 := VKEV(c,0) ! first frequency variable involved in the constraint
f1 := VKEV(c,1) ! second frequency variable involved in the constraint
vke := VKE(c,k) ! minimal distance between frequency variables when polarity are equals
vkd := VKD(c,k) ! minimal distance between frequency variables when polarity are differents
if not cp_post(distance(freq(f0),freq(f1)) >= vkd) then ! posting and propagating the constraint
returned := false ! a contradiction occured proving that no solution of level k exists
break
end-if
implies(pols(f0) = pols(f1),distance(freq(f0),freq(f1)) >= vke) ! posting the distance constraint
end-do
end-if
end-function
!============================================================================
! Evaluating the relaxation level of the current solution
!============================================================================
function evaluate_solution(k:integer): integer
forall (kk in 0..10) nbViols(kk) := 0 ! initialization of the constraints violation table
ksol := -1
forall(f in 0..nbTravel-1) do
settarget(freq(f),getsol(freq(f)))
settarget(pols(f),getsol(pols(f)))
end-do
forall(c in 0..nbCEME-1) do ! for every soft CEM constraints
v0 := VKEV(c,0) ! first frequency variable involved in the constraint
v1 := VKEV(c,1) ! second frequency variable involved in the constraint
fv0 := getsol(freq(v0)) ! value of the first frequency variable in the current solution
fv1 := getsol(freq(v1)) ! value of the second frequency variable in the current solution
pv0 := getsol(pols(v0)) ! value of the first polarity variable in the current solution
pv1 := getsol(pols(v1)) ! value of the second polarity variable in the current solution
kv := 0
forall (kk in 0..k-1) do ! for all level of relaxation 'k'
vke := VKE(c,kk) ! look the
vkd := VKD(c,kk)
if (pv0 = pv1) then ! polarity are equals in solution
if not ((fv0 - fv1 >= vke) or (fv1 - fv0 >= vke)) then ! test of constraint violation
nbViols(kk) += 1 !
kv := kk !
end-if
else ! polarity are different in solution
if not ((fv0 - fv1 >= vkd) or (fv1 - fv0 >= vkd)) then ! test of constraint violation
nbViols(kk) += 1 ! increase number of constraints violation at level 'kk'
kv := kk !
end-if
end-if
end-do
if (kv > ksol) then
ksol := kv ! Relaxation level of the solution equal max of relaxation level of all violated constraints.
end-if
end-do
sumvk1 := 0
forall (kk in 0..ksol-1) sumvk1 := sumvk1 + nbViols(kk) ! compute the sum of violation at level k-1
! Write some statistics about the solution
writeln("" + nbViols(ksol) + " constraints violated at level " + ksol)
writeln("" + sumvk1 + " constraints violated at level < " + ksol)
returned := ksol+1 ! return the relaxation level of the solution
end-function
end-model
|
