Thibault Tricard

Jimmy Etienne

Cédric Zanni

Sylvain Lefebvre

A brick in the wall: Staggered orientable infills for additive manufacturing

Keywords: CCS Concepts:, Computing methodologies → Mesh models; Texturing additive manufacturing, FDM, infill, procedural generation

. Our technique generates dense planar infill trajectories, precisely following an input direction field. Through the anisotropy of the deposition process, the direction field controls the appearance and the physical properties of the 3D printed object. The trajectories are arranged in a staggered layout across layers. Left: A 3D printed gear where the direction of the trajectories is parameterized to adapt to the functionality of the part. Here, trajectories are mostly circular in the rim and hub while being aligned with the spokes. Note the staggered layout in the side views. Right: The anisotropy of the deposition results in anisotropy in the specular reflectance.

Here, this is controlled to pattern the appearance of an otherwise flat part, resulting in a brushed metal effect.

Additive manufacturing is typically conducted in a layer-by-layer fashion. A key step of the process is to define, within each planar layer, the trajectories along which material is deposited to form the final shape. The direction of these trajectories triggers an anisotropy in the fabricated parts, which directly affects their properties, from their mechanical behavior to their appearance. Controlling this anisotropy paves the way to novel applications, from stronger parts to controlled deformations and surface patterning.

This work introduces a method to generate trajectories that precisely follow an input direction field while simultaneously avoiding intra-and inter-layer defects. Our method results in spatially coherent trajectories -all follow the specified direction field throughout the layers -while providing precise control over their inter-layer arrangement. This allows us to generate a staggered layout of trajectories across layers, preventing unavoidable tiny gaps from forming tunnel-shaped voids throughout a part volume.

Our approach is simple, robust, easy to implement, and scales linearly with the input volume. It builds upon recent results in procedural generation of oscillating patterns, generating a signal in the 3D domain that oscillates with a frequency matching the deposition beads width while following the input direction field. Trajectories are extracted with a process akin to a marching square.

INTRODUCTION

Additive manufacturing produces a part layer by layer, accumulating material progressively. As a consequence, the trajectories along which material is deposited or solidified have a direct influence on the final part properties. On the one hand, this provides unprecedented control over the way the volume is formed, enabling novel and unique possibilities to perfectly adjust a part to its final functionality. On the other hand, the trajectories play a crucial role in the part structural soundness, and care must be taken to minimize potential defects.

Often, the deposition strategies have to enforce competing objectives. For instance, the process may be expected to produce dense parts with minimal porosity. However, it is also desirable for the trajectories to follow specific directions within layers, either due to manufacturing constraints or to maximize part functionality, e.g., stiffness under some external loads. In particular, additively manufactured parts tend to be stronger in the direction of deposition, while in the orthogonal direction, the bond between neighboring beads is less strong [START_REF] Rybachuk | Anisotropic mechanical properties of fused deposition modeled parts fabricated by using acrylonitrile butadiene styrene polymer[END_REF]]. In addition, some processes are capable of depositing fibers alongside the trajectories, resulting in even higher ratios of anisotropic mechanical response. Being able to conform the deposition trajectories with specific directions, therefore, paves the way to several applications: stronger parts (e.g. [START_REF] Steuben | Implicit slicing for functionally tailored additive manufacturing[END_REF]), controlled deformations (e.g. [An et al. 2018]), as well as the ability to pattern the surface roughness (light reflection, grip).

In this work, we propose a novel approach to generate trajectories following an input direction field. The field may be user-created, or may be the result of an optimization process (e.g., topology optimization [START_REF] Pantz | A Post-Treatment of the Homogenization Method for Shape Optimization[END_REF]).

Our objective is to produce dense parts, conciliating the contradictory objectives of fully covering an area with evenly spaced material beads while having them follow an arbitrary direction. We seek a result where defects such as gaps have a minimal impact within a layer and across layers.

Even in the absence of a direction field to follow, avoiding defects in infills is a problem that remains a focus of research (we provide details in Section 2). Adding the direction objective makes this problem significantly harder, as following the direction requires paths to move away from one another or to regroup in specific locations, creating branching patterns where tiny gaps are difficult to avoid. Due to spatial coherence of trajectories geometry across layers, such small infill defects within a layer tend to accumulate across layers and degenerate into critical defects: internal tunnel-shaped voids that endanger the structural integrity of the part. This problem is worsened when following a given direction field that is similar from one layer to another.

Our technique generates trajectories for planar layers. However, to address the aforementioned challenges, we generate infill trajectories that not only follow the given directions but also form a coherent structure across layers. In particular, the deposited beads form a staggered structure akin to the layout of a brick wall. As with bricks and mortar [START_REF] Nicholson | The New Practical Builder and Workman's Companion, Containing a Full Display and Elucidation of the Most Recent and Skilful Methods Pursued by Architects and Artificers... Including, Also, New Treatises on Geometry[END_REF]], the intent is to prevent the weaker lateral bond between beads from forming vertical sheets prone to failure. This also prevents tiny gaps in split/merge locations from forming coherent structures across layers within the volume. Please note however that in this paper we focus solely on the modeling aspects of a staggered infill. Measuring to what extent this truly benefits mechanical properties will require an in-depth mechanical study left as future work.

We achieve a staggered layout by building upon the regularized phasor noise technique introduced by Tricard et al. for microstructure design [START_REF] Thibault Tricard | Freely orientable microstructures for designing deformable 3D prints[END_REF]]. This efficient, scalable approach generates oscillating fields within volumes, whose direction can be controlled. We define and extract our trajectories from crossed phasor fields (Section 4.1,4.2 and 4.3), and exploit the exposed argument of the oscillations to obtain a staggered layout (Section 4.5). By controlling the oscillations frequency, we can further adjust the deposition width (Section 4.2).

This results in a simple, robust, and efficient method for generating dense, oriented trajectories. Thanks to the stochastic nature of the underlying synthesis process, the unavoidable tiny gaps in branching locations are evenly distributed within the volume and are further prevented to align and form long tunnel-shaped void by the staggered layout of the beads. We analyze results and discuss applications in Section 5.

PREVIOUS WORK

The processing pipeline for additive manufacturing is responsible for generating the instructions used by the machine to fabricate the final part, layer by layer [START_REF] Livesu | From 3D models to 3D prints: an overview of the processing pipeline[END_REF]]. Our work focuses on the later steps of the pipeline, where the geometry of layers has already been determined, but surface-covering trajectories remain to be computed. In particular, our technique is orthogonal to and compatible with slicing approaches such as adaptive slicing [START_REF] Alexa | Optimal discrete slicing[END_REF][START_REF] Kulkarni | An accurate slicing procedure for layered manufacturing[END_REF]Pandey et al. 2003;[START_REF] Sabourin | Adaptive slicing using stepwise uniform refinement[END_REF].

Techniques computing trajectories for inner parts of a layer are often referred to as infill strategies. For dense infills, the two most common techniques are contour-parallel infills [START_REF] Held | Pocket machining based on contour-parallel tool paths generated by means of proximity maps[END_REF]Yang et al. 2002] and zigzag (raster) infills [START_REF] Mcmains | Layered Manufacturing of Thin-Walled Parts[END_REF]], and their typical combination. Note that we do not review here techniques for sparse infill unless technically related, as our focus is to produce fully filled parts.

The contour parallel and zigzag infills leave gaps within a layer, as it is generally not possible to fully and exactly cover an arbitrarily shaped polygon with fixed-width strokes. These defects are understood as underfill (gaps) and overfill (too much material deposited in a same location). Therefore, several methods proposed to reduce these artifacts, in particular by locally varying the deposition radius [START_REF] Ding | Adaptive path planning for wire-feed additive manufacturing using medial axis transformation[END_REF][START_REF] Hornus | Variable-width contouring for additive manufacturing[END_REF][START_REF] Yuan | Optimization of process planning for reducing material consumption in additive manufacturing[END_REF][START_REF] Kao | Optimal motion planning for deposition in layered manufacturing[END_REF][START_REF] Kuipers | A Framework for Adaptive Width Control of Dense Contour-Parallel Toolpaths in Fused Deposition Modeling[END_REF][START_REF] Xiong | Process planning for adaptive contour parallel toolpath in additive manufacturing with variable bead width[END_REF]]. However, in all these approaches, the direction of the infill patterns is driven by the contour of the shape (contour parallel) or is fixed in space (zigzag), possibly chosen differently in different zones [START_REF] Ding | A tool-path generation strategy for wire and arc additive manufacturing[END_REF].

A recent trend is, however, to consider the direction of the trajectories with respect to the desired functionalitymost often the part structural rigidity, but also to, e.g., trigger controlled deformations under heat [An et al. 2018]. This exploits the natural anisotropy of the mechanical resistance of 3D printed parts [START_REF] Fang | Reinforced FDM: Multi-Axis Filament Alignment with Controlled Anisotropic Strength[END_REF][START_REF] Rybachuk | Anisotropic mechanical properties of fused deposition modeled parts fabricated by using acrylonitrile butadiene styrene polymer[END_REF] which can be further reinforced by the inclusion of fibers. In this work, we do not presume a specific application and assume the input direction field given.

In order to obtain oriented infills, [START_REF] Steuben | Implicit slicing for functionally tailored additive manufacturing[END_REF] optimize a continuous scalar field in which gradient aligns with a principal stress direction. Trajectories are extracted as iso-values in this field. However, there is no specific treatment for maintaining a constant spacing between trajectories. [START_REF] Ezair | Volumetric covering print-paths for additive manufacturing of 3D models[END_REF] propose a curved slicing approach that has the potential to address this issue, as the trajectories are extracted from a given tri-variate field while maintaining their spacing by hierarchical subdivision [START_REF] Elber | Adaptive Isocurve-based Rendering for Freeform Surfaces[END_REF][START_REF] Etienne | Procedural band patterns[END_REF]. [START_REF] Fang | Reinforced FDM: Multi-Axis Filament Alignment with Controlled Anisotropic Strength[END_REF] extract curved trajectories as iso-values of a field that is specially optimized to enforce spacing and direction constraints. The main drawback of these approaches is that computing a tri-variate field under strict direction constraints is a computationally and numerically challenging task. In addition, as we will later show, adjusting spacing through iso-value extraction leads to gaps in coverage and to an alignment of the split/merge locations detrimental to the resulting quality. [START_REF] Boddeti | Optimal design and manufacture of variable stiffness laminated continuous fiber reinforced composites[END_REF] rely on the optimization of 2D stripe patterns in individual layers to produce trajectories. This is based on the global optimization of a periodic parameterization [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF]. This allows a more natural adaptation to changes in direction. Similar methodologies have been applied to generate conforming lattices [START_REF] Wu | Design and Optimization of Conforming Lattice Structures[END_REF]] and flexible microstructures [START_REF] Thibault Tricard | Freely orientable microstructures for designing deformable 3D prints[END_REF], using techniques respectively from hex-meshing [START_REF] Gao | Robust Hex-Dominant Mesh Generation using Field-Guided Polyhedral Agglomeration[END_REF]] and procedural texturing [START_REF] Thibault Tricard | Procedural Phasor Noise[END_REF]]. Our technique is rooted in similar methodologies: we generate periodic 3D signals to extract oriented trajectories across all 2D layers, arranged in a globally coherent staggered layout.

BACKGROUND ON PHASOR NOISE

Our technique relies on the phasor noise introduced by Tricard et al. [START_REF] Thibault Tricard | Procedural Phasor Noise[END_REF]. A phasor noise 𝜙 is a procedural noise allowing the definition of fiber-like patterns. It is defined as the argument of a complex Gabor noise 𝐺 [START_REF] Ares Lagae | Procedural noise using sparse Gabor convolution[END_REF]:

𝜙 (𝑥) = arg (𝐺 (𝑥)) (1)
with :

𝐺 (𝑥) = 𝑗 𝑒 -𝑏 ∥𝑥-𝑥 𝑗 ∥ 2 𝑒 2𝑖 𝑓 𝑗 𝑑 𝑗 .(𝑥-𝑥 𝑗 )+𝑖𝜑 𝑗 (2)
It is an unstructured representation, defined from the convolution of complex Gabor kernels with random seed positions 𝑥 𝑗 . The kernel properties locally define the properties of the noise: its bandwidth 𝑏, frequency 𝑓 𝑗 , the direction of anisotropy 𝑑 𝑗 , and phase shift 𝜑 𝑗 . This generates a scalar field in 3D space: an oscillating signal of precisely controlled direction and scale. In our work, we use the oscillations to define trajectories for material deposition.

Without special treatment, a phasor noise exhibits local defects located around singularities of the phase field. These result in a less coherent pattern, and in our case, would result in less continuous deposition trajectories. In [START_REF] Thibault Tricard | Freely orientable microstructures for designing deformable 3D prints[END_REF]] two methods are combined to significantly decrease the number of defects in the generated patterns. First, the phase 𝜑 𝑗 of neighboring kernels is iteratively adjusted in a simple and efficient optimization process. Second, a closed-form regularization filter is applied when sampling values from the phasor noise. Through these combined methods, highly regular anisotropic patterns are obtained. The high quality and regularity of the fields enable the definition of local parametric patches for microstructure design; for more details, please refer to [START_REF] Thibault Tricard | Freely orientable microstructures for designing deformable 3D prints[END_REF]]. In our work, we extend this methodology to define deposition trajectories for dense infill patterns.

The phasor noise presents several advantages for our application. First, it is a meshless process. In other words, it means that it can be computed in any Euclidean space while following arbitrary input fields. Secondly, it scales to large volumes as to the complexity of the iterative phase alignment, in both space and time, is only proportional to the number of kernels -not the resolution of the generated patterns. Finally, the entire framework maps well onto GPUs, enabling fast, parallel optimization and evaluation of regularized phasor fields. For the remainder of the paper, it is important to remember that the phasor noise is defined in 3D.

OUR APPROACH

Our technique takes as an input a 3D model for which to generate 3D trajectories as well as a direction field 𝑑 𝑈 that is everywhere orthogonal to 𝑧 (the manufacturing direction). The algorithm starts by defining two crossed phasor fields used to generate oscillating patterns everywhere throughout the object volume. The oscillation frequency 𝑓 is chosen to match the target trajectory spacing -typically the nozzle width 𝑤 (with 𝑓 = 1/𝑤). However, by varying the frequency, we naturally support varying deposition widths, see Section 4.2. The direction of the phasor fields closely follows 𝑑 𝑈 .

Once the fields are optimized, we proceed layer by layer to efficiently generate actual trajectories. This step extracts a set of discrete curved trajectories from the continuous definition of two phasor fields. We describe it in Section 4.1.

Note that this is performed in a streaming fashion (we only process a fixed number of slices at a time) and that slices anywhere in space may be extracted in parallel. To prevent trajectories from creating inter-layer artifacts by always being exactly aligned, we introduce a phase shift every other layer. This produces staggered trajectories across layers while still closely following the direction field.

Defining trajectories

In this section, we describe our method to define dense infills -and extract the associated trajectories -following the arbitrary direction field 𝑑 𝑈 given by the user, whose directions must be constrained to the slicing plane. Let us define 𝑑 𝑉 a second field orthogonal to 𝑑 𝑈 , also constrained to the slicing plane (i.e. 𝑑 𝑉 = 𝑑 𝑈 × 𝑧). Note that 𝑑 𝑈 and 𝑑 𝑉 are defined as 3D fields, and their values can evolve along the Z-axis. As a result, the trajectories' direction can evolve across layers. The trajectory extraction is based on two orthogonal periodic variable fields, 𝑈 and 𝑉 , obtained from two phasor fields. To any given point 𝑥 in space we associate a value pair 𝑈𝑉 (𝑥) ∈ [0, 2𝜋] 2 . We note 𝑈 (𝑥) and 𝑉 (𝑥) the two components of 𝑈𝑉 (𝑥) (see Figure 2, left). In practice, 𝑈 (𝑥) (respectively 𝑉 (𝑥)) defines a locally oscillating signal taking the shape of a sawtooth. This defines a periodic variable field following the direction 𝑑 𝑈 (respectively 𝑑 𝑉 ). From this periodic variable field we can construct a set of curves defined by 𝑈 (𝑥) = 𝜋 (respectively 𝑉 (𝑥) = 𝜋). The generated curves 𝐶 𝑈 𝜋 (respectively 𝐶 𝑉 𝜋 ) are orthogonal to 𝑑 𝑈 (respectively 𝑑 𝑉 ). Therefore, we define the trajectories as the set of curves 𝐶 𝑉 𝜋 (see Figure 2, left). The orthogonal set of curves 𝐶 𝑈 𝜋 is used during the computation of inter-trajectories distances and defines the sampling rate of the trajectories (see Section 4.3).

In practice, the extraction of the trajectories is performed in three steps. First, we define and optimize our periodic variable fields (see Section 4.2). Then, for a given printing slice, we extract a planar graph in whose vertices are defined as the intersection between 𝐶 𝑉 𝜋 and 𝐶 𝑈 𝜋 (see Figure 2, right). Edges of the graph are subdivided into two sets, 𝐸 𝑈 and 𝐸 𝑉 .

The edges of 𝐸 𝑈 (resp. 𝐸 𝑉 ) form a tessellation of the 𝐶 𝑉 𝜋 curves (respectively 𝐶 𝑈 𝜋 ). The graph extraction is described in Section 4.3. Finally, we use 𝐸 𝑈 to define the trajectories and use 𝐸 𝑉 to measure the local spacing between them. This information is later used to adjust the deposition flow (see Section 4.4).

Periodic variable fields

We define the two local periodic variable fields 𝑈 and 𝑉 by relying on two orthogonal 3D phasor noises (see Equation 1).

The kernel directions 𝑑 𝑖 are computed by sampling the input direction field 𝑑 𝑈 (resp. 𝑑 𝑉 ) at the kernel position 𝑥 𝑖 . The frequencies 𝑓 𝑖 are chosen as 1/𝑤, with 𝑤 the target deposition bead width. Note that 𝑤 can be made to spatially vary by sampling an additional user input field. Finally, parameter 𝜑 𝑖 is computed by iterative phase alignment [START_REF] Thibault Tricard | Freely orientable microstructures for designing deformable 3D prints[END_REF]].

Graph extraction

As previously stated, we seek to extract two sets of edges 𝐸 𝑈 and 𝐸 𝑉 from the periodic variable field 𝑈𝑉 . Recall that while the fields are defined in 3D, the trajectory extraction is performed in 2D, layer by layer (possibly in parallel). We base our approach on a modified version of the marching square algorithm [An et al. 2018;[START_REF] Wyvill | Data Structure for Soft Objects[END_REF]].

First, we sample our periodic variable field on a regular 2D grid which resolution is chosen to sample at four times the maximum field frequency. This sampling frequency is chosen to help to filter the 2𝜋 jump in the periodic variable fields. Then, for each square cell of the grid, we extract 𝐶 𝑈 𝜋 and 𝐶 𝑉 𝜋 separately. We extract vertices for the intersection between 𝐶 𝑈 𝜋 (respectively 𝐶 𝑉 𝜋 ) along each edge of the cell. More precisely, given 𝑖 and 𝑗 the grid points defining a processed edge, we extract a vertex if 𝑈 (𝑖) -𝜋 and 𝑈 ( 𝑗) -𝜋 have opposite signs. We filter out vertices that do not verify |𝑈 (𝑖) -𝑈 ( 𝑗)| ≤ 𝜋 as they are generated by the 2𝜋 jumps in the field. This condition ensures robustness with respect to the periodicity of the field. The positions of the vertices are chosen based on the linear interpolation of the 𝑈 (resp. 𝑉 ) parameter along the edges.

For a given field, we only extract an edge in a given cell if there are exactly two vertices created. Edges extracted from fields 𝑈 (resp. 𝑉 ) are registered in the set 𝐸 𝑉 (resp. 𝐸 𝑈 ). As this is done for both 𝑈 and 𝑉 , a cell contains from 0 to 2 edges. If the cell contains two intersecting edges (defined from two distinct fields), the latter are subdivided, and a new vertex is created at the intersection point, resulting in four new connected edges.

When each square cell of the grid is processed, we stitch all the partial results together and collapse edges until all vertices either belongs to the two edge sets (vertices of valence 4 verifying 𝑈𝑉 ≈ {𝜋, 𝜋 }) or are curve end-points (vertices of valence one close to singularity or on the border of the domain). This reduces the sampling frequency of 𝐶 𝑈 𝜋 (respectively 𝐶 𝑈 𝜋 ) to the frequency of the 𝑉 (respectively 𝑈 ). Thus the frequency of 𝑈 controls the polygonization error in curved trajectories. In practice, we use the same frequency for both fields as they generate a sampling equal to or close to the nozzle diameter, which is enough for our application. Note that we keep track of the edge origin during the decimation, thus updating the edge sets 𝐸 𝑈 (edge marked as V) and 𝐸 𝑉 (for edge marked as U).

Clipping and ordering trajectories

To extract the final trajectories from the graph, we clip the set of edge 𝐸 𝑈 by the contour of the slice; then, we use a greedy approach to extract the longest uninterrupted paths. For each point, we compute its average distance to all its direct neighbors on the 𝐸 𝑉 set. This average distance -which represents the local bead width -is then clamped in order to remain between half and twice the nozzle size. It is then used to proportionally adjust the deposition volume. Finally, we transform the trajectories into GCode and adjust the material deposition at each point by taking into account the inter-spacing previously computed.

We order the trajectories following a simple nearest-first heuristic (the first trajectory is picked at random, then the next trajectory is the one with closest end-point). This provides reasonably good results at a low cost. An ordering can be visualized in Figure 5, right.

Staggered bead layout

To obtain a staggered layout, we offset the periodic variable fields for every odd layer. This is easily achieved by shifting the phases of the kernels 𝜑 𝑖 when sampling odd layers, adding an offset 𝛿 = 𝜋. Doing so produces a pattern where the beads of one layer are shifted by half the width of the beads from the layer below, forming a coherent staggered layout.

Note that this works naturally across directions and frequency (width) controls, as we directly manipulate the phase of the coherent 3D oscillations. This simple approach produces remarkably regular and aligned patterns, as can be seen in Figure 3.

RESULTS

In this section, we discuss the specific properties of our method and demonstrate the range of possibilities it offers.

All our tests are done with fused-filament 3D printing, the most widespread additive manufacturing technique. Parts were printed with a variety of printers. Blue CuteOcto was printed with a Prusa MK2S, parts in grey with a micro delta, parts in multicolor with a modified Prusa printer equipped with a diamond mixing extruder, other parts were printed using a Creality CR10S pro. Note that all parts shown here are solely composed of infill to illustrate the method.

In practice our infills would likely be combined with a contouring technique for the outer object shell, replacing the internal typical zigzag solid fill. The list of all model sources is provided in Table 1. Our method was implemented in C++ using the Vulkan API, and shaders were programmed with glsl. Trajectories are always extracted on the full bounding box of the object before being clipped to the surface. Therefore we measure runtimes for cubes of increasing volumes and only measure timings related to our method (e.g., alignment of kernel phases and graph extractions, leaving out the actual clipping of trajectories, scheduling as well as GCode generation). Runtimes were measured on an i7-6700HQ CPU (@2.6GHz only one thread used) and GTX1070 graphic card and are given in Table 2. The first thing to notice is that phase alignment runtime is bounded to 28𝑠. This is due to our strategy to compute kernel density and bandwidth such that the phase alignment of kernels can be done in a single pass on GPU while avoiding memory overflow (the maximal number of kernels used is approximately eight millions). If needed, a streaming strategy could be adopted to increase the number of degrees of freedom in larger parts. As expected, the graph extraction (including the sampling of the phasor noises) is linear in terms of processed volume. The grid size for the marching square is based on the target For all tested models, the GPU memory used remains bounded to 1.7𝐺𝐵. Should the memory consumption become problematic, a smaller set of slices and XY tilling could be used.

Comparison to curve hierarchies

When a bi-variate parametric field is available, curve hierarchies may be used to produce trajectories that follow the field while being roughly equally spaced [START_REF] Elber | Adaptive Isocurve-based Rendering for Freeform Surfaces[END_REF][START_REF] Ezair | Volumetric covering print-paths for additive manufacturing of 3D models[END_REF]]. We compare our approach to a curve hierarchy in Figure 4. As can be seen, the curve extraction leaves many gaps. This typically happens where the field exhibits a high divergence, as is the case with the radial field shown in Figure 4 (left). There, steep transitions of local density from 1 to 0.5 occur: the change of hierarchy level introduces a sudden jump. In addition, the gaps occur along specific transition curves of the field, which in turn leads to coherent defects within and across layers, weakening the overall part structure. Note that other curves extraction techniques could alleviate this problem [START_REF] Etienne | Procedural band patterns[END_REF], but this has not been demonstrated for dense infills.

In comparison, our method exhibits fewer transitions with tiny gaps well distributed throughout a slice (Figure 4 (right)) -keeping in mind that the paths are additionally laid out staggered from one layer to the next. Finally, recall that this is achieved directly from the direction field without the need to optimize for a bi-variate parametric field in the first place.

Singularities, staggered layout

Figure 5 shows the CuteOcto model sliced with a radial direction field and constant bead width (nozzle diameter). As can be seen, the branching pattern and split/merge locations are evenly distributed within each slice.

Figure 7 reveals the effect of the staggered layout across layers. It is the same model as Figure 5 with a sectional view of the beads. On the left disabling the phase shifts, on the right our result with the staggered layout. This cut reveals that singularities have a limited extent in the Z direction. It also shows how our staggered layout (right) reduces the influence of singularities by breaking the gap coherence, increasing the bond between beads in these areas. It is important to notice that the staggered layout visibility is reduced on printed parts (see Figure 1(left), 8) due to travel artifacts between trajectories at the boundary of the object (visually similar to larger beads).

Another slicing of CuteOcto presents the combination of direction field and deposition bead width field (see Figure 8 6). Direction field is defined as the gradient of the signed distance field to the surface, while bead width field is radially defined.

Control over the direction field can help for a wide range of applications, from controlling the appearance to the physical properties of the printed object. In the screwdriver handle (see Figure 9), the direction field is used to control the outer texture of the object. A grip-like pattern is created in the locally convex area (by using trajectories locally orthogonal to the surface), while in the concave area, smooth surfaces are used (trajectories parallel to the surface). The direction field can also be used to control the appearance of horizontal parts of an object where filament directions are clearly visible (see Figure 1(right),10). The direction field can also be used to align paths along branches in branch-like structures (see Figure 11). In this case the direction field was created to favor long trajectories across branching regions (see junctions between legs and torso). Finally, the direction field can also be used to adapt fiber direction in mechanical parts (see Figure 1(left),12).

DISCUSSION

We presented a method to stagger the beads forming a slice while orienting them. The staggered layout is achieved by adding a phase shift of 𝜋 on the result of the phasor noise every odd slices. We chose this specific phase shift to offset the beads by half their size at each layer. In practice it is possible to achieve different staggered layouts by varying the phase shift and the sets of layers onto which it is applied.

In our implementation we input a 3D voxel field to define the orientation of the beads. This field is hand painted by the user in a dedicated tool, and then sampled to orient each phasor kernel. In practice a field of 64 × 64 × 64 voxels was used in all our examples. Other definitions of orientation fields could be used -e.g. procedural -provided they can be sampled at the center of each phasor kernel.

When printing with FDM, using dense and strongly oriented patterns can create deformations (shrinkage / warping) while printing. This effect stems from the tendency of the polymer to retract when cooling [Yaman 2018]. This did not create specific difficulties when printing our example parts with PLA material.

Before being printed the paths are ordered using a greedy approach, where the next bead to be printed is the one with an extremity closest to the last printed position. This heuristic encourages nearby paths to be printed successively while being fast to compute. Naturally, other ordering methods could be devised depending on the material and process requirements.

CONCLUSION

We introduced a new approach to extract dense planar trajectories following a direction field. The trajectories are obtained from two periodic variable fields defined in 3D from phasor noises. Interestingly, phasor fields have links with global parametrization techniques for quad-and hex-meshing, as discussed in [START_REF] Thibault Tricard | Freely orientable microstructures for designing deformable 3D prints[END_REF]. Therefore, such approaches could likely be modified for a purpose similar to ours. Our approach however inherits the low setup cost, scalability and simplicity of phasor noises, enabling the efficient extraction of trajectories.

The 3D nature of the field allows to distribute imperfections evenly throughout the volume. It also introduces a way to control spatial coherency across layers, leading to the staggered bead layout. While in-depth mechanical studies are required, we believe this inter-layer layout control to be a promising direction of research, as this degree of freedom has not been available before. Our method presents a few limitations. In comparison to traditional infill patterns, it produces heavier GCode due to the complex curved trajectories. This could be addressed by relying on circular arc GCode capabilities. When generating paths that follow the local curvature of the object, our approach does not guarantee a perfect alignment between the object boundary and the deposition paths. This results in some trajectories being cut when encountering the object boundary. As future work, this could be approached by controlling the phase of the Gabor kernels where a precise alignment is required.

Finally, a promising research direction is to generalize our technique to curved layer printing. The key difficulty remains the collision-free guarantees required, both for the direction field definition and the ordering of paths in 3D.

Fig. 2 .

 2 Fig. 2. Representation of the graph extraction process. First: both phasor fields are sampled; here, their arguments are remapped in [0,1] for visualization. The set of curves 𝐶 𝑈𝜋 is highlighted in green while 𝐶 𝑉 𝜋 is in red. Then intersections of 𝐶 𝑈 𝜋 and 𝐶 𝑉 𝜋 with edges of the grid are computed separately, creating multiple polylines. The intersection between polylines of both sets is computed to define the vertices of the graph (in black here). Finally, we collapse edges until all vertices either belong to the two edge sets (valence equal four) or are curve end-points (valence equal one).

Fig. 3 .

 3 Fig. 3. Trajectories can be offset every other layers by shifting the phases of the kernels. This allows to create staggered beam layouts. Handmade 3D example and sectional view of deposition paths. Left: aligned. Right: staggered.

Fig. 4 .

 4 Fig. 4. Comparison to [Ezair et al. 2018] on a radial direction field. Left : method of [Ezair et al. 2018], steep transition of local density can be observed. Right: our methods, note that unavoidable tiny gaps are well distributed throughout the slice.

Fig. 5 .

 5 Fig. 5. CuteOcto model sliced with a radial direction field, constant target beam width, and staggered beam layout. Left: one layer with colored beads. Right: beads colored by order of printing.

Fig. 6 .

 6 Fig. 6. First slices of the CuteOcto model with constant (left) and varying (right) target bead width.

Fig. 7 .

 7 Fig.7. A sectional view of the beads within the CuteOcto model volume, under a radial direction field. Left: phase shifts disabled, note the limited vertical gaps near singularities. Right: Our staggered layout breaks the gap coherence, increasing the bond between beads in these areas and throughout the volume. Note that singularities are placed differently as the optimization of phasor kernels is non-deterministic.

Fig. 8 .

 8 Fig.8. Bead width can be controlled with an additional user input field. CuteOcto model sliced with a radially defined beam width. Direction field is defined as the gradient of the signed distance field to the surface.

Fig. 9 .

 9 Fig.9. Direction field near the surface of an object can be used to control surface roughness. Here a screwdriver handle have a grip-like surface in convex area (direction orthogonal to the surface) and a smooth surface in concave area (direction in surface tangent plane).

Fig. 10 .

 10 Fig.10. Direction field can be used to texture horizontal area by using the fact that filament deposition direction is clearly visible on 3D printed models.

Fig. 11 .

 11 Fig. 11. Gecko object exhibiting a branching structure. Left: Paths are generated in the direction of branches while allowing long deposition paths in the neighborhood of the junctions. Right: The object was printed with a printer equipped with a diamond mixing extruder. The color mixing ratios were chosen to highlight trajectories direction and continuity.

Fig. 12 .

 12 Fig. 12. A mechanical part. A view from below is provided in the second row to highlight the input direction field.

Table 1 .

 1 Table of all models used in the article, Figures they can be seen in, and their sources

	Model	Figure	Sources
	Gear	1	custom made
	SCF	1	custom made
	Badge	10	thingiverse.com/thing:1670621
	Cute octo says hello 5, 8, 7	thingiverse.com/thing:27053
	Screwdriver handle	9	thingiverse.com/thing:34852
	Gecko	11	thingiverse.com/thing:1363148
	Fandisc	12	github.com/dengwirda/jigsaw-models

Table 2 .

 2 Computation time in function of the volume for a 0.4mm nozzle and a layer height of 0.2mm. the length of an edge is one-quarter of the minimal target bead width), and 64 slices are computed at once.

	Volume (𝑐𝑚 3 ) Total (𝑠) Phase optimisation (𝑠) Graph extraction (𝑠) Time over volume (𝑠/𝑐𝑚 3 )
	1	2.81	1.64	0.79	2.81
	8	17.11	9.27	6.54	2.14
	64	82.05	28.14	48.97	1.28
	512	355.32	28.20	323.12	0.5
	bead width (

ACKNOWLEDGMENTS

This work was supported partly by the french PIA project « Lorraine Université d'Excellence » (ANR-15-IDEX-04-LUE), and the project IMPRIMA(ANR-18-CE46-0004). We thank Pierre-Alexandre Hugron for his help in designing and printing the examples shown in this work .


Xiangjing Lai
email: laixiangjing@gmail.com

Jin-Kao Hao
email: jin-kao.hao@univ-angers.fr

Dong Yue

Two-stage solution-based tabu search for the multidemand multidimensional knapsack problem

Keywords: Metaheuristics, Multidemand multidimensional knapsack problem, Twostage optimization, Solution-based tabu search, Combinatorial optimization

The multidemand multidimensional knapsack problem (MDMKP) is a significant generalization of the popular multidimensional knapsack problem with relevant applications. In this work we investigate for the first time how solution-based tabu search can be used to solve this computationally challenging problem. For this purpose, we propose a two-stage search algorithm, where the first stage aims to locate a promising hyperplane within the whole search space and the second stage tries to find improved solutions by exploring the reduced subspace defined by the hyperplane. Computational experiments on 156 benchmark instances commonly used in the literature show that the proposed algorithm competes favorably with the stateof-the-art results. We analyze several key components of the algorithm to highlight their impacts on the performance of the algorithm.

Introduction

Given a set V = {1, 2, . . . , n} of n items, a set R = {r 1 , r 2 , . . . , r m } of m resources with a capacity upper limit b i for resource r i (1 ≤ i ≤ m), where each item j of V is associated with a profit c j and consumes a given quantity a ij for each resource r i (i ∈ {1, 2, . . . , m}), the popular NP-hard 0-1 multidimensional knapsack problem (MKP) involves selecting a subset of items from V such that the resource consummation of the selected items does not exceed the given capacity upper limit for each resource in R (knapsack constraints), while maximizing the total profit of the selected items. Formally, the MKP can be written as follows:

(M KP ) Maximize z = n j=1 c j x j (1)
s.t. n j=1 a ij x j ≤ b i , ∀i ∈ {1, 2, . . . , m} (2)
x j ∈ {0, 1}, ∀j ∈ {1, 2, . . . , n}

where c j ≥ 0, a ij > 0, b i > 0, ∀i ∈ {1, 2, . . . , m}, ∀j ∈ {1, 2, . . . , n} and Eq.

(3) indicates that the binary decision variable x j (1 ≤ j ≤ n) takes the value of 1 if the item j is selected, 0 otherwise.

The multidemand multidimensional knapsack problem (MDMKP) studied in this work is an important extension of the MKP, where q greater-than-orequal-to constraints are imposed, in addition to m less-than-or-equal-to constraints (Eq. ( 2)). Moreover, unlike the MKP, the profit c j in the MDMKP can take a positive, negative or zero value for each item j (j ∈ V ). Formally, the MDMKP can be formulated as follows [START_REF] Arntzen | Adaptive memory search for multidemand multidimensional knapsack problems[END_REF][START_REF] Cappanera | A local search based heuristic for the demand constrained multidimensional knapsack problem[END_REF]:

(M DM KP ) Maximize z = n j=1 c j x j (4)
s.t. n j=1 a ij x j ≤ b i , ∀i ∈ {1, 2, . . . , m} (5) n j=1
a ij x j ≥ b i , ∀i ∈ {m + 1, m + 2, . . . , m + q} [START_REF] Carlton | A note on hashing functions and tabu search algorithms[END_REF] x j ∈ {0, 1}, ∀j ∈ {1, 2, . . . , n}

where the following conditions are assumed:

b i > 0, a ij ≥ 0 ∀i ∈ {1, 2, . . . , m + q}, ∀j ∈ {1, 2, . . . , n} (8)
n j=1 a ij > b i ∀i ∈ {1, 2, . . . , m + q} (9)

max j {a ij } ≤ b i ∀i ∈ {1, 2, . . . , m} (10
) min j {a ij } < b i ∀i ∈ {m + 1, 2, . . . , m + q} [START_REF] Fréville | The multidimensional 0-1 knapsack problem: An overview[END_REF] In above formulas, the inequalities in Eq. ( 5) are called the knapsack constraints, and those in Eq. ( 6) are called the demand constraints.

Clearly, the classic MKP is a special case of the MDMKP when q equals 0 and the profit c j of item j takes a nonnegative value (i.e., c j ≥ 0, ∀j ∈ V ).

Like the MKP, the MDMKP has a number of practical applications [START_REF] Cappanera | A local search based heuristic for the demand constrained multidimensional knapsack problem[END_REF] like obnoxious and semiobnoxious facility location [START_REF] Cappanera | Discrete facility location and routing of obnoxious activities[END_REF][START_REF] Plastria | Static competitive facility location: An overview of optimisation approaches[END_REF][START_REF] Romero-Moraleses | Semi-obnoxious location models: A global optimization approach[END_REF], capital-budgeting, and portfolio-selection [START_REF] Beaujon | Balancing and optimizing a portfolio of R&D projects[END_REF], among others. On the other hand, the MDMKP is computationally challenging, given that it generalizes the NP-hard multidimensional knapsack problem. Consequently, there is no polynomial-time algorithm for the MDMKP, unless P = NP.

Unlike the MKP that has been subject of intensive studies in the past decades (see e.g., [START_REF] Chu | A genetic algorithm for the multidimensional knapsack problem[END_REF][START_REF] Fréville | The multidimensional 0-1 knapsack problem: An overview[END_REF][START_REF] Glover | Critical event tabu search for multidimensional knapsack problems[END_REF][START_REF] Hanafi | An efficient tabu search approach for the 0-1 multidimensional knapsack problem[END_REF][START_REF] Lai | A two-phase hybrid evolutionary algorithm for the 0-1 multidimensional knapsack problem[END_REF][START_REF] Mansini | CORAL: an exact algorithm for the multidimensional knapsack problem[END_REF][START_REF] Puchinger | The multidimensional knapsack problem: structure and algorithms[END_REF][START_REF] Shih | A branch & bound method for the multiconstraint zero-one knapsack problem[END_REF][START_REF] Vasquez | A hybrid approach for the 0-1 multidimensional knapsack problem[END_REF][START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF][START_REF] Vimont | Reduced costs propagation in an efficient implicit enumeration for the 0-1 multidimensional knapsack problem[END_REF]), the MDMKP receives much less attention until now. Still, there exist several exact and heuristic approaches in the literature. For example, general mixed integer programming solvers like CPLEX can be used to solve instances with n ≤ 100 to optimality within an acceptable time. However, it is usually difficult for the existing exact approaches to find an optimal solution for larger instances. As a result, several heuristic algorithms have been proposed to solve large instances approximately.

Specifically, in 2005, Cappanera and Trubian presented a nested-tabu-search heuristic [START_REF] Cappanera | A local search based heuristic for the demand constrained multidimensional knapsack problem[END_REF], which combines a standard attribute-based tabu search with an oscillation method presented in [START_REF] Glover | Critical event tabu search for multidimensional knapsack problems[END_REF]. In 2006, Arntzen et al. proposed an adaptive memory search method called Almha in [START_REF] Arntzen | Adaptive memory search for multidemand multidimensional knapsack problems[END_REF], which uses a dynamical tabu search mechanism and a weighting scheme to handle infeasible solutions. Their computational results show the Almha algorithm outperforms the previous best MDMKP methods and can be viewed as one of the best performing MDMKP algorithms in the literature. In 2007, Hvattum and Løkketangen investigated the behavior of scatter search on the MDMKP [START_REF] Hvattum | Experiments using scatter search for the multidemand multidimensional knapsack problem[END_REF]. In 2009, Gortázar et al. introduced a black box scatter search method for general classes of binary optimization problems, and assessed their method on the MDMKP and some other binary problems [START_REF] Gortázar | Black box scatter search for general classes of binary optimization problems[END_REF]. In particular, their method uses a static penalty approach proposed in [START_REF] Yeniay | Penalty function methods for constrained optimization with genetic algorithms[END_REF] to handle the constraints of the MDMKP. In 2010, Hvattum et al. proposed an alternating control tree (ACT) search framework for the MDMKP [START_REF] Hvattum | Alternating control tree search for knapsack/covering problems[END_REF], which can lead to an exact algorithm or heuristic algorithm by choosing the routine of solving subproblems. Their computational results show that the associated ACT algorithms have a high performance compared to a previous tabu search algorithm and scatter search algorithm. At the same year, Balachandar proposed a dominance principle based heuristic for the MDMKP [START_REF] Balachandar | On efficient techniques for NP-hard problems[END_REF].

In addition to these studies, there exist some theoretical investigations dedicated to the MDMKP in the literature. For example, Delissa investigated the existence and usefulness of equality cuts for the MDMKP [START_REF] Delissa | The existence and usefulness of equality cuts in the multidemand multidensional knapsack problem[END_REF], while Wishon and Villalobos studied robust efficiency measures for the MDMKP [START_REF] Wishon | Robust efficiency measures for linear knapsack problem variants[END_REF].

To enrich the solution arsenal for the MDMKP, we present in this work the first study of using solution-based tabu search [START_REF] Carlton | A note on hashing functions and tabu search algorithms[END_REF][START_REF] Carlton | Solving the traveling salesman problem with time windows using tabu search[END_REF][START_REF] Woodruff | Hashing vectors for tabu search[END_REF] to effectively solve the MDMKP. Actually, unlike the popular attribute-based tabu search approach [START_REF] Glover | Tabu search[END_REF], solution-based tabu search began to attract attention only very recently. Interestingly, this approach already showed excellent performances on several binary optimization problems as reported in [START_REF] Lai | Solution-based tabu search for the maximum min-sum dispersion problem[END_REF][START_REF] Lai | A two-phase hybrid evolutionary algorithm for the 0-1 multidimensional knapsack problem[END_REF][START_REF] Wang | Effective metaheuristic algorithms for the minimum differential dispersion problem[END_REF]. This work aims thus to investigate the interest of the solution-based tabu search approach for the MDMKP. Compared to attribute-based tabu search, solution-based tabu search has at least two appealing features. First, this approach ensures a stronger intensification ability, which is crucial for locating good local optima. Second, this approach makes the notion of tabu tenure irrelevant, thus simplifying the design of the algorithm and reducing the number of required parameters.

We summarize the contributions of this work as follows. First, based on the solution-based tabu search approach, we introduce an effective two-stage search algorithm for the MDMKP. The first search stage aims to identify a promising hyperplane within the whole search space while the second search stage tries to find improved solutions by examining both feasible and infeasible solutions on the identified hyperplane. For both stages, the solution-based tabu search strategy is employed, which relies on a one-flip and swap neighborhoods and a hash-based mechanism to efficiently determine the tabu status of neighbor solutions. Second, we assess the performance of the proposed algorithm based on 96 benchmark instances commonly used in the literature (n = 100, 250, m = 5, 10 and q = 2, 5, 10). The computational results show that the algorithm improves and matches the best known solutions for 17 and 71 instances respectively. Moreover, we report detailed computational results of the proposed algorithm for 60 additional instances with a large number of constraints (with n = 100, 500, m = q = 30). Third, given that the ideas of the two-stage search framework and solution-based tabu search developed in this work are quite general, they could be applied to solve other related binary optimization problems.

The remaining parts of the paper are structured as follows. In the next section, the proposed two-stage tabu search algorithm is described. In Section 3, we present the computational assessment of the proposed algorithm and report experimental results on the well-known benchmark instances. In Section 4, several essential ingredients of the algorithm are investigated to show how they affect the performance of the algorithm. In the last section, we summarize the present work and provide research perspectives. Specifically, the TSTS algorithm first generates randomly an initial solution by its initialization procedure (Section 2.2). Then the algorithm enters the first search stage where the initial solution is improved by the solution-based tabu search procedure presented in Sections 2.3 (line 4). During this search stage, the algorithm explores both feasible and infeasible solutions within the whole search space to find a high-quality feasible solution. At the end of its search stage, the best (feasible) solution found and the consumed time t are returned. At this point, the second search stage is triggered, which starts from the solution returned by the first stage and uses another solution-based tabu search procedure (Section 2.4) to seek improved solutions (line 5). During this stage, the search is limited to the hyperplane Ω [k] (k being the number of the selected items in the returned solution of the first stage, see Section 2.4.1). Finally, the whole algorithm terminates when a given time limit t max is met, and the best solution found during the search process is returned as the final result of the algorithm (line 6).

Initial Solution

Algorithm 2: Procedure of Generating Initial Solution Function InitialSolution() Input: An instance I, the size of instance (n) Output: An initial solution s = (x 1 , x 2 , . . . , x n ) begin for i ← 1 to n do

x j ← rand() mod 2 /* rand() denotes a random integer */ end s ← (x 1 , x 2 , . . . , x n ) return s end
The initial solution of the TSTS algorithm is generated by a randomized procedure whose pseudo-code is given in Algorithm 2. Specifically, given an instance with n items, the initialization procedure assigns randomly a value from the set {0, 1} to each component x i (i = 1, 2, . . . , n) to obtain an initial solution s = (x 1 , x 2 , . . . , x n ). This random initialization has the advantage of being simple and fast. However, an initial solution generated in this way can be infeasible. If this is the case, its feasibility will be established during the first search stage described below.

Tabu Search Method of the First Search Stage

The first search stage is ensured by a solution-based tabu search algorithm (denoted by T abuSearch 1 ()), whose pseudo-code is given in Algorithm 3. After initiating its tabu lists (lines 3-5), T abuSearch 1 () performs a number of iterations to improve the current solution until 1) a feasible solution is found and no improvement can be observed during α consecutive iterations, where α is a parameter called the depth of tabu search, or 2) the allowed maximum time limit t max is reached (lines 8-17 The ingredients of this tabu search algorithm, including the search space, the neighborhood structures and the tabu strategy, are respectively described in the following subsections.

Search Space and Neighborhood

The search space Ω explored by the T abuSearch 1 () procedure is composed of all feasible and infeasible solutions of the given problem instance, i.e.,

Ω = {(x 1 , x 2 , . . . , x n )|x i ∈ {0, 1}, 1 ≤ i ≤ n} (12)
The neighborhood used by T abuSearch 1 () is a combined neighborhood composed of two basic neighborhoods, namely the one-flip neighborhood N 1 and the swap neighborhood N 2 . The one-flip neighborhood N 1 is defined by the one-flip operator (denoted by F lip(•)). Given a solution s = (x 1 , x 2 , . . . , x n ), a one-flip move F lip(q) consists of changing the value of a variable x q to its complementary value 1 -x q . As such, the neighborhood N 1 (s) of solution s includes all possible solutions that can be obtained by applying the one-flip operator to s. Formally, the N 1 (s) can be written as follows:

N 1 (s) = {s ⊕ F lip(q) : 1 ≤ q ≤ n} (13)
The neighborhood N 2 is defined by the swap operator (denoted by Swap(•, •)). Given a solution s = (x 1 , x 2 , . . . , x n ), let I 1 = {q : x q = 1 in s} and I 0 = {q :

x q = 0 in s}, the swap neighborhood N 2 (s) can be written as follows:

N 2 (s) = {s ⊕ Swap(i, j) : i ∈ I 1 , j ∈ I 0 ; } (14)
The first tabu search algorithm explores the union of these two neighborhoods, i.e., N (s) = N 1 (s)∪N 2 (s), whose size equals to n+|I 1 |×|I 0 |. At each iteration of the algorithm, a best non-tabu solution from N (s) according to the extended evaluation function defined by Eq. ( 17) in Section 2.3.3 and the tabu strategy explained in Section 2.3.2 is selected to replace the current solution s.

Tabu Strategy

In the present tabu search method, we adopt the solution-based tabu strategy to determine the tabu status of neighbor solutions. Specifically, the tabu lists are based on three hash vectors H 1 , H 2 , and H 3 of length of L, where each position of them represents a binary variable, and each hash vector H t is associated with a hash function h t . In particular, the effect of hash functions is to map a candidate solution of the search space Ω to an index of H t , i.e., h t : Ω → {0, 1, 2, . . . , L -1}.

Based on these hash vectors and the associated hash functions, we determine the tabu status of candidate solutions by the following rule. Given a candidate solution s, the hash vectors H t (t = 1, 2, 3) and the associated hash functions h t , s is identified as a tabu solution if

H 1 (h 1 (s)) ∧ H 2 (h 2 (s)) ∧ H 3 (h 3 (s)) = 1.
Otherwise, s is determined as a non-tabu solution.

Following previous studies [START_REF] Carlton | A note on hashing functions and tabu search algorithms[END_REF][START_REF] Wang | Effective metaheuristic algorithms for the minimum differential dispersion problem[END_REF][START_REF] Woodruff | Hashing vectors for tabu search[END_REF], we define our hash functions as follows.

Let s = (x 1 , x 2 , . . . , x n ) be a candidate solution, our hash functions h t (t = 1, 2, 3) are given by:

h t (s) = ( n i=1 ⌊i γt ⌋ × x i ) mod L (15)
where γ t is a parameter that is used to define each hash function and L is the length of hash vectors that is empirically set to 10 8 in this work.

Given a solution s and its hash value h(s), the hash value of its neighbor solutions can be determined in O(1) according to Eq. ( 15). Thus, the time complexity of determining the tabu status of a neighbor solution is O(1).

Extended Evaluation Function

Since the search space explored by the first tabu search algorithm contains both feasible and infeasible solutions, we devise an extended evaluation function F which uses a penalty function P to assess constraint violations.

Let s be a candidate solution in Ω, the penalty value P (s) is defined as the summation of all constraint violations, i.e.,

P (s) = m i=1 M ax{0, n j=1 a ij x j -b i } + q+m i=q+1 M ax{0, b i - n j=1 a ij x j } (16)
Thus, a small (large) function value P (s) means a weak (strong) constraint violation in s. In particular, P (s) = 0 means that s is a feasible solution.

Given this penalty function, the extended evaluation function F (s) is defined as a linear combination of the objective function f (s) in Eq.( 4) and P (s):

F (s) = n j=1 c j x j -λ × P (s) ( 17
)
where λ is a weighting factor that is empirically set to 10 2 in this work.

For any two solutions s ′ and s ′′ in

Ω [k] , s ′ is considered to be better than s ′′ if F (s ′ ) > F (s ′′ ).
As shown in our experimental results (Section 3, Section 4.2 and the Appendix), the first search stage equipped with the extended evaluation function always ends up with a feasible solution for all tested instances, including those with 30 demand constraints and 30 knapsack constraints. In other words, the first solution-based tabu search algorithm is typically able to locate a promising valid hyperplane that is further explored by the second search stage.

Tabu Search Method of the Second Search Stage

The second optimization stage of the TSTS algorithm uses another tabu search algorithm (denoted by T abuSearch 2 (), Algorithm 4) to examine candidate solutions of a given hyperplane Ω [k] (see below). Unlike the first tabu search algorithm, T abuSearch 2 () explores only solutions that contain exactly k selected items. T abuSearch 2 () first initializes the hash vectors (lines 3-5), and then performs a number of iterations to improve the current solution (lines 7-14). At each iteration, the algorithm replaces the current solution by a best non-tabu neighbor solution s ′ in terms of the evaluation function in Eq. ( 17).

During the search, the best feasible solution encountered s * is updated each time a better feasible solution is found, and the hash vectors are accordingly updated by the new solution (line 13). Finally, the algorithm terminates if the time limit t max is reached, and then returns the best feasible solution s * found during the search process.

Search Space, Tabu Strategy and Evaluation Function

The search space Ω [k] explored by T abuSearch 2 () is composed of both feasible and infeasible solutions with a fixed number of k selected items. In other words,

Ω [k] contains all n-dimensional 0-1 vectors with n i=1 x i = k, i.e., Ω [k] = {x ∈ {0, 1} n | n i=1 x i = k}. Ω [k] is also called a hyperplane of the search space Ω defined in Section 2.3.1, i.e., Ω = ∪ n k=1 Ω [k] .
Additionally, like the first tabu search algorithm, T abuSearch 2 () uses the solution-based tabu strategy described in Section 2.3.2 to determine the tabu status of neighbor solutions, and employs the extended evaluation function in Eq. [START_REF] Hvattum | Alternating control tree search for knapsack/covering problems[END_REF]

for i ← 0 to L -1 do H 1 [i] ← 0; H 2 [i] ← 0; H 3 [i] ← 0 end s * ← s /*
H 1 [h 1 (s)] ← 1; H 2 [h 2 (s)] ← 1; H 3 [h 3 (s)] ← 1 end return s * end 2.

Neighborhood Structure

To search effectively the hyperplane Ω [k] , T abuSearch 2 () uses a constrained swap neighborhood N 3 (s). Formally, given a solution s = (x 1 , x 2 , . . . , x n ), the neighborhood N 3 (s) is given by:

N 3 (s) = {s ⊕ Swap(i, j) : i ∈ I 1 , j ∈ I 0 ; f (s ⊕ Swap(i, j)) > f (s * )} (18)
where f (•) is the objective function of the MDMKP in Eq. ( 4), s * is the best feasible solution found so far in the current tabu search run, I 1 and I 0 denote respectively the sets of indices having the value of 1 (selected items) and 0 (non selected items) in s. Clearly, the size of this neighborhood is bounded by

O(|I 1 | × |I 0 |).
It is worth noting that we constraint the neighborhood by f (s ⊕ Swap(i, j)) > f (s * ) to eliminate non-promising neighbor solutions. Similar idea was previously investigated for the related MKP in [START_REF] Vasquez | A hybrid approach for the 0-1 multidimensional knapsack problem[END_REF].

Space and Time Complexities of the Algorithm

At each stage of our TSTS algorithm, in addition to three hash vectors (i.e., H 1 , H 2 and H 3 ) with a length of L, we maintain three solutions (i.e., s, s ′ , and s * ) to follow the search process, where each solution is stored by two vectors (i.e., I1 and I 0 ) with a maximum length of n and a vector W = (w 1 , w 2 , . . . , w m+q ) where w i = n j=1 a ij x j (j ∈ V ) holds. Thus, the space complexity of our TSTS algorithm is bounded by O(n + m + q + L).

In addition, for each neighbor solution in the search space, the time complexity of evaluating its quality is bounded by O(m + q), since there are m + q constraints needed to be checked. Hence, for each iteration, the time complexities of the first and second tabu search stages are respectively bounded by

O((m + q) × (n + |I 1 | × |I 0 |)) and O((m + q) × (|I 1 | × |I 0 |))
according to the size of neighborhoods explored by the algorithm (see Sections 2.3.1 and 2.4.2).

Experimental Results and Comparisons

We evaluate the proposed TSTS algorithm by conducting extensive computational experiments based on four sets of benchmark instances commonly used in the literature and by making a comparison between our results and the state-of-the-art results in the literature.

Benchmark Instances

In this study, we employed four sets of benchmark instances to assess the performance of our TSTS algorithm, where the first two sets of benchmark instances are available at http://www.optsicom.es/binaryss, and the third and fourth sets of benchmark instances are available in OR-Library 1 . The first set contains 48 small instances with n = 100, and the second set contains 48 larger instances with n = 250. In addition, for the instances in the first two sets, the number of knapsack constrains m varies from 5 to 10, and the number of demand constraints q belongs to {2, 5, 10}. The third set includes 30 instances with n = 100, where both the number of knapsack constrains m and the number of demand constraints q equal to 30. The fourth set is composed of 30 large instances with n = 500, m = 30 and q = 30. Our TSTS algorithm employs four parameters, including γ 1 , γ 2 and γ 3 that are used to define the hash functions, and the depth α of tabu search used in the first search stage. The parameters γ 1 , γ 2 and γ 3 are set as in Table 1 according to the analysis shown in Section 4.3 while α is empirically set to 5 × 103 .

In addition, our algorithm was implemented in C++ and compiled by g++ compiler with -O3 option2 . All computational experiments were carried out on a computer with an Intel E5-2670 processor (2.5 GHz and 2G RAM), running the Linux operating system. Moreover, when running the DIMACS machine benchmark procedure dmclique 3 , our processor requires 0.19, 1.17, and 4.54 seconds to solve the graphs r300.5, r400.5, and r500.5, respectively. Finally, due to the stochastic nature of the algorithm, we independently ran the algorithm 30 times to solve each instance, where the time limit t max for each run was set to 60 seconds for instances with n ≤ 250 according to [START_REF] Gortázar | Black box scatter search for general classes of binary optimization problems[END_REF].

For the large instances with n = 500, the time limit t max was set to n seconds, where n is the size of instances.

Computational Results and Comparison

In this section, we report the computational results of our TSTS algorithm on the first two sets of benchmark instances. We provide in the Appendix the computational results of TSTS on the third and fourth sets of benchmark instances for which no detailed results are available for existing algorithms in the literature.

The computational results of our TSTS algorithm on the first set of benchmark instances with n = 100 are summarized in Table 2, together with the results of the Almha algorithm implemented in [START_REF] Gortázar | Black box scatter search for general classes of binary optimization problems[END_REF]. Column 1 gives the names of instances. Columns 2 and 3 indicate respectively the best known objective values (BKV ) and the results of the Almha algorithm, which were reported in [START_REF] Gortázar | Black box scatter search for general classes of binary optimization problems[END_REF] and available at http://www.optsicom.es/binaryss. The results of our TSTS algorithm are reported in columns 4-7, including the best objective value obtained over 30 runs (f best ), the average objective value (f avg ), the worst objective value (f worst ), and the standard deviation (σ f ) of objective values. The row Avg. shows the average result over all instances tested for each column. The rows #Better, #Equal, and #W orse respectively show the number of instances for which our best result f best is better than, equal to, worse than the BKV. Moreover, in terms of f best , our improved results (new lower bounds) are indicated in bold and our worse results are indicated in italic compared to the BKV . Finally, to verify whether there exists a significant difference between our best results (f best ) and the BKV, we provide in the last row the p-values (∈ [0, 1]) from the non-parametric Friedman test, where a p-value less than 0.05 means a significant difference between the compared results. This test was also performed to compare the results of the Almha algorithm and our average results (f avg ) 4 .

Table 2 shows that our TSTS algorithm matches the best known results for 47 out of 48 instances, and improves the best known result for the remaining one instance, leading to an improved Avg. value compared to the averaged BKV (18108.96 vs. 18108.10). When comparing the average objective values f avg of our TSTS algorithm with the results of the Almba algorithm, one can find that both algorithms have a very similar performance, which is confirmed by a large p-value of 0.87. In addition, regarding the average results over all instances (Avg.), the value of f avg of our TSTS algorithm is 18088.27, which is slightly better than that of Almba (i.e., 18083.17). These outcomes show that our TSTS algorithm performs similarly on these small instances n = 100 compared to the state-of-the-art Almba algorithm.

The second experiment aims to evaluate our TSTS algorithm on the set of 48 larger instances with n = 250, and the computational results are summarized in Table 3, where we report the same statistics as in Table 2. We observe from Table 3 that our algorithm matches the best known results for 18 out of 48 instances, improves the best known results for 12 instances, and misses the best known results for the remaining 18 instances. In terms of Avg., the value of f best of TSTS is slightly worse than the value of BKV (i.e., 49821.10 vs. 49836.48), and slightly better than the result of Almha (i.e., 49821.10 vs. 49717.81), which is however slightly better than the value of f avg of TSTS (i.e., 49687.60 vs. 49717.81). This experiment indicates that under the short time limit t max = 60 seconds, TSTS performs globally well on these larger instance with n = 250 especially by finding 12 improved best solutions. Additionally, we observe from Table 3 that TSTS performs particularly well for instances with a large number of constraints and achieves a better result than Almha for most of the 12 instances with 10 demand constraints and 10 knapsack constraints. Finally, we mention that the results of TSTS can be further improved by increasing the time limit (see the detailed results in the Appendix), implying that the current time limit (t max = 60) is too short for the TSTS algorithm on these instances.

Analysis and Discussions

To shed light on the functioning of the proposed algorithm, we now analyze and discuss several essential components of the TSTS algorithm.

Effectiveness and Robustness Analysis for Two-Stage Strategy

To study the effectiveness of the underlying two-stage search strategy, we summarize in Table 4 the computational results about the k values returned by the TSTS algorithm, where the results are based on the experiments in Section 3.3 and the value of k represents the number of selected items in the solution found. It is worth noting that the purpose of the first search stage is just to discover a promising hyperplane Ω [k] that contains high quality solutions, and the second search stage aims to locate improved solutions in the given hyperplane. Hence, the two-stage search strategy can be considered to be relevant and robust if the first search stage is able to reach very stably the identical or close hyperplane to the best known solution (i.e., Ω [k best] ).

The results of small instances with n = 100 are reported in the first 4 columns of Table 4, including the name of instances, the k value of the best solution obtained over 30 runs (k best ), the average k value of solutions obtained, and the standard deviation of k values obtained (σ k ). The results of larger instances with n = 250 are reported in columns 5-8, with the same statistics as in columns 1-4. In addition, the row Avg. shows the average results of standard deviations σ k of k values over all tested instances of each test set.

Table 4 shows that the value of k avg is very close to that of k best for most tested instances, which means that the first search stage of the TSTS algorithm is able to find a hyperplane that is very close to the best hyperplane containing the current best known solution. On the other hand, we observe that the standard deviations σ k of k values obtained are very small for most instances.

In particular, the average standard deviations of k values are respectively 0.22 and 0.66 for the set of small instances with n = 100 and the set of larger instances with n = 250. Hence, this experiment confirms to some extend the effectiveness and robustness of the two-stage search strategy employed by the TSTS algorithm.

Effects of Two Stages on the Performance of Algorithm

Table 5 Comparison between the two search stages on the instances with n = 250. To investigate the respective role of the two stages of our algorithm, we carried out an experiment based on the instances with n = 250. We ran our TSTS algorithm 30 times to solve each instance according to the experimental protocol of Section 3.2. The average results from the first stage and the second stage over 30 independent runs are summarized in Table 5. The first column gives the name of instances, columns 2 and 3 show the objective value (f 1 ) obtained by the first stage and computation time (t 1 ) in seconds needed to reach f 1 . Columns 4 and 5 show the objective value (f 2 ) obtained by the second stage and the computation time (t 2 ) needed to reach f 2 from f 1 . The last two columns indicate the gap between f 2 and f 1 and the improvement ratio (ρ) of f 2 over f 1 , which is calculated as ρ = 100 × (f 2 -f 1 )/f 1 .

Instance f 1 t 1 (s) f 2 t 2 (s) f 2 -f 1 ρ 250-5-2-0-
We observe from Table 5 that the first search stage of the TSTS algorithm is able to obtain a high-quality feasible solution for each instance and the solutions obtained in the first stage can be further improved during the second search stage. Furthermore, the improvements of f 2 over f 1 are significant with an average improvement ratio ρ close to 1%. On the other hand, regarding the computation times needed by the two search stages, we observe that most computational efforts are required by the second search stage and that t 2 is about five times larger than t 1 . Of course, this proportion depends also on the setting of the parameters α and t max . These outcomes indicate that both search stages of the TSTS algorithm are indispensable for the high performance of the algorithm. The first search stage is able to generate high-quality feasible solutions while the second search stage is able to further improve the solutions by performing an intensified search in the reached hyperplane.

Sensitivity Analysis of Hash functions

Now, we investigate the impacts of hash functions on the performance of the algorithm and discuss the sensitivity of the associated parameters. For this purpose, we carried out an additional experiment based on 30 representative instances in terms of the numbers of knapsack and demand constraints. We ran our TSTS algorithm 30 times for each instance and each parameter combination of (γ 1 , γ 2 , γ 3 ), where γ i (i = 1, 2, 3) are the parameters used to define the hash functions h i (see Section 2.3.2 for details). Specifically, we tested 10 different settings, i.e., (γ 1 , γ 2 , γ 3 ) ∈ {(1.1, 1.5, 1.8), (1.3, 1.5, 1.8), (1.3, 1.5, 2.0), (1.5, 2.0, 2.5), (1.6, 1.8, 2.0), (1.6, 1.8, 2.5), (1.8, 2.0, 2.2), (1.8, 2.0, 2.5), (1.9, 2.1, 2.3), (2.0, 2.2, 2.5)}.

The experimental results are summarized in Table 6, where the first column gives the name of instances, the second row shows the settings of parameters, and the average objective values (f avg ) obtained over 30 runs are reported in columns 2-11 for each parameter combination and each instance, respectively. In addition, the rows #Best and Avg. of the table indicate respectively the Table 6. Influence of the hash functions on the performance of algorithm. Each instance was independently solved 30 times for each parameter combination, and the average objective values (f avg ) over 30 runs are reported. number of instances for which the associated setting of parameters yielded the best results and the average results over all instances tested.

favg Instance/(γ 1 , γ 2 , γ 3 ) (1.1, 1.5, 1 .8) (1.3, 1.5, 1.8) (1.3, 1.5, 2.0) (1.5, 2.0, 2.5) (1.6, 1.8, 2.0) (1.6, 1.8, 2.5)
We observe from Table 6 that the algorithm is sensitive to the setting of the parameters γ 1 ,γ 2 , and γ 3 . For the parameter combinations containing two small parameter values, the TSTS algorithm performed badly. For example, the algorithm with the combination (1.1, 1.5, 1.8) yielded the worst results in terms of Avg. However, when all parameters in (γ 1 , γ 2 , γ 3 ) take a relatively large value, the algorithm obtained a much better performance. Taking (1.9, 2.1, 2.3) as an example, the algorithm achieved the best results on 9 instances. In summary, the parameters γ 1 ,γ 2 , and γ 3 have an important impact on the performance of the algorithm, and parameter combinations containing at least two large parameter values lead usually to a good performance of the algorithm.

Effectiveness of Solution-based Tabu Strategy

The solution-based tabu strategy is an essential ingredient of our TSTS algorithm. To show its effectiveness with respect to the popular attribute-based tabu strategy, we created a variant A-TSTS of the TSTS algorithm by replacing the solution-based tabu strategy with the popular attribute-based tabu strategy, while keeping the other TSTS components unchanged. In A-TSTS, the adopted tabu strategy can be simply described as follows. Given an incumbent solution s = (x 1 , x 2 , . . . , x n ), once a 0-1 variable x i (1 ≤ i ≤ n) is flipped as x i ← 1 -x i , x i is forbidden to be flipped again for the next tt iterations (tt is the tabu tenure) and the associated neighbor solutions are excluded for consideration during the period identified by the tabu tenure. We empirically set tt = C + rand[0, 2] where C is a parameter that takes the value of 20 and rand[0, 2] is a random integer in [0, 2]. Finally, the tabu status of a variable is disabled if flipping the variable leads to a solution better than all previously visited solutions (this is the so-called aspiration criterion).

To compare TSTS and A-TSTS, we carried out an experiment based on the set of 48 large instances with n = 250, where both algorithms were run 30 times according to the experimental protocol given in Section 3.2. The computational results are summarized in Table 7 where we show for each algorithm the best results (f best ) obtained over 30 runs, the average results (f avg ), and the worst results obtained (f worst ). In addition, the row 'Avg.' shows the average results for each performance indicator. Finally, to check whether there exists a significant difference between the two algorithms in terms of f best , f avg and f worst , we provide the p-values from the non-parametric Friedman test in the last row, where a p-value smaller than 0.05 implies a significant difference between the compared results. We observe from Table 7 that the solution-based algorithm TSTS dominates the attribute-based algorithm A-TSTS in terms of all indicators, by reporting better results in terms of f best , f avg and f worst on all the instances. Furthermore, the small p-values indicate that the performance differences between the compared results are statistically significant. Therefore, this experiment confirms that under the two-stage framework of this work, the solution-based tabu strategy is much more suitable than the attribute-based tabu strategy for solving the MDMKP.

Discussion about the Solution-based and Attribute-based Tabu Search Approaches

Attribute-based tabu search is a popular approach, whose key idea is to prevent one attribute or a combination of several attributes of a solution from being changed during a number of iterations since their last changes. For such a method, a tabu attribute or a combination of several attributes is usually associated with a number of tabu solutions. However, unlike attribute-based tabu search, solution-based tabu search tries to record all the visited solutions and prevent them from being revisited during the following search process. Thus, solution-based tabu search ensures a stronger intensification search ability.

Recent studies on several binary optimization cases including two dispersion problems (minimum difference dispersion [START_REF] Wang | Effective metaheuristic algorithms for the minimum differential dispersion problem[END_REF] and maximum min-sum dispersion [START_REF] Lai | Solution-based tabu search for the maximum min-sum dispersion problem[END_REF]) and the classic multidimensional knapsack problem [START_REF] Lai | A two-phase hybrid evolutionary algorithm for the 0-1 multidimensional knapsack problem[END_REF] demonstrate that solution-based tabu search is more suitable than attribute-based tabu search. On the other hand, our experience on another dispersion problem (max-mean dispersion [START_REF] Lai | A tabu search based memetic search algorithm for the max-mean dispersion problem[END_REF]) does not confirm the advantage of the solutionbased tabu search approach over the attribute-based tabu search approach. So one interesting question concerning the solution-based and attribute-based tabu search approaches is under what circumstances one approach will be more suitable than the other. Given that solution-based tabu search has been investigated only very recently, little knowledge is currently available, which makes it difficult to provide a meaningful guidance on the choice between these two approaches. To fully characterize these approaches and understand the relationships between these approaches and the optimization problem under consideration, more studies are clearly needed, which constitutes an interesting research perspective.

Finally, the ideas of the present solution-based tabu search algorithm being quite general, they could conveniently be tested on other binary optimization problems, by adjusting the γ parameters of the hash functions (Eq. ( 15), Section 2.3.2) or by increasing the number of the hash vectors and the associated hash functions.

In this work, we investigated the NP-hard multidemand multidimensional knapsack problem, by proposing a two-stage tabu search (TSTS) algorithm that combines two solution-based tabu search procedures and a penalty-based evaluation function to explore different search spaces. Computational results on 156 benchmark instances showed that the proposed algorithm is competitive compared to the state-of-art results in the literature, especially for those instances with a large number of knapsack and demand constraints.

The experimental analysis showed the usefulness of the two-stage search strategy and the respective impacts of two stages on the performance of the algorithm. The first search stage is able to reach a promising hyperplane containing high-quality solutions, and the second search stage is able to find elite solutions by an intensified examination of the given hyperplane. We also showed that the hash functions used by the tabu search algorithms are a key component that significantly influences the performance of the algorithm.

This work enriches the existing tools for effectively solving the MDMKP and invites more research and attention on solution-based tabu search for solving binary optimization problems in the future. Specifically, given that the ideas of the two-stage strategy and the solution-based tabu search procedures developed in this work are quite general, it would be interesting to investigate their effectiveness on other problems, especially those related to subset selection problems for which the search space can be divided into a series of hyperplanes. It is also interesting to study search strategies mixing solution-based and attribute-based approaches. As a more fundamental research perspective, studies on a characterization of both solution-based and attribute-based search approaches are needed, which could ease the choice of one or the other approach to solve additional binary optimization problems.

A Appendix

This Appendix presents the detailed computational results of the proposed TSTS algorithm for the third and fourth sets of benchmark instances for which no detailed results are available in the literature. These instances have 100 or 500 items, 30 knapsack constraints and 30 demand constraints, making them more difficult to solve. For each of these instances, the TSTS algorithm was run 30 times under the condition presented in Section 3.2 and the computational results are summarized in Table A.1, where the symbols have the same meanings as in the previous tables.In addition, we report in Table A.2 the computational results of our TSTS algorithm on instances with n = 250 under a long time limit of t max = 300 (instead of t max = 60 used in Section 3). The results reported in this Appendix can serve as references for future comparative studies of new MDMKP algorithms.

We observe from Table A.1 that TSTS is able to reach the best result (f best ) with a success rate of 100% for 15 out of 30 small instances with n = 100, which means a good robustness of our TSTS algorithm on these instances. Note that for some small instances with n = 100, m = 30, and q = 30, it is difficult for some state-of-the-art algorithms in the literature [START_REF] Hvattum | Alternating control tree search for knapsack/covering problems[END_REF] to obtain a feasible solution. However, it is very easy for our TSTS algorithm to obtain a feasible solution for all these instances. For other instances, the standard deviation (σ f ) of objective values obtained by the TSTS algorithm is relatively small, which indicates a good robustness of the algorithm. Regarding the value of k which represents the number of items selected, it can be seen that the gap between k best and k avg and the standard deviation (σ k ) of k values obtained are very small, implying that the two-stage strategy of the algorithm is very robust and effective.

2

  Two-stage tabu search algorithm for the MDMKP Our two-stage solution-based tabu search (TSTS) algorithm combines two search procedures working on two different search spaces. The first stage of the algorithm performs an exploratory search within the whole search space to find a feasible solution as good as possible. Starting from this solution, the second stage carries out a focused exploitation within the reduced space composed of candidate solutions with exactly k selected items (k being identified by the final solution of the first search stage). To explore both search spaces, TSTS relies on two solution-based tabu search procedures guided by a penaltybased evaluation function. One notices that the two-stage search strategy has been used with success to solve other knapsack problems like the Quadratic Knapsack Problem[START_REF] Chen | An iterated "hyperplane exploration" approach for the Quadratic Knapsack Problem[END_REF] and the classic MKP[START_REF] Vasquez | A hybrid approach for the 0-1 multidimensional knapsack problem[END_REF].The proposed TSTS algorithm is thus composed of two optimization stages, where the first stage identifies a suitable hyperplane Ω [k] (see Section 2.4.2 for the definition of hyperplane) that is exploited intensively during the second optimization stage to locate improved solutions (see Algorithm 1).

	found
	begin
	/* Initialization of solution	*/
	s ← InitialSolution(I)	/* Sections 2.2 */
	/* Optimization of the first stage	*/
	{s, t} ← T abuSearch 1 (s)	/* Sections 2.3 */
	/* Optimization of the second stage	*/
	s ← T abuSearch 2 (s, t, t max ) return s	/* Sections 2.4 */
	end

2.1 General Procedure

Algorithm 1: General procedure of two-stage tabu search algorithm for the MDMKP Function TSTS() Input: Instance I, time limit t max Output: The best solution s *

  Initial solution s, extended evaluation function F , hash vectors H 1 , H 2 , H 3 of length L, hash functions h 1 , h 2 , h 3 , depth of tabu search α, time limit t max Output: The best solution s * found, the time t consumed by the search

	Algorithm 3: Tabu search method used in the first search stage
	Function T abuSearch 1 ()
	Input: begin
	/* Initialization of hash vectors (i.e., tabu lists)	*/
	for i ← 0 to L -1 do
	H 1 [i] ← 0; H 2 [i] ← 0; H 3 [i] ← 0;
	end
	s * ← s
	N oImprove ← 0
	/* Main search procedure		*/
	while (time() < t max ) ∧ ((N oImprove < α) ∨ (s * is infeasible)) do
	Find in N 1 (s) N 2 (s) a best non-tabu solution s ′ in terms of the
	extended evaluation function F in Eq. (17)
	/* N 1 (s) and N 2 (s) are defined in Eqs. (13) and (14), and
	the tabu rule is given in Section 2.3.2	*/
	s ← s ′	/* Update the current solution */
	if (F (s) > F (s * ) then
	s * ← s
	N oImprove ← 0
	end
	/* Update the hash vectors (i.e., tabu lists) with s	*/
	H 1 [h 1 (s)] ← 1; H 2 [h 2 (s)] ← 1; H 3 [h 3 (s)] ← 1
	N oImprove ← N oImprove + 1
	end
	t ← time()
	return {s * ,t}
	end
	2.3.3, a best non-tabu neighbor solution is selected to replace the current solu-
	tion, and then the tabu lists are accordingly updated. Finally, the best feasible
	solution found s * during this search stage and the computation time elapsed
	t are returned as the results of the first search stage. As our experiments in
	Section 3 show, the first search stage always ends up with a feasible solution
	for all tested instances, including those with 30 demand constraints and 30
	knapsack constraints. In other words, the first solution-based tabu search al-
	gorithm is typically able to locate a promising valid hyperplane for the second
	search stage.
		). At each iteration, according to the tabu
	rule and the penalty-based evaluation function defined in Sections 2.3.2 and

  The tabu search method used in the second search stage Function T abuSearch 2 () Input: Initial solution s, extended evaluation function F , penalty function P , hash vectors H 1 , H 2 , H 3 of length L, hash functions h 1 , h 2 , h 3 , time limit t max , and time (t) consumed in the first search stage.

	Output: The best solution s * found
	begin
	/* Initialization of hash vectors (i.e., tabu lists)	*/

to evaluate the solutions in the search space Ω [k] . Algorithm 4:

  Find in N 3 (s) a best non-tabu solution s ′ in terms of the extended evaluation function F in Eq. (17)/* N 3 (s) is defined in Eq. (18), and the tabu rule is given in Section 2.4.1.

	Main search procedure	*/
	while time() < t max -t do

*/ s ← s ′ /* Update the current solution */ if F (s) > F (s * ) ∧ P (s) = 0 then s * ← s end /* Update the hash vectors (i.e., tabu lists) with s */

Table 1

 1 Settings of parameters

	Parameters	Section	Description	Values
	α	2.3	tabu depth of T abuSearch 1 ()	5 × 10 3
	γ 1	2.3	parameter used in the hash function	1.9
	γ 2	2.3	parameter used in the hash function	2.1
	γ 3	2.3	parameter used in the hash function	2.3

Table 2

 2 Computational results and comparison on the 48 instances with n = 100. In terms of f best , the improved results are indicated in bold compared to the best known objective values (BKV ).

	TSTS (this work)

Table 3

 3 Computational results and comparison on the 48 instances with n = 250. In terms of f best , the improved results are indicated in bold and the worse results are indicated in italic compared to the best known objective value (BKV ).

	TSTS (this work)

Table 4

 4 Statistical results over 30 runs in terms of the number k of items in the obtained solutions.

		n=100				n=250
	Instance	k best	kavg	σ k	Instance	k best	kavg	σ k
	100-5-2-0-0	30	29.97	0.30	250-5-2-0-0	79	78.47	0.85
	100-5-2-0-1	31	31.00	0.34	250-5-2-0-1	78	76.00	0.89
	100-5-2-0-2	31	30.83	0.34	250-5-2-0-2	78	77.03	1.08
	100-5-2-0-3	32	31.83	0.42	250-5-2-0-3	79	79.07	0.44
	100-5-2-0-4	31	31.00	0.00	250-5-2-0-4	79	78.33	0.54
	100-5-2-0-5	56	55.83	0.40	250-5-2-0-5	137	134.87	1.50
	100-5-2-1-0	28	27.00	0.30	250-5-2-1-0	71	70.03	0.66
	100-5-2-1-1	29	29.00	0.00	250-5-2-1-1	69	69.23	0.67
	100-5-2-1-2	27	27.27	0.34	250-5-2-1-2	72	72.07	0.63
	100-5-2-1-3	29	29.00	0.00	250-5-2-1-3	69	68.57	0.72
	100-5-2-1-4	29	28.17	0.44	250-5-2-1-4	70	70.67	0.83
	100-5-2-1-5	54	53.00	0.31	250-5-2-1-5	129	127.70	1.07
	100-5-5-0-0	28	28.40	0.46	250-5-5-0-0	76	75.50	0.50
	100-5-5-0-1	30	30.00	0.00	250-5-5-0-1	74	71.87	1.02
	100-5-5-0-2	30	30.00	0.00	250-5-5-0-2	74	74.37	0.55
	100-5-5-0-3	30	30.00	0.00	250-5-5-0-3	76	75.40	0.66
	100-5-5-0-4	29	29.00	0.00	250-5-5-0-4	76	76.00	0.26
	100-5-5-0-5	54	54.00	0.18	250-5-5-0-5	136	134.60	0.80
	100-5-5-1-0	27	26.73	0.48	250-5-5-1-0	68	66.90	0.65
	100-5-5-1-1	28	28.00	0.00	250-5-5-1-1	66	65.40	0.55
	100-5-5-1-2	27	27.37	0.47	250-5-5-1-2	65	65.53	0.62
	100-5-5-1-3	28	27.93	0.34	250-5-5-1-3	65	65.80	0.65
	100-5-5-1-4	27	27.07	0.18	250-5-5-1-4	66	66.90	0.60
	100-5-5-1-5	52	52.00	0.00	250-5-5-1-5	126	126.23	0.80
	100-10-5-0-0	28	28.00	0.00	250-10-5-0-0	69	67.20	1.05
	100-10-5-0-1	28	27.40	0.50	250-10-5-0-1	71	71.07	0.44
	100-10-5-0-2	27	26.97	0.00	250-10-5-0-2	70	69.23	0.56
	100-10-5-0-3	28	27.43	0.42	250-10-5-0-3	68	67.33	0.94
	100-10-5-0-4	28	27.07	0.25	250-10-5-0-4	69	67.53	0.96
	100-10-5-0-5	53	52.60	0.47	250-10-5-0-5	130	128.93	0.81
	100-10-5-1-0	26	26.00	0.18	250-10-5-1-0	65	65.13	0.62
	100-10-5-1-1	26	26.03	0.00	250-10-5-1-1	67	66.73	0.63
	100-10-5-1-2	26	26.03	0.25	250-10-5-1-2	65	63.97	0.71
	100-10-5-1-3	26	26.13	0.30	250-10-5-1-3	65	65.03	0.60
	100-10-5-1-4	27	26.47	0.50	250-10-5-1-4	66	65.43	0.50
	100-10-5-1-5	51	50.93	0.18	250-10-5-1-5	128	127.20	0.40
	100-10-10-0-0	28	28.00	0.00	250-10-10-0-0	68	66.83	0.58
	100-10-10-0-1	27	27.00	0.00	250-10-10-0-1	68	68.47	0.56
	100-10-10-0-2	27	27.00	0.00	250-10-10-0-2	67	66.70	0.64
	100-10-10-0-3	27	27.00	0.26	250-10-10-0-3	68	67.30	0.59
	100-10-10-0-4	28	27.30	0.45	250-10-10-0-4	67	66.13	0.56
	100-10-10-0-5	53	52.97	0.00	250-10-10-0-5	130	129.53	0.56
	100-10-10-1-0	24	24.77	0.42	250-10-10-1-0	64	63.70	0.46
	100-10-10-1-1	25	25.03	0.00	250-10-10-1-1	63	63.03	0.18
	100-10-10-1-2	26	26.00	0.00	250-10-10-1-2	64	63.63	0.48
	100-10-10-1-3	25	25.03	0.37	250-10-10-1-3	63	63.40	0.49
	100-10-10-1-4	25	25.17	0.30	250-10-10-1-4	65	64.00	0.37
	100-10-10-1-5	50	49.83	0.30	250-10-10-1-5	123	123.70	0.46
	Avg.			0.22				0.66

Table 7

 7 Comparison between the solution-based and attribute-based tabu strategies on the instances with n = 250. The better results between the two algorithms are indicated in bold in terms of f best , f avg and f worst .

			f best	favg		fworst
	Instance	A-TSTS	TSTS	A-TSTS	TSTS	A-TSTS	TSTS
	250-5-2-0-0	72278	78486

Table A .

 A 2 shows that our TSTS algorithm improves the best known results for 16 out of 48 instances, matches the best known results for 24 instances, and misses the best known results for only 8 instances. These results indicate Table A.1 Computational results of the TSTS algorithm on the instances with a large number of constraints under the condition presented in Section 3.2. For each of these instances, there are 30 knapsack constraints and 30 demand constraints.Table A.2Computational results of the TSTS algorithm on the instances with n = 250 under a time limit of t max = 300 seconds. In terms of f best , the improved results are indicated in bold and the worse results are indicated in italic compared to the best known objective value (BKV ) reported in the literature.

	Instance	f best	favg	fworst	σ f	k best	kavg	σ k
	100-30-30-0-2-1	11312	11300.00	11252	24.00	25	25.20	0.40
	100-30-30-0-2-2	9945	9945.00	9945	0.00	25	25.00	0.00
	100-30-30-0-2-3 100-30-30-0-2-4 Instance 100-30-30-0-2-5 250-5-2-0-0 100-30-30-0-2-6 250-5-2-0-1 100-30-30-0-2-7 250-5-2-0-2	11195 11324 BKV 9704 78486 23296 75132 22442 71003	11130.33 11290.40 f best 9704.00 78486 23296.00 75132 22224.40 70898	10225 11198 favg 9704 78282.30 23296 74849.07 22126 70676.77	241.96 55.72 fworst 0.00 77986 0.00 74152 125.49 70137	25 25 σ f 25 173.14 50 232.22 51 217.86	25.07 25.00 k best 25.00 79 50.00 78 50.27 78	0.25 kavg 0.00 78.17 0.00 76.10 0.00 0.44 76.87	σ k 0.73 0.94
	100-30-30-0-2-8	23452	23312.90	23182	45.57	51	50.13	0.34
	100-30-30-0-2-9	22756	22756.00	22756	0.00	50	50.00	0.00
	100-30-30-0-2-10	24371	24323.60	24287	28.16	50	50.00	0.00
	100-30-30-0-2-11	33472	33472.00	33472	0.00	75	75.00	0.00
	100-30-30-0-2-12	32670	32670.00	32670	0.00	75	75.00	0.00
	100-30-30-0-2-13	32942	32942.00	32942	0.00	75	75.00	0.00
	100-30-30-0-2-14	35106	35106.00	35106	0.00	75	75.00	0.00
	100-30-30-0-2-15	30930	30788.73	30767	55.41	75	75.00	0.00
	100-30-30-1-5-1	5340	5340.00	5340	0.00	26	26.00	0.00
	100-30-30-1-5-2	4390	4390.00	4390	0.00	25	25.00	0.00
	100-30-30-1-5-3	4227	4227.00	4227	0.00	25	25.00	0.00
	100-30-30-1-5-4	4706	4433.40	4424	50.62	25	25.97	0.18
	100-30-30-1-5-5	2597	2591.90	2546	15.30	25	25.00	0.00
	100-30-30-1-5-6	10808	10647.60	10462	92.10	50	50.00	0.00
	100-30-30-1-5-7	9807	9789.37	9766	5.37	51	50.03	0.18
	100-30-30-1-5-8	10882	10865.47	10753	42.27	50	50.10	0.30
	100-30-30-1-5-9	10595	10595.00	10595	0.00	50	50.00	0.00
	100-30-30-1-5-10	10297	10285.87	9963	59.95	50	50.03	0.18
	100-30-30-1-5-11	11029	11029.00	11029	0.00	75	75.00	0.00
	100-30-30-1-5-12	11884	11823.27	11747	45.74	75	75.00	0.00
	100-30-30-1-5-13	10751	10751.00	10751	0.00	75	75.00	0.00
	100-30-30-1-5-14	11567	11567.00	11567	0.00	75	75.00	0.00
	100-30-30-1-5-15	10671	10368.47	10351	59.75	75	75.00	0.00
	500-30-30-0-2-1	85188	84991.27	84418	185.69	128	127.57	0.62
	500-30-30-0-2-2	82073	81856.17	81504	141.84	129	128.07	0.44
	500-30-30-0-2-3	77393	76995.10	76516	205.93	129	127.40	0.80
	500-30-30-0-2-4	82304	82049.27	81628	158.83	128	127.37	0.75
	500-30-30-0-2-5	83525	83300.07	82779	170.66	129	127.87	0.72
	500-30-30-0-2-6	145967	145705.17	145474	88.16	253	252.73	0.68
	500-30-30-0-2-7	152246	152019.70	151665	132.72	253	252.67	0.79
	500-30-30-0-2-8	157687	157487.13	157087	129.19	254	253.30	0.82
	500-30-30-0-2-9	153751	153548.53	153365	87.25	254	252.33	0.75
	500-30-30-0-2-10	142173	141943.90	141721	110.37	253	252.67	0.75
	500-30-30-0-2-11	185226	184986.67	184781	91.72	377	376.97	0.55
	500-30-30-0-2-12	194614	194444.30	194275	72.88	376	376.87	0.67
	500-30-30-0-2-13	208246	208129.67	207904	76.15	378	377.53	0.72
	500-30-30-0-2-14	215849	215693.30	215381	90.85	378	377.67	0.65
	500-30-30-0-2-15	194224	194037.03	193788	91.74	376	376.57	0.62
	500-30-30-1-5-1	51666	51574.93	51359	65.87	126	126.50	0.56
	500-30-30-1-5-2	50101	49871.50	49566	123.11	126	126.33	0.54
	500-30-30-1-5-3	51226	50979.60	50842	88.53	126	126.13	0.62
	500-30-30-1-5-4	51637	51483.20	51261	91.69	127	126.77	0.42
	500-30-30-1-5-5	52078	51860.97	51655	103.52	128	127.00	0.63
	500-30-30-1-5-6	84052	83834.80	83548	104.57	251	251.03	0.41
	500-30-30-1-5-7	82850	82637.00	82342	116.19	250	250.87	0.56
	500-30-30-1-5-8	82722	82576.90	82419	74.14	250	250.70	0.46
	500-30-30-1-5-9	82825	82451.43	81944	146.82	250	250.23	0.62
	500-30-30-1-5-10	82845	82559.93	82098	143.15	249	249.67	0.60
	500-30-30-1-5-11	88887	88756.90	88556	75.33	374	374.23	0.42
	500-30-30-1-5-12	87254	87136.80	86815	82.34	374	374.20	0.65
	500-30-30-1-5-13	87315	87167.80	87028	61.20	375	374.80	0.54
	500-30-30-1-5-14	87583	87486.57	87261	89.89	374	373.90	0.65
	500-30-30-1-5-15	87956	87814.97	87616	72.05	374	374.10	0.40
	Avg.	62265.52	62139.10	61957.25	70.33	150.88 150.78	0.34

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

The source code of our TSTS algorithm will be available at our website: http: //www.info.univ-angers.fr/pub/hao/mdmkp.html

ftp://dimacs.rutgers.edu/pub/dsj/clique

Since the results of Almha are based on only one run, we mainly use our average results for this comparison.

Acknowledgments

We are grateful to the reviewers for their valuable comments which helped us to improve the paper. This work is partially supported by the National Natural Science Foundation of China (Grant No. 61703213), the Natural Science Foundation of Jiangsu Province of China (Grant No. BK20170904), and NUPTSF (Grant No. NY217154).