Thierry Martinez

François Fages

Abder Aggoun

A Stochastic Continuous Optimization Backend for MiniZinc with Applications to Geometrical Placement Problems

MiniZinc is a solver-independent constraint modeling language which is increasingly used in the constraint programming community. It can be used to compare different solvers which are currently based on either Constraint Programming, Boolean satisfiability, Mixed Integer Linear Programming, and recently Local Search. In this paper we present a stochastic continuous optimization backend for MiniZinc models over real numbers. More specifically, we describe the translation of FlatZinc models into objective functions over the reals, and their use as fitness functions for the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) solver. We illustrate this approach with the declarative modeling and solving of hard geometrical placement problems, motivated by packing applications in logistics involving mixed square-curved shapes and complex shapes defined by Bézier curves.

Introduction

MiniZinc [START_REF] Nethercote | MiniZinc: Towards a standard CP modelling language[END_REF] is a medium-level constraint modeling language which is becoming a standard in the Constraint Programming community. It is high-level enough to express most constraint problems easily, but low-level enough to be mapped onto existing solvers easily and consistently. This mapping is done through a flattening process which takes as input a MiniZinc instance and produces a FlatZinc instance. FlatZinc is a low-level solver input language designed to be easy to translate into the form required by a solver. It is chosen for that reason as target language for MiniZinc.

Currently, there exist FlatZinc backends for Mixed Integer Linear Programming (CPLEX, OR-tools 1 , SCIP, . . . ), Finite Domain Constraint Programming solvers (Choco 2 , Eclipse 3 , Gecode 4 , JaCoP, Opturion-CPX 5 , Oscar, SICStus prolog, . . . ), SAT solvers (MinisatID,. . . ) and recently Local Search (iZplus 6 , Oscar-cbls [START_REF] Björdal | A constraint-based local search backend for minizinc[END_REF]).

Most of FlatZinc implementations are thus dedicated to discrete domains. However, constraint optimization and decision problems over real numbers can be expressed in MiniZinc with high generality. Curently, such continuous constraint problems can be solved either using Linear Programming backends, with restrictions on the linearity of the constraints, or using interval arithmetic backends (e.g. G12ic, Eclipse fzn ic).

In this paper, we study another kind of solver based on stochastic continuous optimization for solving FlatZinc instances over real numbers, using namely the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]. More specifically, we show how a FlatZinc instance over real numbers can be translated into a fitness function which can be directly used by CMA-ES to compute approximate solutions to the problem. The transformation we describe is quite general and applies virtually to any MiniZinc model over real numbers. The choice of CMA-ES among other evolutionary or particle swarm optimization algorithms is motivated by the absence of parameterization for this algorithm and by its performances on hard problems.

For discrete domains, there has been related work on the design of high-level constraint-based modeling languages for local search and genetic algorithms. The seminal work of Van Hentenryck and Michel on Comet [START_REF] Michel | The comet programming language and system[END_REF][START_REF] Hentenryck | Synthesis of constraint-based local search algorithms from high-level models[END_REF] showed how a finite domain constraint model can be compiled into an objective function for local search metaheuristics, such as Tabu search, with default neighborhoods derived from the constraint model. In [START_REF] Björdal | A constraint-based local search backend for minizinc[END_REF], Björdal et al. present a constraint-based local search backend for MiniZinc and show that it produces competitive results on the 2010 to 2014 MiniZinc challenges. In these systems, the local search solver is limited to finite domain constraints and use neighborhoods derived from the finite domains of the variables.

Here in the continuous domain, we illustrate our CMA-ES backend for FlatZinc with the solving of hard geometrical placement problems which, to the best of our knowledge, go beyond the state-of-the-art of declarative constraint modeling and solving. As a matter of fact, the only FlatZinc implementations listed on the MiniZinc web page that parse the FlatZinc instances presented in this paper are those based on exact methods using interval arithmetic (i.e. Eclipse fzn ic and G12ic) but none of them can find solutions in reasonable computation time even for the examples presented here. In [START_REF] Martinez | On solving mixed shapes packing problems by continuous optimization with the cma evolution strategy[END_REF], we have already shown that the non-overlap constraint between squares, cubes, rectangles, boxes, triangles, polygons circles and spheres, can be associated with a measure of overlap between objects which can be used directly as a fitness function in CMA-ES for packing mixed shapes in a bin, with an interesting trade-off between generality and efficiency. The measure of overlap does not need to be the area of the intersection (and should not if one object can be included in another) but can be any measure equal to 0 in case of non-overlap, and capable of guiding the continuous optimization solver by measuring progress toward satisfaction [START_REF] Chernov | Mathematical model and efficient algorithms for object packing problems[END_REF]. On a benchmark of consecutive sizes circle packing problems, we showed that CMA-ES finds solutions at 2% of the best known costs obtained by running the three global optimization methods reported in Castillo et al. [START_REF] Castillo | Solving circle packing problems by global optimization: Numerical results and industrial applications[END_REF]. In [START_REF] Salas | Packing curved objects[END_REF], Salas and Chabert show that the overlap measures which were defined in an ad hoc manner in [START_REF] Martinez | On solving mixed shapes packing problems by continuous optimization with the cma evolution strategy[END_REF], can be computed by interval methods in IBEX7 with a numerical algorithm that automatically measures the penetration depth of two objects of virtually any shape defined by conjunction and disjunction of non-linear inequalities. In this paper, we give general MiniZinc definitions for the penetration depths, or simpler overlap measures, between polygons, circles, and also complex shapes defined by Bézier curves, motivated by packing problems in the cosmetic and automotive industries. This illustrates the performance of MiniZinc-CMAES in terms of both declarative modeling and efficient (yet suboptimal) resolution of very hard geometrical packing problems with complex shapes and continuous rotations.

The rest of the paper is organized as follows. In the next section, we present the translation of a FlatZinc instance over real numbers in a fitness function over the reals, and the interface to the CMA-ES solver. In Section 4 we describe MiniZinc models of continuous packing problems involving continuous rotations, mixed square-curved shapes and complex shapes defined by Bézier curves. There we use some simple distance formulae for circles, Minkowski sums for the penetration depth between polygons [START_REF] Dobkin | Computing the intersection-depth of polyhedra[END_REF] and De Casteljau's numerical algorithm for linearizing Bézier curves. In Section 5, we report on the performance results obtained through the compilation chain from MiniZinc, FlatZinc to CMA-ES, on complex shape packing problems. Finally, we conclude on the general perspective opened by this MiniZinc backend for continuous optimization and novel applications at the intersection of Optimization and Computer-Aided Design.

Compiling FlatZinc Instances over Real Numbers in

Real-valued Fitness Functions

In this section we describe our transformation of a FlatZinc instance containing arithmetic and trigonometric constraints over float variables in a fitness function which aggregates the costs of each constraint violation. This transformation is at the heart of the continuous optimization backend.

Arithmetic expressions

The arithmetic expressions that constitute the constraint satisfaction problem need be rebuilt from the FlatZinc instance, since arithmetic sub-expressions and intermediary variables are introduced by the transformation from MiniZinc to FlatZinc. The constraints that result from this transformation are split in three groups:

1. inequality constraints, which are turned into costs, 2. arithmetic and trigonometric constraints, which always appear to be directed and are turned into functional expressions, 3. and equality constraints, which can either be solved statically if there exists a topological sort of the constraint graph that makes the constraint directed, or turned into a cost otherwise.

Every variable X of the model is associated to an expression [X] defined as [X] = X if X is one of the search variables, or as the arithmetic or trigonometric expression deduced from the constraints on X. The notation is extended to float constants, [f ] = f , and vectors of variables and/or float constants, i.e. for a vector U = (X 1 , . . . , X n ), [U ] denotes the vector ([X 1 ], . . . , [X n ]).

Concerning the first group, the inequality constraints in FlatZinc are either strict or non-strict inequalities between linear expression of the form:

constraint float_lin_lt ( U0 , U1 , K ) ; constraint float_lin_le ( U0 , U1 , K ) ;
where K is a constant, U 0 is a vector of constant coefficient and U 1 is a vector of variables. The semantics is respectively U 0 • U 1 < K and U 0 • U 1 ≤ K where • denotes the scalar product. The cost we associate to such an inequality constraint c is cost

(c) = max(0, [U 0 ] • [U 1 ] -K)
where [U 0 ] and [U 1 ] are the expressions constructed from the arguments. For strict inequality, one could refine the cost to

cost(c) = 1 + [U 0 ] • [U 1 ] -K if [U 0 ] • [U 1 ] ≥ K 0 otherwise
in order to ensure that the cost is null if and only if the constraint is satisfied. However, this is not necessary to guide the search for solutions since the constraints a < b and a ≤ b are equivalent almost everywhere in a continuous setting.

Concerning the second group, every float variable X of the model is considered as a search variable, except if it appears in the result position of one of the following directed constraint: constraint float_min (A , B , X ) ; constraint float_max (A , B , X ) ; constraint float_times (A , B , X ) ; constraint float_sqrt (A , X ) ; constraint float_cos (A , X ) ; constraint float_sin (A , X ) ; If X is in the result position of one of these constraints, then [X] is defined as the expression that computes the associated value.

For the third group, FlatZinc linear equality constraints are of the form:

constraint float_lin_eq ( U0 , U1 , K ) :: defines_var ( X ) :: weight ( w ) ;

The FlatZinc compiler generates the annotation defines_var(X) which directs the constraint from its arguments to its result, in the case where the constraint results from a linear arithmetic expression. In that case, [X] is defined as the expression that computes the linear combination.

For the equality constraints that result from reification, of the form:

constraint float_eq_reif (X ,Y , B ) ; constraint f l o a t _ l i n _ e q _ r e i f (X ,U , B ) ;

we eliminate them statically if enough information is available. Otherwise, the model is currently rejected.

It is worth noting that this minimal handling of reification is needed to cope with the code generated by the MiniZinc compiler for partial functions like sqrt , where the formal argument is unified with the actual argument only if the function is defined for this argument. On the other hand, general reified equality constraints impose integrity constraints on the boolean variables. Turning an integrity constraint into a cost function causes rugged landscapes which may be difficult to explore, although CMA-ES is also a pretty good solver in this case. Since our benchmarks do not use reified constraints nor discrete variables, they are currently out of the scope of our MiniZinc backend.

Cost aggregation

In our backend, the gathering of the cost functions can be tuned by using annotations. First, each constraint can be annotated with a weight which will affect the cost in the aggregation. By default, the weight of a constraint is 1.

annotation weight ( float ) ;

Second, every MiniZinc model contains one and only one solve instruction, which gives the objective. The solve item can be annotated with one of the following annotations which change the definition of the violation cost of the whole constraints. By default, weighted_sum is assumed.

violation cost = 1 - c ( 1 1 + cost(c) ) weight(c)
Third, if the FlatZinc instance is a constraint satisfaction problem, the fitness function is defined to be equal to the violation cost. If the FlatZinc instance requires to minimize an expression e, the fitness function is defined to be equal to α • (1 + violation cost) + e, where α is a coefficient large enough to dominate e. Then α can be set with the following annotation applied to the solve item.

Stochastic Continuous Optimization with CMA-ES

The Covariance Matrix Adaptation Evolution Strategy (CMA-ES8 ) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] is one of the most powerful global optimization strategy for minimizing an objective function over the reals in a "black-box" scenario, i.e. without assuming any property about the objective function. This method is a multi-point method which uses a population of configurations (here valuations of the FlatZinc search variables, e.g. packings defined by the coordinates and orientations of the objects) to sample the search space, estimates the covariance matrix at each sampling, determines the next move in the most promising direction (e.g. translations and rotations of objects), and updates accordingly the multi-variate normal distribution for the next sampling (i.e. mean value and covariance of the variables).

CMA-ES behaves in effect like a second-order method where the landscape is estimated by sampling, according to some multi-variate normal distribution of the variables, which is itself updated during search in the most promising direction to adapt to the landscape, using an estimation of the second-order moment, the covariance matrix. When the objective function does not improve, CMA-ES can be restarted to find different local optima. We refer to [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] for more details on that stochastic optimization algorithm.

One advantage of CMA-ES is that it requires very little parameter tuning. All our benchmarks have been performed with the C implementation of CMA-ES using the same parameter set: a population size of 100, an initial standard deviation of 20 and a stopping criterion based on a difference less than 10 -3 for the fitness function.

CMA-ES thus tries to minimize an arbitrary function f : R n → R, where n is the dimension of the search space (i.e. the number of FlatZinc search variables). The result is a vector x ∈ R n such that f (x) is the smallest value encountered so far. The C implementation of CMA-ES expects that the function f has the following interface: double f(double x[]). The FlatZinc-to-CMA-ES back-end derives such a fitness function from the FlatZinc model according to the transformations described in Section 2.

MiniZinc Models of Geometrical Placement Problems

The problems addressed in this section are taken from the industry of cosmetics packaging. They consist of packing products with various shapes. In this application the objective is to pack a given quantity of a product in a minimum number of bins. The study of different forms in the industry of cosmetics packaging show that convex approximations of objects give poor results but that the objects can be modeled using a combination of Bézier curves. Existing Constraint Programming tools are limited and do not offer capabilities to solve such problems. We show here how to model and solve such problems in MiniZinc.

Overlap Measures Between Objects

The overlap measure should guide the optimization procedure towards a geometrical placement without overlap, i.e., ideally, the cost should decrease as long as the placement gets closer towards a placement without overlap. φ-functions [START_REF] Chernov | Mathematical model and efficient algorithms for object packing problems[END_REF] have been introduced for the same purpose of continuous optimization for geometrical placement problems, using decompositions in half-planes, triangles and circles. In this section, we describe three overlap measures for polygones and complex shapes delimited by Bézier curves using the intersection area, the penetration depth, or the sum of the pairwise distances between the intersection points of the borders.

Intersection area. Object intersection area can be used as an overlap measure, which could seem quite natural. However, this area can be costly to compute and does not guide the optimization well. For example, when an object fully contains another one, the overlap measure remains in a plateau for every possible position where the contained object stays inside the container, giving no direction for getting outside the overlap zone (figure 1). In pratice, we do not use this measure. Penetration depth. The penetration depth is a common measure in computeraided design [START_REF] Cameron | Determining the minimum translational distance between two convex polyhedra[END_REF]: the penetration depth between two objects is the smallest norm such that there exists a translation vector to apply to one of the two objects such to lead to a placement where the two objects do not intersect.

The penetration depth between two circles (C 1 , r 1 ) and (C 2 , r 2 ) is trivial to compute:

pd(C 1 , r 1 ); (C 2 , r 2 )) = max(0, r 1 + r 2 -C 1 C 2 )
that is to say the difference between the sum of their radius and the distance between their centers (figure 2). This distance was already used for circles in [START_REF] Martinez | On solving mixed shapes packing problems by continuous optimization with the cma evolution strategy[END_REF][START_REF] Chernov | Mathematical model and efficient algorithms for object packing problems[END_REF].

C1 C2 r1 r1 r2
Fig. 2. Penetration depth between two circles (C1, r1) and (C2, r2)

More generally, the penetration depth between two objects O 1 and O 2 is equal to the distance between the origin (0, 0) and the complementary of the Minkowski difference

O 1 O 2 = {p 1 -p 2 | p 1 ∈ O 1 , p 2 ∈ O 2 }. pd(O 1 ; O 2 ) = min{ u | u / ∈ O 1 O 2 }
Indeed, this distance is by definition the smallest norm such that there exists a vector u such that u / ∈ O 1 O 2 , that is to say a vector such that for all

p 1 ∈ O 1 , p 2 ∈ O 2 , u = p 1 -p 2 , thus we have O 1 ∩ (O 2 + u) = ∅.
The Minkowski difference of two polygons is a polygon, computable in a time quadratic to the number of edges, and the Minkowski difference of two convex polygons is a convex polygon, computable in a time linear to the number of edges [START_REF] Dobkin | Computing the intersection-depth of polyhedra[END_REF] (figure 3).

Note that the penetration depth only consider translations, whereas the search space we consider for optimization may include rotation angles as additional "dimensions". [START_REF] Salas | Packing curved objects[END_REF] extends the Minkowski difference to consider object rotations as well. This extension is difficult to interpret geometrically and makes overlap measures depend on the choice of the origin for each object. We will restrict ourselves to Minkowski difference in the Euclidean space.

Sum of the pairwise distances between the intersection points of the borders. For the overlap measure between two objects O 1 and O 2 that have non polygonal shapes like those delimited by Bézier curves, or for heterogeneous shapes (for example, when mixing Bézier curves and polygons), we prefer to use another measure simpler to compute: the sum of the pairwise distances between the intersection points of the borders ∂O 1 and ∂O 2 . We suppose that ∂(O 1 ) ∩ ∂(O 2 ) O1 O2 (0, 0) Fig. 3. Minkowski difference between two pentagons O1 and O2. The penetration depth between O1 and O2 is the distance between the origin and the border of the Minkowski difference.

is a finite set, applying infinitesimal offsets if necessary.

spd(O 1 ; O 2 ) = pi,pj ∈∂(O1)∩∂(O2) i<j --→ p 1 p 2
There exist several methods to compute the intersections between two Bézier curves numerically [START_REF] Sederberg | Chapter 7, planar curve intersection[END_REF] (figure 4). We use a dichotomic search by using de Casteljau's algorithm for splitting the curves. The dichotomic search can also be used to compute numerically the intersections between Bézier curves and circles. The intersections between a Bézier curve and a segment can be computed algebraically. Indeed, by changing the frame, we can suppose without loss of generality that the segment lays on the abscissa axis. The Bézier curve (p 0 , p 1 , p 2 ) intersects the axis for every parameter t, 0 ≤ t ≤ 1, such that (1 -t)((1 -t)y p0 + typ 1 ) + t((1 -t)y p1 + typ 2 ) = 0: this is a second-order polynomial in t. For each solution t 0 , it suffices to check that the abscissa (1 -t 0 )((1 -t 0 )x p0 + t 0 xp 1 ) + t 0 ((1 -t 0 )x p1 + t 0 xp 2 ) belongs to the segment. This measure does not fulfill the requirements of an ideal overlap measure: the measure is null when one object is included in the other and is not monotonic with respect to the penetration depth. However, it is locally monotonic in a neighborhood around overlap-free placements: in the context of a local search, by choosing an overlap-free initial placement (spreading the objects enough far ones from the others), this measure experimentally appears to be sufficient to preserve the overlap-freeness during the placement compaction process.

Continuous Packing Model

This section makes use of the overlap measures introduced above to express MiniZinc models for continuous packing of circles, arbitrary polygons with rotation and other complex shapes like rosettes delimited by Bézier curves. The penetration depth is used as measure of overlap between two circles and between two polygons, while the sum of distances between the intersection points is used for every other pair of shapes.

Circles. The following predicate expresses the constraint that the two circles ((x 1 , y 1 ), r 1 ) and ((x 2 , y 2 ), r 2 ) do not overlap. It is worth noticing that the inequality (

x 1 -x 2 ) 2 + (y 1 -y 2 ) 2 > (r 1 + r 2 ) 2 is compiled into the cost function (r 1 +r 2 ) 2 -(x 1 -x 2 ) 2 -(y 1 -y 2 ) 2 ,
which is monotonic with respect to the penetration depth (r

1 +r 2 )-(x 1 -x 2 ) 2 -(y 1 -y 2 ) 2 , introduced in the previous section.
We consider a benchmark of circle placement problems [START_REF] Castillo | Solving circle packing problems by global optimization: Numerical results and industrial applications[END_REF] where there are n circles to place in a circular bin. The usual modelling [START_REF] Martinez | On solving mixed shapes packing problems by continuous optimization with the cma evolution strategy[END_REF] supposes that the circular bin is centered on the origin. The following function compute for each circle the minimum radius for the circular bin to contain the circle (x, y), r). The search variables are the positions of the circle centers.

int : n ; array [1 .. n ] of var float : x ; array [1 .. n ] of var float : y ;

The circle positions are constrained to be non-overlapping.

constraint forall ( i in 1.. n , j in i +1.. n ) ( n o n _ o v e r l a p _ c i r c l e s ( x [ i ] , y [ i ] , radius ( i ) , x [ j ] , y [ j ] , radius ( j ) ) ) ;
The goal is to minimize the radius of the circular bin. Intermediary variables are introduced to store the minimal bounding radius for each circle to circumvent a limitation of the max function for arrays in MiniZinc that needs to know the bounds of the arguments. It is worth noticing that applying min and max functions to overlap measures allows Boolean combinations of geometrical shapes to be expressed. For instance, the figure 6 shows a placement for 20 geometrical rosettes, where each rosette R((x, y), r) is defined as the union of six intersections between pairs of circles:

R((x, y), r) = 6 i=1 C((x + cos(2 • i • π 6 ) • r, y + sin(2 • i • π 6 ) • r), r) ∩ C((x + cos(2 • (i + 2) • π 6 ) • r, y + sin(2 • (i + 2) • π 6 ) • r), r)
Objects with rotations. The placement of each object in the subsequent examples is described by three search variables: the position on the x axis, the position on the y axis, and the rotation angle r.

set of int : position = 1 .. 3; int : x = 1; int : y = 2; int : r = 3;

An object is described by a variable of type array[position] of var float, that is to say an array of three variables. The constants x, y and r are used as projectors: given an array object, components can be accessed as object[x], object[y] and object [r]. (Note that MiniZinc has not yet support for records.)

Points are stored in an array of two coordinates. The points that describe the shapes of an object are expressed in a frame relative to the given object position and orientation. The function image_of_point defined below transforms the coordinates of a point to the global frame. The following function returns the position of an object given its index. Polygons. We consider pentagons with the following vertex coordinates (relative to the object frame). These pentagons approximate the Bézier rosettes that we consider below: the vertices join the ends of the petals (figure 7). It is not exactly the convex hull since one petal goes outside the pentagon but it is close to (and the convex hull of the rosette is not polyhedric). However, this approximation is sufficient to observe the gain obtained in the placements by considering the precise Bézier rosettes instead of such approximations.

We suppose that the following function computes the penetration depth between two (convex) polygons. The Minkowski difference between two polygons can be expressed with arithmetic constraints through reification: the constraint is cumbersome to write directly, but can be automatically generated, for example by a ClpZinc model [START_REF] Martinez | On translating minizinc constraint models into fitness function for evolutionary algorithms: Application to continuous placement problems[END_REF]. Alternatively, the function can be implemented in the back-end as an auxiliary C function, which is the case of our current implementation. We suppose that the following function computes the sum of the distances between the intersections of two sets of curves. This function is implemented in the back-end as an auxiliary C function. Mixing Bézier curves and rectangles. For computing the overlaps between Bézier curves and rectangles, we suppose that the following function computes the sum of the distances between the intersections of a Bézier curve and a rectangle. This function is implemented in the back-end (as an auxiliary C function), but the arithmetic could be expressed in MiniZinc as well. The figure 10 shows a placement found for 16 Bézier rosettes and 16 rectangles. It is worth noticing that even if the optimization procedure has found a non-trivial placement, for instance for the rectangles and the rosettes in the bottom right of the figure, some visually obvious improvements of the placement of the left rosettes are not found in this run of CMA-ES.

Evaluation Results of MiniZinc-CMAES

The following tables 1 and 2 summarize the performance obtained with MiniZinc-CMAES, in terms of computation time, smallest area found, mean area and variance of the area among 50 runs of CMA-ES. For every example, results are averaged over 50 runs. All these results have been obtained with the default parameters of CMA-ES described in Section 3. It is worth noticing that smaller initial standard deviations tend to generate solutions with overlaps that the optimization fails to remove, and bigger standard deviations augment convergence times. Total time is the sum of the computation times for all the 50 restarts: for each problem, all the restarts have been computed in parallel on a cluster, one problem per core. It is worth noticing that Eclipse and G12 with their interval constraint solvers can parse the MiniZinc models of the previous section that do not use predicates defined as auxiliary C functions in the back-end. However the performances are very poor with results obtained only for 3 circles.

Conclusion

We have presented here a stochastic continuous optimization backend for MiniZinc models over real numbers. We have shown the benefits of this approach using the CMA-ES solver for continuous optimization on a series of geometrical placement problems motivated by industrial applications in logistics, involving mixed square-curve shapes, and also complex shapes defined by Bézier curves. Probably because of the novelty of these problems for complex shapes, we have not identified benchmarks for comparing the techniques presented here, but in [START_REF] Martinez | On solving mixed shapes packing problems by continuous optimization with the cma evolution strategy[END_REF] we showed that the solutions found with CMA-ES on circle packing were at just 2% of the best solutions found with dedicated solvers.

The declarative modeling in MiniZinc combined to the solving using the transformation to CMA-ES described in this paper, does not come with any significant overhead and provides fully declarative solutions to very hard geometrical placement problems. The non-overlap constraint has a cost function based on the penetration depths between objects, using Minkowski sums for polygons, and a simpler measure of overlap for Bézier curves. A classical difficulty in the definition of the error function of a conjunction of constraints is the normalization of the error function for each constraint. This has been solved here by letting the modeller specify in MiniZinc the cost function if different from the default cost aggregation function (i.e. the sum of the costs).

The recourse to such a black-box optimization procedure for FlatZinc makes sense especially in presence of non-linear constraints, and in absence of integer variables, but the transformation we have given of a FlatZinc model in a non-negative real-valued cost function is quite general. We have focused on continuous placement problems, but our MiniZinc/CMA-ES can be applied in principle to any constraint model over real numbers. The examples taken here from industrial problems in logistics, including objects defined by Bézier curves, should contribute to open a new domain of application of constraint methods in computational geometry, at the intersection of optimization and computer-aided design.

  ) • cost(c) annotation fuzzy ; violation cost = max c weight(c) • cost(c) annotation probabilistic ;

  annotation alpha ( float ) ; By default, α = 10 10 .

Fig. 1 .

 1 Fig.1. Intersection area between two rectangles O1 and O2 in function of x(O2), with a plateau when O2 is included in O1.

Fig. 4 .

 4 Fig. 4. Intersection points of the borders of two objects delimited by Bézier curves.

  predicate n o n _ o v e r l a p _ c i r c l e s ( var float : x1 , var float : y1 , var float : r1 , var float : x2 , var float : y2 , var float : r2 ) = pow ( x1 -x2 , 2) + pow ( y1 -y2 , 2) > pow ( r1 + r2 , 2) ;

  function var float : b o u n d i n g _ c i r c l e _ r a d i u s ( var float : x , var float : y , var float : r ) = sqrt ( pow (x , 2) + pow (y , 2) ) + r ;

array [ 1 .Fig. 5 .

 15 Fig.5. Example of packing found by MiniZinc-CMAES for 18 circles of radii i -1/2 for 1 ≤ i ≤ 18 (circle packing benchmark of[START_REF] Castillo | Solving circle packing problems by global optimization: Numerical results and industrial applications[END_REF]).

Fig. 6 .

 6 Fig. 6. Placement found by MiniZinc-CMAES for 20 geometrical rosettes, defined as unions of circle intersections.

  function array [ position ] of var float : o bj e ct _p os i ti on ( int : object ) = [ o b j e c t _ p o s i t i o n s [ object , d ] | d in position ];

Fig. 7 .

 7 Fig.7. The pentagon obtained by considering each petal's end as a vertex is not a correct approximation of a Bézier rosette.

Fig. 8 .

 8 Fig. 8. Placement found by MiniZinc-CMAES for 29 pentagons

  set of int : curves = 1 .. 10; set of int : q u a d r a t i c _ b e z i e r _ c o n t r o l _ p o i n t s = 1 ..

Fig. 9 .

 9 Fig. 9. Placement found by MiniZinc-CMAES for 10 Bézier rosettes

  function var float : s u m _ o f _ d i s t a n c e s _ b e t w e e n _ b e z i e r _ a n d _ p o l y g o n _ i n t e r s e c t i o n _ p o i n t s ( array [ int , coordinates ] of var float : curves , array [ int , coordinates ] of var float : vertices ) ;

Fig. 10 .

 10 Fig. 10. Placement found by MiniZinc-CMAES for 16 Bézier rosettes and 16 rectangles, in this run of CMA-ES which stays in a local minimum.

Table 1 .

 1 Computation time for placement of geometrical rosettes. .345 5.499 • 10 5 13 + 13 13 j 23 h 13 min 39 s 6 h 42 min 16 s 85 492.98 1.024 • 10 5 6.985 • 10 5 14 + 14 17 j 9 h 27 min 23 s 8 h 20 min 56 s 1.11 • 10 5 1.398 • 10 5 8.219 • 10 5 15 + 15 23 j 9 h 18 min 63 s 11 h 13 min 35 s 1.078 • 10 5 1.584 • 10 5 1.281 • 10 6

	Roses Total time (50 restarts) Mean time Smallest area found Mean area Variance (area)
	10	17 min 32 s	21 s	25.266	27.303	1.893
	11	31 min 13 s	37 s	27.369	29.692	1.585
	12	41 min 18 s	49 s	28.761	32.511	3.675
	13	1 h 1 min 39 s	1 min 13 s	32.936	34.987	2.648
	14	1 h 17 min 37 s	1 min 33 s	33.816	37.714	2.926
	15	1 h 43 min 19 s	2 min 3 s	37.233	41.085	3.57
	16	2 h 6 min 57 s	2 min 32 s	39.729	43.86	5.891
	17	2 h 26 min 23 s	2 min 55 s	41.883	46.113	5.947
	18	3 h 11 min 20 s	3 min 49 s	43.582	49.123	13.828
	19	3 h 47 min 24 s	4 min 32 s	46.74	52.594	10.769
	20	5 h 8 min 42 s	6 min 10 s	49.006	54.89	8.579

Table 2 .

 2 Computation time for placement of mixed shapes: Bézier's rosettes and rectangles.

http://www.ibex-lib.org

https://www.lri.fr/ ~hansen/cmaesintro.html

Acknowledgements. This work has been funded by the ANR Blanc Simi2 Net-WMS-2 grant ANR-11-BS02-0005. We would like to thank all the partners of this project for fruitful discussions.


Abderrahmane Aggoun
email: abder.aggoun@klsoptim.com

Ahmed Rhiat
email: ahmed.rhiat@icam.fr

François Fages
email: francois.fages@inria.fr

Panorama of real-life applications in logistics embedding bin packing optimization algorithms, robotics and cloud computing technologies

Keywords: computing technologies Optimization, Constraint Programming, Bin Packing, Metaheuristics, Cloud Computing, Operational Logistics, IOT, Robotics, RFID

HAL is

INTRODUCTION

Problems addressed are taken from production, warehousing and distribution. In production, typical problems are the design of optimal packing plans of identical items in containers. In the article we present successful real-word applications using state of the art of advanced optimization algorithms where optimality is often required. The article also presents current state of the art of dealing with complex shapes where the problems to solve are still hard to solve and the exploitation of optimization algorithms is not trivial. In warehousing and distribution, problems addressed concern order preparation of large orders for intermediate platforms or important customers like manufacturers in different sectors. An order is a set of lines. Each line describes the ordered product with its corresponding quantity. The problem consists of packing products of various shapes and sizes into available pallets (large boxes) in a way which optimizes the total number of pallets. The objective is to minimize the total number of pallets whist handling several types of pallets and subject to orientation, fragility, stackability, weight and volume constraints. All these problems are known as bin packing problems [START_REF] Lodi | Two Dimensional Packing Problem: A survey[END_REF][START_REF] Martinez | On Solving Mixed Shapes Packing Problems by Continuous Optimization with the CMA Evolution Strategy[END_REF]. Bin packing algorithms were extended and applied to operational planning of trains. To come up with working solutions, we developed several optimization solvers combining several technologies: Constraint Programming, Metaheuristics like tabu search or specific solvers for handling complex shapes with non-linear inequalities. Solvers are integrated in Full Web applications and Cloud Computing platforms, many of them have been operational for many years. In this article we also present preliminary results of integration of robotics, drones and optimization and the potential applications of such results in production, warehousing and e-commerce.

Constraint Programming

Constraint Programming (CP) [START_REF] Aggoun | The contribution of global constraints for modeling and solving industrial applications[END_REF] has become a reality in the industrial world. Its success represents an important development of Artificial Intelligence and a successful vulgarization of the techniques of Operations Research in decision-making solutions in various domains. Indeed, CP has effectively addressed a large class of combinatorial problems. The Constraint Programming system used is Choco [http://choco.emn.fr] designed to tackle real world "constrained search" problems with a short-term development time and good efficiency. The Choco system provides a number of constraint solvers over different computation domains. These domains have been chosen because of their interest for applications, but also since efficient specialized constraint solving methods exist for these domains. Constraint Programming is based on the idea that many interesting and difficult problems can be expressed declaratively in terms of variables and constraints. The variables range over a (finite) set of values and typically denote alternative decisions to be taken. The constraints are expressed as relations over subsets of variables and restrict feasible value combinations for these variables. A solution is an assignment of values (search) to variables, which satisfies all constraints. It is the ability of CP systems to take into account heterogeneous constraints effectively which, in the 90s, raised a commercial interest in this paradigm.

Bin packing optimization solutions 3.1 Optimization of pallets

The problem addressed is taken from warehousing and distribution. It concerns the order preparation of large orders for intermediate platforms or important customers like manufacturers in different sectors. An order is a set of lines. Each line describes the ordered product with its corresponding quantity. A carton in a pallet may have up to six possible orientations which may increase the complexity in some cases. The problem consists of packing boxes subject to orientation, fragility, stackability, weight and volume constraints of various sizes into available pallets (larges boxes) in a way which optimizes the total number of pallets. The objective is to minimize the total number of pallets.

The type of planning is operational planning which means that some decisions must be taken in real time. The context is the following. The in-house warehouse management system (WMS) is connected to an ERP (SAP). The WMS reads the orders to prepare each lap of time (e.g. 1 minute) from the ERP. The main constraint is that the WMS must compute the optimal number of pallets in real time and the optimal linear meters of the virtual vehicle. This decision making information is sent in real time to the transporters (if they are known) or to the transport department. The objective is to compute this decision making information in just a few seconds. For example, the decision is taken at 09:00 and the preparation will start at 15:00. Usually the exact number of pallets is known once the order preparation is finished. In this case, the consequence is that the transport is informed at the last minute with many disadvantages; e.g. sending a large vehicle by using half of the vehicle capacity. We have noticed that objectives and requirements are subject to many changes depending on the final client. In the same time the number of different sizes of boxes increases each year. Therefore we developed many parametric variants of solvers with different strategies like packing by layers, by piles or a combination of both. A CP program is structured in three parts: declaration of decision variables, statement of constraints and enumeration. Constraint Programming allows to separate the modeling of heuristics and strategies from constraints. The search combines constructive solution procedures and some CP strategies. It takes advantage of the constraint programming propagation (domain of variables is reduced after each step of assignment) and the strategies used to select each step, the best object (container, configuration), the best decision variable and the best value. Search procedures make it possible to find a compromise between optimal solution and quality solutions which are preferred in practice. Different variants of optimizations of pallets are deployed in various sectors: automotive, perfume industries, pharmacy, foodprocessing, chocolate factory, etc.

Vehicle / Container Loading Plans

Optimization of vehicle loading is a well-known problem in TMS (Transport Management Systems). For some sector activities in logistics they are complementary to pallet optimization. For some variants the combinatorial problems to solve are similar. Furthermore, one can find similar constraints. First we distinguish two types of problems: (1) a single container loading problem (SCLP): Boxes of different sizes are stowed in the container as is possible by minimizing the unused. The objective is to maximize the use of a container's volume; (2) a multiple container loading problem (MCLP): Boxes of different sizes are stowed in several containers. The objective is to minimize the total number of required containers. Optim Vehicle / Truck tool provides several engine solvers. Each engine solver is dedicated to a variant of SCLP or MCLP problems. There are different variants of solvers for vehicle loading problems. For example, in some cases the loading is carried out in two steps. First one has to build pallets or piles (grouping compatible items) and then to optimize their loading into vehicles. The number of items to pack may vary from about ten to several thousand. Like the optimization of pallets, the search procedure is a key to compute quality solutions under time constraint.

Bin packing e-commerce 3.4 Applications

This application is a variant of the optimization of pallets where items are products (small cartons) and bins (large cartons). The objective is to attempt to find the smallest number of bins into which to pack all items. The complexity comes from the fact that there is a large number of candidate bins with different sizes (e.g. 4 to 30) with specific constraints. In addition the time allowed to compute a feasible solution is limited to a few milliseconds. This application is powered by a solver with dedicated heuristics. This is a typical application to deploy on the cloud.

Assignment of containers to trains

This application deals with the operational planning of trains [START_REF] Aggoun | Online assignments of containers to trains using constraint programming[END_REF]. Furthermore, assigning a destination of a train is treated beforehand and the results are inputs for component optimization. The number of wagons and their specifications are also known. The problem is known in literature as "bin packing". The objective is to minimize the number of wagons (bins) while placing the maximum number of containers (items). The application is an interactive Full Web application powered by a CP solver featuring the following functionalities: (1) an optimization engine to compute an initial solution and to assist the operator when completing partial solutions;

(2) a fast algorithm to assist the operator in choosing the right container for a wagon.

The constraints taken into account are the maximum weight constraints attached to wagons, the maximum number of containers per wagon, the number of wagons attached to a train and its total weight. Within these limited number of constraints the problem can be solved to optimal degrees. The application handles new additional business constraints, among them: the handling of dangerous goods, the incompatibilities between families of containers, e.g. two containers containing dangerous goods cannot be assigned to the same wagon and the original constraints of wagon-makers which are configuration constraints.

The handling of constraints associated with configurations of the selected wagons make the problem difficult to solve. Therefore bin packing algorithms are enriched and extended to solve variants of cutting stock instances. In addition, one must handle preference constraints to increase the quality of solutions. The use of the cumulative constraint [START_REF] Aggoun | The contribution of global constraints for modeling and solving industrial applications[END_REF] in CP combined with efficient search strategies made it possible to solve the different instances in less than one minute whilst it requires several hours by an operator. One of the interesting features of the application is the possibility to combine CP and Metaheuristics. The first solver which is CP based is launched to compute a first correct solution respecting all constraints. The second solver which is based on Metaheuristics assists interactive capabilities of the Full Web application like moving a container from one wagon to another or handling urgent containers.

Handling Curve shapes with CMA-ES algorithms

This part deals with innovative algorithms in bin packing. The project Net-WMS2, lead by the team Lifeware at INRIA, aims to carry out many studies for complex objects, amongst which solving curve shapes. In addition to packing problems for squares, circles, triangles and polygons [START_REF] Martinez | On Solving Mixed Shapes Packing Problems by Continuous Optimization with the CMA Evolution Strategy[END_REF], we consider continuous packing problems with curve shapes [START_REF] Martinez | A Stochastic Continuous Optimization Backend for MiniZinc with Applications to Geometrical Placement Problems[END_REF]; mixing various shapes and mixing three dimensional packing problems, mixing boxes, spheres and cylinders. On square packing problems, exact methods have been used to find optimal solutions and prove optimality. On curve packing problems, inexact methods like for instance genetic algorithms or hybrid simulated annealing with tabu search methods [START_REF] Talbi | Metaheuristics: from design to implementation[END_REF], are usually used to compute suboptimal solutions. In this study, we evaluate the covariance matrix adaptation evolution strategy (CMA-ES2), one of the most powerful evolutionary algorithms for continuous optimization with arbitrary objective functions, on hard packing problems including square shapes, curve shapes, mixed square and curve shapes, and also continuous rotations. The solver made it possible to model and solve a problem in the perfume industry where items are complex shapes defined as Bézier curves. These packing problems are encountered in automotive and perfume industries. The results are promising and details are in [START_REF] Martinez | A Stochastic Continuous Optimization Backend for MiniZinc with Applications to Geometrical Placement Problems[END_REF][START_REF] Martinez | On Solving Mixed Shapes Packing Problems by Continuous Optimization with the CMA Evolution Strategy[END_REF].

Bin packing and robotics

Robot packing problems is a new class of bin packing problems requiring dedicated solvers to which we refer to as robot packing algorithms. The problem addressed is solved in two phases. In the first phase the solvers computes an optimal or quasi-optimal solution of the problem to solve. The solution is feasible for operators but in many cases infeasible for robots. The first class of constraint is the placement of items in a box (a carton). The robot's movements are limited. The second class of constraints are known as accessibility constraints. Before placing the current item one has to make sure that the placement of remaining items is mechanically feasible. Unfortunately, the handling of these constraints leads to the removal of some conflicting items from the optimal solution.

ARCHITECTURE

A set of bin packing optimization applications are deployed on the cloud. The design of cloud computing applications requires much care. Some particular applications which are used daily are implemented in such a way to allow the maximum time to optimization algorithms whilst minimizing the exchange time between the server on the cloud and the distant client. Therefore the design of such applications requires deep collaboration between optimization and cloud computing engineers, which are complementary skills. The architecture design is J2EE compliant which is built from the ground up to be a service-oriented architecture (SOA), i.e., it provides most of its business logic encapsulated as middleware services thus enabling third-party applications to be integrated.

The chosen architecture enables the possibility to integrate the optimisation components in 2-tier architectures (clientserver) which are still dominant on the market. Packing applications involve, in fact, several distinct business processes and workflows that require multiple modelling rules and dedicated packing solvers. On the other hand we have the need to use such knowledge and solvers wrapped as a business service provided to distinct consumers. Such consumers may vary from integration with commercial systems (Enterprise Resource Planning, Warehouse Management Systems), specific systems in logistics to end users accessing a web interface application, or external applications invoking the same business logic for their internal use. The fact is that these packing services are required to be exposed and integrated with a set of different systems, users and technologies.

Living within a Service Oriented Architecture their interactivity is thus a major concern of the whole decision making application. SOA provides means for composition of services, their interaction and exposure to multiple consumers in a secure and transactional way. The design of a packing application as a composition of services provided by business modules has a consequence on the design and the implementation of optimization algorithms where sharing, common services and data are mandatory principles. Thus, some optimization services exploit capabilities provided by J2EE technologies (AJAX, JQuery, Json, et.); they are tuned making it possible to minimize the data flow (few Kbits) exchanged between the cloud server and the distant client. They take a new turn for the better in efficiency, portability and maintenance. Service-oriented architecture opens new perspectives for agile optimization services to exploit capabilities of RFID and Internet of Things in warehouses and also new applications in robotics, drones, etc.

CONCLUSION

In this article we have shown the contribution of a variety of algorithms and technologies to solve bin packing problems. Industrial application require robustness and flexibility. CP is a good candidate tor modeling and solving such problems making a clear distinction between the statement of constraints and the search procedure which can be tuned to compute good quality of solutions.

In some cases we have shown where combining CP and Metaheuristics [START_REF] Aggoun | Online assignments of containers to trains using constraint programming[END_REF] is more adequate (e.g. assigning containers to wagons). The embedding of optimization solutions in cloud computing environments leads to the organization of optimization algorithms as services opening many perspectives to Metaheuristics and dedicated heuristics to complain with real-time response constraint. Many of the applications described in this article are exploited daily. Items are modeled as regular shapes (1D, 2D and 3D). Over the past few years the demand has been evolving and unfortunately items may have irregular shapes like parts in the automotive industry. Algorithms like grouping items to find bounding boxes (3D) also has its limitation. Solvers based on continuous optimization based on Covariance Matrix Adaptation Evolution Strategy (CMA-ES) opens new perspectives to solve bin packing problems where items are complex shapes which are modeled as Bézier curves and their variants. This approach is more interesting when the solver is integrated as a backend in a declarative modeling language like MiniZinz [START_REF] Nethercote | MiniZinc: Towards a standard CP modelling language[END_REF].

ACKNOWLEDGMENTS

We would like to thank all the partners of the projects Net-WMS-2 (ANR Blanc Simi2) and SOCHART (I-TRANS) for fruitful discussions.