Francis Lazarus
email: francis.lazarus@inria.fr

Anne Verroust
email: anne.verroust@inria.fr

Feature-based shape transformation for polyhedral objects

Keywords: Computer Animation, Computer-Aided Geometric Design, Interpolation, Deformation, Shape transformation

A new technique is presented for computing transformations between polyhedral objects. Here, the correspondence and the interpolation problems are considered jointly. The approach gives the animator high level control of the shape transformation by providing natural speci cation and interaction.

Introduction

We are interested in simulating deformations of geometric shapes. The way the deformation is speci ed is important: it determines both the exibility and the usability of the deformation tool. There exists many approaches to deform a 3D shape from an object using a modelization of a time-dependent deformation process (cf. Bar84, TPBF87, PEF + 90, CJ91, BD93]). Our goal is to provide a general and intuitive tool that computes a shape transformation mainly speci ed by two shapes: the initial and the nal shapes. Such a tool is useful in both animation and design.

In animation, it is well adapted to traditional techniques such as keyframing: the objects are given at each keyframe and the in-between shapes are computed by the shape transformation. Moreover, the visual e ects of metamorphosis between two images or between two shapes are very popular. This e ect is called \morphing" [START_REF] Bishop | There is more than a way to frame a curve[END_REF][START_REF] Chen | Shape averaging and its applications to industrial design[END_REF]) and is often used in movies.

In design, it is used to create new shapes by combining two objects ( CP89, KR91]).

Our concern is concentrated in the 3-D shape transformation problem; thus, we will not describe the 2-D approaches (the main approaches can be found in CP89, BN92, SG92, SGWM93]). The 3-D methods can be classi ed into two families according to the type of information used to compute the intermediate shapes.

The volumic approaches [START_REF] Kent | Shape transformation for polyhedral objects[END_REF][START_REF] Chen | Shape averaging and its applications to industrial design[END_REF].

As the volumic information is used to compute the shape transformation, there is no restriction for the topological correspondence between the initial and the nal shapes: for example, a torus can be transformed into a sphere. Kaul and Rossignac KR91] de ne the internal points of intermediate polyhedral objects using linear interpolation based on Minkowski sums of the internal points of the two original polyhedra. The faces of the deformed shapes have a constant orientation and their vertices move on a straight line between a vertex of the initial shape and a vertex of the nal shape. Moreover, when the two original shapes are convex, these faces form the boundary of the deformed polyhedron. Hughes Hug92] proposes a method for sampled volumetric models. He de nes an interpolation between the Fourier transforms of the two models, separating the processing of high frequencies and low frequencies. The most interesting point in this approach is the underlying idea of dissociating the general shape and the details during the transformation. These two approaches are robust but there is no control over the intermediate shapes: only one shape transformation is computed from two original solids.

The boundary approaches BU89, KPC91, KPC92, Par92].

In all the cases, the original solids are polyhedral objects and they have equivalent surface manifold topologies. For these methods, the shape transformation process is decomposed in two parts:

1. a correspondence process, where the correspondence between the two surfaces is established. As the surfaces are polyhedral, the process is reduced by computing the correspondence between the topological structures of the two objects. In fact, all the approaches build a vertex/edge/face network containing the topological structures of the original objects. 2. an interpolation process: once the correspondence is established, the in-between polyhedra are computed using the corresponding vertices. This is achieved by either linearly interpolating or using a spline interpolating curve between the corresponding vertices (Hermite cubic path with end tangents set equal to the vertex normal in KPC92]). All the approaches consider the two processes separately and propose solutions to the correspondence problem.

First of all, Bethel and Uselton BU89] tackle the correspondence process mainly as a graph correspondence problem. As the geometry of the two polyhedra is not taken into account, the visual e ects produced during the shape transformation are often not intuitive in the general case and are di cult to control. For instance, as the metric relationships between points are ignored, two closed points of the rst object may be associated with two distant points of the second object. If this occurs randomly, it will result in a discontinuous transformation. Kent, Parent and Carlson KPC91, KPC92] use the geometric information of the initial and the nal shapes to build an adjacency graph containing their adjacency graphs. Their approach is very intuitive and convincing for star-shaped objects: in this case, they use a spherical projection to solve the correspondence and the interpolation processes. More precisely, the correspondence process is decomposed in two steps: ? a projection of each adjacency graph is done onto the surface of the unit sphere, ? the common adjacency graph is computed by merging these projected graphs. They generalize their approach in KPC92] mainly by extending the projection method to other classes of polyhedral models. In most of these cases the relationship between the geometry of the object and the mapping method seems less natural, since the correspondence process uses a \physical" approach which is abandoned during the interpolation process. Parent Par92] uses a recursive subdivision process to build the common adjacency graph. In his method, the user can indicate explicitly the correspondence he wants to ensure between areas of the objects (areas composed of connected sets of faces). The common adjacency graph is then obtained by recursively subdividing these areas. Our opinion is that dissociating the correspondence from the interpolation process leads to a di cult control of the deformation of the shape during the interpolation. In a sense, the shape transformation problem has to be taken as a whole, even if two problems have to be solved.

Our method can be classi ed as a boundary approach but it di ers from the previous models in several points:

First, we solve jointly the correspondence and the interpolation processes.

We do not built a common adjacency graph from the initial and the nal adjacency graphs. A common parametrized polyhedral mesh is associated to the initial and nal objects during a sampling process.

We propose to the user the control of the \global" evolution of the deformation. This control is done by specifying two corresponding skeletal structures (here 3D curves or axes) from which the parametrizations of the sampling meshes are built.

In this paper we will focus our attention on objects that are star-shaped around an axis. For these objects there exists a 3D curve inside the object, namely its axis, such that each point of the object can be attached to the axis without ambiguity. These objects can naturally be split into 3 sheets: two hemispherical parts, corresponding to the points that are attached to one extremity of the axis. one cylindrical part, corresponding to the points that are attached orthogonally to the axis, i.e. that belong to one normal cross-section of the axis. Each of these parts admits a natural parametrization, respectively spherical and cylindrical. We shall refer, in this paper, to a more general cylindrical coordinate system, where the z-axis is replaced by a 3D curve (the axis of the object).

The basic scheme

In order to perform the transformation of one object into another, we construct a new polyhedral mesh (a discretization) for each object. Actually, this discretization is one of the main point of our method: using the fact that a three dimensional object is not de ned by a particular polygonal mesh but by any mesh that captures its global geometry and its peculiar features, we choose to associate to the object the mesh corresponding to the parametrization that will be used both during the correspondence and the interpolation processes.

To do so, the user must specify two 3D axes, one for each object, such that the objects are star-shaped around them. The curves are discretized and a 3D frame is associated to each discretization point. As in LCJ93], we have chosen a rotation minimizing orthogonal frame along each 3D curve [START_REF] Bechmann | Animation through space and time based on a space deformation model[END_REF][START_REF] Hughes | Scheduled Fourier volume morphing[END_REF]. Each axis is used to de ne an appropriate parametrization of the associated object, such as described earlier. The shape transformation process will then consist in the interpolation of such parametrization. We propose here an algorithm decomposed in two main steps. Each of these steps will be explained in the following sections.

Using the axes, the sampling of the objects are computed using the natural parametrizations. This sampling process must be followed by an adjustment of the sampling vertices/edges/faces to the singular features of each object. For polyhedral objects these features are sharp edges or conic points. This adjustment is of prime importance for a good approximation of the shapes. This process will be detailed in the next section.

Then we interpolate simultaneously the parameter of each corresponding parametrization (i.e. a 3D frame for a spherical parametrization and a whole 3D curve for a cylindrical parametrization) and the local coordinates of each sampled objects with respect to these parametrizations. This will be described in section 4. At this stage it is useful to mention that the correspondence problem is automatically solved by taking the same discretization in the parameter space of two corresponding parametrizations of the objects.

3 Sampling and adjustment

Sampling

Given an object and its axis, we sample the object by computing the intersections of its boundary with a set of half-lines attached to the axis. These half-lines are built from the axis following two di erent ways: they radiate from one of the two axis extremities, in the half-space determined by the normal crosssection of the axis at this extremity. They intersect the object at two hemispherical sheets. Given one extremity, each half-line is de ned by two angular parameters, as shown in Figure 1. they belong to one normal cross-section of the axis. Here, the collection of intersection points compose the cylindrical sheet of the object. In this case, each half-line is de ned by one angular parameter and an abscissa along the axis (see Figure 1). We will call these half-lines sampling half-lines in the following.

In both cases, the parameters may vary in a rectangular domain. The sampling is thus obtained by mapping square lattices in the parameter spaces onto the object. By connecting the sampling points according to the square grid associated to the square lattice, we de ne a mesh subdivided in three sheets. The mesh is composed of rectangular faces except for the poles of the two hemispheres where triangular faces are joined at the respective pole (see Figure 1). Let us point out that for the hemisphere, one edge of the associated square domain 0; 2 ( 0; 2 ] is mapped on a single point. for each face F of the object, we compute the intersections with the sampling half-lines and F. More precisely, { given the parameters of the vertices of F, the face F is cut in three parts w.r.t. the three sheets (one or two may be empty), { inside each sheet, F is projected in the parameter space, using the parameters of its vertices previously evaluated. More precisely, the sampling parameters appearing in the projected region are located and the corresponding sampling half-lines are intersected with F. This is done by a sweep, the length of the sweep being determined by the the extremal parameters of the vertices of F: ? for the cylindrical part, we cut F with every intermediate cross section and nd the sampling parameters for each intersecting segment (see Figure 2). ? for the hemispherical parts, the process is similar but the intermediate cross sections are the planes determined by the angular parameters (cf. Figure 1).

Adjustment

If we apply this algorithm to a cube with a straight axis as in Figure 3, the resulting object will look very di erent from the original one: the sharp edges of the cube do not appear on the sampling model. One can notice that the salient features (i.e. sharp edges and conic points, for a polyhedron) are of prime importance for the appearance of an object, since they mark (with the contour lines) the discontinuities of the light re ectance function.

Figure 3: The sampling of a cube (center) and the same sampling after adjustment (right).

Hence, if we want to obtain a good reproduction of the shapes, we have to take into account the sharp edges and the conic points of the original model during the sampling process.

In order to minimize the topological change of the sampling mesh previously calculated, we choose to adjust the sampling mesh by locally deforming this mesh, moving its vertices and adding some diagonal edges to the rectangular faces. Thus we proceed in the following two steps:

we nd the sharp edges and conic points of the original object, we deform locally the sampling mesh: for each of these entities, we locate the vertices to be moved, compute their displacement and, if necessary, insert diagonal edges in the sampling mesh.

Sharp edges

A simple angle criterion between the adjacent faces to an edge allows us to select accurately the sharp edges.

A similar criterion can be applied to conic points. An interactive module is also proposed to the user in order to give him greater freedom of choice. We obtain from the above two lists of salient features: a list of sharp edges and a list of conic points.

Deforming the sampling mesh

Deforming the sampling mesh consists of tting the vertices and the edges of the sampling mesh to the sharp edges and conic points and adding diagonal edges in the sampling mesh in such a way that, at the end:

each sharp edge can be decomposed in a set of disjoint edges of the sampling mesh, the conic points belong to the vertices of the sampling mesh.

The displacements of the vertices of the sampling mesh are local and bounded w.r.t. the parametrization. More precisely, a covering of the parameter space in \deformation zones" associated to the vertices of the sampling mesh is done. These zones are rectangular neighborhoods constructed from a parallel mesh passing by medium values (in Figure 4, the deformation zones are the rectangles lled in grey surrounding the selected vertices). This parallel mesh is used to localize the vertices to be moved and adjust them on the sharp edges and the conic points as follows:

for each sharp edge E, a sweep is done (similar to the sweep previously used to sample a face): E is cut with every medium intermediate cross section and the parameters of the intersection points are used to determine the vertices to be moved. One can notice here that the number of selected vertices is not xed: in Figure 4, two vertices have been selected on section D (between intersections of E with sections C' and D') whereas only one has been selected on section B (between intersections of E with sections A' and B'). Figure 4: The selected vertices are marked in grey and the intersection points in black. The black curve is the projection of E in the parameter space. On the left, the vertices have been selected and, on the right, the sampling mesh have been deformed w.r.t. E.

When the set of vertices to be moved have been selected, each vertex V is displaced on a point P belonging to E as follows: { for the cylindrical sheet, we compute the nearest point Q belonging to E of the sampling half-line corresponding to V. If the parameters of Q belong to the deformation zone associated to V then P=Q. In the other case, P is determined by the intersections of E and the boundaries of the deformation zone associated to V.

{ in the hemispherical cases, we use spherical projection. We compute Q belonging to E such that its projection Q' is the nearest point from the projection of V on a sphere centered at the extremity (see Figure 5). We take P=Q when the parameters of Q are inside the deformation zone associated to V. If it is not the case, we proceed as in the cylindrical case, using the intersections of E and the boundaries of the deformation zone associated to V.

The cases of the conic points and of the extremities of the sharp edges are similar: the points are located w.r.t. the mesh of the deformation zones and the vertices containing these points in their deformation zone are moved on them. Once the vertices of the sampling mesh have been moved on the sharp edges, the deformed mesh may not totally cover E. This is the case when the sequence of selected vertices does not form a subpath of the sampling mesh. In this case, a diagonal edge must be added each time there is a \hole" in the sequence of If the neness of the sampling is not su cient, two sharp edges may deform the same vertex of the sampling mesh. As taking into account the two edges to compute the nal position of the vertex has no sense, we move the vertex even though it has been displaced before. Since the position of the conic points are of prime importance, we choose to examine rst the sharp edges (w.r.t. their order of appearance in the list of sharp edges) and then the conic points. The order of the edges inside the list of sharp edges, which determines the relative importance of the salient features, can be changed via an interactive tool. Anyway, the co-occurrence of several constraints over one vertex of the sampling mesh can be a criterion to re ne the mesh.

E Q V Q' axis

Interpolation

We must ensure, before the interpolation process, that the two sampled meshes of the objects have the same topology. This was true before we adjusted these meshes on the sharp edges. But in the previous section we saw that diagonals could be added to the meshes, thus we have to report in each mesh the same diagonals. This is done by splitting the rectangular faces where diagonal edges have been added in two or four triangular faces:

if only one diagonal appears in one of the two sampling meshes, then it cuts the rectangular face in two, if two diagonal edges have been added (one in each sampling mesh), we have to insert a new vertex and split the face into four triangular faces.

As we have already stated, the interpolation is performed at two levels: we interpolate simultaneously the coordinates systems, which are de ned by the axes, and the local coordinates of each object in these coordinates systems.

Axes interpolation

We consider the discretized curve as a polygonal line rather than a di erentiable entity. We discretize the two curves with the same number of points since they correspond to sampled values in the parameter space. Simply it is necessary to nd a good interpolation process between the two polygonal lines. Henceforth, the terms axis, curve or polygonal line will designate the same entity.

As noticed in SGWM93], a linear interpolation of the two polygonal lines seems to be inadequate. For instance, the interpolation of two parallel segments oriented in opposite directions collapses for some interpolation value; this is not acceptable.

Actually we consider the axis as a moving frame, which means that a 3D frame is associated to each point of the axis and we interpolate the whole structure (see Figure 6).

γ i γ i+1 γ i+2 R i R i+1 R i+2 γ
Figure 6: A discretized axis with its frames A moving frame associated to a curve is a couple ( ; R), where R is an orthonormal moving frame (e 1 ; e 2 ; e 3 ) such that e 1 (t) has the same direction as the tangent 0 (t) at the point (t). The basic idea is to interpolate the evolution of these frames rather than the frames themselves. Thus we interpolate the relative coordinates and the relative rotation between two consecutive points and frames of the polygonal line. Let's note ( 1 ; :::; n ) (resp. (R 1 ; :::; R n )) the digitized points (resp. frames) of the curve , then we have:

i+1 = i + R i X i+1 (1) R i+1 = r i+1 R i (2)
where X i+1 are the local coordinates of i+1 in ( i ; R i ), and r i+1 is the rotation between R i and R i+1 . By interpolating the X i and the r i of the two axes, we can construct for any intermediate time t a polygonal line: starting from a point on the curve (for instance an extremity) we recursively compute the next point using equations (1) and (2) where X i+1 and r i+1 have been replaced by the interpolated values. We use linear interpolation for the coordinates X i and quaternion interpolation for the rotation r i (more exactly for the corresponding quaternion q i ). We follow the notations of Sho87] in equation (4). Hence we have:

X i (t) = (1 t)X i;1 + tX i;2 (3) q i (t) = Slerp(q i;1 ; q i;2 ; t) = q i;1 (q 1 i;1 q i;2 ) t (4) where subscripts 1 and 2 account for the initial and the nal curves.

It is clear from the description of the whole process that the interpolation is continuous. Let us point out here that continuity does not only mean that a slight change in the polygonal lines would result in a slight change in their interpolation, but also that the process applied to digitized curves, whose steps of discretization go towards zero, would result in continuous interpolating curves.

Moreover, equation (1) shows that X i has a second and a third components close to zero since it is almost parallel to the tangent curve. One can note from equation (3) that the interpolation is linear at a local level which holds true for the complete length of the curve. For an exact linear interpolation of the length, we could have interpolated separately the length and the orientation of X i . However we found this of no practical use. This curve interpolation leads to good visual results as one can see in Figure 11. Moreover the computation is simple and fast: the whole interpolation process is fully interactive on a Silicon Graphics workstation.

Vertex interpolation

A linear interpolation of the local coordinates is performed between two corresponding vertices of each sampled object. Hence, in the local frame, the tangent to the trajectory of a vertex is radial to the axis at the beginning and the end of the animation. In order to preserve this property for the adjusted vertices, we compute the cubic spline de ned by the vertices and their associated radial directions.

Thus the vertex path is obtained by composing this movement and the movement of the associated 3D a ne frame. This results in a complex movement which could not have been designed directly !

As one can see here, interpolation and correspondence (a consequence of the sampling) are strongly related: in both cases, the main role is played by the sampling half-lines attached to the axis.

A sample of animations

Let us comment our method using a sample of animations.

Figures 7,8 and 9 illustrate the role played by the axes to control the shape transformation. The original objects are a cube and a cone. Two axes are de ned: a vertical line segment and an oblique curve. These axes are discretized using seven frames and, for each planar section, sixteen angular parameters are used to sample the objects. They are drawned in bold for each sampled model associated to an initial or a nal object.

Figure 7: The same axis is used to sample the two original objects.

In Figure 7, the same axis, a vertical line segment, is used to sample the initial and the nal objects. In the second shape transformation the axis has been translated inside the cube. This way to control the shape transformation is similar to what Kent, Parent and Carlson KPC91] showed by moving the center of projection inside one of the star-shaped objects. When taking two di erent axes (here an oblique curve and the vertical line segment used in Figure 7) to sample the initial and the nal objects, the axis interpolation notably in uences the shape transformation. In Figures 8 and9, the two shape transformations di er: on the top, the axis rotate in clockwise order from the left to the right whereas on the bottom, the rotation is made in counterclockwise order. These rotations lead to di erent movements of the salient features and thus to di erent intermediate shapes, as one can see on Figure 8. Three objects have been sampled and three transformations have been constructed in Figure 10. Here, the shape transformations are showned in vertical columns. The parametrization associated to the sh is nearly cylindrical: its associated axis begins and ends near the mouth and the tail. Then, when transformed in a cube, the mouth and the tail form the bases, as it is visible on the intermediate shape of the left column. Figure 11 shows two di erent views of the same shape transformation process. Here, the two original objects are generalized cylinders around two di erent curves. This example illustrates the visual e ect produced by our axis interpolation method.

In Figure 12, the sh has been deformed into a sickle. The axis of the sh is a straight line whereas the axis associated to the sickle is bended on the sharp part of the instrument. One can see the variety of intermediate shapes and the appearance/vanishing of the salient features in the in-between shapes.

Conclusion

We have proposed a method for computing a continuous deformation between two polyhedral shapes. The proposed method achieved: A natural way to solve jointly the correspondence and the interpolation problem.

Interactive speci cation and control of the shape transformation. The use of axes to specify the transformation is intuitive: it can be viewed as the speci cation via a \skeleton" of the global deformation of the shapes. It is a natural extension of the rst approach of Kent, Parent and Carslon KPC91].

Good visual e ects during the transformation, as shown on the examples of the previous section. An e cient way to compute transformation between shapes. The use of an intermediate polyhedral mesh reduces the computational time: for each object, it depends on the number of vertices and faces of the object and on the neness of the discretization. The user may obtain approximate sketches of the deformation process using a rough discretization during the sampling. Moreover, our approach can be well adapted to compute transformations between a sequence of shapes.

In this case, all the objects have to be sampled with the same discretization for the axes and the angular parameter. The axis interpolation and the local interpolation have to be extended into spline interpolation.

A general shape transformation adapted to any objects that are star-shaped around an axis. In particular, our method can be simply extended to implicit surfaces WBB + 90]. For, it is su cient to be able to compute the intersection of a sampling half-line with the surface of the object, which is an easy task in this case. This algorithm has been integrated in ACTION3D, a general interactive modeling system developed jointly by SOGITEC and INRIA. The whole process (sampling and interpolation) is interactive on a Silicon Graphics workstation.

Work in progress include attempts to generalize the method to other families of objects, i.e., to generalize the sampling and parametrization to another kind of skeleton.

Figure 1 :

 1 Figure 1: The three sheets

Figure 2 :

 2 Figure 2: Image of a face in the parameter space

Figure 5 :

 5 Figure 5: Computation of Q in the hemispherical case.

Figure 8 :

 8 Figure 8: The in uence of the axis interpolation in the shape transformations.

Figure 9 :

 9 Figure 9: The same shape transformations in wire-frame.

Figure 10 :

 10 Figure 10: in the horizontal middle row, three intermediate shapes interpolated from the top and bottom original shapes

Figure 12 :

 12 Figure 12: Shape transformation between a sickle and a sh.


L Wu

Philippe De Reffye

François Xavier Le Dimet

Baogang

Optimization of Source-Sink Relationship based on a Plant Functional-Structural Model: A Case study on Maize

published or not. The documents may come   L'archive ouverte pluridisciplinaire


Jakob Puchinger
email: puchinger|raidl@ads.tuwien.ac.at

Günther R Raidl

Bringing Order into the Neighborhoods: Relaxation Guided Variable Neighborhood Search *

Keywords: Variable Neighborhood Search, Integer Linear Programming, Multidimensional Knapsack Problem

In this article we investigate a new variant of Variable Neighborhood Search (VNS):

Relaxation Guided Variable Neighborhood Search. It is based on the general VNS scheme and a new Variable Neighborhood Descent (VND) algorithm. The ordering of the neighborhood structures in this VND is determined dynamically by solving relaxations of them. The objective values of these relaxations are used as indicators for the potential gains of searching the corresponding neighborhoods. We tested this new approach on the well-studied multidimensional knapsack problem. Computational experiments show that our approach is beneficial to the search, improving the obtained results. The concept is, in principle, more generally applicable and seems to be promising for many other combinatorial optimization problems approached by VNS.

Introduction

The general Variable Neighborhood Search (VNS) heuristic including Variable Neighborhood Descent (VND) (Hansen andMladenović, 1999, 2003) is a relatively young metaheuristic concept that has successfully been applied to several combinatorial optimization problems.

It relies on the principle of searching different neighborhood structures in a systematic way.

In [START_REF] Hansen | A Tutorial on Variable Neighborhood Search[END_REF] the following significant question is raised: What is the best order in applying the neighborhoods? Usually, neighborhoods are sorted according to increasing complexity of searching for the best moves, which often, but not always, corresponds to the size of neighborhoods. Sometimes, however, such an order might be hard to estimate or even be misleading. In particular, a prespecified, fixed order might in general not be ideal during all phases of the optimization process. We will introduce a new dynamic approach, in which the neighborhoods are always processed in an order that depends on estimated improvement-potentials.

We call this enhanced variant of general VNS Relaxation Guided Variable Neighborhood Search. It follows the standard VNS scheme but uses a new VND algorithm. As a specific technique for determining improvement-potentials for the neighborhoods considered in VND, we employ (Integer) Linear Programming (ILP/LP) techniques. Note that combining metaheuristics and ILP techniques is a highly promising research area for its own [START_REF] Dumitrescu | Combinations of Local Search and Exact Algorithms[END_REF]Puchinger and Raidl, 2005a), our approach can also be seen as a contribution to this direction.

The VND is guided by always sorting the neighborhoods according to estimations of the improvement-potentials depending on the current solution. For each neighborhood this potential is determined by quickly solving a relaxation. Searching the neighborhoods in this order is expected to increase solution quality and/or to speed up VNS.

In order to evaluate this new approach, we use the Multidimensional Knapsack Problem (MKP). It is a well-studied, strongly NP-hard combinatorial optimization problem occurring in many different application areas. In the last decades a multitude of exact and metaheuristic algorithms were developed for the MKP, some of them being very complex. We examine the benefits of relaxation guided VNS in comparison to standard version of VNS on this problem.

Although it is not our primary goal to compete with the leading algorithms for large and complex MKP instances, we get relatively good results that are in general not far from the best known solutions.

In the next section we present the general scheme of Relaxation Guided VNS. Section 3 introduces the MKP in detail. We then describe the neighborhood structures used for the MKP, together with computational experiments comparing standard VNS and relaxation guided VNS in Section 4. Some extensions of this approach using more neighborhood structures and further experiments are presented in Section 5. We close with some concluding remarks and an outlook on future work.

Relaxation Guided VNS

Relaxation Guided VNS (RGVNS) follows the the general VNS scheme (Hansen andMladenović, 1999, 2003) and incorporates an improved VND, which we call Relaxation Guided Variable Neighborhood Descent (RGVND). Let us assume that the neighborhood structures N 1 , . . . , N kmax are to be used within VND. An important question, which often has a substantial impact on the algorithm's performance, is the order in which the neighborhoods are to be considered. Occasionally, this ordering is straight-forward, since the neighborhoods are nested, such that neighborhood N i always fully contains N -i -1, ∀i = 2, . . . , k max .

An example are classical k-opt neighborhoods for the traveling salesman problem or more general k-exchange neighborhoods. However, a particular strength of VND/VNS often lies in the combination of different types of neighborhoods which do not form such a monotone sequence of extensions.

In such cases, rules of thumb are typically used for defining a meaningful order. Examples are: (a) search smaller neighborhoods first; (b) search neighborhoods first for which faster algorithms exist; and (c) search neighborhoods first which are generally considered to be "more promising". In many situations, however, these strategies do not lead to the overall best ordering. Furthermore, the ordering that is best suited in a particular situation might in general depend on the current solution, i.e. there is not a single ordering which is optimal at any time of the heuristic search. With the exception of our preliminary results we published in Puchinger and Raidl (2005b), we are not aware of any previous work where the ordering of the neighborhoods is determined in an automatic way and, especially, is adapted during the search.

The central point of our extended variant of VND is that we automatically steer the order in which the neighborhood structures are processed in dependence of estimated improvementpotentials. These potentials are devised by quickly solving a relaxation of each neighborhood structure. If the improvement-potentials are able to predict to a certain degree which neighborhoods are more promising and which are less, considering the neighborhoods in this order will yield better results and/or shorter running times.

In the following we will consider maximization problems; minimization problems can be treated analogously.

Assume that we are given a Combinatorial Optimization Problem (COP) defined as

z COP = max{f (x) | x ∈ S}, ( 1
)
with S being a finite set of solutions and f (x) : S → R an objective function. In analogy to the relaxation of an ILP (Wolsey, 1998), we introduce the following formal definition of a relaxation of a COP.

Definition 1 A relaxation R of COP is any maximization problem defined as

z R = max{f R (x) | x ∈ S R } (2)
with the following properties:

(i) S ⊆ S R (ii) f (x) ≤ f R (x), ∀x ∈ S.
The following evident result (Wolsey, 1998) yields a bound for COP.

Proposition 1 If R is a relaxation of COP, z R ≥ z.
If such a relaxation is chosen appropriately, it is substantially faster to calculate its optimal solution than the solution to the original COP. An example is the widely used LP relaxation of an ILP formulation for a COP, which can be solved in polynomial time. Other generic relaxation methodologies from the ILP domain, which are commonly applied with much success, are e.g. the surrogate relaxation, in which constraints are combined, and the Lagrangian relaxation, where constraints are transformed into penalty terms added to the objective function.

It is a prerequisite for our RGVND scheme that the used relaxations can be solved to optimality much faster than their corresponding original neighborhood structures on one side,

Algorithm 1: Relaxation Guided VND (RGVND) Input: A feasible solution x l ← 1 π = DetermineOrderOfNeighborhoods(x) repeat Find the best neighbor x * ∈ N π(k) (x) | f (x * ) ≥ f (x ) ∀x ∈ N π(k) (x) if f (x * ) > f (x) then x ← x * k ← 1 π = DetermineOrderOfNeighborhoods(x) else k ← k + 1 until k = k max return x
and that the obtained bounds are "good" in the sense that they give some indication for the objective values of the best solution in the neighborhood.

As a second precondition we assume that the used neighborhoods are not nested since this would lead to trivial orderings and render our approach meaningless, i.e. N k ⊆ N k and

N k ⊆ N k holds for any N k , N k with k = k .
Algorithm 2: DetermineOrderOfNeighborhoods(x)

for k = 1, . . . , k max do Solve N R k (x) yielding solution value z R k Sort π = (1, . . . , k max ) according to decreasing z R k return π
In Algorithm 1 the pseudocode of RGVND is given. The significant differences to the standard VND scheme, as described in Hansen andMladenović (1999, 2003), are the calls of function DetermineOrderOfNeighborhoods(x) in lines 2 and 8. This function determines the order of the neighborhood structures by first solving their relaxations yielding objective values z R k , and then sorting the neighborhoods according to decreasing z R k . Ties are broken arbitrarily or according to some static heuristic rules.

The Multidimensional Knapsack Problem

The Multidimensional Knapsack Problem (MKP) is a well-studied, strongly NP-hard combinatorial optimization problem occurring in many different applications. It can be defined by the following ILP:

(MKP) maximize z = n j=1 p j x j (3) subject to n j=1 w ij x j ≤ c i , i = 1, . . . , m (4)
x j ∈ {0, 1}, j = 1, . . . , n.

(5)

There are n items with profits p j > 0 and m resources with capacities c i > 0 given. Each item j consumes an amount w ij ≥ 0 from each resource i. The goal according to (3) is to select a subset of items with maximum total profit; chosen items must, however, not exceed resource capacities, see ( 4). The 0-1 decision variables x j indicate which items are selected.

The MKP first appeared in the context of capital budgeting (Lorie and [START_REF] Savage | Three problems in capital rationing[END_REF][START_REF] Manne | On the solution of discrete programming problems[END_REF]. A comprehensive overview of practical and theoretical results for the MKP can be found in the monograph on knapsack problems by [START_REF] Kellerer | Knapsack Problems[END_REF]. A recent review of the MKP was given by [START_REF] Fréville | The multidimensional 0-1 knapsack problem: An overview[END_REF]. Besides exact techniques for solving small to moderately sized instances based on dynamic programming [START_REF] Gilmore | The theory and computation of knapsack functions[END_REF]Weingartner and Ness, 1967) and branch-and-bound [START_REF] Shih | A branch and bound method for the multiconstraint zero-one knapsack problem[END_REF][START_REF] Gavish | Efficient algorithms for solving the multiconstraint zero-one knapsack problem to optimality[END_REF], many kinds of metaheuristics have already been applied to the MKP.

To our knowledge, the method currently yielding the best results, at least for commonly used benchmark instances, was described by [START_REF] Vasquez | A Hybrid Approach for the 0-1 Multidimensional Knapsack Problem[END_REF] and has recently been refined by [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF]. It is a hybrid approach based on tabu search. The search space is reduced and partitioned via additional constraints, thereby fixing the total number of items to be packed. Bounds for these constraints are calculated by solving a modified LP-relaxation. For each remaining part of the search space, tabu-search is independently applied, starting with a solution derived from the LP-relaxation of the partial problem. The improvement described in [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF] lies mainly in an additional variable fixing heuristic.

Various other metaheuristics have been described for the MKP [START_REF] Glover | Critical event tabu search for multidimensional knapsack problems[END_REF][START_REF] Chu | A genetic algorithm for the multiconstrained knapsack problem[END_REF], including several variants of hybrid evolutionary algorithms (EAs); see [START_REF] Raidl | Empirical Analysis of Locality, Heritability and Heuristic Bias in Evolutionary Algorithms: A Case Study for the Multidimensional Knapsack Problem[END_REF] for a recent survey and comparison of EAs for the MKP.

The authors of the current article have developed different approaches, combining exact methods with metaheuristics for solving the MKP (Puchinger et al., 2005[START_REF] Puchinger | The Core Concept for the Multidimensional Knapsack Problem[END_REF][START_REF] Puchinger | Combining Metaheuristics and Integer Programming for Solving Cutting and Packing Problems[END_REF].

4 Relaxation Guided VNS for the MKP

We now focus on the problem-specific details of our RGVNS implementation for the MKP, introducing the used neighborhoods and their relaxations, and present results for indicating the effectiveness of the new approach in comparison to standard VNS.

Representation and Initialization

Solutions are directly represented by binary strings, and all our neighborhoods are defined on the space of feasible solutions only. Previous studies such as [START_REF] Raidl | Empirical Analysis of Locality, Heritability and Heuristic Bias in Evolutionary Algorithms: A Case Study for the Multidimensional Knapsack Problem[END_REF] indicate that for the MKP, this seems to be in general a wise decision, since approaches also considering infeasible solutions to a larger degree turned out to be typically less effective.

We denote by I 1 (x f ) = {j | x f j = 1} the index-set of the items contained in the knapsack of a current solution x f and by

I 0 (x f ) = {j | x f j = 0} its complement.
The initial solution for our VNS is generated using a greedy first-fit heuristic, considering the items in a certain order that is determined by sorting the items according to decreasing values x j of the solution to the MKP's LP-relaxation; ties are broken randomly. See [START_REF] Raidl | Empirical Analysis of Locality, Heritability and Heuristic Bias in Evolutionary Algorithms: A Case Study for the Multidimensional Knapsack Problem[END_REF] for more information on this randomized greedy heuristic.

ILP-Based Neighborhoods

Instead of simple bit-flip or k-exchange neighborhoods that would inherently lead to infeasible solutions, we consider more sophisticated neighborhoods based on the MKP's ILP model.

We want to force a certain number of items of the current feasible solution x f to be removed from or added to the knapsack. This is realized by adding neighborhood-defining constraints depending on x f to the ILP formulation of the MKP.

In the first neighborhood, ILP-Remove-and-Fill IRF (x f , κ), we force precisely κ items from I 1 to be removed from the knapsack and any combination of items from I 0 is allowed to be added to the knapsack as long as the solution remains feasible. This is accomplished by adding the following equation to (3)-( 5):

j∈I 1 (x f ) x j = j∈I 1 (x f ) x f j -κ. ( 6
)
In the second neighborhood, ILP-Add-and-Remove IAR(x f , κ), we force precisely κ items not yet packed, i.e. from I 0 , to be included in the knapsack. To achieve feasibility any combination of items from I 1 may be removed. This is achieved by adding the following equation to (3)-( 5):

j∈I0(x f ) x j = κ. ( 7
)
As relaxations IRF R (x f , κ) and IAR R (x f , κ) for RGVND we use the corresponding LPrelaxations in which the integrality constraints ( 5) are replaced by 0 ≤ x j ≤ 1, j = 1, . . . , n.

Note that depending on the specific instance's characteristics, both neighborhoods may become quite large even for κ = 1. Nevertheless, the LP-relaxations can be solved to optimality very quickly by means of standard LP algorithms. For searching the (integer) neighborhoods we use a general purpose ILP-solver -CPLEX 9.0 in our case -with a certain time limit.

Relaxation Guided RGVND

Our Relaxation Guided Variable Neighborhood Decent (RGVND) for the MKP uses the previously defined neighborhoods IRF (x f , κ) and IAR(x f , κ) with κ = 1, . . . , κ max , where κ max is a prespecified upper limit on the number of items that shall be removed or added.

It otherwise follows the general algorithm already presented in Section 2, i.e. the order in which the neighborhoods are considered is always dynamically determined by solving the LP-relaxations and sorting the neighborhoods according to decreasing solution values. Ties are broken by considering neighborhoods with smaller κ earlier.

Shaking

When RGVND terminates with a solution that is locally optimal with respect to IRF (x f , κ)

and IAR(x f , κ) for κ = 1, . . . , κ max , shaking is performed within the VNS framework. Our shaking flips l ≥ 1 different, randomly selected variables x j of the current solution and applies the following fast repair and local improvement procedures that have originally been proposed

by [START_REF] Chu | A genetic algorithm for the multiconstrained knapsack problem[END_REF].

Both, repair and local improvement, are greedy first-fit strategies and guarantee that any resulting candidate solution lies at the boundary of the feasible region, where optimal solutions are always located. The repair procedure considers all items in a specific order Π and removes selected items (x j = 1 → x j = 0) as long as any capacity constraint is violated.

Local improvement works vice-versa: It considers all items in the reverse order Π and includes

items not yet appearing in the solution as long as no capacity limit is exceeded.

Crucial for these strategies to work well is the choice of the ordering Π. Items that are likely to be selected in an optimal solution should appear near the end of Π. As described by [START_REF] Chu | A genetic algorithm for the multiconstrained knapsack problem[END_REF], Π is determined according to pseudo-utility ratios (efficiency values)

u j = p j m i=1 r i w ij , ( 8
)
where we set the relevance factors r i to the dual variable values of the solution to the LPrelaxation of the MKP. [START_REF] Puchinger | Combining Metaheuristics and Integer Programming for Solving Cutting and Packing Problems[END_REF] and [START_REF] Puchinger | The Core Concept for the Multidimensional Knapsack Problem[END_REF] evaluated in a slightly different context various efficiency measures and concluded that the dual variable values are in general a very good choice.

As usual in general VNS, l runs from 1 to some l max and is reset to 1 if an improved solution is found.

Experimental Comparison

We compare RGVNS to standard VNS with a fixed neighborhood ordering and a VNS variant in which the neighborhoods are always considered in random order (RandVNS).

All three strategies use the same ILP-based neighborhoods IRF (x f , κ) and IAR(x f , κ), for κ = 1, . . . , κ max = 10 for their VND. In standard VNS the neighborhood-ordering is 10), IAR(x f , 10), i.e. the neighborhoods are sorted according to increasing complexity and we always switch between IRF and IAR. For shaking the maximum number of variable flips l max was set to n to not restrict the theoretical possibility of reaching any solution in the search space.

IRF (x f , 1), IAR(x f , 1), IRF (x f , 2), IAR(x f , 2), . . . IRF (x f ,
The algorithms were implemented in C++ using the ILP-solver CPLEX 9.0. The neighborhoods were not always fully explored to provable optimality since this turned out to be less effective. In our experiments, we allowed a maximum of two seconds per neighborhood evaluation. If CPLEX did not terminate within this time, the so far best neighboring solutions was used. The total run-time given to the algorithms was limited to 500 seconds. All experiments were performed on a 2.4GHz Intel Pentium 4 machine.

As in many previous publications for the MKP, we use Chu and Beasley's benchmark library for our experiments, which is available at Beasley's OR-Library1 . In particular, we consider here the largest instances with n = 500 items, m ∈ {5, 10, 30} constraints, and tightnessratios α ∈ {0.25, 0.5, 0.75}. Each instance has been generated randomly such that c i = α • n j=1 w ij for all i = 1, . . . , m. In order to evaluate the results of our experiments statistically, we performed 30 independent runs on the first instances for each of the nine parameter combinations.

For a solution with objective value z, we measure its quality by the %-gap = z LP -z /z LP • 100%, where z LP is the solution value of the MKP's LP-relaxation. Table 1 lists average and median percentage gaps together with standard deviations (in parentheses) for the final solutions of VNS, RandVNS, and RGVNS. We also performed Wilcoxon rank sum tests to see how significant the observed differences in the results of the three algorithms are. Columns p VNS,RGVNS and p RandVNS,RGVNS show obtained error probabilities for the hypotheses that the mean gap of VNS is larger than the one of RGVNS and the mean gap of RandVNS is larger than the one of RGVNS, respectively.

For seven out of the nine instances RGVNS yields significantly better results than the classical VNS approach with fixed neighborhood ordering. For one of the test cases (m = 5, α = 0.75) all obtained results were equal, whereas for the (m = 10, α = 0.75) case the standard approach obtained on average slightly better results than RGVNS, but without statistical significance. When comparing RGVNS to RandVNS, one can observe that in four cases RGVNS was significantly better, in one case results were identical, in two cases RandVNS was better, and in two cases no conclusion can be drawn. As expected the random ordering yields better results than the fixed order, but the relaxation guided approach outperforms both of the naive orderings.

Extending RGVNS for the MKP

The pure RGVNS approach of the previous section yields relatively good results, but it cannot directly compete with today's state-of-the-art metaheuristics for the MKP. We therefore extend the RGVNS described in the previous section by some standard neighborhood structures which can be searched in faster ways. These simpler neighborhood structures are {5, 10, 30} constraints, and tightness-ratios α = {0.25, 0.5, 0.75}. Ten different instances exist for each parameter combination yielding 90 instances in total. As before, the experiments were performed on a 2.4GHz Intel Pentium 4 machine, and each run was terminated after 500s of CPU-time. The neighborhoods are ordered as follows: N 1 :=SWP (x f ), N 2 :=RGF (x f , 1),

N 3 :=AGR(x f , 1), N 4 :=RGF (x f , 2), N 5 :=AGR(x f , 2).
In Table 2 we show results of the following algorithm variants: VNS with neighborhoods N 1 to N 3 only (VNS-N 1-3 ), VNS with neighborhoods N 1 to N 5 only (VNS-N 1-5 ), VNS with the ILP-based neighborhoods (VNS), RGVNS with the ILP-based neighborhoods (RGVNS), RGVNS with additionally N 1 to N 3 , and RGVNS with additionally N 1 to N 5 . Listed are average percentage gaps and numbers #-best indicating for how many instances each approach could find the overall best solutions among the ten instances per group. In the last column (VV05) the percentage gaps of the approach presented in [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF] are shown, note that run-times of more than 80 hours were needed for some instances.

We can observe a clear performance difference between N 1-3 , N 1-5 , and the ILP-based neighborhoods. RGVNS+N 1-3 yields the best total average percentage gap. Furthermore the two RGVNS+ methods yield the highest number of best solutions found.

The fact that the RGVNS variants in general yield the best average percentage gaps over all categories together with the results from the previous section clearly documents the benefits of sorting the neighborhoods according to a dynamically determined potential for improvement.

The RGVNS approach presented here is not able to compete with the approach from [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF] in terms of solution quality. Enhancing RGVNS with other ideas such as the MKP cores or parallel-execution of different algorithms (Puchinger et al., 2005[START_REF] Puchinger | The Core Concept for the Multidimensional Knapsack Problem[END_REF][START_REF] Puchinger | Combining Metaheuristics and Integer Programming for Solving Cutting and Packing Problems[END_REF] can strongly increase its capabilities, and results close to the ones obtained in [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF] have been achieved in significantly shorter running-times.

Conclusion

We presented an extension of traditional variable neighborhood search: Relaxation Guided Variable Neighborhood Search (RGVNS). The order in which the neighborhoods are investigated is dynamically determined by estimating their improvement-potential using quickly determined solutions to relaxations. This idea seems to be particularly useful if the order of the neighborhoods is not obvious and their relaxations can be quickly solved and yield relatively tight bounds. We tested this approach on standard benchmark instances of the multidimensional knapsack problem. The results obtained in our computational experiments show a clear advantage of RGVNS compared to VNS without guidance. In the future we want to apply RGVNS on other problems in order to gain further experience with the presented approach. Also, we will study further possibilities of relaxations, such as surrogate and Lagrangian relaxations.

  Weingartner, H. M. and D. N. Ness (1967). "Methods for the solution of the multidimensional 0/1 knapsack problem." OperationsResearch 15, 83-103. Wolsey, L. A. (1998). Integer Programming. Wiley-Interscience.

Table 1 :

 1 Comparison of VNS, RandVNS, and RGVNS; listed are average and median percentage gaps, standard deviations in parentheses, and error probabilities of Wilcoxon rank sum tests indicating the significance of differences.

			VNS	RandVNS	RGVNS	p VNS,RGVNS p RandVNS,RGVNS
	m	α	mean	median	mean	median	mean	median	W-test	W-test
	5 0.25	0.091	0.096	0.088	0.088	0.082	0.076	0.04	< 0.01
			(0.011)		(0.006)		(0.010)
		0.5	0.042	0.041	0.037	0.036	0.034	0.034	< 0.01	< 0.01
			(0.005)		(0.004)		(0.000)
		0.75 0.023	0.023	0.023	0.023	0.023	0.023	n.a.	n.a.
			(0.000)		(0.000)		(0.000)
	10 0.25	0.251	0.251	0.229	0.236	0.212	0.204	< 0.01	< 0.01
			(0.018)		(0.025)		(0.016)
		0.5	0.115	0.108	0.105	0.108	0.108	0.108	< 0.01	0.42
			(0.009)		(0.009)		(0.007)
		0.75	0.073	0.075	0.071	0.070	0.075	0.079	0.19	< 0.01
			(0.003)		(0.003)		(0.005)
	30 0.25	0.685	0.686	0.639	0.642	0.635	0.614	< 0.01	0.383
			(0.047)		(0.025)		(0.034)
		0.5	0.291	0.304	0.256	0.244	0.272	0.277	< 0.01	< 0.01
			(0.032)		(0.019)		(0.022)
		0.75	0.152	0.154	0.139	0.0136	0.131	0.131	< 0.01	< 0.01
			(0.016)		0.009		(0.000)

Table 2 :

 2 Average results of the different approaches, using the whole Chu and Beasley 500variables instance set.

	VNS-N 1-3	VNS-N 1-5	VNS		RGVNS	RGVNS+N 1-3 RGVNS+N 1-5 VV05
	5 0.25 0.124	0	0.109	0	0.113	2	0.090	2	0.088	4	0.087	6	0.074
	0.5 0.065	0	0.053	0	0.049	0	0.042	7	0.043	5	0.042	4	0.038
	0.75 0.041	1	0.029	2	0.032	1	0.026	7	0.027	6	0.026	8	0.024
	10 0.25 0.357	0	0.293	0	0.271	1	0.234	5	0.230	6	0.232	4	0.174
	0.5 0.180	0	0.137	0	0.131	0	0.108	4	0.108	3	0.103	8	0.082
	0.75 0.103	0	0.083	0	0.084	1	0.069	7	0.069	6	0.072	5	0.057
	30 0.25 0.890	0	0.825	0	0.716	0	0.609	5	0.595	6	0.615	4	0.482
	0.5 0.414	0	0.332	0	0.310	0	0.265	4	0.263	6	0.268	3	0.210
	0.75 0.232	0	0.216	0	0.189	1	0.167	3	0.168	3	0.167	3	0.135
	Average 0.267	0.1	0.231	0.2	0.211	0.7	0.179	4.9	0.177	5.0	0.179	5.0	0.142

m α %-gap #-best %-gap #-best %-gap #-best %-gap #-best %-gap #-best %-gap #-best %-gap

http://people.brunel.ac.uk/˜mastjjb/jeb/info.html

* This work is supported by the European RTN ADONET under grant 504438.

1

not ordered according to relaxations, but explored in a fixed order, before relying on the relaxation based ordering for calling IRF (x f , κ) and IAR(x f , κ) with κ = 1, . . . , κ max . Furthermore, we introduce an additional parameter β max that limits the number of ILP-based neighborhoods to be explored within one call of RGVND. In this way, shaking is performed more frequently and diversification is emphasized. If, for example, we choose β max = 10 and κ max = 10, a total of 20 ILP-based neighborhoods are sorted according to their LPrelaxations, but only the 10 most promising are actually evaluated in the current call of VND.

Swap Neighborhood

The first neighborhood we use is a simple swap SWP (x f ), where a pair of items (x f i , x f j ) | i ∈ I 1 and j ∈ I 0 is exchanged, i.e. x f i := 0 and x f j := 1. Infeasible solutions are discarded. Note that this neighborhood is contained in both, IRF (x f , 1) and IAR(x f , 1). Therefore, including it in the relaxation based ordering does not make sense. Its main advantage is its simplicity and speed.

Greedy Neighborhoods

Based on the ideas of [START_REF] Chu | A genetic algorithm for the multiconstrained knapsack problem[END_REF] and as another simplification of IRF and IAR but an extension of SWP, we define two additional neighborhoods based on greedy concepts.

In the first case, the Remove-and-Greedy-Fill neighborhood RGF (x f , κ), κ items are removed from x f , i.e. a κ-tuple of variables from I 1 (x f ) are inverted. The resulting solution is then locally optimized using the greedy first-fit heuristic from Section 4.4.

In the second case, the Add-and-Greedy-Repair neighborhood AGR(x f , κ), κ items are added to x f , i.e. κ variables from I 0 (x f ) are inverted. The resulting solution, which is usually infeasible, is then repaired and locally improved using the greedy algorithms from Section 4.4.

Computational Experiments

In these experiments, we combined the simpler neighborhoods, which were explored using a best improvement strategy, with the ILP-based neighborhoods. This time, we used the whole set of instances of Chu and Beasley's benchmark library with n = 500 items, m ∈