S Samaddar

I Yudhistira

S Adam

H Courtois

Clemens B Winkelmann
email: clemens.winkelmann@neel.cnrs.fr

Charge Puddles in Graphene near the Dirac Point

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

Electrons in graphene are subjected to a disordered potential created by random charged impurities, either adsorbed on the graphene or buried in the substrate. These lead to inhomogeneities in the local carrier density, that is, charge puddles [1][2][3][4][5][6][7]. Scanning probe microscopies have in particular strongly contributed to unraveling the puddles' spatial properties and have challenged different theories about their origin [1][2][3]. Charge puddles are usually thought of as a limitation to the extent the charge neutrality point can be approached macroscopically, thereby also limiting possible device performances. However, the behavior of the puddles itself unveils the fascinating electronic properties of graphene and, more generally, Dirac materials.

Electrostatic screening in two dimensions (2D) has a counterintuitive behavior. Thomas-Fermi screening entails a characteristic length scale q -1 TF . Unlike in three dimensions, the Thomas-Fermi wave vector q TF in a 2D electron gas (2DEG) is proportional to the density of states at the Fermi level. As a consequence, q TF is energy-and thus carrier density-independent in conventional 2DEGs, while q TF is proportional to k F = √ πn in graphene and other Dirac materials. This has the important consequence that the unscreened potential created by a charged impurity in a medium with dielectric constant κ, V (q) = e 2 /2κ 0 q and the screened potential Ṽ (q) ∝ (q + q TF ) -1 are identical within a multiplicative constant [7]. In other words, near charge neutrality local inhomogeneities in the screened potential can be arbitrarily large. Further, a rough estimate of the lateral extent of charge carrier density correlations is given by q -1 TF , from which a strong growth ∝ 1/ √ n of the puddles size is expected near charge neutrality. The carrier density dependence of both the charge puddles' amplitude and size in a Dirac material has not been reported yet.

In this Letter, we report the direct microscopic observation of the doping disorder landscape in monolayer graphene at different charge carrier densities. The charge inhomogeneities are found to grow, both in spatial extent and in amplitude, as the Fermi level approaches the Dirac point. From transport measurements on the very same graphene sample at study, the microscopic parameters of the disorder potential can be estimated in the frame of the self-consistent screening theory. Calculations of the charge puddles distribution based on these are in very good agreement with the experimental observations. The sample is fabricated on a heavily doped Si substrate covered with thermal oxide. Single layer graphene is prepared by mechanical exfoliation [8]. The number of graphene layers and the absence of surface contamination are confirmed from combined optical, Raman and ex-situ AFM characterization. Using a mechanical mask [9], we deposit the metallic source and drain contacts to form a 4 µm long graphene junction (Fig. 1a). Organic resist is avoided, as to ensure a residue-free surface for scanning probe microscopy. Details of the device fabrication are described in the Supplemental Material file.

The experimental setup is a home-made combined AFM-STM operating at a temperature of 130 mK [10], at which all measurements presented here were obtained. The sample stage allows for in situ multi-terminal transport measurements of the device. AFM is performed by electrical excitation and read-out of a mechanical quartz Length Extension Resonator (LER) [11,12]. This allows to rapidly move the tip to the graphene junction with the help of the position code [13,14]. Scanning tunneling micrographs reveal a clean graphene surface, following the substrate corrugation with a roughness of about 100 pm (see Supplemental Material file for details). Scanning tunneling spectroscopy (STS) is achieved by lock-in measurements of the differential tip-to-sample tunneling conductance G t (x, y) = dI t /dV b , by adding a 12 mV ac modulation at frequency f = 322 Hz to the bias voltage V b , which is uniformly applied to the sample.

Transport measurements are performed with the tip retracted a few hundreds of nm from the sample sur- face. However, approaching the tip to STM contact does not produce a significant effect in the device characteristics. The conductivity of the graphene device shows a perfectly linear behavior at high carrier densities (Fig. 1b), in line with a density independent mobility of about 6000 cm 2 V -1 s -1 . This indicates that carrier transport is dominated by long range disorder, as can be caused by random charge impurities in the substrate [15]. A slight difference between the measured electron and hole mobilities, µ e /µ h = 0.9 ± 0.05 could be related to the difference in their scattering cross-sections off charged impurities [5,16].

The gate voltage at which conductivity is minimized gives the overall charge neutrality point, V 0 g = 29 V. This overall hole doping is consistent with the presence of negatively charged silanol groups on the surface. A residual conductivity σ 0 ≈ 11 e 2 /h is found at the charge neutrality point. Within self-consistent screening theory [7], the above values of residual conductivity and mobility point to a charged impurity distribution with a density n i = 7.5 ± 0.5 × 10 11 cm -2 at a distance 0.1 nm < d < 1 nm below the graphene, in agreement with earlier experiments in similar conditions [5,6,15].

We performed scanning tunneling spectroscopy on the graphene sheet, at distances greater than 1 µm from the metal-graphene interface as to rule out any possible influence of the leads on local properties. Fig. 1c shows the differential tunneling conductance G t (V b ) acquired at a given location, but at different gate voltages V g . A Vshape spectrum, characteristic of graphene, is obtained in every case. A frequently reported gate-independent depression of the tunneling conductance is seen at zero bias [10,17,18]. In addition, the curves display two gatedependent local minima, highlighted by red and black arrows respectively, which move in opposite directions with V g .

The primary minimum V 0 b (red arrows) occurs when the Fermi level of the tip is aligned with the local Dirac point E D (r) of graphene, which can be written as:

V 0 b = -γ sign(V g -V D g ) |V g -V D g |, (1)
where γ = hv F πκ 0 /(e 3 t) employing a plate-capacitor model, with v F = 1.1×10 6 m/s the graphene Fermi velocity, t = 285 nm and κ = 3.9 the SiO 2 thickness and dielectric constant, respectively (details are described in the Supplemental Information file). The local quantity V D g is the gate voltage at which the Fermi levels of both the graphene and the tip are aligned with the Dirac point. It includes the influence of the local gating produced by the tip due to both the tip-sample work function mismatch and the bias voltage [9,19]. In the absence of capacitive coupling to the tip, the spatially averaged V D g would coincide with V 0 g found from transport experiments. In the case of Fig. 1, the experimental gate dependence of the primary minimum V 0 b can be well fitted with Eq. (1), yielding V D g = 28 V (Fig. 1d), the value of γ being determined by known parameters. The nearly exact matching of V D g at that particular position and tip condition with the global value of V 0 g is coincidental, since V D g depends on the position. The secondary minimum (black arrows) occurs when the Fermi level of graphene passes through the Dirac point [9,[18][19][20]. The above analysis provides a detailed understanding of the electron tunneling spectra dependence on the gate potential, at a given location.

Several strategies can be used for mapping the local Dirac point. Performing a complete spectrum with open feedback loop at each position, from which E D is then individually extracted, is the most reliable method but is very time-consuming [2, 10]. Mapping G t with closed feedback loop set to a fixed set-point tunnel current I set t and a single V b that is slightly offset from the average primary minimum V 0 b by δV b , was shown to reproduce qualitatively the E D (x, y) maps [2]. This stems from the fact that, to first approximation, a shift in E D simply shifts the G t (V b ) curves along the V b axis. Complications with this approach arise when one wishes to compare E D maps at different gate voltages because V 0 b itself is a function of V g . Our strategy consists in first determining V 0 b (V g ) at a given position (Fig. 1d) and then mapping G t at a gate dependent bias voltage

V b (V g ) = V 0 b (V g ) ± 100 mV. The sign of the offset is set such that |V b | > |V 0 b |.
We further normalize the differential tunnel conductance maps to the set-point tunnel conductance, writing Gt = (V b /I set t ) G t . This normalization is known to provide a more faithful conversion of the differential tunnel conductance to the density of states when the set-point tunnel conductance is not fixed from map to map [START_REF] Wiesendanger | Scanning Probe Microscopy and Spectroscopy[END_REF]. We have verified the structures observed in Gt maps to be consistent with the E D maps found from full current imaging tunneling spectroscopy (CITS) measurements, which were acquired at selected gate voltages. These full CITS also allow for determining the proportionality factor between the Gt and E D maps (see Supplemental Material file for a more detailed discussion).

Experimental maps of the variations of Gt ∼ E D around their spatially averaged value are shown for several gate voltages in Fig. 2. It is readily seen that not only the lateral extent, but also the amplitude of the doping inhomogeneities, gradually increase as the Fermi level approaches the Dirac point. For proper quantification of the observed inhomogeneities, we introduce the auto-correlation function of the E D maps. Assuming rotational symmetry (which is only approximately true, due to the finite size of the maps), the latter is a function of only r = |r|. The charge puddles' size ξ is determined from fitting the angular average of the auto-correlation function of each E D map to a gaussian decay. The gate dependence of ξ displayed on Fig. 3a shows a strong increase near charge neutrality, which is found at a gate voltage V D g of about 38 V. This value is a spatially averaged property of the map area. Because of the capacitive influence of the tip [START_REF]back gate can be estimated to β ≈ 50, meaning that a tipsample work function difference by 0.2[END_REF], V D g is somewhat larger than the charge neutrality condition V 0 g = 29 V found from transport experiments [9,19].

We further determine the standard deviation of the Dirac point variations δE D over a map. This quantity, which reflects the amplitude of the doping inhomogeneities across the sample, is plotted on Fig. 3b as a function of V g and also shows a marked peak at V D g ≈ 38 V. The error bars on ξ and δE D are mainly associated to the finite size of the maps; a detailed discussion of their determination can be found in the Supplemental Information file. Some asymmetry of the puddles' behavior is observed, which appear somewhat larger and more prominent at large electron doping, than on the hole doped side. As electron doping involves a quite large gate potential of about 60 V, a possible scenario for this asymmetry is that the back-gate eventually influences the substrate impurities distribution itself [6].

55 V 45 V 37.5 V 27.5 V 20 V -10 V (a) (b) (c) (d) (e) (f)
Our main experimental findings are thus that both the amplitude and the spatial extent of the puddles significantly increase as the Fermi level approaches the Dirac point. For a quantitative understanding, we now compare these results to calculations. From Thomas-Fermi theory in 2D, assuming a flat Fermi sea, follows that variations in the local value of E D /e are equal to variations in the screened electrostatic potential Ṽ [START_REF] Ashcroft | Solid State Physics[END_REF]. For a random 2D distribution of charged impurities with density n i at a distance d from the graphene sheet, the auto-correlation function of the screened potential can be written [START_REF] Adam | [END_REF] C(r) = 2πn i e 2 4π 0 κ 2 +∞ 0 q dq 1 (q) e -qd q 2 J 0 (qr),

(2) where J 0 is the zeroth-order Bessel function and (q) is the graphene dielectric function. The latter describes the screening of Dirac fermions which, within Random Phase Approximation (RPA), can be written as [15] (q) = 1 + 4k F r s /q for q < 2k F , 1 + πr s /2 for q > 2k F ,

where r s = e 2 /(4π 0 κhv F ) ≈ 0.8 on SiO 2 is the effective fine structure constant of graphene. The dependence of the correlation function on the mean doping level (and thus on the gate potential) enters here through the dependence of (q) on 2k F . We calculated the auto-correlation function for the screened potential and extracted the correlation length ξ and fluctuation amplitude δ Ṽ = δE D /e. The result for d = 1 nm, shown as dashed curves in Fig. 3a,b, accounts for the overall decrease of both ξ and δE D , that is, stronger screening with increasing charge carrier density. The calculations have no other adjustable parameter than V D g = 38V , the impurity density in the substrate n i being determined from the transport measurements. The puddles size follows in particular the expected trend ξ ∼ q -1 T F ∝ 1/ √ n at high carrier densities, where n ∝ |V g -V D g | is the gate induced charge carrier density. This agreement validates the microscopic picture of random potential fluctuations, for the description of which we call for Thomas-Fermi screening in a Dirac material.

At charge neutrality, for a homogeneous system, there are no excess charges available to screen the impurity potential. Accordingly, Eq. ( 3) predicts that both the amplitude and the correlation length diverge. However, this ignores the fact that the induced charges within the puddles can themselves screen the impurity potential. Accounting for this process self-consistently [15] leads to rewriting the RPA dielectric function of Eq. ( 3) with a corrected charge carrier density. The usual expression of k F = √ πn is then replaced by π(n + n * ) [START_REF] Adam | [END_REF], where n * represents the disorder-induced residual charge carrier density in the graphene sheet which cannot be compensated by a global gate. The self-consistent calculations are plotted in Fig. 3a,b. The gray regions are delimited by the theoretical curves for d = 0.1 nm and 1 nm respectively. The ensuing saturation of both the puddles' size and amplitude at the charge neutrality point is in very good agreement with the experimentally observed trend.

To conclude, this work provides the first microscopic observation of the growth of charge inhomogeneities in graphene near the Dirac point. It further shows that the observed behavior can be very well described with a theory based on a microscopic description of the impurity potential, using parameters found from transport measurements, performed in situ on the very same graphene sample. This observation gives utmost credit to the charged impurity potential scenario as a limiting factor to the exploitability of Dirac physics in graphene on SiO 2 [25].

This work was funded by the European Commission under project no. 264034 (Q-NET Marie Curie Initial Training Network). Work in Singapore was supported by the National Research Foundation of Singapore under its Fellowship programme (NRF-NRFF2012-01). Samples were fabricated at the Nanofab facility at Institut Néel. We thank S. C. Martin, B. Sacépé, A. de Cecco, J. Seidemann and A. K. Gupta for discussions and help with the experiments.

  FIG. 1. (a) Experimental configuration, combining transport measurements with scanning probe microscopy on a monolayer graphene device. A conductive tip mounted on a stiff mechanical resonator serves as the probe for combined AFM-STM. Two isolated electrical contacts (Source and Drain) enable two-probe transport measurements. The atomic force micrograph shows both the topography (vertical scale, varying between 0 and 57 nm) and the phase (varying by 3.8 • from blue to brown). (b) Device conductivity as a function of backgate voltage Vg, measured at a bias voltage of 5 mV. The overall charge neutrality point is found at V 0 g = 29 V. (c) Differential tip-to-sample tunneling conductance Gt as a function of the voltage V b (uniformly applied to the sample), at several values of Vg. The red and black arrows indicate the position of the primary (V 0 b ) and the secondary minimum of Gt respectively (see text). I set t = 50 pA at V b = 0.9 V. The curves are vertically offset for clarity. (d) Variation of the primary minimum with Vg. The black dashed line shows the fit with Eq. (1), yielding the fit parameter V D g = 28 V.

FIG

  FIG. 2. Spatial maps of the variations of the local Dirac point energy ED around its spatially averaged value ĒD, over an area of 100 × 100 nm 2 at different gate voltages (indicated below each figure). The imaging parameters for the original Gt maps (see text) are I set t = 50 pA and bias voltages V b equal to (a) 0.298 V, (b) 0.191 V (c) 0.122 V, (d) -0.145 V, (e) -0.232 V and (f) -0.267 V.

FIG. 3 .

 3 FIG. 3. (a) Puddle size and (b) amplitude of doping inhomogeneities, as a function of the back-gate potential. The dashed lines are calculations without self-consistent correction to carrier density, with d = 1 nm. The gray areas represent the range of values expected from a random impurity distribution at a distance d ranging from 0.1 to 1 nm from the graphene sheet. They are delimited by solid lines corresponding to calculations using the extremal values of d, respectively. For all calculations, the value derived from transport data is used: ni = 7.5 × 10 11 cm -2 .


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

Günther R Raidl

An Evolutionary Algorithm for Column Generation in Integer Programming: an Effective Approach for 2D Bin Packing

We consider the 3-stage two-dimensional bin packing problem, which occurs in real-world problems such as glass cutting. For it, we present a new integer linear programming formulation and a branch and price algorithm. Column generation is performed by applying either a greedy heuristic or an Evolutionary Algorithm (EA). Computational experiments show the benefits of the EA-based approach. The higher computational effort of the EA pays off in terms of better final solutions; furthermore more instances can be solved to provable optimality.

Introduction

The Two-Dimensional Bin Packing (2BP) problem occurs in different variants in important real-world applications such as glass, paper, and steel cutting. A recent survey on 2D packing problems is given in Lodi et al. [START_REF] Lodi | Recent advances on two-dimensional bin packing problems[END_REF]. Among the algorithms for exactly solving the general 2BP problem are the branch and bound algorithm of Martello and Vigo [START_REF] Martello | Exact solutions of the two-dimensional finite bin packing problem[END_REF] and the hybrid Branch and Price / Constraint Programming algorithm presented by Pisinger and Sigurd [START_REF] Pisinger | Using decomposition techniques and constraint programming for solving the two-dimensional bin packing problem[END_REF].

In many cases there is a special requirement on the cutting patterns: only orthogonal guillotine cuts are allowed, i.e., pieces may only be cut horizontally or vertically from one border to the one opposite. Furthermore, the number of stages of such cuts, i.e., the height of the cutting tree of each bin, is often limited in real-world applications. The case of two-stage cutting was first considered by Gilmore and Gomory [START_REF] Gilmore | Multistage cutting-stock problems of two and more dimensions[END_REF]. More recently two-stage 2BP was considered in Lodi et al. [START_REF] Lodi | Models and bounds for two-dimensional level packing problems[END_REF] and Belov and Scheithauer [START_REF] Belov | A branch-and-cut-and-price algorithm for onedimensional stock cutting and two-dimensional two-stage cutting[END_REF]. Three-stage cutting problems were treated in Vanderbeck [START_REF] Vanderbeck | A nested decomposition approach to a 3-stage 2-dimensional cutting stock problem[END_REF] and Puchinger et al. [START_REF] Puchinger | Solving a real-world glass cutting problem[END_REF], where particular real-world problems with specific additional properties were considered.

In Sec. 2 we present an Integer Linear Programming (ILP) model for classical 3-stage 2BP, based on the model of [START_REF] Lodi | Models and bounds for two-dimensional level packing problems[END_REF]. In Sec. 3 a column generation formulation and a Branch and Price (B&P) framework, based on [START_REF] Pisinger | Using decomposition techniques and constraint programming for solving the two-dimensional bin packing problem[END_REF], are proposed. We describe a greedy heuristic in Sec. [START_REF] Lodi | Models and bounds for two-dimensional level packing problems[END_REF] and an evolutionary algorithm in Sec. 5 for solving the pricing problem within the B&P approach, i.e., for generating new columns. In Sec. 6 experimental results are given and analyzed.

...void space

Three-Stage Two-Dimensional Bin Packing

The 2BP problem consists of a set of n rectangular items, each having a height h i and a width w i , i = 1, ..., n. The objective is to pack them into a minimum number of rectangular bins, each having height H and width W . Items may not overlap and we do not consider rotation.

A feasible layout for 3-stage 2BP consists of a set of bins, each bin consists of a set of stripes, each stripe consists of a set of stacks, and each stack consists of items having equal width. Every such pattern can be reduced into its so-called normal form by moving each item to its uppermost and leftmost position, so that void space appears only at the bottom of stacks, to the right of the last stack in each stripe and below the last stripe, see Fig. 1. In the sequel we consider only patterns in normal form.

In [START_REF] Lodi | Models and bounds for two-dimensional level packing problems[END_REF] a polynomial-sized ILP model for 2-stage 2BP has been proposed. We extend this model in order to get a polynomial-sized ILP formulation for 3-stage 2BP.

The items are sorted so that h 1 ≥ h 2 ≥ . . . ≥ h n . The order of the items within each stack is not relevant, so they can always be ordered according to their indices. A solution may contain at most n stacks. We label each stack with the index of the highest item it contains, i.e., the smallest item index. Similarly, a solution has at most n stripes, and a stripe's label is the label of its highest stack. Finally at most n bins are needed and we label each of them with the smallest index of the stripes it contains. The model uses the following 0/1-variables:

-α j,i , j = 1, . . . , n, i = j, . . . , n: rectangle i is contained in stack j; -β k,j , k = 1, . . . , n, j = 1, . . . , n: stack j is contained in stripe k; -γ l,k , l = 1, . . . , n, k = l, . . . , n: stripe k is contained in bin l; -δ l,i,j , l = 1, . . . , n -1, i = l + 1, . . . ,
n, and j = l, . . . , i -1: item i contributes to the total height of all stripes in bin l; i.e., item i appears in stack j, stack j appears in stripe j, and stripe j appears in bin l.

The 3-stage 2BP problem can now be stated as the following ILP:

minimize n l=1 γ l,l (1)
subject to

i j=1 α j,i = 1, ∀i = 1, . . . , n (2)
n i=j+1 α j,i ≤ (n -j)α j,j , ∀j = 1, . . . , n -1 (3) α j,i = 0, ∀j = 1, . . . , n -1 ∀i > j | w i = w j ∧ h i + h j > H (4) n k=1 β k,j = α j,j ∀j = 1, . . . , n (5)
n i=j h i α j,i < n i=k h i α k,i + (H + 1)(1 -β k,j ), ∀k = 2, . . . , n, ∀j = 1, . . . , k -1 (6) n i=j h i α j,i ≤ n i=k h i α k,i + H(1 -β k,j ), ∀k = 1, . . . , n -1, ∀j = k + 1, . . . , n (7)
n j=1 w j β k,j ≤ W β k,k , ∀k = 1, . . . , n (8)
k l=1 γ l,k = β k,k , ∀k = 1, . . . , n (9)
n i=l h i   γ l,i + i-1 j=l δ l,i,j   ≤ Hγ l,l , ∀l = 1, . . . , n -1 ( 10
)
α j,i + γ l,j -1 ≤ δ l,i,j ≤ (α j,i + γ l,j )/2, ∀l = 1, . . . , n -1, ∀i = l + 1, . . . , n, ∀j = l, . . . , i -1 (11) n k=l+1 γ l,k ≤ (n -l)γ l,l , ∀l = 1, . . . , n -1 (12)
The objective function (1) minimizes the number of used bins. Equations [START_REF] Berkey | Two-dimensional finite bin packing algorithms[END_REF] state that each item has to be packed once. In (3) it is ensured that items are only assigned to a used stack j, i.e., the stack contains item j. The fact that the items packed into the same stack must have identical width is guaranteed by [START_REF] Lodi | Models and bounds for two-dimensional level packing problems[END_REF]. In an implementation, it is not necessary to use all of the α j,i , j = i, but only those for which w i = w j and h i + h j ≤ H. However, we keep them in our model for the sake of clarity. Each used stack is packed exactly once according to [START_REF] Lodi | Recent advances on two-dimensional bin packing problems[END_REF]. Constraints ( 6) and (7) ensure for each stripe k and each contained stack j that the stack's height does not exceed the stripe's height (= stack k's height). Constraints [START_REF] Pisinger | Using decomposition techniques and constraint programming for solving the two-dimensional bin packing problem[END_REF] guarantee that width W is not exceeded and no stacks are packed into unused stripes (β k,k = 0). Equations ( 9) ensure each used stripe being packed into a bin. In [START_REF] Vanderbeck | A nested decomposition approach to a 3-stage 2-dimensional cutting stock problem[END_REF] it is guaranteed that height H is not exceeded and no stripes are packed into unused bins (γ l,l = 0). Constraints (11) force variables δ l,i,j to be set to one iff rectangle i appears in stack j, stripe j is used, and stripe j appears in bin l. In (12) it is ensured that no stripes are packed into an unused bin.

Column Generation Formulation

A general introduction to integer linear programming, B&P, and column generation can be found in Wolsey [START_REF] Wolsey | Integer Programming[END_REF]. Our column generation formulation for 3-stage 2BP is based on the set covering model from [START_REF] Pisinger | Using decomposition techniques and constraint programming for solving the two-dimensional bin packing problem[END_REF] and the ILP from the last section. Let P be the set of all feasible 3-stage packings of a single bin. The 0/1variable x p indicates whether packing p ∈ P appears in the solution. For every rectangle i = 1, . . . , n and every packing p ∈ P, let A p i = 1 iff packing p contains rectangle i; otherwise A p i = 0. The 3-stage 2BP problem can now be formulated as:

minimize p∈P x p (13) subject to p∈P x p A p i ≥ 1 ∀i = 1, . . . , n, (14)
x p ∈ {0, 1} ∀p ∈ P. ( 15
)
In general, P is too huge for explicitly considering all variables x p , p ∈ P . Fortunately, we can use delayed column generation to solve the Linear Programming (LP) relaxation of the problem without explicitly considering the majority of the variables [START_REF] Wolsey | Integer Programming[END_REF]. We start with a small set of initial patterns P ⊂ P taken from an initial feasible solution, and solve the LP relaxation of the problem restricted to P . Based on the obtained solution, we search for a new variable/pattern whose inclusion in the restricted problem might improve the result. The extended LP is resolved and the whole process repeated until no further improvements are possible.

In addition, every M -th iteration the restricted problem is solved to integrality by using branch and cut, possibly providing a new incumbent solution (M = 100 turned out to be a reasonable choice).

The reduced costs of a packing p ∈ P are

c π p = 1 - n i=1 A p i π i , ( 16
)
where π i are the dual variables of the restricted LP-relaxed problem. Only variables with negative reduced costs can improve the current solution of the master problem leading us to the challenge of finding such a variable/pattern.

The Pricing Problem consists of finding a packing p with the smallest reduced costs c π p . It is a 3-Stage 2D Knapsack Packing (2DKP) problem with respect to the profits π i . It can be modeled as follows:

maximize n i=1 π i i j=1 α j,i (17) subject to i j=1 α j,i ≤ 1, ∀i = 1, . . . , n (18) n i=1 h i i j=1 δ i,j ≤ H (19) α j,i + β j,j -1 ≤ δ i,j ≤ α j,i + β j,j 2 ∀i = 1, . . . , n, ∀j = 1, . . . , i
(20) and the constraints: (3), ( 4), ( 5), ( 6), [START_REF] Martello | Exact solutions of the two-dimensional finite bin packing problem[END_REF], and (8).

Variables α j,i and β k,j have the same meaning as in the ILP of Sec. 2. The 0/1variables δ i,j are set to one iff item i contributes to the total height of all used stripes, i.e., iff item i appears in stack j, and stack j appears in stripe j.

Branching If no further variables with negative reduced costs can be determined, and the difference between the solution value of the LP-relaxed restricted problem and the value of the so-far best integer solution is greater than or equal to one, branching becomes necessary. We use a branching rule similar to the one described in [START_REF] Pisinger | Using decomposition techniques and constraint programming for solving the two-dimensional bin packing problem[END_REF]. The solution space is divided into two parts, where two different items i 1 and i 2 have not to be or have to be in the same bin. We always choose the two highest possible items from a pattern/variable with an LP solution value closest to 0.5.

The first branch corresponds to adding the constraint

p∈P x p A p i 1 A p i 2 = 0, (21)
the second branch corresponds to adding the two constraints

p∈P x p A p i1 (1 -A p i2 ) = 0 and p∈P x p (1 -A p i1 )A p i2 = 0. ( 22
)
In the actual implementation we do not explicitly add the constraints (21) and ( 22), but the variables violating them are fixed to zero.

The following constraints have to be added to the pricing problem in order to guarantee that patterns violating the branching constraints cannot be generated. In the first branch

i 1 j=1 α j,i 1 + i 2 j=1 α j,i 2 ≤ 1, (23)
and in the second branch

i1 j=1 α j,i 1 = i2 j=1 α j,i 2 . ( 24
)
In the sequel, we call i 1 and i 2 conflicting if constraint ( 23) is active and say that i 1 induces i 2 and vice-versa if equation ( 24) is active.

Initial feasible solution

In order to initialize the column generation algorithm, a feasible solution is needed. The packing patterns of its bins are used as initial P . This solution is generated by the order-based Finite First Fit heuristic from [START_REF] Puchinger | Solving a real-world glass cutting problem[END_REF] without considering the additional application-specific constraints described there. The heuristic is called 20n times for different item orders, and the best obtained solution is used. The first five item orders are determined by sorting the items according to decreasing height, width, area, 2h i + w i , and h i + 2w i ; all further orders are random permutations.

Solving the Pricing Problem

The pricing problem is solved by a greedy First Fit heuristic respecting the Branching Constraints (FFBC); see Fig. 2.

Similarly to the initialization heuristic, FFBC considers the items in a given order. One item after the other is packed into the first stack it fits. If the item does not fit into any existing stack, a new stack is created in the first stripe it fits. If no such stripe exists and there is enough space left in the bin, a new stack is created and packed into a new stripe. Otherwise, the algorithm proceeds with the next item. If the addition of an item to a stack would increase the corresponding stripe's height, we check if enough vertical space is left in the bin and actually add the item with a probability of 0.5.

The constraints resulting from branching are handled as follows. If an item is considered for packing, we first check whether there are any conflicts with already packed items in the bin (checkNoConflicts in Fig. 2); if there are, the item is skipped. Otherwise we recursively check if other items are induced, and if any of those stay in conflict with any other induced or already packed item (recursiveCheck in Fig. 2). If such a conflict occurs, none of these items can be packed. Otherwise, we immediately try to also pack all the induced items. If this is impossible we skip the whole chain of items.

FFBC is iteratively applied to up to 100 different item orders until a solution with negative reduced costs is found. The first five orders are determined by sorting the items according to decreasing π i , π i h i •w i , π i h i +w i , π i h i , and π i w i ; all further orders are random permutations.

If the heuristic does not find a packing pattern with negative reduced costs, the general purpose ILP-solver CPLEX is finally called on the pricing problem in order to perform an exact search, eventually proving that no such patterns exist anymore.

An Evolutionary Algorithm for the Pricing Problem

Since the 2DKP problem is strongly NP-hard, calling the ILP-solver may be very time-consuming. A more sophisticated metaheuristic, performed when FFBC did not find a variable with negative reduced costs and before solving the problem in an exact way, could lead to a faster overall column generation since significantly fewer calls of the ILP-solver may be needed. We decided to apply an Evolutionary Algorithm (EA) operating directly on stripes, stacks, and items.

Structure of the EA We use a standard steady-state algorithm with binary tournament selection and duplicate elimination. In each iteration, one new candidate solution is created by always applying recombination, and applying mutation with a certain probability. The new solution replaces the worst solution in the population if it is not identical to an already existing solution.

Representation and initialization

The chosen representation is direct: Each chromosome represents a bin as a set of stripes, each stripe as a set of stacks, and each stack as a set of item references. Using such a hierarchy of sets makes it easy to ignore the order of items, stacks, and stripes and to avoid symmetries.

Initial solutions are created via the FFBC heuristic using randomly generated item orders. These orders are created in a biased way by assigning each item i a random value r i ∈ [0, 1) and sorting the items according to decreasing r i π i .

Recombination This operator first assigns a random value r s ∈ [0, 1) to each stripe s in the two parent solutions. All these stripes are then sorted according to decreasing r s p s , with p s being the sum of the π i of the items contained in stripe s. The stripes are then considered in this order and packed into the offspring's bin when they fit into it (i.e., their height is smaller than the remaining unused height of the bin). Identical stripes of both parents appear twice in the ordered list, but they are considered at their first appearance only.

When all stripes have been processed, repairing is usually necessary in order to guarantee feasibility. First, the bin is traversed in order to check if items appear twice, the first of these items is deleted. Then, the branching constraints are considered: Items conflicting with others are removed. Afterwards, we try to pack induced items; if this is not possible the corresponding original items are also removed from the bin. Finally, FFBC is applied to the remaining items, possibly improving the solution.

Mutation The mutation operator removes a randomly chosen item i from the bin. If the branching constraints induce other items for i, they are also deleted. Finally, FFBC is applied to the remaining items for local improvement.

Experimental Results

We performed experiments on the benchmark instances from Berkey and Wang [START_REF] Berkey | Two-dimensional finite bin packing algorithms[END_REF] (classes 1 to 6) and Martello and Vigo [START_REF] Martello | Exact solutions of the two-dimensional finite bin packing problem[END_REF] (classes 7 to 10).

We compare CPLEX 8.1 directly applied to the ILP model ( 1) to ( 12) and the two variants of the B&P approach with and without the EA for solving the pricing problem. The B&P algorithm was implemented using the opensource framework COIN/Bcp (version 2004/04) [START_REF] Lougee-Heimer | The Common Optimization INterface for Operations Research: Promoting open-source software in the operations research community[END_REF], the LPs were solved using COIN/Clp. The computational experiments were performed on a Pentium 4 PC with 2.8 GHz. The EA's population size was 100, the mutation was performed with probability 0.75, and the EA terminated when either 1 000 iterations were performed without an improvement of the best solution or after a total of 100 000 iterations. Each of the experiments had a time limit of 1 000 seconds, which was occasionally exceeded because CPLEX is given the same time limit of 1 000 seconds.

Table 1 shows results obtained for the 10 problem classes; in each class there are 50 instances divided into 5 subclasses with n = 20, . . . , 100 items. For each of the considered algorithms (CPLEX, B&P, B&P with EA), average objective values z of finally best integer solutions, numbers of instances solved to provable optimality Opt (out of 10), and average times t in seconds are given. The last rows show totals and averages over all instances.

When CPLEX is directly applied to the ILP model, 335 out of 500 instances could be solved to provable optimality. This is not bad, but substantially less than B&P's 402 instances and in particular the 409 completely solved instances of the EA-enhanced B&P. The differences in the objective values of the finally best integer solutions found by the three algorithms for instances that could not be solved to optimality are in general relatively small. Nevertheless, B&P's solution values are in several cases significantly better than those of CPLEX, and the EAenhanced B&P performs best on average. The two variants of B&P with and without the EA exhibit approximately the same total running times. Applying CPLEX directly was significantly slower in most cases. Thus, the application of the EA within the B&P framework is worth the additional effort.

Conclusions and Future Work

For 3-stage 2BP, we presented a compact ILP model having only O(n 3 ) variables.

In practice, however, the proposed column generation approach having a number of potential variables that grows exponentially with n turns out to be more efficient. Using the described EA as an additional strategy for solving the pricing problem pays off in terms of a higher capability of solving instances to provable optimality, but also slightly better average solution values. The combination of B&P and an EA in this form is also highly promising for other combinatorial optimization problems. Research on more sophisticated interaction and a parallel execution of these algorithms will be done next.

Fig. 1 .

 1 Fig. 1. A three-stage cutting pattern for one bin in normal form.

Fig. 2 .

 2 Fig. 2. First fit heuristic respecting the branching constraints.

Table 1 .

 1 Experimental results of the presented algorithms.

				CPLEX			B&P		B&P with EA
	Class	n	z	Opt	t [s]	z	Opt	t [s]	z	Opt	t [s]
		20	7.2 10.0	0.0	7.2 10.0	0.3	7.2 10.0	6.3
		40	13.6	6.0	404.8 13.6	8.0	206.1 13.6	8.0	204.0
	1	60	20.2	3.0	748.9 20.2	8.0	256.9 20.1	8.0	221.6
		80	27.8	0.0 1000.6 27.6	9.0	182.1 27.6	9.0	183.9
		100 32.7	0.0 1001.3 32.3	4.0	845.6 32.0	6.0	590.4
		20	1.0 10.0	0.0	1.0 10.0	0.1	1.0 10.0	0.1
		40	2.0	9.0	100.2	2.0	9.0	113.0	2.0	9.0	118.3
	2	60	2.8	7.0	300.8	2.8	7.0	317.9	2.8	7.0	411.2
		80	3.4	7.0	302.4	3.4	7.0	377.9	3.4	7.0	410.1
		100	4.1	8.0	206.7	4.1	8.0	344.0	4.1	8.0	220.8
			5.4 10.0	0.0	5.4 10.0	0.2	5.4 10.0	0.3
		40	9.7	8.0	307.2	9.8	9.0	124.3	9.7 10.0	7.0
	3	60	14.2	5.0	704.3 14.2	7.0	386.9 14.1	9.0	184.4
		80	20.3	0.0 1000.4 19.5	7.0	423.9 19.3	8.0	280.2
		100 23.9	0.0 1000.8 23.2	2.0	950.1 22.9	3.0	946.9
		20	1.0 10.0	0.0	1.0 10.0	0.1	1.0 10.0	0.1
		40	2.0	9.0	100.1	2.0	9.0	100.7	2.0	9.0	166.6
	4	60	2.6	7.0	353.7	2.7	6.0	413.3	2.7	6.0	736.5
		80	3.3	7.0	300.9	3.3	7.0	402.3	3.3	7.0	303.3
		100	4.0	7.0	302.0	4.0	7.0	378.2	4.0	7.0	306.9
		20	6.6 10.0	0.0	6.6 10.0	0.2	6.6 10.0	0.8
		40	12.3 10.0	24.2 12.3 10.0	23.0 12.3 10.0	3.1
	5	60	18.3 10.0	10.3 18.3 10.0	19.0 18.3 10.0	89.0
		80	25.0	5.0	530.1 24.8	9.0	199.4 24.8	9.0	213.7
		100 29.4	1.0	901.5 28.9	5.0	621.6 28.9	7.0	536.3
		20	1.0 10.0	0.0	1.0 10.0	0.1	1.0 10.0	0.1
		40	1.9 10.0	14.2	1.9	6.0	401.7	1.9	6.0	405.7
	6	60	2.3	8.0	200.2	2.3	8.0	202.2	2.3	8.0	201.3
		80	3.0 10.0	0.4	3.0 10.0	3.1	3.0 10.0	3.1
		100	3.6	6.0	400.8	3.6	6.0	405.5	3.6	6.0	405.6
		20	5.7 10.0	0.0	5.7 10.0	0.3	5.7 10.0	0.7
		40	11.5	6.0	566.3 11.5 10.0	4.9 11.5 10.0	39.7
	7	60	16.2	0.0 1000.2 16.2	9.0	128.5 16.1 10.0	24.7
		80	23.5	0.0 1000.4 23.2 10.0	60.2 23.3	9.0	157.3
		100 28.0	0.0 1000.7 27.1 10.0	292.1 27.1 10.0	269.6
		20	6.1 10.0	0.0	6.1 10.0	0.8	6.1 10.0	0.9
		40	11.4 10.0	0.6 11.5	9.0	133.7 11.4 10.0	120.1
	8	60	16.4 10.0	8.6 16.5	9.0	116.0 16.5	9.0	118.6
		80	22.6	8.0	288.6 22.7	9.0	177.8 22.9	7.0	346.5
		100 28.2	8.0	410.0 28.3	7.0	506.4 28.4	6.0	508.9
		20	14.3 10.0	0.0 14.3 10.0	0.1 14.3 10.0	0.2
		40	27.8 10.0	0.0 27.8 10.0	0.4 27.8 10.0	0.6
	9	60	43.7 10.0	0.1 43.7 10.0	1.4 43.7 10.0	1.5
		80	57.7 10.0	0.2 57.7 10.0	3.6 57.7 10.0	3.8
		100 69.5	8.0	200.3 69.5 10.0	7.9 69.5 10.0	8.5
		20	4.5 10.0	0.0	4.5 10.0	0.7	4.5 10.0	0.5
		40	7.7	9.0	185.8	7.7	9.0	166.1	7.8	8.0	217.4
	10	60	10.7	3.0	822.1 10.7	2.0	857.2 10.5	3.0	825.1
		80	14.0	0.0 1000.4 13.9	0.0 1048.8 13.6	0.0 1114.5
		100 16.9	0.0 1000.7 16.9	0.0 1048.5 16.6	0.0 1151.4
	Total		741.0 335.0 17702.0 737.5 402.0 12254.9 735.9 409.0 12068.2
	Average		14.82 6.70 354.04 14.75 8.04 245.10 14.72 8.18 241.36

This work is supported by the Austrian Science Fund (FWF) under grant P16263-N04.


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

Günther R Raidl

Combining Metaheuristics and Exact Algorithms in Combinatorial Optimization: A Survey and Classification

In this survey we discuss different state-of-the-art approaches of combining exact algorithms and metaheuristics to solve combinatorial optimization problems. Some of these hybrids mainly aim at providing optimal solutions in shorter time, while others primarily focus on getting better heuristic solutions. The two main categories in which we divide the approaches are collaborative versus integrative combinations. We further classify the different techniques in a hierarchical way. Altogether, the surveyed work on combinations of exact algorithms and metaheuristics documents the usefulness and strong potential of this research direction.

Introduction

Hard combinatorial optimization problems (COPs) appear in a multitude of realworld applications, such as routing, assignment, scheduling, cutting and packing, network design, protein alignment, and many other fields of utmost economic, industrial and scientific importance. The available techniques for COPs can roughly be classified into two main categories: exact and heuristic methods. Exact algorithms are guaranteed to find an optimal solution and to prove its optimality for every instance of a COP. The run-time, however, often increases dramatically with the instance size, and often only small or moderately-sized instances can be practically solved to provable optimality. In this case, the only possibility for larger instances is to trade optimality for run-time, yielding heuristic algorithms. In other words, the guarantee of finding optimal solutions is sacrificed for the sake of getting good solutions in a limited time.

Two independent heterogeneous streams, coming from very different scientific communities, had significant success in solving COPs:

-Integer Programming (IP) as an exact approach, coming from the operations research community and based on the concepts of linear programming [START_REF] Dantzig | Linear Programming and Extensions[END_REF]. -Local search with various extensions and independently developed variants, in the following called metaheuristics, as a heuristic approach.

Among the exact methods are branch-and-bound (B&B), dynamic programming, Lagrangian relaxation based methods, and linear and integer programming based methods, such as branch-and-cut, branch-and-price, and branch-and-cut-andprice [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF].

Metaheuristics include, among others, simulated annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], tabu search [START_REF] Glover | Tabu Search[END_REF], iterated local search [START_REF] Lourenço | Iterated local search[END_REF], variable neighborhood search [START_REF] Hansen | An introduction to variable neighborhood search[END_REF], and various populationbased models such as evolutionary algorithms [START_REF] Bäck | Handbook of Evolutionary Computation[END_REF], scatter search [START_REF] Glover | Fundamentals of scatter search and path relinking[END_REF], memetic algorithms [START_REF] Moscato | A gentle introduction to memetic algorithms[END_REF], and various estimation of distribution algorithms [START_REF] Larrañaga | Estimation of Distribution Algorithms[END_REF]. Recently there have been very different attempts to combine ideas and methods from these two scientific streams. Dumitrescu and Stützle [START_REF] Dumitrescu | Combinations of local search and exact algorithms[END_REF] describe existing combinations, focusing on local search approaches that are strengthened by the use of exact algorithms. In their survey they concentrate on integration and exclude obvious combinations such as preprocessing.

Here, we present a more general classification of existing approaches combining exact and metaheuristic algorithms for combinatorial optimization. We distinguish the following two main categories:

-Collaborative Combinations: By collaboration we mean that the algorithms exchange information, but are not part of each other. Exact and heuristic algorithms may be executed sequentially, intertwined or in parallel. -Integrative Combinations: By integration we mean that one technique is a subordinate embedded component of another technique. Thus, there is a distinguished master algorithm, which can be either an exact or a metaheuristic algorithm, and at least one integrated slave.

In the following sections this classification is further refined and examples from the literature are presented, reflecting the current state-of-the-art. Figure 1 gives an overview of this classification.

Collaborative Combinations

The different algorithms and approaches described in this section have in common that they are top-level combinations of metaheuristics and exact techniques; no algorithm is contained in another. We further distinguish whether the algorithms are executed sequentially or in an intertwined or even parallel way.

Sequential Execution

Either the exact method is executed as a kind of preprocessing before the metaheuristic, or vice-versa. Sometimes, it is difficult to say if the first technique is used as initialization of the second, or if the second is a postprocessing of the solution(s) generated by the first.

Clements et al. [START_REF] Clements | Heuristic optimization: A hybrid AI/OR approach[END_REF] propose a column generation approach in order to solve a production-line scheduling problem. Each feasible solution of the problem consists of a line-schedule for each production line. First, the squeaky wheel optimization (SWO) heuristic is used to generate feasible solutions to the problem. SWO is a heuristic using a greedy algorithm to construct a solution, which is then analyzed in order to find the problematic elements. Higher priorities, such that these elements are considered earlier by the greedy algorithm, are assigned to them, and the process restarts until a termination condition is reached. SWO is called several times in a randomized way in order to generate a set of diverse solutions. In the second phase, the line-schedules contained in these solutions are used as columns of a set partitioning formulation for the problem, which is solved using MINTO. This process always provides a solution which is at least as good as, but usually better than the best solution devised by SWO. Reported results indicate that SWO performs better than a tabu-search algorithm.

Applegate et al. [START_REF] Applegate | On the solution of the traveling salesman problem[END_REF] propose an approach for finding near-optimal solutions to the traveling salesman problem. They derive a set of diverse solutions by multiple runs of an iterated local search algorithm. The edge-sets of these solutions are merged and the traveling salesman problem is finally solved to optimality on this strongly restricted graph. In this way a solution is achieved that is typically superior to the best solution of the iterated local search.

Klau et al. [START_REF] Klau | Combining a memetic algorithm with integer programming to solve the prize-collecting Steiner tree problem[END_REF] follow a similar idea and combine a memetic algorithm with integer programming to heuristically solve the prize-collecting Steiner tree problem. The proposed algorithmic framework consists of three parts: extensive preprocessing, a memetic algorithm, and an exact branch-and-cut algorithm applied as post-optimization procedure to the merged final solutions of the memetic algorithm.

Plateau et al. [START_REF] Plateau | A hybrid search combining interior point methods and metaheuristics for 0-1 programming[END_REF] combine interior point methods and metaheuristics for solving the multiconstrained knapsack problem. The first part is an interior point method with early termination. By rounding and applying several different ascent heuristics, a population of different feasible candidate solutions is generated. This set of solutions is then used as initial population for a path-relinking (scatter search) algorithm. Extensive computational experiments are performed on standard multiconstrained knapsack benchmark instances. Obtained results show that the presented combination is a promising research direction.

Sometimes, a relaxation of the original problem is solved to optimality and the obtained solution is repaired to act as a promising starting point for a subsequent metaheuristic. Often, the linear programming (LP) relaxation is used for this purpose, and only a simple rounding scheme is needed. For example, Feltl and Raidl [START_REF] Raidl | An improved hybrid genetic algorithm for the generalized assignment problem[END_REF] solve the generalized assignment problem using a hybrid genetic algorithm (GA). The LP-relaxation of the problem is solved using CPLEX and its solution is used by a randomized rounding procedure to create a population of promising integral solutions. These solutions are, however, often infeasible; therefore, randomized repair and improvement operators are additionally applied, yielding an even more meaningful initial population for the GA. Reported computational experiments suggest that this type of LP-based initialization is effective.

Vasquez and Hao [START_REF] Vasquez | A hybrid approach for the 0-1 multidimensional knapsack problem[END_REF] heuristically solve the multiconstrained knapsack problem by reducing and partitioning the search space via additional constraints that fix the total number of items to be packed. The bounds for these constraints are calculated by solving a modified LP-relaxation of the multiconstrained knapsack problem. For each remaining part of the search space, parallel tabu-search is finally performed starting with a solution derived from the LP-relaxation of the partial problem. This hybrid algorithm yields excellent results also for large benchmark instances with up to 2 500 items and 100 constraints.

Lin et al. [START_REF] -Z. Lin | A hybrid genetic/optimization algorithm for finite horizon partially observed markov decision processes[END_REF] describe an exact algorithm for generating the minimal set of affine functions that describes the value function of the finite horizon partially observed Markov decision process. In the first step a GA is used to generate a set Γ of witness points, which is as large as possible. In the second step a component-wise domination procedure is performed in order to eliminate redundant points in Γ . The set generated so far does, in general, not fully describe the value function. Therefore, a Mixed Integer Program (MIP) is solved to generate the missing points in the final third step of the algorithm. Reported results indicate that this approach requires less time than some other numerical procedures.

Another kind of sequential combination of B&B and a GA is described by Nagar et al. [START_REF] Nagar | A meta-heuristic algorithm for a bi-criteria scheduling problem[END_REF] for a two-machine flowshop scheduling problem in which solution candidates are represented as permutations of jobs. Prior to running the GA B&B is executed down to a predetermined depth k and suitable bounds are calculated and recorded at each node of the explicitly stored B&B tree. During the execution of the GA the partial solutions up to position k are mapped onto the correct tree node. If the bounds indicate that no path below this node can lead to an optimal solution, the permutation is subjected to a mutation operator that has been specifically designed to change the early part of the permutation in a favorable way.

Tamura et al. [START_REF] Tamura | An approximate solution method for combinatorial optimisation[END_REF] tackle a job-shop scheduling problem and start from its IP formulation. For each variable, they take the range of possible values and partition it into a set of subranges, which are then indexed. The chromosomes of the GA are defined so that each position represents a variable, and its value corresponds to the index of one of the subranges. The fitness of a chromosome is calculated using Lagrangian relaxation to obtain a bound on the optimal solution subject to the constraints that the values of the variables fall within the correct ranges. When the GA terminates, an exhaustive search of the region identified as the most promising is carried out to produce the final solution.

Parallel or Intertwined Execution

Instead of a strictly sequential batch approach, exact and heuristic algorithms may also be executed in a parallel or intertwined way. Such peer-to-peer combinations of exact/heuristic techniques are less frequent. An interesting framework for this purpose was proposed by Talukdar et al. [START_REF] Talukdar | Asynchronous teams: Cooperation schemes for autonomous agents[END_REF][START_REF] Talukdar | Asynchronous teams[END_REF] with the so-called asynchronous teams (A-Teams). An A-Team is a problem solving architecture consisting of a collection of agents and memories connected into a strongly cyclic directed network. Each of these agents is an optimization algorithm and can work on the target problem, on a relaxation-i.e., a superclass-of it, or on a subclass of the problem. The basic idea of A-Teams is having these agents work asynchronously and autonomously on a set of shared memories. These shared memories consist of trial solutions for some problem (the target problem, a superclass, or a subclass as mentioned before), and the action of an agent consists of modifying the memory by adding a solution, deleting a solution, or altering a solution. A-Teams have been successfully utilized in a variety of combinatorial optimization problems, see e.g. [START_REF] Chen | Job-shop-scheduling by a team of asynchronous agents[END_REF][START_REF] Talukdar | Asynchronous teams[END_REF].

Denzinger and Offerman [START_REF] Denzinger | On cooperation between evolutionary algorithms and other search paradigms[END_REF] present a similar multi-agent based approach for achieving cooperation between search-systems with different search paradigms. The TECHS (TEams for Cooperative Heterogenous Search) approach consists of teams of one or more agents using the same search paradigm. The communication between the agents is controlled by so-called send-and receive-referees, in order to filter the exchanged data. Each agent is in a cycle between searching and processing received information. In order to demonstrate the usefulness of TECHS, a GA and a B&B based system for job-shop scheduling is described. The GA and B&B agents exchange only positive information (solutions), whereas the B&B agents can also exchange negative information (closed subtrees). Computational experiments show that the cooperation results in finding better solutions given a fixed time-limit and in finding solutions comparable to the ones of the best individual system alone in less time.

Integrative Combinations

In this section we discuss approaches of combining exact algorithms and metaheuristics in an integrative way such that one technique is a subordinate embedded component of another technique.

Incorporating Exact Algorithms in Metaheuristics

We start by considering techniques where exact algorithms are incorporated into metaheuristics.

Exactly Solving Relaxed Problems The usefulness of solutions to relaxations of an original problem has already been mentioned in Section 2.1. Besides exploiting them to derive promising initial solutions for a subsequent algorithm, they can be of great benefit for heuristically guiding neighborhood search, recombination, mutation, repair and/or local improvement. Examples where the solution of the LPrelaxation and its dual were exploited in such ways are the hybrid genetic algorithms for the multiconstrained knapsack problem from Chu and Beasley [START_REF] Chu | A genetic algorithm for the multidimensional knapsack problem[END_REF] and Raidl [START_REF] Raidl | An improved genetic algorithm for the multiconstrained 0-1 knapsack problem[END_REF].

Exactly Searching Large Neighborhoods A common approach is to search neighborhoods in local search based metaheuristics by means of exact algorithms. If the neighborhoods are chosen appropriately, they can be relatively large and nevertheless an efficient search for the best neighbor is still reasonable. Such techniques are known as Very Large-Scale Neighborhood (VLSN) search [START_REF] Ahuja | A survey of very large-scale neighborhood search techniques[END_REF].

Burke et al. [START_REF] Burke | Effective local and guided variable neighborhood search methods for the asymmetric travelling salesman problem[END_REF] present an effective local and variable neighborhood search heuristic for the asymmetric traveling salesman problem in which they have embedded an exact algorithm in the local search part, called HyperOpt, in order to exhaustively search relatively large promising regions of the solution space. Moreover, they propose a hybrid of HyperOpt and 3-opt which allows to benefit from the advantages of both approaches and gain better tours overall. Using this hybrid within the variable neighborhood search metaheuristic framework also allows to overcome local optima and to create tours of high quality.

Dynasearch [START_REF] Congram | Polynomially Searchable Exponential Neighbourhoods for Sequencing Problems in Combinatorial Optimisation[END_REF] is another example where exponentially large neighborhoods are explored. The neighborhood where the search is performed consists of all possible combinations of mutually independent simple search steps and one Dynasearch move consists of a set of independent moves that are executed in parallel in a single local search iteration. Independence in the context of Dynasearch means that the individual moves do not interfere with each other; in this case, dynamic programming can be used to find the best combination of independent moves. Dynasearch is restricted to problems where the single search steps are independent, and it has so far only been applied to problems, where solutions are represented as permutations.

For the class of partitioning problems, Thompson et al. [START_REF] Thompson | The theory of cycle transfers[END_REF][START_REF] Thompson | Cycle transfer algorithm for multivehicle routing and scheduling problems[END_REF] defined the concept of a cyclic exchange neighborhood, which is the transfer of single elements between several subsets in a cyclic manner; for example, a 2-exchange move can be seen as a cyclic exchange of length two. Thompson et al. showed that for any current solution to a partitioning problem a new, edge-weighted graph can be constructed, where the set of nodes is split into subsets according to a partition induced by the current solution of the partitioning problem. A cyclic exchange for the original problem corresponds to a cycle in this new graph that uses at most one node of each subset. Exact and heuristic methods that solve the problem of finding the most negative-cost subset-disjoint cycle (which corresponds to the best improving neighbor of the current solution) have been developed. Puchinger et al. [START_REF] Puchinger | Solving a real-world glass cutting problem[END_REF] describe a combined GA/B&B approach for solving a realworld glass cutting problem. The GA uses an order-based representation, which is decoded using a greedy heuristic. The B&B algorithm is applied with a certain probability enhancing the decoding phase by generating locally optimal subpatterns. Reported results indicate that the approach of occasionally solving subpatterns to optimality may increase the overall solution quality.

The work of Klau et al. [START_REF] Klau | Combining a memetic algorithm with integer programming to solve the prize-collecting Steiner tree problem[END_REF] has already been mentioned in Section 2.1 in the context of collaborative sequential combinations. When looking at the memetic algorithm we encounter another kind of exact/heuristic algorithm combination. An exact subroutine for the price-collecting Steiner tree problem on trees is used to locally improve candidate solutions.

Merging Solutions Subspaces defined by the merged attributes of two or more solutions can, like the neighborhoods of single solutions, also be searched by exact techniques. The algorithms by Clements et al. [START_REF] Clements | Heuristic optimization: A hybrid AI/OR approach[END_REF], Applegate et al. [START_REF] Applegate | On the solution of the traveling salesman problem[END_REF], and Klau et al. [START_REF] Klau | Combining a memetic algorithm with integer programming to solve the prize-collecting Steiner tree problem[END_REF], which were already discussed in Section 2.1, also follow this idea, but are of sequential collaborative nature. Here, we consider approaches where merging is iteratively applied within a metaheuristic.

Cotta and Troya [START_REF] Cotta | Embedding branch and bound within evolutionary algorithms[END_REF] present a framework for hybridizing B&B with evolutionary algorithms. B&B is used as an operator embedded in the evolutionary algorithm. The authors recall the necessary theoretical concepts on forma analysis (formae are generalized schemata), such as the dynastic potential of two chromosomes x and y, which is the set of individuals that only carry information contained in x and y. Based on these concepts the idea of dynastically optimal recombination is developed. This results in an operator exploring the potential of the recombined solutions using B&B, providing the best possible combination of the ancestors' features that can be attained without introducing implicit mutation. Extensive computational experiments on different benchmark sets comparing different crossover operators with the new hybrid one show the usefulness of the presented approach.

Marino et al. [START_REF] Marino | Improving graph colouring with linear programming and genetic algorithms[END_REF] present an approach where a GA is combined with an exact method for the Linear Assignment Problem (LAP) to solve the graph coloring problem. The LAP algorithm is incorporated into the crossover operator and generates the optimal permutation of colors within a cluster of nodes, hereby preventing the offspring to be less fit than its parents. The algorithm does not outperform other approaches, but provides comparable results. The main conclusion is that solving the LAP in the crossover operator strongly improves the performance of the GA compared to the GA using crossover without LAP.

Exact Algorithms as Decoders In evolutionary algorithms, candidate solutions are sometimes only incompletely represented in the chromosome, and an exact algorithm is used as decoder for determining the missing parts in an optimal way.

Staggemeier et al. [START_REF] Staggemeier | A hybrid genetic algorithm to solve a lot-sizing and scheduling problem[END_REF], for example, present a hybrid genetic algorithm to solve a lot-sizing and scheduling problem minimizing inventory and backlog costs of multiple products on parallel machines. Solutions are represented as product subsets for each machine at each period. Corresponding optimal lot sizes are determined when the solution is decoded by solving a linear program. The approach outperforms a MIP formulation of the problem solved using CPLEX.

Incorporating Metaheuristics in Exact Algorithms

We now turn to techniques where metaheuristics are embedded within exact algorithms.

Metaheuristics for Obtaining Incumbent Solutions and Bounds

In general, heuristics and metaheuristics are often used to determine bounds and incumbent solutions in B&B approaches. For example, Woodruff [START_REF] Woodruff | A chunking based selection strategy for integrating meta-heuristics with branch and bound[END_REF] describes a chunking-based selection strategy to decide at each node of the B&B tree whether or not reactive tabu search is called in order to eventually find a better incumbent solution. The chunking-based strategy measures a distance between the current node and nodes already explored by the metaheuristic in order to bias the selection toward distant points. Reported computational results indicate that adding the metaheuristic improves the B&B performance.

Metaheuristics for Column and Cut Generation

In branch-and-cut and branch-and-price algorithms, the dynamic separation of cutting-planes and the pricing of columns, respectively, is sometimes done by means of heuristics including metaheuristics in order to speed up the whole optimization process.

Filho and Lorena [START_REF] Filho | Constructive genetic algorithm and column generation: an application to graph coloring[END_REF] apply a heuristic column generation approach to graph coloring. They describe the principles of their constructive genetic algorithm and give a column generation formulation of the problem. The GA is used to generate the initial columns and to solve the slave problem (the weighted maximum independent set problem) at every iteration. Column generation is performed as long as the GA finds columns with negative reduced costs. The master problem is solved using CPLEX. Some encouraging results are presented.

Puchinger and Raidl [START_REF] Puchinger | An evolutionary algorithm for column generation in integer programming: an effective approach for 2D bin packing[END_REF][START_REF] Puchinger | Models and algorithms for three-stage two-dimensional bin packing[END_REF] propose new integer linear programming formulations for the three-stage two-dimensional bin packing problem. Based on these formulations, a branch-and-price algorithm was developed in which fast column generation is performed by applying a hierarchy of four methods: (a) a greedy heuristic, (b) an evolutionary algorithm, (c) solving a restricted form of the pricing problem using CPLEX, and finally (d) solving the complete pricing problem using CPLEX. Computational experiments on standard benchmark instances document the benefits of the new approach. The combination of all four pricing algorithms in the proposed branch-and-price framework yields the best results in terms of the average objective value, the average run-time, and the number of instances solved to proven optimality.

Metaheuristics for Strategic Guidance of Exact Search French et al. [START_REF] French | Using a hybrid geneticalgorithm/branch and bound approach to solve feasibility and optimization integer programming problems[END_REF] present a GA/B&B hybrid to solve feasibility and optimization IP problems. Their hybrid algorithm combines the generic B&B of the MIP-solver XPRESS-MP with a steady-state GA. It starts by traversing the B&B tree. During this phase, information from nodes is collected in order to suggest chromosomes to be added to the originally randomly initialized GA-population. When a certain criterion is fulfilled, the GA is started using the augmented initial population. When the GA terminates, its fittest solution is passed back and grafted onto the B&B tree. Full control is given back to the B&B-engine, after the newly added nodes were examined to a certain degree. Reported results on MAX-SAT instances show that this hybrid approach yields better solutions than B&B or the GA alone.

Kotsikas and Fragakis [START_REF] Kostikas | Genetic programming applied to mixed integer programming[END_REF] determine improved node selection strategies within B&B for solving MIPs by using genetic programming (GP). After running B&B for a certain amount of time, information is collected from the B&B tree and used as training set for GP, which is performed to find a node selection strategy more appropriate for the specific problem at hand. The following second B&B phase then uses this new node selection strategy. Reported results show that this approach has potential, but needs to be enhanced in order to be able to compete with today's state-of-the-art node selection strategies.

Applying the Spirit of Metaheuristics Last but not least, there are a few approaches where it is tried to bring the spirit of local search based techniques into B&B. The main idea is to first search some neighborhood of incumbent solutions more intensively before turning to a classical node selection strategy. However, there is no explicit metaheuristic, but B&B itself is used for doing the local search. The metaheuristic may also be seen to be executed in a "virtual" way.

Fischetti and Lodi [START_REF] Fischetti | Local Branching[END_REF] introduced local branching, an exact approach combining the spirit of local search metaheuristics with a generic MIP-solver (CPLEX). They consider general MIPs with 0-1 variables. The idea is to iteratively solve a local subproblem corresponding to a classical k-OPT neighborhood using the MIP-solver. This is achieved by introducing a local branching constraint based on an incumbent solution x, which partitions the search space into the k-OPT neighborhood and the rest: ∆(x, x) ≤ k and ∆(x, x) ≥ k + 1, respectively, with ∆ being the Hamming distance of the 0-1 variables. The first subproblem is solved, and if an improved solution could be found, a new subproblem is devised and solved; this is repeated as long as an improved solution is found. If the process stops, the rest of the problem is solved in a standard way. This basic mechanism is extended by introducing time limits, automatically modifying the neighborhood size k and adding diversification strategies in order to improve the performance. Reported results are promising.

Danna et al. [START_REF] Danna | Exploring relaxation induced neighbourhoods to improve MIP solutions[END_REF] present an approach called Relaxation Induced Neighborhood Search (RINS) in order to explore the neighborhoods of promising MIP solutions more intensively. The main idea is to occasionally devise a sub-MIP at a node of the B&B tree that corresponds to a certain neighborhood of an incumbent solution: First, variables having the same values in the incumbent and in the current solution of the LP-relaxation are fixed. Second, an objective cutoff based on the objective value of the incumbent is set. Third, a sub-MIP is solved on the remaining variables. The time for solving this sub-MIP is limited. If a better incumbent could be found during this process, it is passed to the global MIP-search which is resumed after the sub-MIP termination. CPLEX is used as MIP-solver. The authors experimentally compare RINS to standard CPLEX, local branching, combinations of RINS and local branching, and guided dives. Results indicate that RINS often performs best.

Conclusions

We gave a survey on very different, existing approaches for combining exact algorithms and metaheuristics. The two main categories in which we divided these techniques are collaborative and integrative combinations. Some of the combinations are dedicated to very specific combinatorial optimization problems, whereas others are designed to be more generally useful. Altogether, the existing work documents that both, exact optimization techniques and metaheuristics, have specific advantages which complement each other. Suitable combinations of exact algorithms and metaheuristics can benefit much from synergy and often exhibit significantly higher performance with respect to solution quality and time. Some of the presented techniques are mature, whereas others are still in their infancy and need substantial further research in order to make them fully developed. Future work on such hybrid systems is highly promising.

Fig. 1 .

 1 Fig. 1. Major classification of exact/metaheuristic combinations.

This work is supported by the Austrian Science Fund (FWF) under grant P16263-N04.

http://www.isye.gatech.edu/faculty/Martin Savelsbergh/software http://www.ilog.com

http://www.dashoptimization.com/