Solution procedures for block selection and sequencing in flat-bedded potash underground mines

Phosphates, and especially potash, play an essential role in the increase in crop yields. Potash is mined in Germany in underground mines using a conventional drill-and-blast technique. The most commercially valuable mineral contained in potash is the potassium chloride that is separated from the potash in aboveground processing plants. The processing plants perform economically best if the amount of potassium contained in the output is equal to a specific value, the so-called optimal operating point. Therefore, quality-oriented extraction plays a decisive role in reducing processing costs. In this paper, we mathematically formulate a block selection and sequencing problem with a quality-oriented objective function that aims at an even extraction of potash regarding the potassium content. We, thereby, have to observe some precedence relations, maximum and minimum limits of the output, and a quality tolerance range within a given planning horizon. We model the problem as a mixed-integer nonlinear program which is then linearized. We show that our problem is NP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {NP}}$$\end{document}-hard in the strong sense with the result that a MILP-solver cannot find feasible solutions for the most challenging problem instances at hand. Accordingly, we develop a problem-specific constructive heuristic that finds feasible solutions for each of our test instances. A comprehensive experimental performance analysis shows that a sophisticated combination of the proposed heuristic with the mathematical program improves the feasible solutions achieved by the heuristic, on average, by 92.5%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$92.5\%$$\end{document}.


Introduction
This paper considers one of the biggest German potash mines. In Germany, the potash ores are generally found in deep deposits. Hence, potash mines are typically underground mines with an area-wide expansion of the deposit, a so-called flat-bedded deposit. According to Musingwini (2016), optimization in underground mine planning is not as well developed and widely applied as in open pit mine planning, although the logic of the planning is the same. O' Sullivan et al. (2015) state by comparing the common mathematical formulations for both open-pit and underground mining that scheduling in underground mines is more complex than the scheduling in the open-pit mines based on the complex structure of precedence constraints, the characteristics of the operations and activities, and the irregularity of the size and shape of the blocks. Musingwini indicates that the main reason for the complexity difference between open-pit and underground mine planning is that the direction of mining in open-pit mines is essentially down and outward to the pit limits. However, in underground mines, there are numerous permutations of the direction of mining depending on the mining method chosen.
The mostly applied extraction method for flat-bedded deposits of limited thickness is a drill-and-blast technique using the room-and-pillar mining method. By the use of the room-and-pillar mining method, the material is extracted across a horizontal plane, and pillars, arranged in regular patterns, are left for support purposes. Thus, a grid-like structure is formed, as demonstrated in Fig. 1 (Hamrin 2001; Schulze et al. 2016).
By employing a conventional drill-and-blast technique, the mining activities are conducted at the salt faces (cf. Fig. 1). At each salt face, the following discrete steps (mining operations) must be processed in chronological order (K+S 2013): 1. scaling the mine roof and sidewalls, 2. removing the scaled material, 3. bolting the roof with anchors, 4. drilling large diameter boreholes, 5. removing drilling dust, Fig. 1 Grid structure caused by the room-and-pillar mining method. Reprinted from Schulze et al. (2016) 1 3 Solution procedures for block selection and sequencing in… 6. drilling blast holes, 7. filling blast holes with explosive substances, 8. blasting, 9. transporting broken material to a feeder breaker.
During mining operation 4 (see the enumerated list above), large drill jumbos drill three adjacent horizontal boreholes with a diameter of 0.28 m and a length of 7 m at each salt face. The large boreholes act in particular as a direction guideline for blasting. After the detonation (mining operation 8), chambers, so-called rooms, are created in the direction of the mining activity (cf. Fig. 1). The height and the width of the exploded area are then scaled based on a given plan. Thus, after each detonation, a three-dimensional block of potash is removed, and a new room of the same size arises. During mining operation 9, the raw material is delivered using loaders from each salt face to a feeder breaker, where the lumps are broken. After completing mining operations 1 to 9, the position of the current salt face is shifted by the respective block's length, exposing a new salt face, which can be operated on directly afterward.
The onward transportation of crude materials from feeder breakers takes place using a conveyor belt system to a bunker near the shaft. The excavated crushed potash is transported from the bunker through the shaft to the surface. The potash ores are rich in potassium chloride, sodium chloride, and different associated minerals. Potassium chloride is a valuable salt that is mainly used as a fertilizer. Furthermore, it is an integral additive in the chemical, medical as well as human and animal food-processing industry (Chesworth 2008;USGS 2011;Schulze et al. 2016). After transporting the crude salt to the surface, potassium chloride is separated from the extracted potash by flotation, recrystallization, or electrostatic separation in aboveground processing plants. For each technical device, there is an optimal operating point at which the device has the best performance. This point can be determined based on various factors. The processing plants above ground can be most cost-effectively utilized if the amount of potassium contained in the extracted material is equal to a specific value. The equivalent content of potassium oxide is often reported to indicate the percentage of potassium by weight in the potassium chloride, where 100% potassium chloride is precisely equal to 63.17% potassium oxide (cf. Heinz and von der Osten 1982, p. 147). As mentioned, at a salt face, a block of potash can be unearthed. For each block, the amount of potash is measured according to the dimensions of the excavation. Moreover, the potassium contained in each block of potash in percent is determined based on geological investigations. The percentage of potassium contained in the extracted potash from a block is defined as the quality value of the corresponding block. Accordingly, the quality value of a block multiplied by the amount of potash obtained from that block determines the amount of potassium contained. In general, the blocks are different regarding the amount and the quality value of the potash contained. Moreover, not all the available salt faces and the corresponding blocks can be mined within a given planning horizon. Therefore, the quality of the amount of potash extracted within different time intervals 1 3 strongly depends on the selected blocks and the sequence of the extraction. Vast fluctuations in the quality value result in non-homogeneous output that leads to high costs for aboveground processing. Because quality fluctuations occur frequently owing to the way in which the potash is extracted, quality-oriented mining of blocks plays a decisive role in reducing the processing costs. Newman et al. (2010) classify the existing approaches and models in mining companies according to the planning horizon or the hierarchy level into long-term (strategic), medium-term (tactical), and medium-and short-term (tactical-operative) problems. In this paper, we consider a medium-term planning horizon. Within the planning horizon, we want to have a homogeneous output regarding the quality value of the extracted potash. More precisely, the quality value of the extracted material should deviate as little as possible from a given quality target value so that the mineral processing above ground is conducted economically. In this regard, we want to answer two questions: 1. which block should be excavated (i.e., block selection), and 2. if a block is extracted, the time at which it is extracted (i.e., block sequencing). Taken together, we solve a block selection and sequencing problem with a qualityoriented objective function at a tactical planning level. For our problem, precedence relations, maximum and minimum limits of the output, and a quality tolerance range have to be taken into account.
The remainder of the paper is organized as follows: In Sect. 2, the characterization of the problem at hand and a literature review are given. Sect. 3 introduces a mathematical formulation for our problem. Since some constraints are not linear, more decision variables and constraints are introduced to linearize the mathematical program. Subsequently, we show that our problem is NP-hard in the strong sense. Accordingly, in Sect. 4, we devise a constructive heuristic, which is embedded in a multi-start environment and provides good, feasible solutions for the problem through a sophisticated time-saving procedure. In Sect. 5, based on some generated realistic problem instances, we compare the introduced mathematical program and the proposed constructive heuristic. We also show if we solve our problem heuristically and use the feasible solution found as an initial solution for the mathematical program, we obtain much better results, especially for large and challenging problem instances. The paper concludes with a summary of the achieved results and an outlook on future research in Sect. 6.

Problem specification and related literature
For a proper operation, from a geographical point of view, underground mines are usually divided into smaller spatial areas, so-called mining districts. Accordingly, an underground potash mine has, on average, up to 5 mining districts. The area of a mining district may be several square kilometers. Due to this spatial expansion, several tipple areas are constructed for a mining district to divide the area into smaller parts avoiding long transportation routes. In each mining district, depending on its area, 4 to 6 tipple areas are involved. In a tipple area, a feeder breaker is installed, where the lumps of the extracted potash in that tipple area are delivered to and crushed. As mentioned, the mining operations are conducted at a salt face.
Subsequently, the detonation occurs in the mining direction, in which a block of potash is extracted. After mining operation 8 (blasting), a three-dimensional block of potash with the known dimensions is removed, and a room is created. A chain of consecutive blocks that are extracted one after another in a certain mining direction is defined as an underground location. Since a tipple area can have an extent of up to a few 100 meters, several underground locations (usually between 5 and 11) are assigned to the single tipple areas. That means the blocks of material extracted from each underground location are transported to the feeder breaker installed in the corresponding tipple area. Figure 2 illustrates a tipple area with three underground locations (UL 1 to UL 3) that are assigned to a feeder breaker, i.e., the potash unearthed from UL 1 to UL 3 is transported to the illustrated feeder breaker. In UL 1, three blocks are already removed in the mining direction. The associated rooms are designated by squares with solid lines. After mining a block, a new salt face (a potential block) in the mining direction becomes available. Geological sampling and mining investigations determine how many of the consecutive blocks in an underground location can be removed within a considered planning horizon. In Fig. 2 within the planning horizon. For better clarity, we summarize in Table 1 the aforementioned mining terms. The excavation of a block can only then be started if the previous blocks in the same underground location have been fully excavated. A block has, at most, one direct predecessor and, at most, one direct successor block in the corresponding underground location. The time needed to remove a block is the sum of the processing times of the required mining operations 1 to 9. Each mining operation must be processed by one machine and by one skilled worker. According to the speed of the assigned machine and the skill level of the assigned worker, different processing times are needed to accomplish a mining operation. Machines and workers are scheduled at the beginning of each work shift at an operational planning level (see, e.g., Schulze and Zimmermann 2017;Seifi et al. 2019, andSeifi et al. 2020). Since we deal with a block sequencing problem at a tactical planning level, the processing times required for the extraction of blocks must be estimated. In doing so, the dimensions or shapes of blocks are crucial factors. Moreover, the current status of a block at the beginning of the planning horizon affects the needed processing time. For example, a block for which mining operations 3-9 must be carried out has a shorter processing time than a block for which all mining operations 1-9 must still be processed.
In every tipple area, the extracted material is initially carried out via the loaders from underground locations to the assigned feeder breaker. The conveyance of the extracted potash from the particular tipple areas of each mining district takes place through a conveyor belt system to a central bunker system close to the shaft, from where the material ultimately reaches the surface. The capacities of the conveyor system, the bunker, and the processing plants above ground are limited. Hence, in each work shift, an upper limit of the output for each tipple area and each mining district must be observed. On the other hand, the primary task of mining companies is the extraction of raw minerals. Accordingly, a lower Underground location A chain of consecutive blocks that can be removed in the mining direction limit for the total output over the planning horizon must be considered. Moreover, a minimum output in each work shift and every mining district must be satisfied. In our problem, we consider a planning horizon (e.g., a month) that is a union of some smaller sub-intervals (e.g., weeks). Within those sub-intervals, the quality value of the entire extracted potash must be within a given quality tolerance range. The assumptions regarding the length of the planning horizon and the corresponding sub-intervals can be customized according to the current situation of the mine at hand. Assume that a set of blocks is excavated within a specific time interval. The quality value of the entire extracted material within that time interval is the weighted average of the quality values of the removed blocks. Due to our objective, the quality value of the extracted material should deviate as little as possible from a given quality target value. The quality target value in percent typically represents the optimal operating point of the processing plant above ground. We speak of a negative deviation if the weighted average of the quality of the extracted material is less than the quality target value. Analogously, there is a positive deviation if the weighted average of the quality of the extracted material is greater than the quality target value. For formulating the objective function, we first determine the absolute deviation's value of the weighted average of the output quality from the predetermined quality target value in each considered sub-interval in the planning horizon. The value of the weighted average of the output quality and, thus, the values of negative and positive deviations are determined based on the amount of material extracted within the considered time interval. Since the amount of material removed is not known a priori, determining the deviations from the given target value requires some nonlinear constraints in the mathematical formulation that must be linearized. The aim is then to minimize the average of the calculated deviations over the number of subintervals considered in the planning horizon.
Altogether, we minimize a quality-oriented objective function observing the following groups of constraints: Precedence relations between the blocks in an underground location; Minimum limit of the output over the entire planning horizon, in each work shift, and for every mining district; Maximum limit of the output for each tipple area and every mining district within every single work shift; and Quality tolerance range over each certain sub-interval in the planning horizon. Newman et al. (2010), Kozan and Liu (2011), as well as Blom et al. (2019) published survey articles on the application of operations research methods in the field of mining. Lately, Leite et al. (2020) gave a review on state-of-the-art applications of operational research techniques to mining problems taking the mentioned surveys into account. Leite et al. (2020) consider (1) layout and design problems, (2) production and scheduling problems, and (3) operational equipment allocation problems at strategic, tactical, and operational planning levels, respectively; consequently, they review the published articles in both open-pit and underground mines.
For more convenience, in Table 2, we list the most significant works published in the previous decade that introduce a mathematical formulation for a block scheduling problem in the field of mining. Under columns "Constraints," we observe the three types of constraints we deal with in the problem at hand (Quality, Quantity, and Precedence). The "OF" column specifies whether the objective function is quantity-("Q") or monetary-related ("M").
Regarding our quality-related objective function, we first have to calculate the weighted average of the output quality. Assume that a set of blocks B is available to be extracted. Let A b and Q b be the amount of material (in tonnes) and the  Solution procedures for block selection and sequencing in… quality value (in percent) of block b ∈ B , respectively. Moreover, let x b be the binary decision variable that is 1 if block b is excavated. Then, the quality value of the entire extracted material is the weighted average of the quality values of the removed blocks, denoted by q . The value of q (in % ) is calculated as follows: In the literature, the average proportion of the most valuable mineral contained in the ore is indicated as "ore grade." Martinez and Newman (2011) formulate a mixed-integer linear program to minimize the deviations from a given target demand for three different ore grades in LKAB's Kiruna iron ore mine. However, the target demand for a particular ore grade is given in tonnes. Accordingly, expressed in our notation, they only consider the linear term (the amount of extracted iron in tonnes) to measure the deviations. Thus, we categorized the proposed objective function as quantity-related. Azzamouri et al. (2018) minimize the deviations from the demand for a certain grade of the extracted material, where the objective function and the related constraints are formulated quantity-oriented, too. Analogously, Clausen (2013) determines the amount of potassium oxide contained in the extracted potash in tonnes to minimize the deviations from a given target value in an underground potash mine.
In the literature, the quality-related constraints are mostly known as "grade control constraints" or "grade blending constraints." Those constraints ensure that q is within a permitted tolerance range [Q − , Q + ] . Thus, the following inequalities must apply: If we multiply both sides of the above inequalities with the denominator of q , we obtain linear constraints (see, e.g., Rivera Letelier et al. 2020). Except for Montiel and Dimitrakopoulos (2015) and Blom et al. (2016), all the quality-related constraints observed in the papers from Table 2 are linear grade control constraints. Montiel and Dimitrakopoulos (2015) maximize discounted profits in an open-pit copper mine. They penalize the deviations regarding mining, processing, transportation, and blending targets and consider a penalty cost in the objective function. In the problem they consider, multiple material types are sent to the available processes or to stockpiles where they are blended to meet the quality requirements. The grade of the material handled in a given period corresponds to the grade of the stockpiles. Montiel and Dimitrakopoulos emphasize in their work that the corresponding blending constraint is nonlinear; consequently, the authors propose a risk-based heuristic approach to tackle the problem without suggesting any linear mathematical formulation. Blom et al. (2016) consider a multiple mine, multiple time-period, open-pit production scheduling problem. The authors define the productivity of a mine in terms of the desirable utilization of dig and trucking resources, i.e., dig and trucking resources should be fully utilized. In each period, ore produced at each mine is transported by rail to a set of ports and blended into products for shipping. The objective function minimizes the deviation between the composition of port products and desired bounds and maximizes the productivity achieved at each mine. The observed objective function is denoted by "Other" in Table 2. The authors propose a nonlinear mathematical program to ensure blending constraints at each port while generating schedules for each mine. To tackle the problem, Blom et al. suggest a decomposition-based algorithm that can find high-quality solutions.
Our literature review suggests that no one, to the best of our knowledge, has introduced a linearization of the nonlinear quality-related constraints taking a qualityoriented objective function into account. It is straightforward to show that the scheduling problems observing upper and lower limits for the operational resources are NP-hard. The computational time required for solving NP-hard problems may be affected by the numerical parameters of the input data. A NP-hard problem in the strong sense is still NP-hard even when all of its numerical parameters are bounded by a polynomial in the length of the input. The contributions of this paper are: • to introduce a linear mathematical program for the problem described; • to show that the problem at hand is NP-hard in the strong sense; and • to introduce a solution approach that provides high-quality solutions for realistically sized problem instances.
In the next section, we first introduce the original mathematical program in its nonlinear structure and then linearize the corresponding nonlinear constraints. Lastly, we show that our problem is NP-hard in the strong sense.

Mathematical model
In this section, we introduce a linear mathematical program for the block selection and sequencing problem described in Sect. 2. From an operational point of view, we consider max sub-intervals in the given planning horizon. Each sub-interval consists of several periods, where each period represents one work shift in the corresponding sub-interval. Each work shift is represented by time interval (t − 1, t] . Let T max denote the number of work shifts in the planning horizon under consideration. Thus, a planning horizon of T max work shifts is a union of time intervals (t − 1, t] for t = 1, 2, … , T max and the point in time 0. Our mathematical formulation is based on the discrete completion times for the extraction of the blocks that are selected and mined in the considered planning horizon. Accordingly, we only focus on the discrete points in time that represent the end times of work shifts in the planning horizon. Note that a completion time at point 0 is not relevant since the point in time 0 is not the end time of any work shift in the planning horizon under consideration. Let T be the set of positive, discrete points in time in the entire planning horizon. Moreover, let T 1 , T 2 , … , T max be the disjoint subsets of T , where ⋃ ∈{1,…, max } T = T . The elements of T are the positive, discrete points in time that represent the end times of the work shifts contained in sub-interval . In Fig. 3, a planning horizon Solution procedures for block selection and sequencing in… with T max = 12 work shifts and max = 3 sub-intervals is depicted. The whole planning horizon is a union of disjoint time intervals [0, 4], (4, 6], and (6, 12]. Since the point in time 0 is not the end time of any work shift in T , we have Table 3 demonstrates the sets, parameters, and decision variables used in the mathematical program. Note that the auxiliary decision variable q is used for better clarity in the mathematical formulation and is calculated by definition as follows: The mathematical program for our problem is formulated as follows: The planning horizon and the associated sub-intervals Set of blocks assigned to tipple area k ∈ K r B r Set of blocks assigned to mining district r ∈ R K r Set of tipple areas in mining district r ∈ R Positive continuous decision variable that substitutes the product of decision variables x bt and + q Weighted average of the quality value of the extracted blocks during sub-interval 1 3

Solution procedures for block selection and sequencing in…
We minimize the average of the positive and negative deviations of the quality of the output from a given quality target value over the predetermined sub-intervals (cf. objective function (1)). A block can be excavated at most once within the entire planning horizon (constraint set (2)). By definition of decision variable x bt , if x bt = 1 , the point in time t represents the completion time of block b. Constraint set (3) determines the completion times of blocks. Note that if a block is not excavated, the completion time is 0. Constraint set (4) guarantees that the completion time of a block must be greater than or equal to its processing time, i.e., the mining of a block must start during the considered planning horizon. Constraint set (4) is active only if a block is excavated. In the corresponding big-M formulation, we can choose M equal to T max . Constraint sets (5) and (6) observe a precedence relation (5)). Moreover, constraint set (6) ensures that the completion time of block l must be at least by Z l time units greater than the completion time of block b if block l is ever processed. Constraint sets (7) to (11) ensure the minimum and maximum limits of the output. Constraint sets (12) and (13) guarantee that the weighted average of the quality value of the excavated blocks during every sub-interval ( q ) is within a permitted tolerance range. Those constraints are like "grade control constraints" or "grade blending constraints" explained in Sect. 2.
Constraint sets (14) and (15) determine the lower bounds for + and − , respectively. Since objective function (1) must be minimized, + and − take the value of the positive and negative deviations of the quality of the output from Q . However, the constraint sets are inherently nonlinear. We write the mathematical formula of q in the inequalities and reduce the left-hand side of each inequality to a common denominator. If we multiply both sides of the inequalities by the common denominator, we obtain the following nonlinear inequalities: In order to formulate a mixed-integer linear program, we introduce positive continuous decision variables + bt and − bt to substitute the products x bt ⋅ + and x bt ⋅ − on the right-hand sides of the above inequalities, respectively. We, therefore, have the following constraint sets: By substituting a product of a binary decision variable and a continuous decision variable, the upper bound for the substitution variable (here + bt and − bt ) must be determined. The upper bound is related to the maximum value that the continuous decision variable can take if the binary decision variable is 1. Furthermore, the substitution variable takes the value 0 if the binary decision variable is 0. Constraint sets (14-2) and (15-2) guarantee that decision variables + bt and − bt take the value of zero for all t ∈ T if block b is not excavated within sub-interval . Note that each block can be mined only once within the entire planning horizon. Thus, it is sufficient if we consider the summation of decision variables + bt ( − bt ) and x bt over the points in time t ∈ T . Otherwise, if a block is removed at any point in time t ∈ T , the left-hand side of constraint set (14-2) (constraint set (15-2)) can at most have the value of Q + − Q ( Q − Q − ): Solution procedures for block selection and sequencing in… Note that the maximum values of the positive deviation ( + ) and the negative deviation ( − ) are bounded by Q + − Q and Q − Q − , respectively. Subsequently, the following constraint sets determine the values of + and − that are used in (1): If we replace constraint sets (14), as well as (15) by constraint sets (14-1), (14-2), and (14-3) as well as constraint sets (15-1), (15-2), and (15-3), respectively, we have a mixed-integer linear program, where constraint set (16) indicates that decision variables x bt are binary, and in addition to non-negativity constraint sets (17) and (18), the following non-negativity constraint set must be considered: We denote the proposed mixed-integer linear program for our block selection and sequencing problem with BSSP.
In what follows, we show that our block selection and sequencing problem is NP -hard in the strong sense. We introduce a restricted case of BSSP (RBSSP) to which a pseudo-polynomial transformation from the well-known 3-PARTITION problem can be easily constructed.
We consider an underground mine that has one mining district, over a planning horizon of T max = T work shifts with max = 1 . Let the number of blocks be n = 3T . We assume that there is only one block in each underground location ( V = � ). Let the processing time of all blocks be equal to 1 work shift. Hence, constraint sets (3) to (6) do not have to be observed. Moreover, we assume that the maximum output for each tipple area is a very large number with the result that constraint set (11) is always satisfied. Since there is only one mining district in the restricted problem, constraint sets (7) and (9) are the same. By assuming the minimum output for the entire planning horizon equal to ∑ t∈T P − t , constraint sets (7) and (9) are redundant. Furthermore, let the quality value of all blocks be equal to Q with Q − < Q < Q + . Thus, the deviation from the quality target value is a constant number regardless of which blocks have been excavated. Consequently, all of the quality-related constraint sets are met. We assume ∕T . Hence, we can denote P + 1t = P − t with P. Consequently, the restricted problem, RBSSP, can be formulated as follows: (2)) (cf. constraint set (8)) (cf. constraint set (10)) RBSSP is not an optimization, but a so-called feasibility problem, i.e., a specific decision problem. In this restricted problem, we intend to determine T subsets B 1 , B 2 , … , B T of blocks that are removed at points in time 1, 2, … , T , respectively, subject to the constraint sets of the problem. Garey and Johnson (1979) showed that 3-PARTITION is NP-complete in the strong sense.
• 3-PARTITION Input: a finite set U = {u 1 , u 2 , … , u 3m } , a bound X ∈ ℕ , and a size s(u b ) ∈ ℕ for each b = 1, … , 3m , such that X 4 < s(u b ) < X 2 for each b and ∑ 3m b=1 s(u b ) = mX. Question: are there m disjoint subsets U 1 , U 2 , … , U m of U such that: If we consider each element u b of the given set U in 3-PARTITION as a block b ∈ B in RBSSP, a transformation from an arbitrary instance of 3-PARTITION to an instance of RBSSP is given by T = m , A b = s(u b ) , and P = X . This transformation

Min. Constant number
can be performed in time polynomially in the input length. The length of the constructed RBSSP-instance is polynomially related to the length of the given 3-PAR-TITION instance. Furthermore, the largest number in the constructed instance is equal to the largest number of the given 3-PARTITION instance; consequently, the conditions of a pseudo-polynomial transformation are met. In a solution of RBSSP, we have T subsets of B . Those subsets play the same role as the sets U 1 , … , U m in the desired partition of U in the 3-PARTITION. As a result, a solution for RBSSP exists if and only if the desired partition exists for the given 3-PARTITION instance. Thus, RBSSP is NP-complete in the strong sense, and optimization problem BSSP is NP-hard in the strong sense.
For NP-hard problems in the strong sense, an optimal solution, especially for large and challenging problem instances, cannot generally be found using a MILPsolver within a reasonable amount of time. In the next section, we propose a heuristic approach that solely seeks a subset of the feasible region using some rules to provide good, feasible solutions.

Heuristic solution procedure
In this section, we introduce a constructive heuristic that is embedded in a so-called multi-start algorithm. Constructive heuristics gradually generate a complete solution based on a partial solution. In our heuristic algorithm, at each point in time t, eligible blocks are determined according to two different factors. On the one hand, a block is eligible if it can be completed at the considered point in time according to its processing time and the completion status of its predecessors. On the other hand, a block is not eligible if its extraction results in an overrun of the upper limit of the output in the associated tipple area and the related mining district. Based on a specific priority rule, the elements of the eligible set are prioritized. Then, a roulettewheel selection procedure is applied to randomly select a block that is scheduled next in the planning horizon (its completion time is set to t). After that, the status of the mine and, accordingly, the eligible set are updated. The selection procedure continues until the eligible set at the considered point in time is empty. By the use of the roulette-wheel selection procedure, the next block to be scheduled is randomly selected from the eligible set. That means, if the method is carried out several times in a multi-start environment, it is very likely that a large number of feasible solutions are generated. Table 4 outlines the sets, parameters, and variables used in our heuristic algorithm.
Algorithm 1 describes the developed multi-start heuristic in detail.
In our minimization problem, all of the decision variables are non-negative. Thus, the objective function cannot take a value smaller than 0. With the prescribed value of > 0 , we denote a tiny quality tolerance value of the production process. Accordingly, a feasible solution is called an optimal solution if the associated value of the objective function lies in the narrow interval [0, ] . For a given problem instance, feasible solutions are generated using priority rule Ψ until the objective value of a feasible solution is within the predefined interval or a given time limit is exceeded (while-loop at line 2). Within an initialization step (line 3 in Alg. 1), we store for each block b ∈ B the associated mining district R b and tipple area K b . Moreover, we store the blocks with no predecessor block in set E . For all b ∈ E , the earliest start time is 0 ( es b = 0 ). Furthermore, the residual capacities res rt and res kt are set to the corresponding maximum limits of the output P + rt and P + kt , respectively. After the initialization step, we pass the entire planning horizon in a for-loop (line 4). For the current point in time t, we determine the related sub-interval and put those blocks from E into set E T that can be completed until t (if-condition at line 7). At the beginning of the planning horizon ( t = 1 ), the earliest start times of all blocks b ∈ E are 0. Accordingly, only the blocks b ∈ E with processing times Z b = 1 can be added to E T . If no block can be inserted into E T , the output at the point in time t is 0. Hence, constraint set (8) is violated, and no feasible solution can be found. Thus, the algorithm terminates (line 10). Otherwise, blocks from E T are reconsidered according to the residual capacities of the corresponding mining districts and tipple areas. For block b ∈ E T , if A b does not exceed the associated residual capacities (if-condition at line 13), b is inserted into E CT . In a realistic problem instance, we always have: A tiny quality tolerance value of the production process New value of q if block b is mined next in sub-interval in % Therefore, at least one block from E T (if there is any) can always be inserted into E CT at each point in time. The block that must be scheduled next is selected from E CT according to the given priority rule Ψ (line 16). For our heuristic, we consider four different priority rules. Based on priority rule Ψ ∈ {1, 2, 3, 4} , we determine for each block b in E CT a priority value b . In the following, we explain the way how b is calculated according to a specific priority rule.
For this rule, we put emphasis on the quality value of a block. A block with a smaller deviation from the quality target value takes a larger priority value. This priority rule gives a block that may contribute to a better objective function value a greater probability to be extracted next. Since some blocks may have no deviation from the quality target value, we add a tiny number > 0 to the denominator of the fraction. In our work, we set = 0.0001.

priority rule 2:
For priority rule 2, we additionally consider the amount of material extracted from a block. If quality values of blocks are the same, a block with a larger amount of material is given a larger priority value, or-in other words-for the same amount of material, a block that has a smaller deviation from the quality target receives a larger priority value (like priority rule 1). Using priority rule 2, the lower limits of output are also considered to avoid generating infeasible solutions. priority rule 3: Here, we calculate the value of q b that represents the change of q if block b is excavated next. A block takes a larger priority value if its excavation results in a smaller deviation from the quality target value. We can, therefore, give a more considerable priority value to the blocks that allow the best possible improvement in the objective function value in each step.
Priority rule 4 helps to avoid infeasible solutions in terms of constraint set (9). For block b, if the lower limit of the output for the corresponding mining district R b is achieved, we set the priority value of block b equal to a tiny number > 0 . Otherwise, block b receives a priority value equal to A b . Using priority rule 4, we focus only on finding feasible solutions regarding the lower limits of output.
According to the priority values, all blocks b ∈ E CT receive a probability value b as follows: 1 3

Solution procedures for block selection and sequencing in…
We apply a roulette-wheel selection procedure, where each block occupies an area on the roulette-wheel proportional to its b -value (Michalewicz and Fogel 2004, Sect. 6.1). That is equivalent to partitioning the interval [0, 1] into |E CT | parts, where the b-th sub-interval (part) has the width b representing block b. Subsequently, a random number between 0 and 1 is generated. The sub-interval that contains the random number determines block b * . At lines 17 and 18, the associated residual capacities and the value of q are updated. If b * has a successor l (line 19), we insert l in E . The earliest start time of l is set equal to the completion time of b * , i.e., es l = t (line 21). Block b * is then removed from E CT . Moreover, all blocks b ∈ E CT with have to be removed from E CT (lines 22 and 23). The repeat-until-loop in Alg. 1 will be executed until E CT is empty.
After violating the condition of the repeat-until-loop, t is incremented by one. The entire for-loop at line 4 will be executed for all t ∈ T . We then check the feasibility of the solution found according to constraint sets (7), (8), (9), (12), and (13). Feasible solutions are stored, and infeasible solutions are discarded.
We observe the loops in Alg. 1 to determine the time complexity of the presented heuristic approach in every run. The for-loop at line 4 is executed T max times. Moreover, the for-loop at line 6 and the repeat-until-loop are conducted, at most, |B| times. Line 11, the for-loop at line 12, as well as line 16, if-condition at line 19, and lines 22 and 23 within the repeat-until-loop have, at most, |B| steps. Accordingly, the algorithm has a time complexity of O(T max ⋅ |B| 2 ) that implies a pseudo-polynomial algorithm. Note that in practical applications, T max is given as a constant number. Therefore, the time complexity of the algorithm is O(|B| 2 ) , and the heuristic proposed is polynomial.
In the next section, we compare the results achieved by a MILP-solver with the solutions provided by the proposed heuristic approach. Moreover, we introduce a sophisticated combination of the heuristic and the mathematical program to improve the results obtained by the heuristic procedure for the most challenging problem instances, for which our MILP-solver cannot find a feasible solution within a reasonable amount of time.

Computational study
In this section, we perform an experimental performance analysis for the proposed mixed-integer linear program and the constructive heuristic. The computational study is executed on randomly generated problem instances that are based on realworld data derived from Clausen (2013). Table 5 shows some parameters that are typical for a potash mine and used to generate realistic problem instances.
To normalize the problem instances, we set the number of blocks in each mining district r equal to 325. Thus, the problem instances with more underground locations have fewer blocks in each underground location and vice versa.
The allowed maximum output in each mining district r for time interval (t − 1, t] depends on the number of loaders available in the corresponding mining district during the associated time interval. A loader can typically transport 750 tonnes of crude material within a work shift. In each work shift, 1 to 4 loader(s) are available so that P + rt ∈ {750, 1500, 2250, 3000} . The minimum output in each mining district depends on the capacity of loaders, too. The parameter P − r can be determined as a certain percentage of the capacity of all the available loaders in mining district r within the entire time horizon: The parameter varies between 70% and 90% , with a step size of 5% . Furthermore, the lower limit of the total output is the sum of P − r over the mining districts contained in the underground mine: The upper limits of the output for tipple areas P + kt are randomly chosen from the set {1000, 2000, … , 5000} . Those numbers are given due to the characteristics of typical feeder breakers. Even though the parameters P + kt have different values in terms of tipple areas k, they are the same for all time intervals (t − 1, t] in the planning horizon. In each time interval (t − 1, t] , an output greater than 0 has to be achieved. We set the parameter P − t as follows: We consider a planning horizon of one month that without loss of generality, is supposed to have only four weeks. Each sub-interval represents a time interval of one Number of underground locations in each tipple area 5-11 Number of blocks in each underground location 5-10 Amount of material contained in each block in tonnes A b 700-1200 Quality value of each block in % Q b 9.8-16.2 Lower limit of the quality value of the output in % Q − 11.1 Upper limit of the quality value of the output in % Q + 14.1 Quality target value for each sub-interval in % Q 12.6 Processing time of the first blocks in each underground location measured in work shifts Processing time of other (not the first) blocks in each underground location measured in work shifts 6-12 1 3 Solution procedures for block selection and sequencing in… week ( max = 4 ). A week contains 18 work shifts, so that the planning horizon has T max = 72 work shifts.
The quality values of blocks are stated in percentage and usually have decimals. We multiply the quality values of blocks and, accordingly, the quality parameters from Table 5 by 100 to avoid rounding errors. For example, a quality value of 13.7% is 1370 in a problem instance. As mentioned in Sect. 4, from a practical point of view, a solution is called an optimal solution if the associated value of the objective function lies within the small interval of [0, ] with > 0 . Parameter is the quality tolerance value of the production process and is chosen to be 0.1% . Like the quality values of blocks, we multiply by 100 as well; is, therefore, equal to 10.
We distinguish between 5 different test sets according to the number of mining districts. Hence, there is 1 mining district in the first test set, 2 mining districts in the second test set, and so on (5 mining districts in the fifth test set). We randomly generated 20 problem instances for each test set using the introduced parameter ranges. Additionally, we consider five different levels of , i.e., ∈ {70%, 75%, … , 90%} , concerning the lower limit of the output ( P − r and accordingly P − , see above). Thus, we have 100 problem instances for each -level. Test sets with 5 mining districts represent the largest problem instances. On the other hand, the problem instances with a -level of 90% are the most challenging problem instances to solve. We used GAMS 24.9 to implement our mixed-integer linear program and solved the problem instances using CPLEX 12.7.1. We set the solver parameters in the way that the solution procedure terminates in the following cases: 1. if a solution found is optimal; 2. if a solution found is in the predefined interval [0, 10]; or 3. if a time limit of 1800 s is exceeded.
The corresponding solution procedure is called CPLEX. According to the cases mentioned above, if we say that CPLEX finds an optimal solution for a problem instance, case 1. or case 2. holds.
The heuristic algorithm is implemented in the programming language ++ and executed with the compiler Microsoft Visual Studio 2010. Since we use four different priority rules in the heuristic, we set the upper limit of the solution time equal to 450 s for each priority rule. Thus, for each problem instance, we run the multi-start heuristic approach for 1800 ( 4 × 450 ) s to have a relatively fair comparison with the results of CPLEX, which has a time limit of 1800 s, too. The resulting multi-start heuristic approach with four randomized priority rules is called MSH-4R.
All tests are executed on an Intel(R) i7-7700K@4.20GHz machine with 64 GB RAM under Windows 10.
We can compare the results achieved by CPLEX and MSH-4R from different aspects. On the one hand, the number of problem instances for which no feasible solution can be provided is important. On the other hand, we want to know whether a feasible solution obtained is optimal, and if not, what can be said about the solution quality. We discussed that the problem instances with a -level of 90% are the most challenging problem instances to solve, and the ones with 5 mining districts are the largest problem instances. Accordingly, the results of each solution procedure for each -level and regarding the size of the problem instances are of most interest. Our computational study suggests the following results: • In general, the greater the problem instances, the more time CPLEX needs to find an optimal solution. • For -levels of 70% to 85% , CPLEX finds an optimal solution for almost all problem instances. However, for the most challenging problem instances with = 90% , CPLEX can hardly find any feasible solution. • MSH-4R can find for all of the problem instances at least one feasible solution that is not far from the best solutions found. • MSH-4R using priority rule 3 provides the most high-quality solutions in comparison to the other priority rules. • If we use the solutions found by MSH-4R using priority rule 3 as initial solutions for CPLEX, the results found for the most challenging problem instances are improved by, on average, 92.5%.
A detailed explanation of the results achieved by the solution approaches is given in the following. Table 6 compares the results achieved by CPLEX with the solutions found by MSH-4R. Each row in Table 6 shows the information for a test set that consists of 20 problem instances. Columns # Optimal depict the number of test instances for which CPLEX and MSH-4R could find an optimal solution (see the explanation for an optimal solution above). The number of test instances for which a feasible solution could be found but the optimality of the solution could not be proven or shown is given under columns # Feasible for both solution procedures. Columns # Unknown present the number of test instances for which no feasible solution is found within the considered time limit (1800 s). Note that the numbers appearing under each column cannot be greater than 20. Columns Solution time report the average solution time in seconds for the considered 20 problem instances in each test set. According to the setting, the average solution time cannot be greater than 1800 s. Analogous to CPLEX, MSH-4R terminates if a solution found by a given priority rule is in the predefined interval [0,10]. In that case, we say that MSH-4R finds an optimal solution for a problem instance. Note that MSH-4R is not able to state whether a feasible solution with an objective function value greater than 10 is optimal or not. For a problem instance solved by MSH-4R, an optimal solution may be found by using any priority rule. In those cases, the procedure terminates before the time limit (450 s) is exceeded for the respective priority rule. The time in seconds under Solution time MSH-4R is the sum of solution times taken by every priority rule and cannot be greater than 1800 s ( 4 × 450).
We see in Table 6 that CPLEX could find for the whole 100 problem instances with = 70% an optimal solution. The greater the number of mining districts in a test set, the more time CPLEX needs to find an optimal solution. For = 70% , the solution time is, on average, 14 s for the 20 problem instances with 1 mining district and 132 s for problem instances with 5 mining districts. For the whole 100 problem instances with = 70% , CPLEX needed, on average, 53 s to prove the optimality of the solutions found. For = 75% , there is the same trend. All of the test instances are solved to optimality by CPLEX; the solution time is, on average, 26 s for problem instances with 1 mining district and 330 s for problem instances with 5 mining districts. Notably, the computational times are higher (for the 100 problem instances with = 75% , the computational time is, on average, 187.6 s) since the problem instances are more challenging. For -levels of 80% and 85% , the same trend can be seen. For = 80% ( = 85% ), CPLEX could optimally solve 98 (90) problem instances. At the -level of 80% ( 85% ), the solution time is, on average, 202 (324) s for the 20 problem instances with 1 mining district and 676 (1020) s for the 20 problem instances with 5 mining districts. For the other problem instances with = 80% and = 85% that could not be solved  to optimality, CPLEX could find a feasible solution (cf. columns # Feasible). For the most challenging problem instances with = 90% , CPLEX could find for only seven problem instances with 1 mining district an optimal solution within, on average, 1241 s. Furthermore, for only two other problem instances, a feasible solution could be found by CPLEX. In other words, for 91 problem instances at the -level of 90% , CPLEX cannot find any feasible solutions within a reasonable amount of time.
On the contrary, MSH-4R was able to find an optimal solution for 12 problem instances each at the -levels of 70% , 75% , 80% , and 85% as well as for nine problem instances with = 90% (altogether 57 of the whole 500 problem instances). Hence, it cannot be said at which level MSH-4R works best. MSH-4R could find for every problem instance at every level and with every size at least one feasible solution (the numbers under # Unknown MSH-4R are all 0), even for the -level of 90% . To evaluate the solution quality of the results achieved by MSH-4R, we calculate the average of the absolute deviations from the best solutions found (in % 100 ), which is stated in the last columns Ave. deviation from best of Table 6. The following example illustrates how the deviation from the best solution found is calculated.
Example 1 Let 6.64 be the objective function value of the best solution found by CPLEX ( S C ), and 24.34 the objective function value of the best solution found by MSH-4R ( S H ). Since S C ≤ = 10 is true, we consider S C as an optimal solution. Accordingly, the deviation of S C from the best solution found is 0, and the deviation of S H from an optimal solution is S H − = 14.34 % 100 . Now, let S C = 16.34 and S H = 30.71 . The best solution found is S C and, consequently, the deviation of S C from the best solution found is 0. Independent of whether we can prove the optimality of S C or not, the deviation of S H from the best solution found is S H − S C = 14.37 % 100 .
For the problem instances with = 70% , the solutions found by MSH-4R deviate, on average, 11.30 % 100 ( 0.113% ) from the optimal solution found by CPLEX. The minimum value is for the test set with 1 mining district (7.42 % 100 ), and the maximum value belongs to the test set with 3 mining districts (14.60 % 100 ). For = 70% , MSH-4R could find for only 12 problem instances an optimal solution. Hereby, seven problem instances that are solved to optimality by MSH-4R have 1 mining district. On the other hand, for the test set with 3 mining districts, no optimal solution is found by MSH-4R. Since the average is calculated over the whole 20 problem instances for each test set, the differences regarding the average deviations from the best solutions found can be explained by the different number of optimal solutions found for each test set. The same trend exists at -levels of 75% , 80% , and 85% , where the solutions found by MSH-4R deviate, on average, 11.3, 11.12, and 10.45 % 100 from the best solution found, respectively. We can conclude that the quality of the solutions found by MSH-4R is quite promising for the -levels of 70% , 75% , 80% , and 85% , where optimal solutions for 388 of 400 problem instances are known, which is a strong measure to evaluate the results. In addition, we see that for = 90% , where CPLEX could find for only nine problem instances a feasible solution, MSH-4R found for all of the 100 problem instances at least one feasible solution. Therefore, it is reasonable to use MSH-4R to solve the problem instances at the -level of 90%.
In the next step, we compare the applied priority rules in MSH-4R (cf. Sect. 4) in Table 7. The first columns on the left-hand side of Table 7 show the number of problem instances for which a specific priority rule exclusively found the best solution. For example, in the first row of Table 7, the sum of the numbers under Number of best-known solutions exclusively found for rules 1 to 4 is 14. That means there are six problem instances in the test set for which the best solution could be found by means of at least two different priority rules. At the same line, we see that, e.g., priority rule 2 could exclusively find the best solutions for two problem instances, which means that the other priority rules could not find the best solution for those problem instances. We see in Table 7 that for 407 of the whole 500 problem instances, the best solutions were exclusively obtained using priority rule 3,  which suggests priority rule 3 as the most potent rule. Remember that priority rule 3 enables the best possible improvement regarding the objective function in each step. Let t bs and t fs be the times in seconds that are taken to find the best (bs) and the first feasible solution (fs), respectively. Those times are given on the right-hand side of Table 7. MSH-4R using priority rules 1, 2, 3, and 4 needed, on average, 225, 202, 213, and 220 s, respectively, to find the best feasible solutions. The first feasible solution could be quickly found by MSH-4R regardless of the applied priority rules (e.g., priority rules 2 and 4 find the first feasible solution by no longer than, on average, 0.05 s).
We saw that CPLEX performs for = 90% worse, in particular, with respect to the number of feasible solutions found (cf. Table 6). At the -level of 90% , where we have the most challenging problem instances, and where it is tough to find a feasible solution by a MILP-solver, MSH-4R can find at least one feasible solution for each problem instance. In the following, we introduce another solution approach to tackle the problem instances with = 90% . The idea is to support CPLEX by starting with an initial feasible solution that is found by MSH-4R. We distinguish between the following kinds of initial solutions: first solution (fs): MSH-4R using priority rules 2 and 4 is the fastest approach to find a feasible solution (cf. Table 7). That may be the case because both priority rules consider the amount of material of a block for determining the priority values to observe the lower limits of output. If we observe the quality of the solutions found, MSH-4R using priority rule 2 (priority rule 4) could exclusively find the best solution for 18 problem instances (one problem instance) with = 90% . Remember that priority rule 2 considers both block's amount of material and the block's quality deviation from the quality target. By contrast, priority rule 4 focuses only on the amount of material, which justifies the higher quality of the solutions found using priority rule 2 than using priority rule 4. Thus, we solve the problem instances at the -level of 90% by MSH-4R using priority rule 2 and store the time t fs . Then, we give the first feasible solution found as an initial solution to CPLEX and set the upper limit of time equal to 1800 − t fs seconds. The corresponding solution approach is called CPLEX with fs. best solution (bs): Table 7 shows that MSH-4R using priority rule 3 performs best among all four priority rules. MSH-4R using priority rule 3 could exclusively find for 67 problem instances with = 90% the best solution within, on average, 221.2 s. Hence, we solve the problem instances at the -level of 90% by MSH-4R using priority rule 3, set the upper limit of the solution time equal to 250 s, and store the best feasible solution found. Then, we solve the problem instances by CPLEX using the best feasible solutions found by MSH-4R as an initial solution with a time limit of 1550 s. The corresponding solution approach is called CPLEX with bs. Table 8 depicts that using an initial solution leads to a much better performance of CPLEX at the -level of 90% . Without using an initial solution, CPLEX could find for only nine of 100 problem instances a feasible solution (thereof seven optimal solutions) within, on average, 1688 s (cf. Table 6, last 5 rows). If we apply CPLEX with fs, 42 problem instances can be solved to optimality within, on average, 1414.2 s. Typically, the larger the problem instances, the harder it is to solve. For the problem instances with 1 mining district, CPLEX with fs could find an optimal solution for 12 of 20 problem instances. For the problem instances with 5 mining districts, only five of 20 problem instances can be solved to optimality by CPLEX with fs. CPLEX with bs performs even better, where 45 problem instances can be optimally solved within, on average, 1403 s (some more problem instances within some less time in comparison to CPLEX with fs). The same trend regarding the size of the instances is recognizable. We can see that 14 of 20 problem instances with 1 mining district and five of 20 problem instances with 5 mining districts can be solved to optimality by CPLEX with bs.
Finally, we compare the results achieved by MSH-4R with the combination of MSH-4R and CPLEX. For this purpose, the average of the absolute deviations from the best solution known is depicted in Table 9 for each test set with = 90% . Let S H i be the objective function value of the best solution found by MSH-4R for problem instance i. Moreover, let S fs C i and S bs C i denote the objective function values of the solutions found by CPLEX with fs and bs for problem instance i, respectively. The numbers in parentheses are the average of (S H i − S fs C i )∕(S H i ) as well as (S H i − S bs C i )∕(S H i ) over the problem instances in every corresponding test set. We see that the solutions found by CPLEX with an initial solution are much better than the solutions found by MSH-4R. In particular, if we first solve a problem instance heuristically using priority rule 3 and then give the best feasible solution found as an initial solution to CPLEX (CPLEX with bs), the solutions are improved by, on average, 92.5%.

Conclusion
In this paper, we introduced a mixed-integer linear program and a heuristic approach to solve a block selection and sequencing problem occurring in underground potash mines. In the problem under consideration, we minimized the deviations of the output quality value from a prescribed quality target value. The deviations were calculated for certain sub-intervals within a given planning horizon.
To evaluate the solution approaches, we randomly generated 100 problem instances (5 test sets of 20 problem instances each) based on realistic data. We solved the problem instances on 5 different levels in terms of the lower limit of the output using the proposed mixed-integer linear program and the suggested constructive heuristic. We can conclude that if the MILP-solver starts with an initial solution, an optimal solution for a problem instance can be found more quickly. We can also say, if we solve a problem instance within a reasonable amount of time heuristically using introduced priority rule 3 and then give the best feasible solution found as an input to the MILP-solver, we obtain quite promising results. Hence, practice-relevant problems in potash mines can be solved with an acceptable quality within an adequate time frame. Consequently, the results achieved could support the decisions of underground mining operators to provide the aboveground processing plants with homogenous output by a systematic comparison between the actual and target performance.
Further research will apply metaheuristics to improve the results achieved by the proposed constructive heuristic. Uncertainty regarding the processing times needed to excavate the blocks must also be investigated. On the other hand, the material excavated in a work shift can be stored in a bunker during a specific time interval. That can lead to another value of the output quality in comparison to having all the material conveyed directly to the surface. Ongoing research can distinguish those cases to make the generated solutions more realistic.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.