An analytics-based heuristic decomposition of a bilevel multiple-follower cutting stock problem

This paper presents a new class of multiple-follower bilevel problems and a heuristic approach to solving them. In this new class of problems, the followers may be nonlinear, do not share constraints or variables, and are at most weakly constrained. This allows the leader variables to be partitioned among the followers. We show that current approaches for solving multiple-follower problems are unsuitable for our new class of problems and instead we propose a novel analytics-based heuristic decomposition approach. This approach uses Monte Carlo simulation and k-medoids clustering to reduce the bilevel problem to a single level, which can then be solved using integer programming techniques. The examples presented show that our approach produces better solutions and scales up better than the other approaches in the literature. Furthermore, for large problems, we combine our approach with the use of self-organising maps in place of k-medoids clustering, which significantly reduces the clustering times. Finally, we apply our approach to a real-life cutting stock problem. Here a forest harvesting problem is reformulated as a multiple-follower bilevel problem and solved using our approach.


Introduction
Mathematical Optimisation has been used to solve practical problems in areas such as scheduling, planning, cutting stock, transportation as well as other fields ( Vanderbei 2007;Kantorovich 1960;Bazaraa et al. 2010). In practical applications, however, optimisation problems are rarely straightforward since challenges such as several levels of optimisation and/or nonlinearities arise.
A bilevel optimisation problem is one in which one optimisation problem is nested within another. There is a leader (or outer/upper-level) problem, and a follower (or inner/lower-level) problem. The variables are split into leader variables and follower variables. The leaders solution must optimise an objective, with the constraint that the follower optimises a different objective. Both the leader and follower influence each other: the leader makes a decision while taking the followers objective into account, and the follower reacts to the leader's decisions in an optimal manner. Due to this interaction, the solution of bilevel problems is difficult in the general case (Colson et al. 2007). Besides, bilevel problems become more difficult to solve when they contain nonlinearities (such as the problems analysed in this paper). A general mathematical formulation of the bilevel problems is as follows: where represents the leader decision vector and the follower decision vector. The objective functions at the upper-and lower-level are represented by F and f, respectively. G and g represent the inequality constraints at the upper-and lower-levels, respectively, and equality constraints may also exist. X and Y are the bound constraints for the upper-level decision vector and lower-level decision vector, respectively. The condition ∈ argmin forces to be an optimal solution to the follower problem, and can be considered as a computationally expensive constraint.
One type of bilevel problem is the single-leader, multiple-follower problem, which has been used in applications like toll-setting, resource management, conflict resolution and many others. While several classical and evolutionary solution approaches exist for solving these multiple-follower problems, they are either not applicable in cases in which the follower problems are not traditional optimisation problems, or do not scale up appropriately. For this reason, in this paper, we present an alternative heuristic decomposition approach that decomposes the bilevel lower-level problem into smaller sub-problems (each sub-problem is composed of a follower) which are solved to optimality. Then, these solutions are used for computing an approximate solution to the upper-level bilevel problem (composed by the leader).
Furthermore, intending to reduce the computation time of solving this type of bilevel problems (which typically are large-scale when associated with real applications), our approach uses data-driven/analytics-based methods (specifically, we use clustering techniques). In data-driven optimisation, the goal is to use insights (1) from any available data to improve the quality of the decisions made or the solution approaches taken. Analytics techniques can therefore be used to enhance decision making by learning useful features from the data or using data to simplify mathematical models, thus reducing the complexity and/or solution times of the models. This paper is an extension of the work originally presented in the 2019 World Congress on Global Optimization (WCGO) Fajemisin et al. (2020). In this extension, we have focussed on providing a detailed literature review, the use of selforganising maps as an alternative to clustering and applying our approach to a reallife cutting stock problem: the forest harvesting problem. The main contributions of the work presented are: (i) a new class of multiple-follower bilevel problems is proposed; (ii) a novel analytics-based heuristic decomposition approach for solving this class of large-scale bilevel multiple-follower problems is given; (iii) the forest harvesting problem is reformulated as a bilevel optimisation problem to take into account the cutting operations of the harvesters. This paper is organised as follows: a literature review of bilevel decomposition approaches is given in Sect. 2. The new class of bilevel multiple-follower problems we consider is given in Sect. 3. The heuristic decomposition method is described in Sect. 4, while numerical examples demonstrating the advantage of the analyticsbased approach over evolutionary approaches are given in Sect. 5. Our approach is applied to the bilevel cutting stock problem in Sect. 6, and general conclusions are given in Sect. 7.

Related work
Heuristic decomposition approaches divide the problem into subproblems, which are solved individually, and the results assembled into a solution to the original problem. None of these steps needs to be optimal, for example in Elkamel et al. (1997) the subproblems are solved optimally then combined into a near-optimal solution, while in Krüger et al. (1995) the subproblems are solved approximately. In these papers, scheduling problems are addressed. As an example of another application domain, in Lin et al. (2011) the authors present a decomposition heuristic for a network design problem. Our approach decomposes the problem into follower subproblems, and it is an approximate approach because we use sampling and clustering.
The following subsections provide an overview of the current solution approaches for bilevel problems and the motivation for an analytics-based approach.

Decomposition approaches in the literature
Bilevel problems are known to be non-convex, non-differentiable and strongly NP -hard even in the simplest cases (Colson et al. 2007), and so most approaches for solving them involve some form of decomposition of the problem into more tractable forms. The most common method for single-level reduction is to replace the lower-level problem with its Karush-Kuhn-Tucker (KKT) (Bertsekas 1999) conditions (Bard 1984;Dempe and Zemkoho 2012). The authors in Visweswaran et al.

3
(1996) present a decomposition approach for solving linear and quadratic bilevel problems. They transform the problem into a single-level one by replacing the inner problem with its KKT conditions. The problem is then decomposed into a series of primal and relaxed-dual sub-problems, whose solutions are used as lower and upper bounds. This procedure is run iteratively until a global optimum is found. This method is similar to that in Zeng and An (2014), where the problem is first reformulated into a single-level problem using KKT conditions and strong duality. Then, a similar iterative process is carried out until an optimal solution is found.
More complex approaches combine the use of the KKT conditions with other techniques. For example, in Kristianto et al. (2013) the stochastic bilevel problem is reduced to a single level using KKT conditions and is then solved using Benders' decomposition. Related approaches include (Saharidis and Ierapetritou 2009), where an algorithm for solving mixed-integer bilevel linear problems based on Benders' decomposition is presented. Similarly, (Raidl et al. 2014) uses logic-based Benders' decomposition to solve a bilevel vehicle routing problem, which is combined with a variable neighbourhood search heuristic to speed up search time and improve scalability. To address the problem of weak Benders' cuts, (Nishi et al. 2011) uses a Lagrangian relaxation method to generate stronger cuts for simultaneous scheduling and routing problem for automated guided vehicles.
In Dempe and Franke (2016) lower-level problems are replaced with their Fritz-John conditions (John 1948), and an algorithm is presented for solving problems with fully convex lower-levels. This method is applied in Dempe and Franke (2014) to solve a bilevel road pricing problem. Nogales Martín and Miguel (2004) show a relationship between one bilevel decomposition algorithm and a direct interior-point method based on Newton's method. The authors in Iyer and Grossmann (1998) present a decomposition algorithm for solving a network planning problem. The upper level is solved to get an upper bound, which is then used to get a solution for the lower level problem which provides a lower bound. The process occurs iteratively, adding integer cuts along the way until a small enough gap between the bounds is achieved.
In Sugiyama et al. (2012), the authors solve a railway crew rostering problem. Their decomposition is in the form of cuts to reduce the feasible region of the master problem. Local search is also incorporated to improve the upper bound generated by solving the sub-problems. In Caramia and Mari (2016) the decomposition algorithm proposed consists of solving the single-level relaxation (SLR) of the Bilevel Facility Location (BFL) problem, solving the slave problem (SVP) which is the BFL for a given fixed set of open facilities, generating cuts based on the structure of the problem, and repeating until a stopping criterion is reached.
As with single-follower bilevel problems, both classical and evolutionary approaches have been used for solving Bilevel Multiple-Follower (BLMF) problems. Lu et al. (2006) presents a general framework and solutions for these problems. Nine classes of multiple-follower problems are presented (none of which include the problem class proposed in Sect. 3) with corresponding models presented for each class. Also, an extended Kuhn-Tucker approach is presented for solving the uncooperative model to optimality. A practical example in the form of a road network problem is given. Similarly, ) uses a Kuhn-Tucker (KT) approach for BLMF problems in which the followers may or may not have shared variables.  and  use extended KT approaches. In  a branch-and-bound algorithm is used to solve the problems with referential-uncooperative followers. Calvete and Galé (2007) reformulates a problem with multiple followers into one with one leader and one follower by replacing the lower levels with an equivalent objective and constraint region. This method cannot be applied to the BPMSIF, as neither its objectives nor its inducible region are equivalent to those of the problem class addressed in Calvete and Galé (2007).
There is also literature on applying the Kth-best approach to problems with multiple followers.  presents the theoretical properties of BLMFs and Kth-best approach for solving such problems, while  uses the approach for problems with shared variables among followers. Similarly, (Zhang et al. 2008b) presents an extended Kth-best approach for solving referential-uncooperative BLMF decision problems, and provides an application in the form of a logistics planning problem.
Fuzzy approaches to solving BLMFs include ) which uses a fuzzy interactive algorithm to solve problems with partially shared variables among the followers. Zhang and Lu (2010) combines fuzzy models and a Kthbest algorithm to solve cooperative multiple-follower problems. Fuzzy models combined with a branch-and-bound algorithm have also been used in Zhang et al. (2008aZhang et al. ( , b, 2007 to solve problems with shared decision variables among the followers. Decomposition approaches that also incorporate evolutionary approaches exist. A Co-Evolutionary Decomposition-based Bilevel Algorithm (CODBA) is presented in , in which an algorithm is first used to generate a set of points from a discrete solution space. This allows them to generate a population of solutions for the lower-level problem. Several sub-populations of the lower-level problem are generated, and the best individuals in the subpopulations are allowed to co-evolve. CODBA II  is an improvement in which parallelism and co-evolution are implemented at both levels of the bilevel problem. Evolutionary approaches also exist where an iterative approximation of the reaction set is used to approximate the lower-level problem (Sinha et al. 2013. Additionally, decomposition approaches involving evolutionary approaches are given in Li and Wang (2010); Li et al. (2016). In terms of multiple-follower problems, the literature includes (Angelo and Barbosa 2015), where a differential evolution method is used to solve cases in which there is information shared between the followers. Liu (1998) presents a genetic algorithm for solving nonlinear multilevel problems with multiple followers. Also, (Islam et al. 2016) extend their bilevel memetic algorithm to solve problems with multiple followers using a combination of global and local search. The authors in Ke et al. (2016) combine fuzzy programming with an evolutionary algorithm, as well as neural networks to solve a multi-follower problem with non-cooperative followers. Reviews on the use of metaheuristic approaches in bilevel optimisation are given in Talbi (2013); .

Motivation for the analytics-based heuristic decomposition approach
In summary, both classical and evolutionary approaches have been applied to single-and multi-follower problems. Lu et al. (2006) presents a general framework and solutions for nine classes of multi-follower problem, but none are applicable to the new class of problems we consider in Sect. 3. The authors in Calvete and Galé (2007) reformulate a problem with multiple followers into one with one leader and one follower, by replacing the lower levels with an equivalent objective and constraint region. This method also cannot be applied to our problem, as neither its objectives nor its inducible region are equivalent to the problem class of Calvete and Galé (2007). Additionally, the methods proposed in Calvete and Galé (2007); Lu et al. (2006) assume that the followers are linear, which is not the case with our class of problems. Most classical methods for handling bilevel problems require assumptions of smoothness, linearity or convexity, while we make no such assumptions. Evolutionary and meta-heuristic techniques also do not make these assumptions (Angelo and Barbosa 2015;Islam et al. 2016;Liu 1998) but most are computationally intensive nested strategies. They are efficient for smaller problems but do not scale up well to large-scale problems. In contrast, our analytics-based approach scales up well as the number of followers increases (see Sect. 5).

A new class of bilevel problems
In Bilevel Multi-Follower (BLMF) problems there may be several followers, and multi-leader problems are also known.
For a BLMF problem with Q followers, let represent the leader decision vector, and q the decision vector for follower q ( q = 1 … Q ). The leader chooses a strategy , following which each follower selects its own strategy in response to . The BLMF problem may be either cooperative, partially cooperative or uncooperative, depending on the leader and follower objectives. Based on the type of interaction between followers, nine classes of linear BLMF problems are identified in Lu et al. (2006). Problems in which the followers do not share objectives or constraints are known as "independent" and take the form: Several researchers have worked on bilevel optimisation with multiple independent followers (Lu et al. 2006;Zhang and Lu 2010). However, we strengthen this independence condition to one we call strong independence.
An analytics-based heuristic decomposition of a bilevel…

Definition 1 A Bilevel Problem with Multiple Strongly Independent Followers (BPMSIF) is one in which:
(i) the followers do not share each others' follower or leader variables, so that can be partitioned into q parts: q ( q = 1 … Q); (ii) follower functions f q ( q , q ) are allowed to be nonlinear; (iii) variables from different follower problems are not tightly mutually constrained.
In (iii) weak constraints such as a single linear inequality are allowed. Thus the BPMSIF has the form: where F, f q may be any (possibly nonlinear) objective functions, G, g q may be any set of (possibly nonlinear) constraints, the G constraints are weak, and X q , Y q may be vectors of any variable domains (real, integer, binary, or richer Constraint Programming domains such as set variables). Examples of bilevel real applications with nonlinearities are (Herskovits et al. 2000), (John 1948) and (Calvete et al. 2011). The last two are applied to planning problems. In this paper, in Sect. 6 we address a cutting stock problem associated with a forestry harvesting real application, which is a BPMSIF, as stated in Definition 1.
Problem (2) satisfies the features of a BPMSIF. Firstly, each follower function here can be seen to be a function of only its variables q and a portion of the leader's variables q . Secondly, G( 1 , … , Q , 1 , … , Q ) ≤ 0 is weak and may, for example, take the form of a simple weighted sum such as ∑ Q q B q y q ≤ b , where the B q and b are constants. The BPMSIF is different from multi-leader problems such as those of (DeMiguel and Xu 2009; Lu et al. 2016;Ramos et al. 2016). Its constraint region is: The projection of Ω onto the leader's decision space is: The feasible set for follower q is affected by a corresponding part q of a given leader decision vector so that: Each follower's rational reaction set is given as: Finally, the inducible region ( IR ) is: As in standard bilevel programming min and argmin have been used without loss of generality: each subproblem may involve maximisation.

Analytics-based heuristic decomposition approach
In this section, we first explain the Monte Carlo simulation and then we analyse two alternatives: k-medoids clustering and self-organising maps.

Monte Carlo simulation
For each follower q a large number S of feasible solutions for the leader vector q associated with that follower are generated, using Monte Carlo simulation. To avoid bias the q are generated using Hypersphere Point Picking (Marsaglia 1972;Muller 1959), which uniformly samples from a vector space. This results in a set sq ( s = 1 … S , q = 1 … Q ). The associated follower problems are then solved using the sq to obtain a corresponding set of follower vectors sq . We now have multiple potential leader solutions, together with their corresponding follower solutions for each follower problem f q ( q , q ).

k-medoids clustering
In order to model and solve the BPMSIF as an Integer Linear Program (ILP), the large number of potential solutions sq must be reduced to a manageable size, which we do via k-medoids clustering (de Amorim and Fenner 2012). Unlike k-means where centroids are means of data points, k-medoids selects actual data points. Once the large set sq has been clustered using K clusters, the medoids from each cluster are selected to represent the large generated set. The corresponding kq are then selected from sq so that we now have a smaller representative set of assignments kq and kq . The most common algorithm for k-medoid clustering is the Partitioning Around Medoids (PAM) algorithm (de Amorim and Fenner 2012). This is inefficient on very large datasets, so instead we use the CLARA algorithm which is a combination of PAM and random sampling (Kaufman and Rousseeuw 2009;Wei et al. 2000). The BPMSIF can now be transformed into a single-level optimisation problem: An analytics-based heuristic decomposition of a bilevel… The constraint q = kq → q = kq ensures that if q is assigned a value in kq then q is assigned the corresponding value in kq . This constraint can either be linearised using the big-M approach, or implemented directly using CPLEX indicator constraints (IBM 2017).

Self-organising maps: an alternative to clustering
As an alternative to k-medoids clustering (which can be time consuming), a neural network approach is also considered. Self-organising maps (SOMs) are a type of artificial neural network in which, given a set of input data vectors, individual neurons in the map compete to align themselves with the vectors they are best matched with (Kohonen 1982). SOMs operate in two phases-a training phase in which the map is built using sample data, and a mapping phase in which new data is automatically classified based on the map built in the previous phase, with the number of neurons of the SOM determining the number of classes. In this paper however, the intention is not to use SOMs for classifying new data. Instead, the map generated by the SOM trained over the large data set is used to select a smaller but highly representative subset of the data. (Note that the size of the reduced set is equal to the number of neurons.) This is analogous to clustering (Bullinaria 2004; Lampinen and Oja 1992). This reduced dataset is then used to solve the bilevel optimisation problem. Figures 1, 2 and 3 give an illustration of this approach.
Starting with a 2-D lattice of neurons, as shown in Fig. 1, the map is trained on the underlying data. After several iterations, the map will have aligned itself with  Fig. 2. Once the map is aligned, the data points closest to each neuron is selected, one data point for each neuron. Figure 3 shows a zoomed-in section of Fig. 2, with the neurons as octagonal structures and the selected data points in red. This reduced dataset can then be used to solve the transformed multiple-follower bilevel problem.
To use the SOM for clustering, i.e. reducing sq to kq , the number of neurons in the lattice is first set to K. The SOM is then trained on the sq over several iterations. Once the training phase is complete, the closest follower vector to each An analytics-based heuristic decomposition of a bilevel… neuron is selected, resulting in a set kq . Finally, the leader vectors kq corresponding to the kq are selected, resulting in a clustered set.

Numerical examples
To illustrate and evaluate our approach, we use two example problems. Monte Carlo simulation and clustering were done in Java and R (using the CLARA package Maechler et al. (2017)), respectively. The CPLEX 12.6 solver was also used on a 3.0 GHz Intel Xeon Processor with 8 GB of RAM.

A benchmark problem
The first problem considered, a two-follower problem, is Example 2 from (Bard 1988), and is used simply to illustrate our approach. This toy problem can be solved optimally with complete methods, but it is given here simply to demonstrate our technique (which is meant to be used in much larger problems). This is a BPMSIF as the followers are strongly independent: the followers do not share each others' follower or leader variables, and the follower problem variables are not mutually constrained. The leader vector = (x 1 , x 2 , x 3 , x 4 ) is partitioned among the followers with variables (x 1 , x 2 ) occurring in follower 1 and (x 3 , x 4 ) in follower 2. The variables = (y 11 , y 12 ) and = (y 21 , y 22 ) are also computed by followers 1 and 2, respectively.
To solve this problem using the analytics-based decomposition method, denote (x 1 , x 2 ) by a vector 1 and (x 3 , x 4 ) by a vector 2 . A large number S of assignments for 1 and 2 which satisfy the bounds of the 's are generated (Marsaglia 1972;Muller 1959), and denoted by s1 and s2 ( s = 1 … S ), respectively. For each 1 in s1 the corresponding follower problem is solved as an ILP, obtaining assignments s1 ; similarly for s2 . Next, the s1 vectors are clustered using k-medoids to get the most diverse set of assignments k1 , ( k = 1 … K ). The k1 vectors that correspond to the k1 are then selected. The same is done for s2 to obtain k2 along with its corresponding k2 . Using this decomposition, problem (3) can now be rewritten as a standard optimisation problem: (3) maxF( , 1 , 2 ) = (200 − y 11 − y 21 )(y 11 + y 21 ) + (160 − y 12 − y 22 )(y 12 + y 22 ) s.t.
In experiments, as K increases better solutions were found, with the highest value of 6594.05 obtained when K = 160 giving = (8.13, 3.80, 11.23, 16.82) , 1 = (0.74, 11.20) and 2 = (28.04, 0.00) (rounded to 2 decimal places). The clustering time when K = 160 is 234.53 s. The solution is 0.09% less than optimal, but the strength of our approach is in its ability to handle large-scale problems, as demonstrated next.

A large-scale problem
In this experiment, a problem with arbitrarily many followers is evaluated. The problem is also evaluated for the optimistic case in which the followers' solutions lead to the best objective function value for the leader.
where ∑ q q q = ∑ q ∑ n a qn x qn , and are the variables controlled by the leader and followers, respectively, and Q is the total number of followers. Both the and are vectors of real numbers. The leader variables are partitioned among the followers such that each follower contains one q each, and each q is of size n. Each component of the vector x qn is constrained to be ≤ a given upper bound x max qn . q , q , q , q and q are vectors of constants.

An analytics-based heuristic decomposition of a bilevel…
where M is a sufficiently large constant.

Evaluation
The values used for the problem are N = 6 , x min qn = 0 , x max qn = 10 , y max qn = 10 , ( ∀q, n ). . The number of followers Q was varied between 10 and 1000, and the problem was solved using both our analytics-based heuristic decomposition approach (using S = 1000 , K = 30 for each follower) and two genetic algorithms, and the results are shown in Figs. 4 and 5. Genetic algorithms are frequently used in bilevel optimisation, so this evaluation looks at the performance of this approach for the example problem. The first genetic algorithm is the Nested Bilevel Evolutionary Algorithm (N-BLEA) used in  and has been well-used for solving bilevel problems. The second is the Multi-Follower Genetic Algorithm (MFGA) described in Algorithm 1 and specifically designed for this problem.

Fig. 4
Comparing approaches: objectives and timings for Q = 10 to 100 1 3 An analytics-based heuristic decomposition of a bilevel…

Fig. 5
Comparing approaches: objectives and timings for Q = 100 to 1000

N-BLEA parameters
In order to select the parameters to use, the problem with 100 followers ( Q = 100 ) was first solved while varying some algorithm parameters. The number of parents and number of offspring ( = ) were varied from 3-8. For each of these values, the number of generations (maxGens) was also varied from 50-200 in steps of 50. This operation was run 10 times for each value of , and maxGens, and the average objective value was recorded.
It was seen that the following settings produced the best solutions on average: = = 8 , number of generations maxGens = 150 , tournamentSize = 5 , number of random individuals added to pool r = 2 , crossoverProbability = 0.9 and mutationProbability = 0.1 . The constraint handling method used by the algorithm is given in Deb (2000), and the variance-based termination criteria were set to 0.000001.

MFGA parameters
These were also varied using 100 followers. The population size popSize was varied from 30-90, and the maximum number of generations maxGens from 50-500. The MFGA parameters selected were therefore: maxGens = 500 , popSize = 50 . This population size was selected because, although there is little difference between its objective value and the best objective at popSize = 80 , the difference in time taken is almost 50% less. Uniform crossover with a crossover rate of 0.5 ( 50% ) was used. Other parameters 1 3 An analytics-based heuristic decomposition of a bilevel… are elitePercentage = 0.20 , tournamentSize = 5 , mutationRate = 0.015 and fitnessFunction = ∑ Q q ∑ N n a qn x qn + ∑ Q q ∑ N n b qn y qn .

Comparing all 3 approaches
For both N-BLEA and MFGA, each problem size was solved 10 times, and the average objective values and solution times were recorded. It should be noted that the poor performance of N-BLEA is due to the operation of its crossover operator which is additive in nature, and frequently violates the bounds of the vectors. This crossover operator results in offspring which are frequently infeasible, and are thus heavily penalised by the constraint handling scheme. MFGA was designed to avoid this problem: since vector generation is done using Hypersphere Point Picking with the appropriate boundaries, it always produces feasible offspring. For 10-100 followers, the solution found by the MFGA was better in 7 out of 10 of the cases, though our approach finds a close solution in a fraction of the time (Fig. 4). However, as the problems get larger, (Q from 100 to 1000) our approach is much better in terms of both the solution quality and the runtime (Fig. 5), especially as Q gets larger. This demonstrates the scalability of our approach. Reduction in a very large set of potential solutions to a much smaller (but highly representative) set using medoids allows the ILP to choose the best solution from a vast number of possibilities.

A large-scale bilevel cutting stock problem
In this section we first explain the real application of forestry harvesting. Then, we provide its bilevel reformulation and we evaluate it.

A multiple stock size cutting stock problem
In the classic CSP, all stock items have the same known dimensions, which makes it easier to compute the possible cutting patterns. In certain problems however, stock items come in several different dimensions and these types of problems are known as Multiple Stock Size Cutting Stock Problems (MSSCSP) (Wäscher et al. 2007).
An example of such a problem is the forest harvesting problem, as the trees differ in size from each other, sometimes significantly. In this problem, a forest is subdivided into areas called "blocks", with each block having a number of trees to harvest. This partitioning is illustrated in Fig. 6.
There are Q blocks each with market value V q ( q = 1 … Q ). Each block has R trees, with tree r ( r = 1 … R ) in block q having dimensions qr . Each tree can be cut into L different log types, with each log type having a particular monetary value.
This problem was classified in Climent et al. (2016) as * ∕V∕D∕R (see Dyckoff's typology Dyckhoff (1990)), where * means any dimensionality, V means that the total amount of items in stock (i.e. the total number of trees) is sufficient to accommodate all the demanded products (hence, only some of the stock will be cut), D means that all large pieces (stock items) are different (in terms of shape) and R indicates many products of few different types are demanded. The feature V (any demand can be fulfilled) requires that the stock items to be cut need to first be selected. Using the more recent typology presented in Wäscher et al. (2007), the problem can be classified as a variant of the MSSCSP in which there is a heterogeneous assortment of large pieces.
In practice, a semi-autonomous harvesting machine (e.g. Komatsu forest machines) cuts a tree into logs in order to maximise its total value using an algorithm A , which is typically based on a Dynamic Programming (DP) algorithm. DP is an approach that enables the solution of complex problems by dividing them into a collection of simpler sub-problems (see Anderson et al. 2015; Bellman and Dreyfus 2015 for example). The sub-problems must be sequential and independent, and the problem of cutting a tree stem satisfies these properties since it is recursive (maximise by cutting the first product and then maximising the cutting of the remainder). Let L be the length of a section of a tree stem measured from the base of the stem, and L max be the total usable length of the stem. If y k is the length of a short log of type k cut at a distance L − y k from the base of the stem, and (r(y k , L) is the associated product value, this recursive relationship can be represented as where f (0) = 0 and 0 ≤ L ≤ L max (Eng et al. 1986).
The semi-autonomous nature of the harvesting machines is due to their construction, as they are hard-coded to produce log types with the highest possible monetary value wherever possible. Thus, we only have indirect control over the cutting (6) f (L) = maximum k,y k ∈Y(L) (r(y k , L) + f (L − y k )) of trees via a set of L continuous variables called a weight vector ∈ ℕ + L . Each w l represents the weight (usually the price in € ∕m 3 ) associated with product type l.
Blocks are sold wholesale, i.e. either a block's trees are completely cut or none of its trees are cut. For each block q, a set of product types L , a vector of tree dimensions qr in the block, and a weight vector q is passed to A . The result is a product vector q ∈ ℕ + L showing the total amount of each log type obtained from the block whose trees are cut, denoted by q = ⟨a 1 , … , a L ⟩ where a l represents the volume in m 3 of units of log type l, (l = 1 … L) obtained. Consequently, A can be represented by the following mapping function for block q: Solving A with fixed q assigns values to the q variables, and the value obtained from the piece is q ⋅ q . In some applications actual market prices are used as the weights. For this problem however, the q are manipulated to obtain the desired product yields. This is the only control we have over how tree stems are cut due to the hard-coding of the harvesting machines.
Given a demand vector denoting the desired yield D l for each product type l ∈ L , the problem is to determine which blocks should be selected for harvesting, as well as the weight vector to use for each such block, in order to meet the demands while minimising the total value of the harvested blocks to conserve natural resources. There are no restrictions on which subsets of blocks can be chosen, and the trees in a block are either all cut, or none of them are.
There are several approaches that aim to achieve the desired yield D l in the literature. The authors in Kivinen (2004) and (Kivinen 2006) use genetic algorithms to try to improve the fit between the yields obtained by the harvester and the demand with varying levels of success. In Divvela and Sinha (2012), the authors use a priceweighted apportionment degree (AD) index to try to improve the fit between output and demand. This approach still prioritises logs with higher value and may not fulfil demand, leading to overproduction of unwanted logs and consequently waste. The paper (Malinen and Palander 2004) uses flexible variations on the AD to improve the fit between demand and supply; however their approach is not guaranteed to be optimal. In Marshall et al. (2006), the authors provide three mathematical models for bucking to order using a small set of market prices, targeting certain cutting patterns, and using the AD index, respectively. The paper (Kivinen et al. 2005) compares four different measures to determine the similarity between the demand and output log distributions. None of the four are shown to be superior, even though they can be used in the industry to some extent. The authors in Dems et al. (2015Dems et al. ( , 2017 use the priority list approach where higher value log types are prioritised. This approach also only considers a few cutting patterns which are assumed to be sufficient, although this is not always the case. The paper (Sessions et al. 1989) adjusts the price iteratively, but using only a small set of prices. In Duffner (1980), the authors also vary price but how they do this is not stated.
The analytics-based heuristic decomposition approach used in this paper is a good fit for this problem since a much larger number of prices can be evaluated, thus creating a good approximation of the distribution relating the prices to the products.
Also, separating the harvester operation ( A ) from the rest of the linear program using bilevel formulation allows for the more efficient solution of problems with a large number of blocks. Additionally, the use of analytics approaches presents a new way of solving a real-world bilevel problem.

Bilevel reformulation
The above problem can be naturally modelled as a multiple-follower bilevel optimisation problem. This reformulation of the forest harvesting problem as a bilevel problem is novel, and is one of the contributions of this paper. Here, the leader's objective is to select a set of blocks to harvest to fulfil demand, while each follower q seeks to harvest its block to get the optimal product vector q , given a price input q . Define binary variables h q = 1 if block q cuts its raw stock of trees, and product vectors of integers qr to describe the product yields from raw r in block q's stock. V q is the monetary value of block q estimated by the forest providers. qr are the dimensions of an uncut tree stem r in block q. The problem is thus: This is a nonlinear, mixed-integer, bilevel optimisation problem with multiple followers which we call the Bilevel Cutting Stock Problem with Multiple stock sizes (BCSPM). It is also large: there might be hundreds of blocks and hundreds of (sampled) trees per block, hence tens of thousands or more follower problems (since A is evaluated for each r), as well as a large number of product types. It cannot be solved by classical bilevel methods but it could be tackled by evolutionary methods. Metaheuristic approaches (popularly used in industry) such as (Murphy et al. 2004) and (Dueck and Scheuer 1990) have been tried, with very poor results obtained.
The model above does not have the strong independence property because all the follower problems corresponding to a block share the same variables. It can however be transformed so that it does, by grouping each block's follower problems into a single problem via new vectors of integer variables q , which model the total yield from each block: where each qr (q = 1, … , Q) solves: An analytics-based heuristic decomposition of a bilevel… Now the followers are strongly independent: each uses a unique set of variables q , q and none of the follower variables are mutually constrained. The decomposition method detailed in Sect. 4 can now be applied.
For each q a number of feasible solutions ( sq ) are generated. Each follower problem is then solved for the sq using the cutting simulator A , to get corresponding sq . Next, k-medoids clustering is applied for each follower q, resulting in the selection of a diverse set of assignments kq , together with the corresponding kq . The problem can now be formulated as an ILP: where h q = 1 indicates that all block q's trees are cut, and x qk = 1 indicates that they are cut using weights with index k. If block q is not selected then h q = 0 which forces x qk = 0 for k = 1 … K.

Evaluation of the bilevel cutting stock problem
To empirically study the performance of our approach, real data from an industrial partner was used. A smaller evaluation using 8 blocks (Q = 8) with each block's trees partitioned into a maximum of 4 different types of products was done in Prestwich et al. (2015). Results showed that the approach came to within 0.4% of optimality bound (see Fig. 7). It was also seen that the total clustering time increased exponentially with k.
A larger evaluation was carried out on a different forest with 263 blocks ( Q = 263 ) and the stock partitioned into a maximum of 11 product types. The total volume of raw material was 6149.781 m 3 , with the data obtained from the industrial partner. 1000 random weight vectors sq , ( S = 1000 ) were generated for each block, giving a total of 263,000 cutting pattern vectors sq . The total time for generating these was where each q (q = 1, … , Q) solves: approximately 18 h and 17 min. The total clustering time for all blocks with k = 125 was 31 h. 12 different instances of random demands were solved. In 6 of these instances, the demand for product types was in the range [0, 300] m 3 (low demand), while the remaining 6 had demands in the range [300, 600] m 3 (high demand). ILP solution times were nominal, taking less than 5 s for all instances evaluated. Due to the high clustering times, self-organising maps (SOM) were used an alternative to k-medoids. The SOM experiments were done using the Java Kohonen Neural Network Library (JKNNL) Rybarski and Habdank-Wojewódzki (2006). For the small problem ( Q = 8 ), the SOM was trained on the 80,000 cutting pattern vectors A chart comparing the increase in clustering times for both k-medoids and SOM approaches is shown in Fig. 8.
In terms of scalability, it can be seen in Fig. 8 that using the SOM approach offers a significant improvement in data reduction time, which makes it more useful than k-medoids for large problems. In Fig. 9, when k = N = 125 for Q = 8 , the total clustering times are 788.29 and 45.51 s for the k-medoids and SOM, respectively. When Q = 263 , the difference in clustering and map training times is even more striking and is better visualised using a logarithmic scale (Fig. 10). Here, the total times taken are 111898.35 s (31 days) and 226.72 s (4 min) for the k-medoids and SOM, respectively.

Conclusions
In this paper, a novel analytics-based heuristic decomposition approach for new class of bilevel multiple-follower problems is presented. Our approach can be used to solve large-scale multiple-follower bilevel problems more efficiently than standard approaches, as complete approaches are not practicable due to excessive computational times.
Two numerical examples were solved using this approach, and the results compared to those obtained by using evolutionary algorithms (which is a standard approach for large-scale bilevel problems). For the first example, a toy problem is solved for demonstration purposes to within 0.09% of the optimal. This shows that even for small-scale problems, the analytics approach is competitive as it is able to cover the space of the follower problems adequately.
The second example was an arbitrarily large problem evaluated for up to a thousand followers. The results were compared with those from two evolutionary approaches and it was seen that as the size of the problem increased, the heuristic decomposition approach produced significantly better results than the standard approaches. This shows that the decomposition approach is much more scalable as the number of followers increases, in terms of both runtime and solution quality.
Besides, a large-scale MSSCSP with applications in the forestry industry was evaluated. This problem was reformulated as a large-scale multiple-follower bilevel problem and solved using the heuristic decomposition method. Evaluation on a small-scale problem showed that up to a point, increasing the number of clusters got the solution to about 0.4% from the optimal. To reduce the clustering time required for large problems, self-organising maps were used as an alternative to k-medoids clustering with significant speed-up seen.
A possible future line of work would be the extension of the presented analyticsbased heuristic decomposition method to non-independent followers. Furthermore, our approach could be applied to other complex industrial applications that can be modelled as multi-follower bilevel problems