A branch and bound algorithm for numerical Max-CSP

Click here to download Manuscript: Constraints-SI.tex 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Click here to view linked References A Branch and Bound Algorithm for Numerical Max-CSP. J.-M. Normand∗† A. Goldsztejn‡ M. Christie§ F. Benhamou¶ October 13, 2009 Abstract The Constraint Satisfaction Problem (CSP) framework allows users to define problems in a declarative way, quite independently from the solving process. However, when the problem is over-constrained, the answer “no solution” is generally unsatisfactory. A Max-CSP Pm = hV, D, Ci is a triple defining a constraint problem whose solutions maximize the number of satisfied constraints. In this paper, we focus on numerical CSPs, which are defined on real variables represented as floating point intervals and which constraints are numerical relations defined in intension. Solving such a problem (i.e., exactly characterizing its solution set) is generally undecidable and thus consists in providing approximations. We propose a Branch and Bound algorithm that provides under and over approximations of a solution set that maximize the maximum number mP of satisfied constraints. The technique is applied on three numeric applications and provides promising results. 1 Introduction CSP provides a powerful and efficient framework for modeling and solving problems that are well defined. However, in the modeling of real-life applications, under or over-constrained systems often occur. In embodiment design, for example, when searching for different concepts, engineers often need to tune both the constraints and the domains. In such cases, the classical CSP framework is not well adapted (a no solution answer is certainly not satisfactory) and there is a need for specific frameworks such as Max-CSP. Yet in a large number of cases, the nature of the model (over-, well or under-constrained) is a priori ∗ University of Nantes, LINA, France. E-Mail: Jean-Marie.Normand@univ-nantes.fr Lab, Universitat de Barcelona, Spain. ‡ CNRS, LINA, France. E-Mail: Alexandre.Goldsztejn@univ-nantes.fr § INRIA/IRISA, LINA, France. E-Mail: Marc.Christie@univ-nantes.fr ¶ University of Nantes, LINA, France. E-Mail: Frederic.Benhamou@univ-nantes.fr † EVENT 1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 not known and thus the choice of the solving technique is critical. Furthermore, when tackling over-constrained problems in the Max-CSP framework, most techniques focus on discrete domains. In this paper, we propose a branch and bound algorithm that approximate the solutions of a numerical Max-CSP. The algorithm computes both an under (often called inner in continuous domains) and an over (often called outer) approximation of the solution set of a Max-CSP along with the sets of constraints that maximize constraint satisfaction. These approximations ensure that the inner approximation contains only solutions of the Max-CSP while the outer approximation guarantees that no solution is lost during the solving process. As a consequence, since we aim at computing the sets of inner and outer approximations of a numerical Max-CSP, it has to be composed of constraints that give rise to continuum of solutions, namely inequalities. This paper is organized as follows: Section 2 motivates the research and presents a numerical Max-CSP application. Section 3 presents the definitions associated to the Max-CSP framework. A brief introduction to Interval analysis and its related concepts is drawn in Section 4, while our branch and bound algorithm is described in Section 5. Section 6 presents the work related to MaxCSP in both continuous and discrete domains while Section 7 is dedicated to experimental results showing the relevance of our approach. 2 A Motivating Example We motivate the usefulness of the numerical Max-CSP approach on a problem of virtual camera placement (VCP), which consists in positioning a camera in a 2D or a 3D environment w.r.t. constraints, namely cinematographic properties, defined on the objects of a virtual scene. For example, one wishes to see a character from a cinematographic point of view (e.g., profile, front, back, etc.) or from a typical shot distance (close-up, long shot, etc.). Within a classical CSP framework, these cinematographic properties are transformed into numerical constraints in order to find an appropriate camera position and orientation. One of the counterparts of the expressiveness offered by the constraint programming framework lies in the fact that it is easy for the user to specify an over-constrained problem. In such cases, classical solvers will output a no solution statement and does not provide any information neither on why the problem was over-constrained nor on how it could be softened. The Max-CSP framework will, on the other hand, explore the search space in order to compute the set of configurations that satisfy as many constraints as possible. This provides a classification of the solutions according to satisfied constraints and enables the user to select a particular solution, or to remodel the problem. For the purpose of illustration, we present a small over-constrained problem. The scene consists of three characters defined by their positions, and by an orientation vector (that denotes the front of the character). VCP frameworks generally provide the angle property that specifies the orientation the character should have on the screen (i.e., profile, etc.), and the shot property that spec- 2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 ifies the way a character is shot on the screen. The first properties obviously constrains both the location and the orientation of the camera, whereas the second constrains the distance to the character. In Figure 1, a sample scene is displayed with three characters A, B and C, each of which must be viewed from the front and with a long shot. The regions corresponding to each property are represented (wedges for the angle property and circles for the shot distance). No region of the space fully satisfies all the constraints since the intersection of the wedges and the circles is empty, meaning that the original problem is over-constrained. Conversely, the MaxCSP can be solved and its solutions are represented as dark areas. The solving process associated to this simple example is provided in section 7.2. For clarity sake, we need to emphasize the fact that Figure 1 represents a simplified version of the properties generally offered in VCP frameworks. Indeed, the shot constraint is represented here by a circle which does not take into account the actual shot properties. A shot property aims at constraining the size of an object on the screen. Classical cinematographic shot distances are taken into account by those properties, it is thus possible to constrain the size of an object (or a character) by using for example close-up, medium close-up (also called plan américain), or long shot distances. In order to achieve those kinds of shot, the camera has to be situated within a certain distance to the object involved. The correct modeling of such constraints should thus be done by placing a 2D ring (for a top-view scheme such as Figure 1 or by a 3D torus in 3D environments) around the object. For additional information on VCP frameworks and the modeling of the corresponding properties into a constraint satisfaction problem, the reader is invited to refer to [8]. 3 The Max-CSP framework The current section is dedicated to the Max-CSP framework and its related concepts. According to [28, 11], a Max-CSP P is a triplet hV, D, Ci where V = {x1 , . . . , xn }, D = {D1 , . . . , Dn } and C = {c1 , . . . , cm } are respectively a set of variables, a set of domains assigned to each variable and a set of constraints. For clarity sake, we use vectorial notations where vectors are denoted by boldface characters: the vector x = (x1 , . . . , xn ) represents the variables and D = D1 × . . . × Dn is interpreted as the Cartesian product of the variables domains. To define solutions of a Max-CSP, a function satc is first defined as follows, for each constraint c ∈ C and for each n-dimensional x in D:  1 if c(x) is true satc (x) = (1) 0 otherwise Once the satisfaction function satc is defined for a constraint c, we can introduce a satisfaction function, satC , defined over a set of constraints C and ∀x ∈ D as follows: X satC (x) = satc (x) (2) c∈C 3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Figure 1: Top-view of a 3D scene representing the camera’s position search space for a VCP problem. Requested orientations of objects are represented by dotted wedges areas centered on each object of interest, shot distances are represented as circles. Light gray areas represent intersection of constraints while dark gray areas represent the max-solution regions of the over-constrained problem, i.e., where the maximum number of constraints are satisfied. Obviously, we have satC (x) ≤ #C since at most all the constraints of C are satisfied. In the sequel, the maximum number of satisfied constraints of a Max-CSP P is denoted by mP and formally defined as follows: mP := max satC (x). x∈D (3) In the context of Max-CSP solving, one is usually more interested in the maximum number mP of satisfied constraints than in the variables assignments that actually satisfy the mP constraints. However, for our purpose, the set of these Max-CSP solutions is of central importance, since we want to know the assignments. We call it the solution set of the Max-CSP, and denote it by MaxSol(P). Formally, the following set of solution vectors (or variables assignments in the discrete case): X  MaxSol(P) = x ∈ D : (4) satc (x) = mP . c∈C 4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 4 Interval Analysis for Numerical CSPs Numerical CSPs present the additional difficulty of dealing with continuous domains, not exactly representable in practice. As a consequence, it is impossible to enumerate all the values that can take variables in their domains. Interval analysis (cf. [26, 15]) is a key framework in this context: soundness and completeness properties of real values belonging to some given intervals can be ensured by using correctly rounded computations over floating point numbers. We now present basic concepts allowing us to ensure that an interval vector contains no solution (resp. only solutions) of an inequality constraint. 4.1 Intervals, interval vectors and interval arithmetic An interval [x] with floating-point bounds (floating points are numbers representable in practice on a computer, cf. [14]) is called a floating-point interval. It is of the form: [x] = [x, x] = {x ∈ R : x ≤ x ≤ x} where x, x ∈ F. Given a floating-point value x, we denote by x⊕ (resp. x⊖ ) the smallest floating point value strictly greater than x (respectively, the greatest floating point value smaller than x). An interval vector (also called a box) [x] = ([x1 ], . . . , [xn ]) represents the Cartesian product of the n intervals [xi ]. The set of interval vectors with representable bounds is denoted by IFn . Operations and functions over reals are also replaced by an interval extension enjoying the containment property (see [24]). Definition 1 (Interval extension [26]) Given a real function f : Rn → R, an interval extension [f ] : IFn → IF of f is an interval function verifying [f ]([x]) ⊇ {f (x) : x ∈ [x]}. (5) Interval extensions of elementary functions can be computed formally thanks to their simplicity, e.g., for f (x, y) = x + y we have [f ]([x], [y]) = [x + y, x + y]. To obtain an efficient framework, these elementary interval extensions are used to define an interval arithmetic. For example, the interval extension of f (x, y) = x + y is used to define the interval addition: [x]+[y] = [x+y, x+y]. Therefrom, an arithmetical expression can be evaluated for interval arguments. The fundamental theorem of interval analysis [24] states that the interval evaluation of an expression gives rise to an interval extension of the corresponding real function. For example: [x] + exp([x] + [y]) ⊇ {x + exp(x + y) : x ∈ [x], y ∈ [y]}. 4.2 (6) Interval Contractors Discarding all inconsistent values from a box composed of variables domains in a CSP is an intractable problem when the constraints are real ones. There exist many techniques to obtain a contraction, like 2B consistency (also called hull 5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Figure 2: Diagrams (a) shows the constraint c (satisfied in dark gray regions) in an original box [x] that represents the search space. Diagram (b) shows how applying an interval contracting operator over the constraint c and the box [x] gives rise to a new box [y] in which as many inconsistent values of [x] as possible have been removed. consistency)[20, 4], box consistency[3] or interval Newton operator[26]. Those methods rely on outer contracting operators that discard values according to a given consistency. These methods are instances of the local consistency framework in artificial intelligence [21]. Definition 2 (Interval contracting operator) Let c be an n-ary constraint, a contractor for c is a function Contractc : IFn −→ IFn which satisfies x ∈ [x] ∧ c(x) =⇒ x ∈ Contractc ([x]). (7) Figure 2 illustrates the application of an interval contracting operator in a given search space [x] over a constraint c. The contractor aims at discarding from the original box [x] as many inconsistent values as possible according to a given consistency. This process gives rise to a new box [y] that encompasses the constraint c as tightly as possible. Experiments of Section 7 are performed with two different contractors dedicated to inequality constraints f (x) ≤ 0. The first is called the evaluation contractor and is defined as follows: Given an interval extension [f ] of a real function f  ∅ if [f ]([x]) > 0 Contractf (x)≤0 (x) = (8) [x] otherwise. The evaluation contractor is thus very simple while evaluating a box [x], either every value in [x] are consistent and the box is kept, otherwise, [x] is rejected and a null box is returned. The second, called the hull contractor, is based on the hull consistency, cf. [20, 4] for details1 . Given an interval extension [f ] of a real function f , the hull 1 A contractor based on box consistency [4] could also be used but the optimality between hull and box consistency is problem dependent. 6 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 contractor is applied as follows: Contractf (x)≤0 (x) = [y] (9) where [y] is the smallest box encompassing the solutions of the real function f , in practice, only a polynomial approximation of the hull consistency is enforced using HC3 or HC4 like algorithms (cf. [20, 1] for more details). Interval contractors can be used to partition a box into smaller boxes where the constraint can be decided. To this end, we introduce the inflated set difference in the following subsection. 4.3 Inflated Set Difference When an interval contractor is applied in a search space on a constraint, it outputs a tighter box which encompasses the constraint more precisely. As for our branch and bound algorithm for numerical Max-CSP, we want to apply in sequence such a procedure on each constraint c ∈ C of the problem in order to get a set of boxes in which the constraints are satisfied. In order to compute the regions of the search space that maximize constraint satisfaction for a numerical Max-CSP P, we want to compute the inflated set difference between every box obtained by inner and outer contractions over a constraint c. This is achieved by introducing the operator idiff([x], [y]) which computes a set of boxes {[d1 ], . . . , [dk ]} included in [x] such that [di ] ∩ [y] = ∅ and [d1 ] ∪ · · · ∪ [dk ] ∪ [y]⊕ = [x], where [y]⊕ is the smallest box that strictly contains [y]. The inflated set difference is illustrated in two dimensions for clarity sake, in Figure 3 located on page 11. When combining the results obtained by interval contraction over a constraint and the inflated set difference operator, we can compute a set of boxes that does not contain any solution of the constraint involved in the outer contraction process. Indeed, we can observe thanks to diagrams (a) and (c) of Figure 3 that the computation of: idiff([x], Contractc ([x])) gives rise to a set of boxes ⊗ that do not contain any potential solution to the constraint c. Moreover, following [9, 2] this technique can also be used to compute boxes that contain only solutions, applying the contractor to the negation of the constraint: idiff([x], Contract¬c ([x])) outputs a set of boxes that contain only solutions. Indeed, applying an outer contractor operator on the negation of a constraint encompasses the region of the search space that cannot satisfy the constraint. Thus, the set difference between the search space [x] and the box [y] (that does not contain any solution to the constraint c) gives rise to the set of boxes ⊗ that contains only solutions to the constraint c, (cf. diagrams (b) and (c) of Figure 3). Let us illustrate this process over a simple example. Consider the CSP P = hx, [−1, 1], x ≤ 0i, being composed of only one variable x that takes its 7 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 values in the floating-point interval [−1, 1] and that involves a single inequality: x ≤ 0. The first process is to apply an outer contracting operator on the constraint in order to remove as many inconsistent values as possible from the original search space ([−1, 1]). Then, we observe that: Contractx≤0 ([−1, 1]) = [−1, 0] meaning that the interval [−1, 0] contains all the solutions of the constraint x ≤ 0 within the domain [−1, 1]. What we want to do now is to compute the interval that does not contain any solution to the CSP P. This is achieved by computing a set difference (denoted by the operator \) between the original search space [−1, 1] and the result of the outer contraction [−1, 0]: [−1, 1]\[−1, 0] = ]0, 1]. This set difference computation outputs an semi-open interval ]0, 1] that is not easily dealt with in the context of interval analysis. In order to bypass this limitation, we introduce the operator idiff that computes the inflated set difference between two intervals. This operator is thus able to compute a set difference without using any open intervals. Indeed, if we compute the inflated set difference on the intervals used in our previous example, the result is the following: idiff([−1, 1], [−1, 0]) = [0⊕ , 1]. where 0⊕ is the smallest representable number greater than 0. As a consequence, the constraint x ≤ 0 is proved to be false inside [0⊕ , 1]. Nonetheless, in order to provide a full covering of the initial domains of the search space, the interval [−1, 0⊕] is kept for further processing instead of [−1, 0]. To sum up, the outer contraction has proved that the interval [−1, 0] contains all the solutions (but not only solutions) to the CSP P because the outer contracting step has removed as many inconsistent values as possible. As a consequence and thanks to the inflated set difference operator idiff, we also proved that the interval [0⊕ , 1] does not contain any solution of P. Finally, the domain [−1, 0⊕ ] remains to be processed, and we look for solutions therein by computing Contractx>0 ([−1, 0⊕ ]) = [0, 0⊕ ]. This time, the outer contracting operator is applied to the negation of the only constraint of P. We use again idiff([−1, 0⊕ ], [0, 0⊕ ]) = [−1, 0⊖ ] (where 0⊖ is the largest representable number lower than 0) to infer that [−1, 0⊖] contains only solutions. Finally, the constraint cannot be decided for the values of the interval [0⊖ , 0⊕ ], it is called a boundary box. 5 The Algorithm The main idea of the algorithm is to attach to each box, computed by a tree search, the set of constraints that were proved to be satisfied for all vectors of 8 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 this box (S) and the set of constraints that have not been decided yet, i.e., unknown constraints (U). The pseudo-code of our branch and bound algorithm for numerical Max-CSP is presented in Algorithm 2 and detailed in the last subsection. Before presenting the concepts and procedures required by our algorithm, we draw a global sketch of the algorithm. First, for each constraint c of the problem, the algorithm computes an outer contraction of the current box [x] (i.e., the current search space). This step removes as much inconsistent values as possible for constraint c in [x]. In order to ensure the exhaustive exploration of the search space, and thus to enforce the completeness property of the algorithm, a set difference operator is applied in a second step between the original search space [x] and the result of the outer contraction [y], giving rise to a set of boxes {[zi ]}. From this set, only the boxes that can satisfy more constraints than the current lower bound (m)of the algorithm are kept for further study. In a third step, an inner contraction is applied for the constraint c in the original search space [x]. This is achieved by applying a classical outer contraction procedure on the negation ¬c of the constraint. This step outputs a box [w] that removes all inconsistent values for ¬c, as shown in Figure 3. As a consequence the set difference between boxes [x] and [w] (i.e., [x]\[w]) only contain values consistent with constraint c. However, since the original search space [x] has already given rise to the set of boxes {[zi ]} in the outer contraction step, in order to speed up the solving process, we do not compute the set difference between [x] and [w] but between {[zi ]} and [w], as shown in Figure 4. This gives rise to a new set of boxes {[ui ]} that allows us to update the lower bound of the algorithm as well as to further reject uninteresting boxes. Finally, we iterate on the previous steps until no more new boxes are created or all the constraints have been decided. In order to converge, bisection steps might be applied in order to obtain smaller boxes where more constraints can be decided. The remainder of this section is constructed as follows: first we propose the definition of a SU-box, which represents a triplet consisting of a box [x] and the two sets of constraints, S and U, described above as the sets of satisfied and unknown constraints for [x]. Then, we present the specification and usage of outer and inner contractions in the following three subsections. Finally, as stated above, we describe the core algorithm in the last subsection. 5.1 SU-boxes The main idea of our algorithm is to associate to a box [x] the sets of constraints S and U which respectively contain the constraints proved to be satisfied for all vectors in [x] and the constraints which remain to be treated. Such a box, altogether with the satisfaction information of each constraint is called a SUbox 2 and is denoted by h[x], S, Ui. The following definition naturally links a SU-box with a Max-CSP P: 2 SU-boxes have been introduced in [27] in a slightly different but equivalent way. 9 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Definition 3 Let P be a CSP whose set of constraints is C. Then, a SU-box h[x], S, Ui is consistent with P if and only if the three following conditions hold: 1. (S ∪ U) ⊆ C and (S ∩ U) = ∅ 2. ∀x ∈ [x] , ∀c ∈ S , c(x) (the constraints of S are true everywhere in [x]) 3. ∀x ∈ [x] , ∀c ∈ (C\(S ∪U)) , ¬ c(x) (the constraints of C which are neither in S nor in U are false everywhere in [x]). A SU-box h[x], S, Ui will be denoted by hxi, without any explicit definition when there is no confusion. Given a SU-box hxi := h[x], S, Ui, [x] is denoted by box(hxi) and is called the domain of the SU-box. The interest of the two sets of constraints S and U is to update the bounds of the branch and bound algorithm. Indeed, the cardinality of S (denoted by #S) represents the lower bound on the number of satisfied constraint by a SU-box hxi, while the cardinality of both S and U (i.e., #S + #U) will bound the maximum number of constraints that can be satisfied by hxi. 5.2 The Function InferSUouter Let us consider a SU-box h[x], S, Ui, a box [y] ⊆ [x] and a constraint c, where the box [y] is obtained by computing [y] = Contractc ([x]). Therefore, as described in Section 4.3, every box of idiff([x], [y]) does not contain any solution of c. This is illustrated by Figure 3. The SU-boxes created with the boxes from the inflated set difference can thus be updated by discarding c from their set of constraints U. This is formalized by the following definition: Definition 4 The function InferSUouter (h[x], S, Ui, [y], c) returns the following set of SU-boxes: {h[x′ ], S, U ′ i : [x′ ] ∈ idiff([x], [y])} ∪ {h[y]⊕ ∩ [x], S, Ui}, (10) where U ′ = U\{c}. Here, the box [x] is partitioned to idiff([x], [y]) and [y]⊕ ∩ [x]. Each of these boxes is used to form new SU-boxes, with a reduced set of constraints U ′ for the SU-boxes coming from idiff([x], [y]). This is formalized by the following proposition, which obviously holds: Proposition 1 Let hxi := h[x], S, Ui be a SU-box consistent with a CSP P, c ∈ U and [y] satisfying: ∀x ∈ idiff([x], [y]) ⇒ ¬ c(x). We define T := InferSUouter (hxi, [y], c). As a consequence, every SU-box of T is consistent with P. Furthermore, ∪ {box(hx′ i) : hx′ i ∈ T } is equal to [x]. 10 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Figure 3: Diagrams (a) and (b) show the constraint c (satisfied in dark gray regions and unsatisfied in light gray regions), the initial SU-box h[x], S, Ui and the contracted box [y] for two different situations: (a) outer contraction and (b) inner contraction for constraint c. Diagram (c): Five new SU-boxes are inferred from the contraction of [x] to [y]. The contraction does not provide any information on the SU-box ⊕ = h[x] ∩ [y], S, Ui, hence its sets of constraints remain unchanged. The SU-boxes ⊗ = h[x′ ], S ′ , U ′ i are proved to contain no solution of c (if outer pruning was applied) or only solutions (if inner pruning was applied), their sets of constraint can thus be updated: U ′ = U\{c}, and, S ′ = S for outer pruning, while S ′ = S ∪ {c} and U ′ = U for inner pruning. 5.3 The Function InferSUinner This function is equivalent to InferSUouter , but deals with inner contraction. As a consequence, the only difference is that, instead of discarding the constraint c from U for the SU-boxes outside [y], this constraint is stored in S. Figure 3 also illustrates these computations. The following definition and proposition mirror Definition 4 and Proposition 1. Definition 5 The function InferSUinner (h[x], S, Ui, [y], c) returns the following set of SU-boxes: {h[x′ ], S ′ , U ′ i : [x′ ] ∈ idiff([x], [y])} ∪ {h[y]⊕ ∩ [x], S, Ui}, (11) where S ′ = S ∪ {c} and U ′ = U\{c}. Proposition 2 Let hxi := h[x], S, Ui be a SU-box consistent with a CSP P, c ∈ U and [y] satisfying: ∀x ∈ idiff([x], [y]) ⇒ c(x). We define T := InferSUinner (hxi, [y], c). As a consequence, every SU-box of T is consistent with P. Furthermore, ∪ {box(hx′ i) : hx′ i ∈ T } is equal to [x]. Following principle of InferSUouter , the procedure InferSUinner is applied to a unique constraint. It computes a set difference between the current search space [x] and [y] which is obtained by inner contraction on ¬c. The set of resulting SUboxes {h[x′ ], S ′ , U ′ i have both their sets of Satisfied and U updated accordingly, i.e., the constraint c is added to the former while being removed from the latter. 11 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Figure 4: Left hand side illustrates the [xk ] boxes obtained from Figure 3. Moreover, the box [y] represents the contraction (inner or outer) of a box. Right hand side: The process illustrated by Figure 3 is repeated for each box [xk ], leading to twelve new SU-boxes, among which ⊗k have their sets of constraints updated. To sum up, merging information on constraints satisfaction obtained on the SU-boxes computed through the outer and inner contractions steps is done by applying a set difference. Indeed, when computing these intersections, we can update the sets S and U and thus maintain the bounds of the algorithm. Uninteresting SU-boxes will be automatically rejected when processed since they cannot satisfy at least m constraints (current lower bound of the Branch and Bound algorithm). 5.4 The Function InferLSUtype Finally, the functions InferSUouter and InferSUinner are naturally applied to a set of SU-boxes T in the following way: [ InferSUtype (hxi, [y], c), (12) InferLSUtype (T , [y], c) := hxi∈T where type is either outer or inner. This process is illustrated by Figure 4. As a direct consequence of propositions 1 and 2, the following proposition holds: Proposition 3 Let T be a set of SU-boxes, each of them being consistent with a CSP P, and [y] satisfying ∀x ∈ idiff([x], [y]) ⇒ ¬ c(x) (respectively ∀x ∈ idiff([x], [y]) ⇒ c(x)). We define T ′ := InferLSUouter (hxi, [y], c) (respectively T ′ := InferLSUinner (hxi, [y], c)). As a consequence, the SU-boxes of T ′ are all consistent with P. Furthermore, ∪ {box(hx′ i) : hx′ i ∈ T ′ } = ∪ {box(hx′ i) : hx′ i ∈ T }. 5.5 (13) The Function MultiStartRandomSearch Like every Branch and Bound method, our algorithm greatly benefits from having a good lower bound on the number of satisfied constraints of the problem. 12 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Indeed, the sooner one gets a good lower bound m, the more uninteresting boxes will be rejected during the search. An uninteresting box being a box that cannot satisfy at least m constraints. This remark applied to a SU-box meaning that a SU-box hxi can be rejected as soon as (#S + #U) < m, which means that hxi will never be able to satisfy m constraints. As a direct consequence, those regions of the search space can be rejected because they will not be part of the solution set MaxSol(P) of the numerical Max-CSP. The procedure presented in this section is very simple but will greatly improve the performance of our algorithm by allowing the computation of an early good lower bound for the algorithm. Indeed, MultiStartRandomSearch(C, D, n) is a function that inputs the set of constraints C of the numerical Max-CSP P, the search space D and a number n of samples that will be generated in order to approximate a first lower bound on the maximum number of satisfied constraints of P. Thus, the function outputs the maximum number of satisfied constraints from the n samples, giving an insight of a lower bound on the number of constraints that can be satisfied. Algorithm 1 gives the pseudo-code of the MultiStartRandomSearch function as implemented in our method. Each vector p is obtained thanks to a function GenerateRandomPoint(D) that generates a random value in each of the dimensions of the search space D. Let us recall here that the present algorithm only aims at computing the outer and inner solution sets of a numerical Max-CSP, the former set ensuring that no solution is lost during the solving process while the latter will be composed of boxes that contain only solutions of the Max-CSP. As a consequence, our algorithm requires the CSPs to be composed of numerical inequalities only. Algorithm 1: The MultiStartRandomSearch function computes a lower bound m for a set of constraints C in a search space D. Input: C, D, n Output: m 1 m ← 0; 2 for i ← 0 to n do 3 p ← GenerateRandomPoint(D); 4 b ← satC (p); 5 if ( m < b ) then m ← b; 6 end 7 return m; Remark 1 When dealing with inequality constraints, satC (p) is rigorously computed using the interval extension of the constraint, i.e., f (p) ≤ 0 ⇐=sup [f ]([p, p]) ≤ 0. Section 7.1 highlights the impact of this preprocessing step on the performance of the branch and bound algorithm. 13 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 5.6 The Branch and Bound Algorithm The pseudo-code of our technique is presented in Algorithm 2. Our algorithm computes three sets of SU-boxes L (SU-boxes to be processed), D (SU-boxes where all constraints are decided, i.e., where (#U) = 0) and E (SU-boxes that will not be processed anymore because considered as too small). The mainline of the while-loop aims at maintaining the following property: [ box(hxi) : hxi ∈ L ∪ D ∪ E , (14) MaxSol(P) ⊆ while using inner and outer contractions as well as the function InferLSUtype in order to decide more and more constraints (cf. lines 11–17), and hence moving SU-boxes from L to D (cf. Line 21). Note that a lower bound m for mP is updated (cf. Line 19) and used to drop SU-boxes that are proved to satisfy less constraints than m (cf. Line 14). When L is finally empty, (14) still holds and hence: MaxSol(P) ⊆ ∪{box(hxi) : hxi ∈ D ∪ E}. (15) Eventually, the sets of SU-boxes D and E are used to obtain new sets of SUboxes M and B, which describe the max-solution set (i.e., both the inner and outer approximation of the solution set of the numerical CSP), and an interval [m, m] for mP . In other words the sets M and B encompass the solution set of the numerical Max-CSP P. M is the set of SU-boxes that have been proved to contain only solutions (i.e., for which (#S) ≥ m)) while B ensures not to lose any solutions of the P by keeping boundary boxes that encompass the border of the solution set, namely SU-boxes such that (#S) + (#U) ≥ m. At this stage of the algorithm, we already have a lower bound m on mP . Now, as (14) holds, an upper bound can be computed in the following way: m = max{(#S) + (#U) : h[x], S, Ui ∈ D ∪ E}. (16) Two cases can then happen, which are handled at lines 29–36: 1. If m = m then the algorithm has succeeded in computing mP , which is equal to m = m. In this case the domains of the SU-boxes of D are proved to contain only max-solutions of P: ∀hxi ∈ D , box(hxi) ⊆ MaxSol(P). (17) Furthermore, the max-solutions which are not in D obviously have to be in the SU-boxes of: B = {h[x], S, Ui ∈ E : (#S) + (#U) ≥ m}. (18) As a consequence, the algorithm outputs both an inner and an outer approximation of MaxSol(P). 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Algorithm 2: Branch and Bound Algorithm for Numerical MaxCSP. Input: P = hx, C, [x]i, ε Output: ([m, m], M, B) 1 L ← {h[x], ∅, Ui}; /* SU-boxes to be processed 2 E ← {}; /* Epsilon SU-boxes 3 D ← {}; /* SU-boxes where each constraint is decided 3 4 m ← MultiStartLocalSearch(P, [x]) ; /* Lower bound on m 5 while ( ¬empty(L) ) do 6 h[x], S, Ui ← extract(L)4 ; 7 if ( wid ([x]) < ε ) then 8 E ← E ∪ {h[x], S, Ui}; 9 else 10 T ← {h[x], S, Ui}; 11 foreach ( c ∈ U ) do 12 [y] ← Contractc ([x]); 13 T ← InferLSUouter (T , [y], c); 14 T ← {h[x], S, Ui ∈ T : (#S) + (#U) ≥ m}; 15 [y] ← Contract¬c ([x]); 16 T ← InferLSUinner(T , [y], c); 17 end 18 foreach ( h[x′ ], S ′ , U ′ i ∈ T ) do 19 if ( (#S ′ ) > m ) then m ← (#S ′ ); 20 if ( (#U ′ ) = 0 ) then 21 D ← D ∪ {h[x]′ , S ′ , ∅i}; 22 else 23 {[x1 ], [x2 ]} ← Bisect([x′ ])5 ; 24 L ← L ∪ {h[x1 ], S ′ , U ′ i, h[x2 ], S ′ , U ′ i}; 25 end 26 end 27 end 28 end 29 m ← max{(#S) + (#U) : h[x], S, Ui ∈ D ∪ E}; 30 if ( m = m ) then 31 M ← {h[x], S, Ui ∈ D : (#S) ≥ m}; 32 B ← {h[x], S, Ui ∈ E : (#S) + (#U) ≥ m}; 33 else 34 M ← ∅; 35 B ← {h[x], S, Ui ∈ D ∪ E : (#S) + (#U) ≥ m}; 36 return ([m, m], M, B); */ */ */ */ 3 The constraints satisfaction is tested on uniformly distributed random points, the lower bound m being updated accordingly. 4 SU-boxes are stored in L thanks to (#S) + (#U ), most interesting first. 5 Bisect splits the box [x′ ] into 2 equally sized smaller boxes: [x1 ] and [x2 ]. 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 2. If m < m, mP is still proved to belong to the interval [m, m]. However, the algorithm has not been able to compute any solution of the Max-CSP P, hence M = ∅. In this case, all the max-solutions are obviously in the domains of the SU-boxes of: B = {h[x], S, Ui ∈ D ∪ E : (#S) + (#U) ≥ m}. (19) As a consequence of Proposition 3, constraints are correctly removed from U and possibly added to S while no solution is lost. This ensures the correctness of the algorithm (a formal proof is not presented here due to its burdensomeness). 5.7 Discussion on the Complexity and on the Success of the Algorithm The worst case complexity of the foreach-loop at lines 11–17 is exponential with respect to the number of constraints in U. Indeed, the worth case complexity of the loop running through lines 11–17, is (2n)q , where n is the dimension and q is the number of constraints. This upper bound is very pessimistic, and not met in any typical situation. Nonetheless, we expect a possibly high number of boxes generated during the inflated set difference process when the number of constraints is important. In order to bypass this complexity issue, that could lead the algorithm to generate too many boxes within the outer and inner contractions steps, we have chosen to switch between different contracting operators. As soon as the number of constraints becomes greater than a threshold (tc ), contractions are performed with the evaluation contractor instead of the hull contractor. This will reduce creation of boxes within the inner and outer contracting steps, thus preventing the algorithm from falling into a complexity pit. Experiments on relatively high number of constraints (cf. Section 7.3) have shown that tc should be set between 5 and 10 in order to speed up the solving process. The question of the success of the algorithm, i.e., m = m, naturally arises. Even if ε is taken arbitrarily small, there are situations where the algorithm does not succeed. There are typically two cases where this will happen: 1. First, if the max-solution set has no interior, e.g. P = x, {f (x) ≤ 0, f (x) ≥ 0}, Rn , with f (x) = x1 2 + x2 2 − 1; for which mP = 2 and MaxSol(P) = {x ∈ Rn : f (x) = 0}. 2. Second, it can happen that even when the solution set has a non empty interior, interval contractors are unable to prove the existence of solutions, e.g. P = x, {|x| + x ≤ 0}, R for which mP = 1 and MaxSol(P) = Sol(P) = (−∞, 0 ]. Both cases are considered as untypical in the sense that any arbitrarily small perturbation of the constraints make the max-solution set empty. Very different algorithms, like formal techniques, have to be used to solve these untypical situations. 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 6 Related Work A large literature is related to Max-CSP frameworks and solving techniques. Most formulate the problem as a combinatorial optimization one by minimizing the number of violated constraints. Approaches are then shared between complete and incomplete techniques depending on the optimality criterion and the size of the problem. For large problems, heuristic-based techniques such as Min-conflict [23] and Random-walk [29] have proved to be efficient as well as meta-heuristic approaches (Tabu, GRASP, ACO). However all are restricted to the study of binary or discrete domains, and their extension to continuous values requires to redefine the primitive operations (e.g., move, evaluate, compare, etc.). Though some contributions have explored such extensions for local optimization problems (cf. for example Continuous-GRASP[16]), these are not dedicated to manage numerical Max-CSPs. More general techniques built upon classical optimization schemes such as Genetic Algorithms, Simulated Annealing or Generalized Gradient can be employed to express and solve Max-CSPs, but do not guarantee the optimality of the solution, nor compute a paving of the search space. Regardless the fact that applications of CSP techniques might be either in continuous or in discrete domains, both present the same need for MaxCSP models. Interestingly very few contributions have explored this class of problems when optimality is required. Here, we report the contribution of Jaulin et al. [18, 17] that looks at a problem of parameter estimation with bounded error measurements. The paper considers the number of error measurements as known (i.e., mP is given beforehand) and the technique computes a MaxCSP approximation over a continuous search space by following an incremental process that runs the solving for each possible value of mP and decides which is the correct one. At each step the technique relies on a classical evaluationbisection process. 7 Experiments In order to illustrate the pertinence of our approach, we propose a set of three different benchmarks, ranging from a virtual camera placement (VCP) problem, as described in Section 2 to parameter estimation with bounded error measurements. Furthermore, we wanted to emphasize the flexibility offered by our method, by applying a slightly modified version of our algorithm to a classical Operational Research problem, namely the facility location problem. Before presenting the results we have obtained, we first detail the usefulness of the simple function MultiStartRandomSearch that aims at giving a first good lower bound on the number of satisfied constraints in a numerical Max-CSP, thus boosting the performance of our method. 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 7.1 Impact of the MultiStartRandomSearch function Before presenting the benchmarks of our algorithm, and in order to emphasize the benefits of the function MultiStartRandomSearch presented in Section 5.5, we have realized some experiments on a simple toy problem. The problem consists in generating random 3D balls in a search space D = [−20, 20]× [−20, 20]× [−20, 20] and trying to find the regions of the search space that maximize ball intersections. Table 1 presents the results obtained when the MultiStartRandomSearch was used in order to get as soon as possible a good lower bound, and second when the first lower bound was set to 0 and updated classically by the algorithm (cf. Line 19 of the Algorithm 2). Updating the lower bound will allow rejecting as many uninteresting boxes as possible during the search, thus increasing the performance of the algorithm. Different values of random points have been chosen in order to evaluate the impact on the global performance of the algorithm. The timings are represented in seconds as well as the maximum number of satisfied constraints for each problem. The number of balls generated in D is denoted after the name of the benchmark, e.g., Balls3D-150 represents a benchmark where 150 random 3D balls have been generated in D. The number of random sample points generated is represented in the #smp column. If this number is set to 0, then the function is not called and the classical Branch and Bound algorithm is applied without preliminary searching for a good lower bound. Finally, the lower bound found by the MultiStartRandomSearch method is given in the mMSRS column, while the interval [m, m] of satisfied constraints computed by the algorithm is represented in the eponym column. Timings are listed in the last two columns, total and time spent in MultiStartRandomSearch (TMSRS ) are represented in seconds. Results (cf. Table 1) show significant improvements when using the function MultiStartRandomSearch. Indeed while the original algorithm fails at solving a number of benchmarks (T.O. stands for time out, which results from a saturation of the memory of the computer), one can observe that the preprocessing step greatly improves the results. Of course, the number of sample points (#smp) that should be generated in order to get a good lower bound is problem dependent. Results only aim at emphasizing the usefulness of such a preprocessing step. As a consequence, while the first two problems are solved with our algorithm (i.e., m = m), the benchmark Random-Balls3D-500 should be tackled with a lower minimum bisection size ε in order to solve it. Indeed Table 1 shows that our algorithm has only computed an outer approximation of the solution set of the Random-Balls3D-500 benchmark. The simple procedure presented here (evaluation of the constraints on random points in the search space) could obviously be improved by adding a local search procedure. However, such random search for Max-CSP in continuous domains is, to the best authors knowledge, yet to be proposed. The authors are aware that generating random points in the search space is certainly not the best heuristic when tackling continuous problems, nevertheless this straightforward improvement has shown interesting results. As future work, 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Benchmark Balls3D-50 Balls3D-50 Balls3D-50 Balls3D-50 Balls3D-150 Balls3D-150 Balls3D-150 Balls3D-150 Balls3D-500 Balls3D-500 Balls3D-500 Balls3D-500 Balls3D-500 Balls3D-500 ε 0.5 0.5 0.1 0.1 0.5 0.5 0.1 0.1 0.5 0.5 0.5 0.1 0.1 0.1 #smp 0 50 0 50 0 50 0 50 0 50 200 0 50 200 mMSRS 0 29 0 29 0 85 0 85 0 279 283 0 279 283 [m, m] [30, 31] [30, 31] [0, 50] [31, 31] [85, 87] [85, 87] [0, 150] [86, 86] [0, 500] [280, 291] [283, 291] [0, 500] [284, 288] [284, 288] TMSRS (s) 0 0.007 0 0.006 0 0.016 0 0.012 0 0.047 0.2 0 0.045 0.167 Ttotal (s) 256 1.297 T.O. 11.272 1412.2 1.219 T.O. 2.478 T.O. 7.719 4.053 T.O. 22.522 8.551 Table 1: Impact of the MultiStartRandomSearch preprocessing step on the performance of our branch and bound algorithm for numerical Max-CSPs. Random 3D balls are generated in D = [−20, 20] × [−20, 20] × [−20, 20] with a minimum bisection size set to ε. Timings obtained on an IntelCore 2 Duo 2.4GHz laptop 4GB RAM. the authors will consider replacing this procedure by state of the art continuous local optimizers such as MINOS [25] or KNITROS [6] for example. Those solvers will certainly prove to be much efficient than our simple procedure, however it could be interesting to study the impact of those solvers on the global performance of our algorithm. The overhead due to their use could reduce the benefits provided by the MultiStartRandomSearch procedure. Indeed, the aim of this procedure is not to find an optimum value for the number of satisfied constraints, but rather to be able to find a relatively good bound on this number in order to reject as many boxes as possible during the exhaustive search performed by our algorithm. 7.2 Camera Control Problem This subsection is dedicated to the virtual camera placement (VCP) example presented in Section 2. Let us recall briefly that this motivating example involves 2 variables and 6 constraints expressed as dot products and Euclidean distances in 2D. The constraints involved in the numerical Max-CSP P are defined over three objects A, B and C abstracted as 2D points representing (xA , yA , xB , yB , xC , yC ) their positions in the scene and three 2D vectors (xOA , yOA , xOB , yOB ,etc.) representing their orientations. The variables of P are the camera 2D position, namely x and y. Finally, the acceptable distances are denoted by rA , rB and rC while accepted orientations are represented by 19 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Figure 5: Pavings generated by our algorithm on the VCP problem. Grey boxes represent the max solution set (cf. Figure 1), while black boxes illustrate boundary boxes encompassing the solution set. the values of the dot products αA , αB , αC . The constraints of P are thus the following: p c1 : (x − xA )2 + (y − yA )2 ≤ rA p (x − xB )2 + (y − yB )2 ≤ rB c2 : p (x − xC )2 + (y − yC )2 ≤ rC c3 : c4 : (x−xA )∗xOA +(y−yA )∗yOA c5 : (x−xB )∗xOB +(y−yB )∗yOB √ (x−xB )2 +(y−yB )2 ) ≥ αB c6 : (x−xC )∗xOC +(y−yC )∗yOC √ (x−xC )2 +(y−yC )2 ) ≥ αC √ (x−xA )2 +(y−yA )2 ) ≥ αA The corresponding paving is shown in Figure 5 while results are presented in Table 2 within the search space {[−30; 30], [1; 30]} and with ε = 0.01. The algorithm outputs m = m = 5 meaning that the set of inner boxes satisfy five constraints out of six, let us recall that the problem is over-constrained as shown in Figure 1 located on page 1. Although Christie et al. [7] do not compute an outer and and inner approximation of the solution set of the Max-CSP (they only aim at computing 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Method Numerical Max-CSP Christie et al. ε 0.01 0.01 Toriginal (seconds) 3.2 4.8 Tapprox (seconds) 3.2 ≈ 15 Table 2: Approximated timings thanks to [12] for the VCP example between our approach (PentiumM750 laptop, 1GB RAM) and Christie et al. [7](Pentium T7600 2.33Ghz 2GB RAM). an inner approximation), their approach is similar enough to compare timings obtained on the VCP example. As discussed in Section 5.7, this example illustrates the powerfulness of the gathering of the solution sets in connected components. Indeed Figure 5 shows that the search space is divided into three connected components each encompassing a different set of equivalent solutions. Presenting the user such useful information on the solution sets might help him improving the modeling of his problem or be able to chose solutions that fulfill specific constraints interesting for the user. 7.3 The facility location problem This example is based on a well-known geometric problem in Operational Research: the facility location problem [13]. The problem is defined by a number of customer locations that must be delivered by some facilities. We consider that a customer c can be served by a facility f if the distance between f and c is smaller than a delivery distance d. The locations of the customers are inputs of the facility location problem, which thus consists in positioning the smallest number of facilities such that each customer can be served by at least one facility. Strictly speaking, the facility location problem in Operational Research is not exactly posed as this benchmark. Indeed, the original problem is defined as a combinatorial one, that is knowing a set of facilities locations and a set of customers locations, which facilities should be ”opened” in order that at least every customer is delivered by a facility. Nevertheless, the authors believe that this slightly modified version of the facility location problem is realistic enough to be used in our set of experiments. Although the facility location problem is not strictly speaking a Max-CSP, we can use the Max-CSP framework to build a greedy algorithm to find an estimation of the number of facilities. The idea consists in computing the MaxCSP areas of the problem which corresponds to the areas of the search space that maximizes customers delivery. Instead of computing the whole MaxSol paving, we stop the algorithm as soon as we have found a MaxSol SU-box h[xm ], Sm , Um i. We then remove from the original problem all the constraints of Sm (which corresponds to all the customer that have been delivered during the iteration) and re-run the algorithm again. We stop when the problem is 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Method Greedy Numerical Max-CSP Basic OR technique Time (seconds) 69.25 ≈ 140 Number of facilities 39 45 Table 3: Timings for the facility location problem for 500 customers randomly generated in {[−5; 5], [−5; 5], [−5; 5]} on a PentiumM750 laptop, 1GB RAM. empty, i.e., when all the constraints have been solved. The numbers of runs of the algorithm then correspond to a number of facilities that solves the problem. Of course this is a quite naive approach but it is very easy to implement thanks to our algorithm. Our greedy implementation might not for now be as effective as OR approaches but it could be improved in order to reduce the number of runs needed to solve the problem. In order to test our greedy algorithm, we have implemented a basic OR technique to solve the same problem (i.e., discretization of the search space and resolution of a linear integer programming problem) in Wolfram Mathematica6 . Indeed, classical OR algorithms compute the intersection areas of all deliverable areas around the customers (e.g., circles in 2D) before using a linear programming algorithm to solve the problem. We believe that the higher the dimensions of the problem, the more difficult the intersection computation will be and the more our algorithm could compete the OR approach. The classical facility location problem is, to the best of the authors knowledge, only defined in 2D in Operational Research. Nonetheless, in order to challenge our hypothesis that the greater the number of dimensions, the better our algorithm will perform, we propose a 3D version of the facility location benchmark. The definition of the problem remains unchanged, only instead of trying to position the facilities in a 2D space, we propose to try to find their coordinates in a 3D space. In 3D this problem boils down to modeling each customer by a random 3D position in the search space {[−5; 5], [−5; 5], [−5; 5]} and trying to position some facilities such that the distance between the facility and the customer is inferior to d = 2. Results are shown in Table 3. The solution found by our algorithm involves 3 × 39 variables (3 coordinates per computed facility). Each iteration of our algorithm removes on average 13 constraints (ranging from 33 to 1 constraint removed at each iteration) of the 500 original ones. Although, this basic OR technique seems less efficient than our greedy algorithm, this is at least encouraging and thorough experiments need to be carried out to compare our greedy approach with more specific OR techniques. 6 The corresponding Mathematica notebook is available for download at: http://www.goldsztejn.com/src/FacilityLocation3D.nb. 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 7.4 Parameter Estimation with Bounded Error Measurements This experiment is based on a parameter estimation problem presented in [17]. The problem consists in determining unknown parameters p = (p1 , p2 ) ∈ ([−0.1, 1.5], [−0.1, 1.5]) of a physical model thanks to a set of experimental data obtained from a system. The physical model is defined as follows and has been taken from [18, 17, 22]. The vector of original available data y is: y = (7.39, 4.09, 1.74, 0.097, −2.57, −2.71, −2.07, −1.44, −0.98, −0.66)T It corresponds to ten scalar values taken at times: t = (0.75, 1.5, 2.25, 3, 6, 9, 13, 17, 21, 25)T But here in order to simulate outliers, Jaulinet al. [17] arbitrarily replaced two data points by zero values and thus the array of data becomes: y = (7.39, 0, 1.74, 0.097, −2.57, −2.71, −2.07, 0, −0.98, −0.66)T We assume that these data are to be described by a vector ym ∈ Rny described by a model with a fixed structure but unknown parameter vector p ∈ Rnp . The purpose of parameter estimation, is to find p such that ym (p) fits y best. In this experiment, we are talking of bounded-error estimation, thus we consider that the parameters p are considered admissible if the error e(p), defined as: e(p) = y − ym (p) belongs to some set of admissible errors. Here we consider the physical system being defined (cf. [17]) such that the ith component of ym (p) is given by: ym,i (p) = 20 exp(−p1 ti ) − 8 exp(−p2 ti ) The set of all feasible model outputs is the box defined by: y = [y − emax , y + emax ] with emax = (4.695, 1., 1.87, 1.0485, 2.285, 2.355, 2.035, 1., 1.49, 1.33)T Here the model used is a two-parameter estimation problem taken from Jaulinet al. [18, 17] which is a 2D extension of a problem presented by Milanese and Vicino [22]. In fact, the solution set of the Max-CSP does not contain the true parameters values. Indeed, the maximum number of satisfied constraints is 9, while there are 2 outliers and the model contains 10 constraints (corresponding to the vectors given above). This means that an outlier is compatible 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Method Numerical Max-CSP Jaulin et al. ε bisection 0.005 0.005 Toriginal (seconds) 0.46 ≈ 180 Tapprox (seconds) 0.46 ≈ 1.9 Table 4: Original and approximate timings obtained from a Pentium M750 laptop processor (Numerical Max-CSP) and from a 486 DX4-100 processor (Jaulin et al.) thanks to [12]. with the error measurements. Hence, the Max-CSP cannot directly help solving this problem. From an application point of view, it is possible to provide an upper bound on the number of outliers (this assumption is realistic since manufacturers are able to ensure a percentage of maximum failures of their probes). Algorithm 2 is easily modified to output the set of vectors which satisfy at least a fixed number of constraints instead of the solution set of the Max-CSP. In order to compare our approach with [17], we set a maximum of 3 outliers, Figure 6 shows the corresponding paving. Although Jaulin et al. do not compute an inner and outer approximation of the solution set, Table 4 compares our results. Approximate timings have been computed using [12] in order to compare the results that were obtained on totally different computers. Time comparisons are presented to give an idea of the gain offered by our algorithm, which is approximately 4 times more efficient (although no exact comparison of the pavings computed by the two algorithms is possible, the pavings presented in [17] are clearly less accurate than the pavings computed by our approach). Moreover, we want to emphasize the fact that our algorithm is much more scalable since Jaulin et al. have to re-run their algorithm for each possible number of outliers whereas we only need one run to compute the solution set of the Max-CSP. Moreover our approach fully characterize this solution set w.r.t. the satisfied constraints of the algorithm. 8 Conclusion In this paper we have presented a branch and bound algorithm for numerical Max-CSP that computes both an inner and an outer approximation of the solution set of a numerical Max-CSP. The algorithm is based on the notion of SU-box that stores the sets of constraints Satisfied or Unknown for each box of the search space. The algorithm outputs a set of boxes that encompass the inner and outer approximation of the solution set of a numerical Max-CSP. For each box, it computes the sets of constraints that are satisfied, thus giving the user additional information on the classes of solutions of the problem. An algorithm of connected component identification could then be applied on the solution SU-boxes in order to create regions of the search space that share similar characteristics (i.e., that satisfies the same constraints). Users could thus benefit from this feature when modeling conceptual design problems for example. In continuous domains, a CSP is generally defined as over-constrained when 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Figure 6: Paving generated for 3 possible outliers. Grey boxes represent inner approximation (i.e., satisfying at least 7 of the 10 constraints) while black boxes illustrate boundary boxes encompassing the solution set. no solution can be found for it. Conversely, a CSP is said to be under-constrained when the number of solutions to the constraint set is infinite. Finally, a CSP is well constrained when the number of solutions of the constraint set is finite. The Branch and Bound algorithm presented in this paper has been designed to solve over-constrained problem. Nonetheless, our method can address any kind of CSP: well-constrained, under-constrained and over-constrained. Indeed, it is often a problem in itself to determine whether a CSP is well-constrained, over-constrained or under-constrained, this is the reason why our algorithm has been designed to be able to solve any kind of CSP, of course at the price of some overhead induced by additional computations. Indeed, if the CSP is well-constrained, the only difference for our B&B algorithm is that the maximum number of satisfied constraints is equal to the number of constraints of the CSP. The whole solving process remains the same. As for the under-constrained case, since we never reject boxes as soon as we are certain that a box is better (i.e., that is satisfies at least one more constraint), we are able to output the set of solutions to the CSP. A drawback of our algorithm is the possible combinatorial explosion created within the inner and outer contractions steps. Indeed a set difference is computed between all the boxes created by inner and outer contractions applied to 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 each constraint of the problem. This leads to poor performance when the number of constraints of the problem becomes too important. In order to bypass this limitation, we have modified our algorithm in order to perform the inner and outer contractions steps only when the number of constraints is not too high. This is achieved by testing the number of constraints contained within the studied SU-boxes. If this number is higher than an arbitrary threshold (tc ) we choose to use a simple box evaluation procedure instead of inner and outer contractors. A problem naturally arises on how to define this threshold, experiments have shown that for a problem consisting of 100 constraints this threshold could be set around 5-10 in order to speed the solving process. Our algorithm could be extended to manage hierarchical constraints in a predicate locally-better way [5]. Hierarchical constraints could be modeled as different “layers” of numerical Max-CSP problems that could be solved in sequence, starting from the top-priority constraints down to the lowest ones and maximizing constraints satisfaction of each “layer” of hierarchical constraints. Finally, an advantage of our algorithm is that each SU-box contains the set of constraints that are fulfilled within this box. We are thus able to give the user information on which constraint is fulfilled in each part of the search space, allowing him to choose between the different solutions of the over-constrained problem. An automatic connected component identification algorithm (cf. [10]) could be processed after solving a problem with our algorithm in order to present the user the sets of equivalent SU-boxes, allowing him to choose which constraints he wants to relax into the original problem. Another perspective we are currently working on lies in the use of the qintersection algorithm introduced by Jaulin and Walter [19], which is able to compute a set of boxes (Cartesian products of domains) that lies in the intersection of q boxes. Acknowledgments The authors are grateful to the anonymous referees whose comments significantly improved the paper. References [1] F. Benhamou, F. Goualard, L. Granvilliers, and J.-F. Puget. Revising hull and box consistency. In Proceedings of the 1999 International Conference on Logic Programming, ICLP’99, pages 230–244, Cambridge, MA, USA, 1999. Massachusetts Institute of Technology. [2] F. Benhamou, F. Goualard, E. Languenou, and M. Christie. Interval Constraint Solving for Camera Control and Motion Planning. ACM Transactions on Computational Logic, 5(4):732–767, 2004. 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [3] F. Benhamou, D. McAllester, and P. van Hentenryck. CLP(intervals) revisited. In International Logic Programming Symposium, ILPS ’94, pages 124–138, Cambridge, MA, USA, 1994. MIT Press. [4] F. Benhamou and William J. Older. Applying Interval Arithmetic to Real, Integer and Boolean Constraints. Journal of Logic Programming, 32(1):1– 24, 1997. [5] A. Borning, B. Freeman-Benson, and M. Wilson. Constraint Hierarchies. Lisp and Symbolic Computation, 5(3):223–270, 1992. [6] R. H. Byrd, J. Nocedal, and R. A. Waltz. KNITRO: An Integrated Package for Nonlinear Optimization. In Large-Scale Nonlinear Optimization, pages 35–59. Springer Verlag, 2006. [7] M. Christie, J.-M. Normand, and C. Truchet. Computing inner approximations of numerical MaxCSP. In Interval analysis, constraint propagation, applications, IntCP 2006, 2006. [8] M. Christie and J-M. Normand. A Semantic Space Partitionning Approach to Virtual Camera Control. In In Proceedings of the Annual Eurographics Conference, volume 24, pages 247–256, 2005. [9] H. Collavizza, F. Delobel, and M. Rueher. Extending Consistent Domains of Numeric CSP. In IJCAI ’99: Proceedings of the Sixteenth International Joint Conference on Artificial Intelligence, pages 406–413, 1999. [10] N. Delanoue, L. Jaulin, and B. Cottenceau. Counting the Number of Connected Components of a Set and Its Application to Robotics. In PARA, pages 93–101, 2004. [11] S. de Givry, J. Larrosa, P. Meseguer, and T. Schiex. Solving Max-SAT as Weighted CSP. In Principles and Practice of Constraint Programming, CP 2003, pages 363–376. Springer-Verlag, 2003. [12] J. Dongarra. Performance of Various Computers Using Standard Linear Equations Software. Technical Report CS-89-85, University of Tennessee, 2007. [13] Z. Drezner and H.W. Hamacher, editors. Facility Location. Applications and Theory. Springer, New-York, 2002. [14] D. Goldberg. What Every Computer Scientist Should Know About Floating-Point Arithmetic. Computing Surveys, 23(1):5–48, 1991. [15] B. Hayes. A Lucid Interval. American Scientist, 91(6):484–488, 2003. [16] M. J. Hirsch, C. N. Meneses, P. M. Pardalos, and M. G. C. Resende. Global optimization by continuous grasp. Optimization Letters, 1:201–212, 2007. 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [17] L. Jaulin, E. Walter, and O. Didrit. Guaranteed robust nonlinear parameter bounding. In CESA’96, IMACS Multiconference (Symposium on Modelling, Analysis and Simulation), pages 1156 – 1161, 1996. [18] L. Jaulin and E. Walter. Guaranteed nonlinear parameter estimation from bounded-error data via interval analysis. Mathematics and Computers in Simulation, 35(2):123–137, 1993. [19] L. Jaulin and E. Walter. Guaranteed robust nonlinear minimax estimation. IEEE Transaction on Automatic Control, 47(11):1857–1864, 2002. [20] O. Lhomme. Consistency Techniques for Numeric CSPs. In IJCAI ’93: Proceedings of the Thirteenth International Joint Conference on Artificial Intelligence, pages 232–238, 1993. [21] A.K. Mackworth. Consistency in Networks of Relations. Artificial Intelligence, 8(1):99–118, 1977. [22] M. Milanese and A. Vicino. Estimation theory for nonlinear models and set membership uncertainty. Automatica, 27(2):403–408, 1991. [23] S. Minton, M.D. Johnston, A.B. Philips, and P. Laird. Minimizing conflicts: a heuristic repair method for constraint satisfaction and scheduling problems. Artificial Intelligence, 58(1-3):161–205, 1992. [24] R. E. Moore. Interval Analysis. Prentice-Hall, Englewood Cliffs, N.J., 1966. [25] Bruce A. Murtagh and Michael A. Saunders. Minos 5.5 user’s guide. Technical report, Systems Optimization Laboratory, Department of Operations Research, Stanford University, 1998. [26] A. Neumaier. Interval Methods for Systems of Equations. Cambridge University Press, 1990. [27] J.-M. Normand. Placement de caméra en environnements virtuels. PhD thesis, Université de Nantes, 2008. [28] T. Petit, J.-C. Régin, and C. Bessière. Range-Based Algorithm for MaxCSP. In Principles and Practice of Constraint Programming, CP 2002, pages 280–294. Springer-Verlag, 2002. [29] R.J. Wallace. Analysis of Heuristic Methods for Partial Constraint Satisfaction Problems. In Principles and Practice of Constraint Programming, CP 1996, pages 482–496. Springer-Verlag, 1996. 28