Diskussionsbeiträge des Fachbereichs Wirtschaftswissenschaft der FernUniversität in Hagen
Herausgegeben vom Dekan des Fachbereichs
A parallel hybrid local search algorithm for the container loading problem
Daniel Mack, Andreas Bortfeldt und Hermann Gehring Diskussionsbeitrag Nr. 353
Dezember 2003
Daniel Mack, Andreas Bortfeldt and Hermann Gehring
Abstract:
In this contribution, a parallel hybrid local search algorithm for the three-dimensional contai- ner loading problem (CLP) is proposed. First a simulated annealing method for the CLP is developed, which is then combined with an existing tabu search algorithm to form a hybrid metaheuristic. Finally, parallel versions are introduced for these algorithms. The emphasis is on CLP instances with a weakly heterogeneous load. Numerical tests based on the well- known 700 test instances from Bischoff and Ratcliff are performed, and the outcome is com- pared to methods from other authors. The results show a high solution quality obtained with reasonable computing time.
Key words:
Cutting and packing, container loading problem, local search, metaheuristic, simulated annea- ling, tabu search, hybrid algorithm, distributed-parallel search.
Fachbereich Wirtschaftswissenschaft, FernUniversität in Hagen Profilstr. 8, D-58084 Hagen, BRD
Tel.: 02331/987–4433 Fax: 02331/987–4447
E-Mail: andreas.bortfeldt@fernuni-hagen.de
A parallel hybrid local search algorithm for the container loading problem
Daniel Mack, Andreas Bortfeldt and Hermann Gehring
1 Introduction
The efficient use of transport capacities is of high economic relevance in real goods processes.
Suboptimal loading of containers, palettes, etc. leads to high follow-up costs all over the distri- bution system. Other side effects are increased traffic activity and negative consequences on environmental resources. Today, the recent research activities in modern optimization tech- niques and artificial intelligence, as well as the development of more powerful computers, are allowing the implementation of elaborated search algorithms for packing problems. In this con- tribution, we present a parallel hybrid algorithm based on the well known tabu search (TS) and simulated annealing (SA) search strategies to solve the container loading problem, with an em- phasis on weakly heterogeneous rectangular box sets.
The container loading problem can be characterized as follows: let a single rectangular con- tainer and a set of rectangular packing pieces (referred to as boxes) be given. In general, the total volume of these boxes exceeds the container volume. The objective is to determine a fea- sible arrangement of a subset of boxes, which maximizes the stowed box volume and meets given constraints, despite some boxes remaining unloaded. An arrangement is called feasible, if:
- each box is placed completely within the container;
- each box is not overlapping with another box;
- each box is placed parallel to the side walls of the container.
The stock of boxes consists of a relatively large number of individual items, grouped by a rela- tively small number of distinct box types. This kind of box set can be described as weakly he- terogeneous, as opposed to homogeneous box sets (i.e. a single box type is used) and strongly heterogeneous sets with a greater number of box types and only a few items per type.
Additionally, two restrictions have to be respected. In order to guarantee the stability of the load, the so-called stability constraint specifies that both horizontal dimensions of each box are to be supported according to a pre-defined percentage. In any case the centre of gravity of each box must be supported in order to avoid the box tipping over. It is assumed that the centre of gravity and the geometric centre of each packing piece coincide. Furthermore, in order to avoid unstable piles of boxes, the supported part of each box is calculated in relation to the container floor. In figure 1, the centre of gravity of the marked box is supported by underlying boxes, but it is not supported (indirectly) by the container floor. As a consequence, this box is not loaded.
x z centre of gravity z
Figure 1. Stability constraint: avoiding unstable box layers.
A second restriction, the orientation constraint, limits the possible spatial orientation variants for certain box types. Depending on the goods stored within the boxes, it may be specified that one or two side dimensions may not be used as the height.
The problem described here belongs to the wider class of cutting and packing problems (C&P).
These can generally be described as a geometric combinatorial exercise (cf. DYCKHOFF and FINKE 1992):
- small geometric objects must be attributed to bigger objects (combinatorial element of the problem);
- geometric limitations are to be respected (geometric element of the problem).
Since packing problems are known to be NP-hard (cf. SCHEITHAUER 1992), in general they re- main too large and too complex to be solved exactly and a heuristic approach is often the only viable option. Over recent years, many intelligent tentative search-strategies have been intro- duced for three-dimensional packing problems. Tree-search algorithms are proposed, e.g., by MORABITO and ARENALES (1994), ELEY (2002), PISINGER (2002) and HIFI (2002). GEHRING
and BORTFELDT (1997) and BORTFELDT and GEHRING (2001) present genetic algorithms (GA) while tabu search algorithms (TS) are introduced by SIXT (1996) and by BORTFELDT and GE- HRING (1998). A simulated annealing adaptation for the three-dimensional packing problem is given by SIXT (1996).
Parallelization and hybridization of intelligent search-strategies present some opportunities for improving the solution quality and efficiency. Parallelization is especially pro-mising for ex- tremely hard problems with large solution spaces as the considered multi-dimensional C&P problems (cf. BORTFELDT 1995). Only a few parallel approaches exist for three-dimensional packing problems, e.g. a parallel GA from GEHRING and BORTFELDT (2002) and a parallel TS from BORTFELDT, GEHRING and MACK (2002). This parallel TS algorithm (referred to in the following as the “PTS”) serves as a starting point for the contribution at hand.
The PTS is hierarchically structured in three modules:
- the lowest module (named “basic heuristic”) is designed to generate feasible solutions to the problem;
- the middle module uses the principles of tabu search in order to direct the search through the solution space;
- the uppermost module, a parallelization module, supervises the concurrent execution of sev- eral instances of the sequential TS, which can cooperate through the exchange of best solu- tions.
In this contribution, the PTS will be extended to a parallel hybrid system in several steps:
- the simulated annealing strategy is implemented as an alternative to tabu search;
- the SA and TS approaches are combined to form a sequential hybrid system (HYB);
- finally, a parallel SA approach (PSA) and a parallel hybrid method (PHYB) are introduced.
Figure 2 gives an overview of the different algorithms and their interdependencies. As figure 2 also shows, the algorithms, PSA and PHYB, consist of modules that are analogous to the PTS.
uppermost level
middle level
lowest level
parallel tabu search (PTS)
parallel hybrid search (PHYB)
parallel simulated annealing (PSA)
tabu search (TS)
hybrid search (HYB)
simulated annealing (SA)
basic heuristic Figure 2. Proposed algorithms for the CLP and their interdependencies.
The primary goal of the extensions introduced in this paper is to improve the solution quality.
Additionally, the following methodological aspects will be considered:
- General suitability of TS and SA for concrete adaptation to packing problems; since identical problem-specific components, namely the basic heuristic and the neighborhood structures, are used for the TS and the SA approach, a fair comparison between both search techniques is possible;
- synergetic effects of combining the two strategies to form a hybrid approach;
- benefits of parallelizing the different strategies (TS, SA and hybrid) and of an exchange of good solutions between the cooperating instances.
TS and SA have some points in common. Both strategies are distinct, qualified, intelligent lo- cal search techniques. Both have to avoid the problem of a stagnation of the search process in local optima of inferior quality, a typical problem observed in most hill climbing or Greedy algorithms. Common to both local search strategies is also the general way they are imple- mented. On the one hand, some problem-specific decisions regarding the representation and generation of solutions and the neighborhood structures employed are to be made. On the other hand, some generic considerations are made to adapt and configure the general search strategy to the concrete problem. For a comprehensive introduction of the search strategies TS and SA the reader is referred to REEVES (1993).
The paper is structured as follows: in section 2 and 3 the problem-specific and the generic as- pects of the SA approach are described. In section 4 the TS and SA algorithms are combined to form a hybrid search technique. The SA approach and the hybrid algorithm are extended to form parallel algorithms in section 5. In section 6 the results of extensive numerical tests are reported and evaluated. Finally, in section 7, the contribution is summarized.
2 Problem-specific aspects of the simulated annealing approach
In general, problem-specific decisions necessary for the design of a simulated annealing proce- dure involve the following components:
- The generating process, generally a heuristic, to compose feasible solutions;
- a representation of feasible solutions determining the solution space;
- a neighborhood structure representing, for each solution, the set of feasible solutions that can be reached in one single step;
- an objective function to evaluate the quality of a solution in order to compare it to other solu- tions.
In the following these topics are addressed in relation to the problem at hand. In addition, a concept for the diversification of the search is described.
Generation of feasible solutions
The basic heuristic of the original PTS, responsible for determining (single) feasible solutions, will be maintained. It has proven its efficiency for weakly heterogeneous box sets. In the inter- est of readability, the basic heuristic will be briefly described. A more detailed description can be found in BORTFELDT,GEHRING and MACK (2002).
By means of the basic heuristic, the container is loaded in several iterations. Within a single iteration, a so-called packing space is filled with one or more boxes as shown in figure 3. A packing space is an empty rectangular space within the container with defined side dimensions.
Initially the complete interior of the container is used as the packing space.
container
Block 2 packing space
block 1
block 2
Local arrangement relative to a packing space consisting of 2 blocks
Figure 3. Components of the basic heuristic: packing space, local arrangement and blocks.
Only arrangements with a predefined structure are considered. These are called local arrange- ments. A local arrangement consists of one or two so-called blocks placed as neighbors in the packing space. Each block is formed by boxes of the same type, arranged in an identical spatial orientation variant. Thus, the space of a block is completely filled by boxes (cf. figure 3). In the case of two-blocks-arrangements, three variants “beside”, “in front of” and “above” are possi- ble depending on the relative position of the blocks.
For a given packing space, the possible arrangements are generated and evaluated according various criteria. Evaluation includes the residual space obtained by realizing the specific ar- rangement so as to generate new promising packing spaces from the unused space. Subse- quently one specific arrangement is chosen. The unused part of the packing space is divided into new residual packing spaces to be filled later in the iterative process. In each iterative step, the smallest available packing space is examined. The steps of the basic heuristic are summa- rized in figure 4.
(1) Input(container data, box data, parameter, packing sequence Pseq[jsp], jsp = 1,…,maxPseq).
(2) Initialize:
the set of residual boxes BRes:= set of all boxes;
the packing space list PsList:= { Container };
the packing space index ips:= 0;
the stowing list StList:= {}.
(3) Determine the current packing space pcurr as the packing space from PsList with mini- mum volume and delete pcurr from PsList.
(4) Initialize the arrangement list for pcurr: ArrList:= {};
generate and evaluate all local arrangements for pcurr;
insert the local arrangements in descending order with respect to the evaluation into ArrList.
(5) If ArrList is empty, go to step (9).
(6) Update the packing space index ips:= ips+1;
determine an arrangement index iar for the current packing space by means of the given packing sequence Pseq and the space index ips: iar:= Pseq[ips];
insert the pair (pcurr, ArrList[iar]) to the stowing list StList.
(7) Determine the residual packing spaces for pcurr after filling pcurr with arrangement ArrList[iar], and insert these new packing spaces into PsList.
(8) Update the set of residual boxes BRes.
(9) If the packing space list PsList is not empty, go to step (3).
(10) Output stowing list StList.
Figure 4. Overall algorithm for the basic heuristic.
The stowing list StList records the position of every packing space and the placement of each loaded box within its packing space. In this way, StList describes unequivocally a complete feasible solution.
Representation of feasible solutions
The basic heuristic is called once to generate one single feasible solution to the packing prob- lem. A vector Pseq, called packing sequence, is handed over to control the generating process (see step 6 in figure 4). Each entry in Pseq represents one packing space and determines the arrangement required to realize in that space. For instance, if Pseq contains the entry “2” in po- sition 3, then the second best local arrangement generated for the third packing space is real- ized. The transformation process, leading from a packing sequence to an evaluated solution, is illustrated in figure 5.
{ 13, 2, 7, 9 } transformation
v = f(s) evaluation
packing sequence Pseq valid solution s volume
utilization Every entry in Pseq encodes
one specific local arrangement relative to a specific packing space.
Figure 5. Transformation of a packing sequence into a valid solution.
Some transformation rules are introduced to ensure that a packing sequence is transformed un- ambiguously into a feasible solution. If the packing sequence contains more entries than filla- ble packing spaces, the surplus entries are not relevant. If it turns out during transformation that there are more packing spaces to fill than entries in the packing sequence, the entry “1” is added several times to the packing sequence in order to realize local optima for the lacking packing spaces. If the entry at a specific vector position is higher than the total number of de- termined local arrangements, the arrangement index is interpreted in a cyclical way. Since packing sequences represent feasible solutions in an unambiguous manner, they can now be used to introduce a convenient neighborhood structure.
Neighborhood structure
According to the initial PTS, a large and a small neighborhood are distinguished. In both cases, a solution is considered as a neighbor when the corresponding packing sequences differ in ex- actly one significant index position common to both packing sequences. At this vector position, a different local arrangement for the relevant packing space is realized. For the large neighbor- hood, all index positions may be modified within an iteration, while in the case of a small neighborhood, this index is pre-determined and cyclically incremented once per iterative step.
In order to enable the algorithm to find the global optimum, a neighborhood type should guar- antee that one or more global optimal solutions are reachable from any feasible solution by a sequence of moves based on the neighborhood type. This so-called reachability condition is met by each of the neighborhood types introduced above. In the case of the small neighbor- hood, several cycles through the sequence of packing spaces may be required.
The type of neighborhood is specified by the parameter nbhType. For the initial PTS, an addi- tional parameter, nbhDecrease, is introduced in order to limit the number of neighbors evalu- ated for both types of neighborhood. The lowest value nbhDecrease = 1 corresponds to the maximum number of neighbors.
Objective function
The quality of a solution s is determined by the volume utilization of the load f(s). f is ex- pressed as the percentage relative to the total container volume. Since high volume utilization is aimed for, f is a benefit rather than a cost function.
Diversification concept
The search process is divided into several search phases. For diversification purposes, an inde- pendent set of parameters for the generating basic heuristic is applied in each of the search phases. The form of local arrangements can be influenced by the parameters phMinSupport and aboveArr. The number of local arrangements generated for a given packing space is limited by the parameter maxArr. The evaluation of local arrangements is controlled by arrEvalMode, while cutEvalMode affects the principles of generating residual packing spaces. For a detailed description of the five phase-parameters (and the parameter, nbhDecrease, mentioned above) the reader is referred to BORTFELDT,GEHRING and MACK (2002).
3 Generic aspects of the simulated annealing approach
Simulated annealing was first introduced by Metropolis et al. (1953) for modelling the physical cooling process of a solid in a heat bath. KIRKPATRICK et al. (1983) proposed to adapt the SA paradigm to combinatorial problems. SA is essentially local search, in which a move to an infe- rior solution is made with an acceptance probability that decreases as the search progresses.
The level of accepting a downhill move (when searching for a maximum) depends on:
- the magnitude of the decrease in the benefit function f, and
- a dynamic parameter, called temperature t, in analogy to its physical origin.
The passage to inferior solutions should allow the algorithm to escape from local optima and penetrate formerly unexplored regions of the search space. The decrease of the acceptance probability is motivated by the empirical evidence that the need to leave a local optimum is high at the beginning and decreases during the search process. Convergence results for SA have been proven using Markov chains representing the sequences of random neighbors deter- mined during the search process (cf. AARTS and KORST 1997). The usefulness of these theo- retical results in practice is however limited, since the annealing has to be very slow to guaran- tee convergence, the algorithm often requiring more iterations than an exhaustive enumeration.
Nevertheless, SA was implemented successfully in many applications, producing good results in acceptable computing time. Figure 6 shows the overall algorithm of the SA strategy for a maximization problem.
Select an initial solution: s:= s0;
select an initial temperature: t:= tStart (tStart > 0);
repeat repeat
select solution s' in neighborhood N(s) at random;
df := f(s') - f(s);
if df > 0 then s:= s';
else
generate random number x in (0,1);
if x < exp( df / t ) then
s := s';
endif endif
until predefined number solPerTemp of solutions per temperature t is examined;
reduce temperature t by reduction function r(t): t:= r(t);
until stopping condition is true.
Figure 6. Overall algorithm of simulated annealing for maximization problems.
Generic decisions involve the adaptation of the general components of SA to the concrete prob- lem. In the following, these general components are adapted to the given packing problem, making use of the problem-specific components of the original PTS described in section 2.
Initial solution
The search starts with the Greedy solution, s0, described by the packing sequence {1,1,1,..}, i.e.
for every packing space identified during the loading process, the highest ranked element ArrList[1] is realized. However, the starting solution does not play an important role, since at
high temperatures, the rate of accepted moves to inferior solutions is high, destroying the initial solution rapidly.
Selection of random neighbors
To perform a single iteration, i.e. to examine a further solution, exactly one member of the cur- rent neighborhood is selected at random. Since a solution is fully described by its packing se- quence, two decisions are made to determine a valid neighbor: First, a changing index in the packing sequence is chosen. Then, the new entry for this position is defined.
The parameter nbhType specifies how to choose the index depending on the nature of the neighborhood. The large neighborhood encloses all packing sequences differing in one arbi- trary significant common position ipr. Hence, the changing position is chosen by a random process among all available positions. In line with the original PTS, in the case of a small neighborhood, the differing position is not chosen randomly, but it is increased in each iterative step.
For the retained vector position representing a specific packing space, a random integer j be- tween 1 and n is generated, n representing the total number of local arrangements created for the given packing space. Subsequently the jth local arrangement for the given packing space is realized. The new packing sequence, s', is then transformed into a feasible solution in order to determine its volume utilization, f(s'), and to compare it to the best volume utilization, f(s), found so far.
The restriction of the neighborhood size by the parameter nbhDecrease, initially designed for the PTS to limit the number of solutions to be transformed, is unnecessary here because, unlike the TS strategy, only one neighbor is visited per iteration. Therefore, nbhDecrease is set to its minimum value 1.
Acceptance of neighborhood solutions
Since a selected neighborhood solution is accepted in accordance with the overall algorithm (cf. figure 6), only the calculation of the acceptance probability is explained in detail.
The acceptance probability p is calculated by p = exp(1000*df/t), where t represents the tem- perature parameter, and df is the difference in the relative volume utilization of the load, ex- pressed as a percentage of the container volume. In our case of a maximization problem, df is negative for inferior moves and lies in the interval [-100,0]. As an example, suppose that the volume utilization of an initial solution corresponds to f1 = 92% of the container volume, while the inferior neighbor has f2 = 91%. The difference df = f2-f1 is then -1 and the probabil- ity of accepting this inferior solution at a temperature t = 900 is p = exp( (-1) * 1000/900) = 0.33, while p = 0.14 when t = 500. The coefficient 1000 in the formula for calculating p was determined empirically and it enables the use of low temperature values up to 1000 “degrees”.
Cooling schedule
A cooling schedule is specified by (cf. AARTS and KORST 1997):
- a reduction function r(t) for decreasing the temperature t during the search,
- a parameter solPerTemp, that determines the number of examined solutions for a given value of temperature, t, and
- a temperature range given by the initial value tStart and the final value tEnd of the tempera- ture parameter, t; tStart is, at the same time, the maximum value, tEnd is the minimum value of t.
Above all, a cooling schedule is characterized by its reduction function. In the literature differ- ent types of reduction functions or schedules are recommended (cf., e.g., DOWSLAND 1993).
Three schedule variants, namely geometric, linear and adaptive cooling have been taken into account. As to the adaptive cooling, a temperature decrease should only be processed after the cost function has stabilized within certain limits according to the floating average and standard deviation. However, a sufficient stabilization of the cost function could not be observed and the results did not justify the calculation costs. The best results in pre-tests were gained by the geometric schedule, which has been proven to perform well in many optimization problems (cf. DOWSLAND 1993). The reduction function of the retained geometric cooling is given by r(t) = m* t, where the constant coefficient m is chosen from the interval (0,1).
Some experiments have been carried out to decide whether the parameter, solPerTemp, should be designed as a constant or as a dynamic parameter depending on temperature, t. Since in- creasing values for solPerTemp during the annealing process did not result in a higher level of solution quality, the more simple alternative of a constant parameter was finally preferred.
The determination of appropriate values for the cooling parameters m, solPerTemp, tStart and tEnd, remains the subject of experimental fine-tuning, which is postponed until section 6. Note please that the cooling parameters also specify the termination criterion of the search. The total number of iterations or examined solutions corresponds to:
nTotalIter = solPerTemp * ln( tEnd/tStart ) / ln(m).
Diversification of the search
As already stated in section 2, the search process is divided into a number of nPhases search phases. A restart is performed at the beginning of each phase. Each phase starts again with the Greedy solution s0 and works with the same cooling schedule. For diversification purposes, the configuration of the phases should be different. However, since SA is a randomised process, identical configuration of all phases does not necessarily produce redundant results, since each search run pursues its individual random walk. The respective phase-parameters, phMinSup- port, aboveArr, maxArr, arrEvalMode and cutEvalMode, were described in section 2.
As in the TS procedure, the termination of the overall method is caused by two criteria. The search is terminated, if all nPhases search phases have been carried out or if the calculation time exceeds a given time limit, maxTime, or if an optimal solution is identified. The best solu- tion over all search phases is finally determined as the solution to the problem.
Modifications to the standard SA paradigm
Some modifications of the standard SA algorithm have been tried out in order to improve effi- ciency and solution quality. According to different recommendations in DOWSLAND (1993), the following variants were envisaged:
- Constant temperature annealing: The temperature is set to a promising value, no cooling during the search process is performed. This measure is motivated by the empirical result that most useful work of the SA search process usually takes place over a limited tempera- ture range in the middle of the run.
- Geometric reheating: After a specified number of consecutive rejected solutions, the system is considered to be frozen, and reheating seems a good way to avoid getting stuck in a local optimum. Alternatively, a special version is used, in which each accepted solution results in a decrease in temperature and each rejected solution in an increase. A promising pre- determined ratio of accepted to rejected moves was applied.
- Multiple consecutive runs of the schedule were compared to a single slowly cooled run. To perform this test, the total number of evaluated solutions and the factor m were maintained, while the number solPerTemp of generated solutions per temperature and the number of con- secutive runs were varied.
- Threshold accepting: In this variant, the decision whether to accept inferior moves only de- pends on the deterioration of the benefit function which should not exceed a static parameter thresh (threshold). The stochastic element relative to the acceptance of inferior moves is omitted.
Some effort was invested to find promising parameter settings for each of the modifications above. Nevertheless, none of the tested SA variants was able to reproduce the results of the simple geometric schedule. Therefore, all the tested variants were finally rejected. However, it cannot be excluded that there are good settings for the proposed modifications, especially for an adaptive temperature setting during the search. This is due to the fact that, even if SA is easy to implement, fine-tuning is a long, complicated and sensitive process, and that in most cases, it is impossible to verify all conceivable settings of the parameters for a representative number of problems.
4 Hybrid local search approach
General considerations
A further opportunity for improving the solution quality could lie in an appropriate combina- tion of the local search strategies TS and SA. This new hybrid approach should integrate the advantages of both metaheuristics and avoid their respective weaknesses.
In the case of SA a neighborhood is left immediately after a randomly selected neighbor has been accepted. Therefore, SA has the ability, at least at high temperatures, to move quickly from any current solution to any other solution of the feasible search space. For the same rea- son SA risks omitting promising neighbors. TS, on the other hand, can perform an intensive examination of limited regions of the search space, since all neighbors of a given solution are visited. At the same time TS risks wasting computing time. This will occur if the currently in- vestigated region does not include any global optimal (or near optimal) solution. In brief, a hy- brid approach should make use of the different features of both strategies according to the bal- ance between extensive and intensive search.
DOWSLAND (1993, p. 51) suggests that approaches of combining SA with other heuristic meth- ods can be “confined to some form of pre-processing before the annealing process is invoked, or post-processing in which attempts are made to improve solutions encountered by the anneal- ing process”. A variant in which SA is preceded by several TS-phases in order to produce a good starting solution for SA turned out to be inappropriate in the given case, since starting the run at high temperatures will effectively destroy the characteristics of the initial configuration rapidly, while starting it at lower temperatures risks getting stuck in this starting solution.
We prefer, therefore, to apply SA first to determine a good starting solution for TS. This an- nealing solution is post-processed by one or more TS-iterations. The first post-processing TS-
iteration should ensure that at least the nearest local optimum to the final annealing solution has been found. The subsequent TS-iterations are intended to an intensive examination of the region around the final annealing solution.
Hybrid approach and diversification concept
The next challenge consists of combining the hybridization described above with the diversifi- cation concept introduced with the definition of SA as well as TS. The following conventions apply to this issue:
- Within the hybrid approach, both SA and TS are processed in distinct search phases with specific parameter sets.
- A final solution determined by SA is post-processed by a single TS-phase consisting of a pre-determined number of TS-iterations, with the aim of improving the solution quality.
- In this post-processing TS-phase the annealing solution is used as a starting point and the phase parameter set underlying this SA-solution is adopted. This measure avoids problems that can occur when the transferred packing sequence is decoded. To ensure that all neighbors are evaluated, nbhDecrease is set to one (cf. section 2).
A supplementary phase-parameter hybNat is introduced to implement different hybridization variants according to the frequency of alternating SA/TS-phases (cf. figure 8):
- hybNat = 1: a post processing TS-phase is provided for every single SA phase in order to improve the best SA phase solution;
- hybNat = 2: a single post-processing TS-phase is applied to the best annealing solution of all applied SA-phases.
hybNat=1 SA1
SA2
SA3 TS1
TS2
TS3
hybNat=2 SA1
SA2
SA3
TS2
time
restarts restarts
time
Figure 7. Variants of hybridization for different values of the parameter hybNat.
In figure 7, for the setting hybNat = 2, phase parameter set 2 is assumed to produce the best annealing solution for the three sets, and as a consequence, only this result is post-processed by a TS-phase using the corresponding phase parameter set 2.
From figure 7, it can also be stated that by applying hybNat = 2, the calculation time can be reduced compared to hybNat = 1, since only one TS-phase is processed. On the other hand, it cannot be excluded that for hybNat = 2, the results are inferior to the results of hybNat = 1.
The attempt to apply a TS-phase more frequently during the annealing process within a particu- lar SA-phase, after a fixed number of SA-iterations or after a certain number of rejected moves, and to reinject the determined TS-solution into the ongoing SA-phase, increased the calculation time significantly without any notable positive effect on the solution quality. As a conse- quence, this variant was rejected.
5 Parallel implementation
General considerations
According to TOULOUSE et al. (1996), three types of parallelization strategies seem appropriate for the methods used in combinatorial optimization: (1) parallelization of operations within an iteration of the solution methods, (2) decomposition of the problem domain or search space, and (3) multi-search threads with various degrees of synchronization and cooperation. Which of these types is suited for the parallelization of an optimization method depends mainly on the goal pursued by the parallelization.
Since an enhancement of the solution quality is the focus here, an approach of type (3) is cho- sen for the parallelization. The container-loading problem is subjected to several processes.
Each process is an instance of the sequential version of the TS, SA or hybrid algorithm, respec- tively, and solves the complete problem. Each of the autonomous processes is assigned to a workstation of a local area network (LAN). Hence, the parallel approach is a distributed- parallel method.
Three parallel approaches are considered and tested alternatively: parallel TS (PTS), parallel SA (PSA) and parallel hybrid (PHYB). Any forms of mixing the search strategies by assigning distinct search strategies to the different cooperating instances are not taken into account. This is motivated by the initial objective to assess the net improvement due to the parallelization for each of the three search strategies (cf. figure 2). Furthermore, the determination of adequate combinations of the search strategies in the above mentioned sense would require some con- siderable extra fine-tuning effort.
Parallelization concepts for SA algorithms
The parallelization of TS is dealt with extensively in BORTFELDT,GEHRING andMACK (2002).
Therefore, the following comments are limited to the parallelization options for SA (and for the hybrid strategy as a derivative of the SA approach).
AZENCOTT (1992) presents five possible ways of parallelizing the general SA algorithm. Two basic variants, “simultaneous independent searches” and “simultaneous periodically interacting searches” are of interest with regard to the concrete SA-implementation. As the name suggests, the independent searches use several processes that do not exchange information during the search. The final result is determined by the best outcome of the individual processes. By con- trast, simultaneous periodically interacting searches cooperate by a periodical transfer of solu- tions, which serve as starting points for further search of the individual processes.
AARTS AND KORST (1997) distinguish generic parallel SA algorithms, alternatively with or without communication. An improvement in the speed is the aim of AARTS and KORST’swork, leading to a division of the serial processes and distribution of the tasks among the available processors. This approach is not appropriate for the retained implementation of type (3) ex- posed above, and as a consequence, this variant of parallelization is not considered.
Communication of the parallel instances
The original PTS includes the possibility for the processes to communicate by the exchange of high quality solutions at the end of each search phase. This transmission of solutions should serve to intensify the search within the regions of best solutions. The improvements due to this communication feature turned out to be very small.
As to the general parallel SA models presented above, AZENCOTT (1992) shows that, in theory and given some premises for the cooling schedule, p periodically interacting annealing proc- esses are less efficient than p simultaneous independent annealing schemes. Hence, periodic interaction between p processors performing identical annealing searches would be “a waste of computing time and might just as well be replaced by a single interaction at the end of the available computing time to select the best terminal configuration” (AZENCOTT, 1992, p. 40).
Experimental comparisons between sequential annealing and p simultaneous periodically inter- acting annealing conducted by GRAFFIGNE (1992) seem to confirm the superiority of independ- ent annealing searches.
In spite of these negative prognoses, considerable effort was invested in order to configure communicating hybrid instances, in analogy to the parallelization techniques employed in BORTFELDT,GEHRING andMACK (2002):
- by using a blackboard or ring model, alternatively;
- by using unconditional versus conditional adoption of foreign solutions;
- by applying different combinations of parameter sets to the instances;
- by trying out different sequels of alternating SA and TS phases;
- by reducing the starting temperatures in the subsequent phases.
The results were fairly disappointing, and as a consequence, the idea of communicating SA or hybrid instances was finally abandoned. The implemented configuration with p processes without communication can be summarized as follows:
- p instances of the serial SA algorithm or hybrid algorithm, respectively, are run simul- taneously and independently;
- the parallel search ends if one of the following conditions is met: all instances have finished the calculation, a predetermined CPU-time bound maxTime has been exceeded, or a global optimal solution has been found by one of the instances (i.e. all boxes are loaded or the con- tainer is 100% full);
- the parallel result is determined by the best output of all p instances at the end of the search.
In the case of the PTS, the phase parameters were different for every parallel instance, in order to guarantee a certain diversification of the search. In the case of PSA and PHYB, the inherited random character of SA prevents the search attaining similar solutions, even if all processors are configured identically, since:
- the neighbors are selected randomly, and during the search, a large number of neighbors are determined, so that the probability of pursuing identical walks for distinct test runs is very small;
- the acceptance of inferior solutions is subject to a random decision process.
The details of the parameter settings are discussed in the next section.
6 Numerical tests
The following tests are based on the well-known reference problems from BISCHOFF and RATCLIFF (1995), consisting of 700 problem instances grouped in seven test cases. Within each test case, the problems present a constant number of distinct box types. The nature of the load varies from weakly to strongly heterogeneous box sets. The 700 BISCHOFF and RATCLIFF prob- lem instances, called hereafter ‘BR-problems’, include an orientation constraint, and the mini- mum stability coefficient minSupport is set to 55%.
The test environment consists of several Intel Pentium computers with a core-frequency of 2GHz combined in a local area network (LAN). The method was encoded in C using MS Vis- ual C++ 6.0.
In the following the procedure and the results of a configuration pre-test are reported before the different algorithms (cf. figure 2) are compared. Finally, a comparison test including some methods from other authors will be presented.
Configuration pre-test
The pre-test implies several steps. After determining the most promising basic heuristic pa- rameter sets, the fine-tuning of the cooling parameters is discussed. Finally, the settings for hy- bridization and parallelization are fixed.
The phase-parameter sets for the SA algorithm are determined in a limited series of experi- ments. These parameters steer the process of finding valid solutions and may influence the overall solution quality. Initially the six best phase-parameter sets were adopted, in analogy to the PTS; they are presented in table 1.
Phase-set-Id Parameter
1 2 3 4 5 6
arrEvalMode 1 1 0 1 1 1
cutEvalMode 0 1 1 0 0 1
maxArr 50 50 50 50 100 100
aboveArr 1 1 1 0 1 1
phMinSupport 55 55 55 100 55 55
nbhType large large large large large large
nbhDecrease 1 1 1 1 1 1
Table 1. Phase-parameter sets for SA and hybrid search.
By monitoring the results of the individual phases for all instances during the pre-tests, it was possible to simulate the effects of reducing the number of phases. It turned out that a reduced number of four phases results in a great saving in calculation time of up to 42%, while the ob- jective function decreased insignificantly by about 0.1% of the container volume, on average.
In particular, the application of parameter set 3 causes a sophisticated evaluation of the con- ceivable residual space structures and consumes a great part of the computing time. As a con- sequence, only the parameter sets 1, 2, 5 and 6 are applied to the SA algorithm, and only four SA-phases are carried out.
The parameter sets in table 1 differ from the original PTS settings in several points, especially relative to the neighborhood parameters (cf. BORTFELDT, GEHRING and MACK 2002, table 2).
Moreover, the TS algorithm has been run with six instead of four search phases. This should, however, not affect the comparability of the different approaches. On the one hand, it seems to be legitimate to use the best settings empirically found for every search strategy. On the other hand, the reduction of the number of search phases in the case of the SA algorithm helps to keep the average calculation times comparable.
Next, the fine-tuning of the cooling parameters for the geometric schedule is discussed. Since it is very hard to set all four parameters simultaneously, the parameters are grouped to distinct categories:
- m and solPerTemp determine the speed of the annealing and are therefore considered to be
“rate parameters” of the cooling schedule;
- tStart and tEnd set the covered “temperature” range.
One cooling-rate- and one temperature-range-parameter are fixed to empirically good values.
Thus, the fine-tuning is limited to the process of synchronizing the remaining parameters, tEnd (range parameter) and solPerTemp (rate parameter). Good starting temperatures of 900 and 700 degrees, representing a probability p of 0.33 and 0.24, respectively, of accepting a downhill move of 1% of the container volume, were rapidly found. A cooling coefficient, m = 0.999, produced good results in the validation tests. Some effort was then invested in order to find ap- propriate combinations of the remaining parameters, solPerTemp and tEnd. To achieve this fine-tuning, the parameters are combined in such a way as to fix the total number of calculated solutions, nTot. The parameter solPerTemp was set to 4, 8, 12 and 16, alternatively and it turned out that solPerTemp = 8 is a promising setting.
As to the hybridization, its type is determined by hybNat = 1, accepting longer calculation times (cf. figure 7), since the optimization of the volume utilization is the focus of this test. The neighborhood parameters, nbhDecrease and nbhType, for the post-processing TS-phase are set to values leading to maximum neighborhoods in order to allow the most intensive search (cf.
table 2). The number of iterations, maxIter, performed within a TS-phase is set to a unique value.
Basic heu- ristic
SA cooling parameters Parameters for the post- processing TS phase Serial
phase
Phase-
set-Id tStart tEnd m solPer- Temp
nbh- Type
nbh- Decrease
max- Iter First 1 700 300 0.999 8 large 1 10 Second 2 700 300 0.999 8 large 1 10
Third 5 700 300 0.999 8 large 1 10
Fourth 6 700 300 0.999 8 large 1 10 Table 2. Parameters for the tests.
Remaining decisions involve the configuration of the parallelization. Since during the pre- tests, the application of differing parameter settings for all parallel instances presented no fur- ther opportunity for an amelioration of the solution quality, every process is set up with identi- cal phase-, cooling- and TS-parameters. These parameters are directly derived from the respec- tive serial settings. The parallel tests are based on four instances or processors.
Results for all variants of the proposed metaheuristics
The results of extensive tests applying the different search strategies to the 700 BR-problems are reported in table 3.
While the deterministic TS results were found within a single test run, SA and hybrid search require several runs to guarantee representative results. Consequently, the serial results in table 3 represent average values over twelve test runs, the parallel results are calculated from three runs. By aggregating the results of all the tests performed, the result for a single test run of 64 instances could be derived (PHYB.XL).
The combination of parallelization and hybridization yields the best solution quality. It can be stated that the volume utilization decreases for problems with more than eight distinct box
types (i.e. BR4-BR7), reflecting the specific design of the basic heuristic for weakly heteroge- neous box sets.
Average volume utilization as a percentage of the container volume Serial methods Parallel methods
Test case Number of box
types
Average number of boxes
per type TS SA HYB PTS PSA PHYB PHYB.XL BR1 3 50 93.23 93.04 93.26 93.52 93.24 93.41 93.70 BR2 5 27 93.27 93.38 93.56 93.77 93.61 93.82 94.30 BR3 8 17 92.86 93.42 93.71 93.58 93.78 94.02 94.54 BR4 10 13 92.40 92.98 93.30 93.05 93.40 93.68 94.27 BR5 12 11 91.61 92.43 92.78 92.34 92.86 93.18 93.83 BR6 15 9 90.86 91.76 92.20 91.72 92.27 92.64 93.34 BR7 20 7 89.65 90.67 91.20 90.55 91.22 91.68 92.50 Average – – 92.00 92.53 92.86 92.70 92.91 93.20 93.78 Table 3. Test results for the Bischoff and Ratcliff test cases.
The standard deviation for the 100 results within a particular test case is given by the mean value of 1.34% for the serial algorithm HYB and 1.28% for the parallel algorithm PHYB, on average over the seven test cases.
On average, over 700 test instances, the calculation time determined on 2GHz-Pentium-PCs for serial SA is 72s (81s for the corresponding parallel version), while hybrid search required 205s (222s for the parallel version). Note please that the computing time for a parallel run is given by the maximum time needed by the individual serial instances employed. The results in table 3 are illustrated in figure 8.
89 90 91 92 93 94 95
BR1 BR2 BR3 BR4 BR5 BR6 BR7
TS SA HYB PTS PSA PHYB PHYB.XL
Figure 8. Mean values of the volume utilization (as a percentage of the container volume) for the 7 test cases of Bischoff and Ratcliff.
test cases mean volume utilization
per test case
The following detailed observations regarding the improvement of the solution quality, the cal- culation costs and the number of evaluated solutions can be stated.
SA versus TS
The application of the SA strategy yields a higher solution quality than the application of the TS strategy. The average improvement obtained by SA for the 700 test instances compared to TS is 0.53% of the container volume for the serial, and 0.21% for the corresponding parallel version. Note that the parallel TS profits from the additional effect of communicating in- stances, while this is not the case for the independent SA processes. From the results, it can also be concluded that the difference between SA and TS is increasing with the complexity of the problem, given by the number of distinct box types. Hence, SA exhibits its conceptual ad- vantages with larger solution spaces. For small solution spaces (3-5 box types), TS produces better results, reflecting the ability of TS to explore such small search spaces efficiently.
The average calculation time of SA increases significantly compared to TS. In fact, the time to transform a single packing sequence into a valid solution and calculate the volume utilization function is higher for SA by a factor 2.5. This can be explained by some disadvantages of SA regarding the transformation process: TS is able to group the large number of solutions to be evaluated per iteration in order to reuse a certain part of formerly transformed packing se- quences in future transformations. On the other hand, the random walk of SA prevents this fea- ture. As a consequence, SA has to operate an expensive initialization of the data structures when evaluating a solution. Additionally, the process of the basic heuristic can be stopped in TS if it turns out that no new optimum will be reached. In SA, as any inferior solution possibly constitutes the starting point of the further search, the volume utilization has to be calculated for any solution, preventing an anticipated stopping of the calculation process.
The total number of evaluated solutions is increasing with the application of SA. Therefore, the results have to be interpreted carefully with regard to a fair comparison with TS. However, separate tests show that the net improvement of SA is 0.2 - 0.3% of the container volume for the sequential version when synchronizing the corresponding number of evaluated solutions in SA and TS and using absolutely identical basic heuristic parameters. This result remains re- markable, because in most publications, TS is considered to be the more promising strategy.
Hybridization and parallelization
Figure 9 summarizes the net improvements (as a percentage of the container volume) found in the tests due to hybridization and parallelization, compared to the ‘pure’ search methods, TS and SA.
TS (92.00%) PTS (92.70%)
HYB (92.86%)
PHYB (93.20%) PSA (92.91%)
SA (92.53%) +0.86% +0.33%
+0.70% +0.34% +0.38%
Figure 9. Net improvements due to hybridization and parallelization, for 700 Bischoff and Ratcliff instances.
The parallel hybrid version generates the best results within the tests. The isolated effects of the hybridization and parallelization are comparable. Starting from the ‘pure’ serial SA implemen- tation, hybridization yields 0.33% on average, against 0.38% for the parallelization. This ob- servation is also valid considering the hybridization and parallelization of TS (see figure 9), where the improvements were 0.86% and 0.70%, respectively.
Comparing the results from other authors
Table 4 compares the previous results with the outcome of the works of BISCHOFF et al. (1995), BISCHOFF andRATCLIFF (1995) and ELEY (2002). The detailed test case results reported for the method from ELEY were provided by the author.
Author’s calculations Test case BISCHOFF
et al.
(1995)
BISCHOFF and RATCLIFF
(1995)
ELEY
(2002) Serial HYB
Parallel HYB
Parallel HYB.XL BR1 81.76 83.79 88.05 93.26 93.41 93.70 BR2 81.70 84.44 88.44 93.56 93.82 94.30 BR3 82.98 83.94 89.23 93.71 94.02 94.54 BR4 82.60 83.71 89.24 93.30 93.68 94.27 BR5 82.76 83.80 88.99 92.78 93.18 93.83 BR6 81.50 82.44 88.91 92.20 92.64 93.34 BR7 80.51 82.01 88.36 91.20 91.68 92.50 Average 82.00 83.50 88.75 92.86 93.20 93.78 Table 4. Comparing the results to methods from other authors.
The results of the metaheuristic proposed in this paper for the 700 BR-problems clearly domi- nate the results of other authors which have been considered in table 4. However, only the vol- ume utilization of the load is taken into account. In particular, specific stability aspects are not considered. The results of the other authors indicated in table 4 were determined with fully supported boxes (minSupport = 100%), while the algorithm proposed here uses minSupport = 55%. Nevertheless the average support of all boxes loaded in the results is 95% for the reported solutions. Moreover, when setting minSupport to 100%, the presented metaheuristic remains superior, with an average volume utilization of 92.41% (variant serial hybrid HYB). The calcu- lation time remains reasonable compared to the times reported by the other authors.
Increasing the number of parallel processors
In general, the parallel strategy multiplies the number of evaluated solutions and achieves an improvement compared to the serial algorithm, by exploring a greater part of the solution space. The solution quality is increasing continuously with the number of processors. In figure 10, the results of a progressive increase of the number of parallel instances are summarized.
The indicated values represent averages of five distinct test runs obtained by random combina- tions of the distinct test runs performed with four processes and four phases per process. The volume utilization is increasing as a logarithmic function, leading to a maximum volume of 93.69% with 44 processors relying on the four phases of table 1.
93,59 93,62
93,68 93,69
93,2
93,51
93,64
93,55
93,66
93,37 93,46
93,00 93,10 93,20 93,30 93,40 93,50 93,60 93,70
4 8 12 16 20 24 28 32 36 40 44
Figure 10. Average hybrid result over 700 BR-problem instances, using a larger number of parallel processes configured with four search phases.
A calculation with 64 PCs, including 20 instances relying on the six-phases-parameter-sets of table 1, leads to a mean utilization volume of 93.78%. In this way, the residual volume could be reduced, on average, by 20% compared to the serial TS, i.e. from 8.0% to 6.22%. The aver- age calculation time for this test increased significantly to 596s.
7 Summary
This paper presents a simulated annealing algorithm for solving the container loading problem with a single container to be loaded. The algorithm constitutes a further development of the parallel tabu search algorithm from BORTFELDT, GEHRING and MACK (2002). The SA algo- rithm has been extended to a hybrid metaheuristic by post-processing the final annealing solu- tion with multiple tabu search iterations. The parallelization follows the concepts of TOULOUSE
et al. (1996) and AZENCOTT (1992), realizing simultaneous independent runs of the cooling schedule. Some forms of communication between the instances were tested, but all had to be rejected because of the disappointing results.
In the application of both search strategies, SA and TS, simulated annealing clearly wins over tabu search, especially for problems with larger solution spaces. The comparison is valid from a methodical point of view, since both strategies use identical problem-specific components, including the kernel heuristic, to generate solutions and the corresponding neighborhood defi- nition. However, the fine-tuning of the steering parameters of SA was a long and complicated process compared to the parameterization of TS.
Hybridization has been shown to be an effective and powerful extension of SA and TS, inte- grating the core advantages of both search strategies. In the parallel version, the improvements due to parallelization were comparable to those found with the hybridization. Parallelization showed an additional opportunity for improving the solution quality, especially when a large number of processors was used.
The presented approach constitutes one of the first SA algorithms proposed for three- dimensional packing problems. Additionally, the attempt of a fair comparison of the different metaheuristics, SA and TS, is not very often reported in the literature. Moreover, the volume utilization could be increased significantly for a large set of standard problem instances.
no. of processes mean volume
utilization [%]
processes
Acknowledgement
The authors wish to thank Mr. Attila Zadanji for his valuable assistance with the parallel tests of the heuristic presented.
References
AARTS, E. and KORST, J. (1997). Simulated Annealing and Boltzmann Machines. New York etc., John Wiley.
AZENCOTT, R. (1992). Parallel Simulated Annealing: An Overview of Basic Techniques. In:
Azencott, R. (ed.) (1992), Simulated Annealing: Parallelization Techniques, New York etc., John Wiley, pp. 37–46.
BISCHOFF, E.E. and RATCLIFF, M.S.W (1995). Issues in the Development of Approaches to Container Loading. Omega, Vol. 23, pp. 377–390.
BISCHOFF, E.E., JANETZ, F. and RATCLIFF, M.S.W. (1995). Loading Pallets with Non-identical Items. European Journal of Operational Research, Vol. 84, pp. 681–692.
BORTFELDT,A.(1995). Informierte Graphensuchverfahren und genetische Algorithmen zur Lö- sung von Containerbeladeproblemen (Dissertation). Verlag Köster, Berlin.
BORTFELDT,A.andGEHRING,H.(1998).Ein Tabu Search-Verfahren für Containerbela- deprobleme mit schwach heterogenem Kistenvorrat. OR Spektrum, Vol. 20, pp. 237–
250.
BORTFELDT, A. and GEHRING, H. (2001). A Hybrid Genetic Algorithm for the Container Load- ing Problem. European Journal of Operational Research, Vol. 131, pp. 143–161.
BORTFELDT, A., GEHRING, H. and MACK, D. (2002). A Parallel Tabu Search Algorithm for Solving the Container Loading Problem. Parallel Computing, Vol. 29, pp. 641 – 662.
DOWSLAND, K.A. (1993). Simulated Annealing. In: Reeves, C.R. (ed.) (1993), Modern Heuris- tic Techniques for Combinatorial Problems, Blackwell Scientific Publications, Oxford, pp. 20–69.
DYCKHOFF, H. and FINKE, U. (1992). Cutting and Packing in Production and Distribution. Hei- delberg, Physica-Verlag.
ELEY,M. (2002). Solving Container Loading Problems by Block Arrangements. European Journal of Operational Research, Vol. 141, pp. 393–409.
GEHRING, H. and BORTFELDT, A. (1997). A Genetic Algorithm for Solving the Container Loading Problem. International Transactions in Operational Research, Vol. 4, pp. 401–
418.
GEHRING, H. and BORTFELDT, A. (2002). A Parallel Genetic Algorithm for Solving the Con- tainer Loading Problem. International Transactions in Operational Research, Vol. 9, pp. 497–511.
GRAFFIGNE, C. (1992). Parallel Annealing by Periodically Interacting Multiple Searches: An Experimental Study. In: Azencott, R. (ed.) (1992), Simulated Annealing: Parallelization Techniques, New York etc., John Wiley, pp. 47–79.
HIFI, M. (2002). Approximate Algorithms for the Container Loading Problem. International Transactions in Operational Research, Vol. 9, pp. 747–774.
KIRKPATRICK, S., GELLAT, C.D. and VECCHI, M.P. (1983). Optimization by Simulated Anneal- ing. Science, Vol. 220, pp. 671–680.
METROPOLIS, N., ROSENBLUTH, A., ROSENBLUTH, M., TELLER, A. and TELLER,E.(1953).
Equation of State Calculations by Fast Computing Machines. Journal of Chemical Physics, Vol. 21, pp. 1087–1092.
MORABITO, R. and ARENALES, M. (1994). An AND/OR-graph Approach to the Container Loading Problem. International Transactions in Operational Research, Vol. 1, pp. 59–
73.
PISINGER, D. (2002). Heuristics for the Container Loading Problem. European Journal of Op- erational Research, Vol. 141, pp. 143–153.
REEVES, C.R. (ed.) (1993). Modern Heuristic Techniques for Combinatorial Problems. Oxford, Blackwell Scientific Publications.
SCHEITHAUER, G. (1992). Algorithms for the Container Loading Problem. Operational Re- search Proceedings, 1991, Berlin, Springer, pp. 445–452.
SIXT, M. (1996). Dreidimensionale Packprobleme. Lösungsverfahren basierend auf den Meta- heuristiken Simulated Annealing und Tabu-Suche (Dissertation). Frankfurt am Main, Peter Lang, Europäischer Verlag der Wissenschaften.
TOULOUSE, M., CRAINIC, T.G. and GENDREAU, M. (1996). Issues in Designing Parallel and Distributed Search Algorithms for Discrete Optimization Problems. Publications CRT- 96-36, Centre de recherche sur les transports, Université de Montréal, Canada.
