Bi-objective optimization of the tactical allocation of job types to machines: mathematical modeling, theoretical analysis, and numerical tests

We introduce a tactical resource allocation model for a large aerospace engine system manufacturer aimed at long-term production planning. Our model identiﬁes the routings a product takes through the factory, and which machines should be qualiﬁed for a balanced resource loading, to reduce product lead times. We prove some important mathematical properties of the model that are used to develop a heuristic providing a good initial feasible solution. We propose a tailored approach for our class of problems combining two well-known criterion space search algorithms, the bi-directional (cid:2) -constraint method and the augmented weighted Tchebycheﬀ method. A computational investigation comparing solution times for several solution methods is presented for 60 numerical instances.


Introduction
Allocating resources among competing activities is one of the most popular type of optimization problems considered in various application areas, such as manufacturing, hospitals, and finance.Manufacturing, Planning, and Control (MPC) is a subject of interest to many practitioners as well as operations research experts due to its positive impact on efficiency for many companies globally.Some of the traditional MPC systems are integrated MPS (Master Production Schedule;Blackstone Jr., 2013) and MRP (Materials Requirement Planning; Plossl and Orlicky, 1999) models that optimize simultaneously the production and the purchase of raw materials over a shorter time horizon, taking the allocation of products to resources as fixed.It is well understood that having more options of verified resources could be useful for the planners to tackle any demand variations in the short term.In some industries, there is, however, a need to perform time-consuming or costly verification before re-allocating a product to a new resource.Hence, it is beneficial to identify such verifications (we use the term qualifications) before they are required.
Our proposed model (for case company GKN Aerospace, a leading manufacturer of aerospace engine systems) identifies the routes that a product may take through the factory over a fairly long time horizon (a quarter of a year in this work).The need of our allocation model is due to the requirement of costly as well as time and resource demanding qualification of production processes when new products are introduced or significant changes arise in the machining capacity (in this work, resources implies machines).The model considers multiple products, each having different part types that require machining operations (e.g., milling, turning, and grinding), and the capacity (in hours) of machines as well as the time horizon are considered to be finite.There is a tendency (at GKN) to utilize only a few capable machines (with respect to speed, tolerances, and multiple functionality) at high loading levels (Blackstone Jr., 2013, p. 94).This causes a resource loading imbalance in the production system that leads to variations in the delivery times of products, resulting in long queues of jobs (operations done on products/parts) in front of some machines.We relieve imbalance by identifying multiple routes by means of qualifying alternative machines for a given type of job.In summary, our model allocates jobs (operations performed on products/parts) to machines with the bi-objective goal to minimize both qualification cost and loading imbalance, while demand in each time period and various other side constraints are satisfied.

Qualifications of machines for jobs
A qualification of a machine for a job type involves verifying that all encoded tasks can be performed by the machine within acceptable tolerances.It involves writing a computer program for the control system, buying new fixtures, training the machining staff, and performing simulations to validate expected output.Once a machine is qualified for a given job type, it can be used for orders of the same job type in all subsequent time periods; hence, the associated costs are one-time costs, as opposed to set-up/start-up (changeover) costs commonly considered in capacitated lot-sizing models (Pochet and Wolsey, 2006, p. 137). 1 Table 1 illustrates the routings in a dummy production system with three machines and the two operations milling and turning, performed on a single product/part.For simplicity, we refer to job type as the combination of a product/part type and the operation performed on it.In the first time period, milling is done in machines 1 and 3, in the second time period, the same operation is done in machines 1 and 2, and in the third time period, it is done only in machine 2. For machine 2, the box around the M indicates that the milling operation is to be qualified and performed in time period 2. For machine 3, the turning operation is to be qualified in time period 2 and performed in time period 3.The blue and green colors indicate the time periods when a given machine is qualified to Table 1 Routings for a single part/product: M (milling) and T (turning) indicate time periods (t) when machines (k) are used for the respective purpose; indicates time period and machine qualification for milling and turning, respectively perform milling and turning, respectively.As a few machines are pre-qualified (used previously) for some job types, not all machines need to be qualified.For example, machine 1 is pre-qualified for both M and T.

Contribution
The contribution of this work is twofold.We introduce a tactical resource allocation problem (TRAP) with unrelated parallel machines (with several company related side constraints and assumptions), and a corresponding mathematical model.We consider the trade-off between resource loading imbalance and qualification costs, which is to the best of our knowledge not considered in the literature.Some mathematical properties of the model are utilized in a starting heuristic for reducing the computing time.We then suggest a tailored algorithm for the TRAP, starting with a bidirectional -constraint method, followed by the Augmented weighted Tchebycheff (AWT) method.
The approach tackles some well-known shortcomings of the -constraint method and is shown to shorten computing times, in particular for some difficult numerical instances.Our model can be used by most organizations where resources are shared among products/parts/customers and costly and/or time-demanding preparations are required the first time a product/part/customer needs to be (re)allocated to a new resource.

Outline
The scope of the TRAP is described in Section 2.1 along with some background, while Section 2.2 provides a brief introduction to multi-objective optimization methods.Section 3.1 presents a mathematical model for the TRAP; it is proven to be NP-hard in Section 3.2.Section 3.3 incorporates decision makers' preferences in the TRAP and mathematical properties of the model are discussed in Section 3.4; these are then utilized in a starting heuristic for reducing computation time; see Section 4. Section 5 presents criterion space search methods and the numerical tests performed.In Section 5.1, our suggested approach is motivated.Section 5.2 introduces a modified bi-directional -constraint method and implementation details.Section 6 presents a method to generate instances for our problem that closely represent the industrial setting, then followed by a comparison of our proposed approach with several state-of-the-art criterion space search methods on the generated instances.

Background and scope
The research field of production planning is broad.We provide a brief orientation of the field before diving into the specific variant of the production planning problem studied in this work.One popular way of classifying production planning models is by acknowledging the time horizons considered.This simplifies to some extent the decision variables and model parameters.Several authors (e.g., Gupta and Maranas, 1999;Min and Zhou, 2002) classify production planning problems as strategic, tactical, or operational.A strategic planning affects the design, the product structure, and the configuration of the factory.It is usually done two to five years in advance (the range varies among industries).Tactical planning involves determining material flows, inventory levels, and capacity utilization; it is usually done one to four years in advance and acts as an input to operational models, such as machine scheduling.

Tactical planning literature and scope of the TRAP
According to Díaz-Madroñero et al. (2014), there are five main categorizations of tactical planning problems: (a) number of products/items and structure of bill-of-materials/levels (e.g., series, assembly, general, and arborescence) (Pochet and Wolsey, 2006, Chapter 13); (b) distribution (stochastic or deterministic) of the demand; (c) time discretization; (d) capacity constraints; and (e) types of objective functions.Further, most industrial problems are multi-item and some are also multi-level.
Mainly two types of time discretizations are considered: short-time buckets, in which their is enough time to manufacture one part/item; long-time buckets, in which multiple parts or an entire product is produced.A classic example in which small time buckets are used is the discrete lot-sizing and scheduling problem (DLSP; Lasdon and Terjung, 1971), in which the integration of the lot-sizing and scheduling problems is enhanced by the short time buckets.Hybrid models refer to combining large and small time buckets; a classic example of a hybrid model is the general lot-sizing and scheduling problem (GLSP; Fleischmann and Meyr, 1997).Transchel et al. (2011) consider a hybrid model using macro (discrete) and micro (continuous event-based) time scales.Although in most companies resources are shared by multiple products, parallel machines are considered in few articles (Wörbelauer et al., 2018; for parallel production lines).It is computationally hard to perform both scheduling and lot-sizing in a monolithic model for industrial instances.Furthermore, most companies prepare the schedule of a job only when a customer order has been created in the MRP system that is generally few hours/days (varying between companies) before the delivery time/date is fixed; so-called rolling horizon.Hence, long-term scheduling is generally not practically useful as modeling all the uncertainties in advance is a hard task.Therefore, there is merit in separating lot-sizing and scheduling problems.Some articles pertain to bi-objective lot sizing problems, for example, Ammar et al. (2019) (objectives inventory cost and set-up time) and Rezaei and Davoodi (2011) (objectives service level and total cost).Many articles highlight demand uncertainty, which is often modeled using stochastic programming (Lan et al., 2011;Nourelfath, 2011), less commonly by fuzzy sets (Chen and Huang, 2010;Lan et al., 2011), and by robust approaches (Wei et al., 2011;Genin et al., 2008).Most capacity constraints in tactical level models regard machines, manhours, fixtures, tools, and inventory levels.The most common type of objective function minimizes cost/time (processing time, set-up time, and fixture costs; Bradley and Glynn, 2002;Mieghem, 2003).The drawback of using such functions is, however, that most of the cost measures rely heavily on the used accounting principles, which are generally misleading (Myrelid and Olhager, 2019).
Our model focuses on long-range resource (machine) loading, a time frame in which reliable and detailed (weekly/daily) demand predictions do not exist.Consequently, short-time buckets are not relevant.Instead our time discretization uses long-time buckets (a quarter of a year for our case) wherein it is reasonable to assume a constant material flow (i.e., precedence between jobs can be ignored) and a deterministic demand due to various fixed contracts.The manufacturing of products involves activities such as cutting, welding, heat treatment, and quality control on either the final product or parts that are later assembled into the final product.We focus only on cutting operations, which make up the majority of the total lead time of products.We propose a capacity allocation plan that promotes a balanced resource loading enabled by new qualifications.This is done by means of a bi-objective optimization model, which is to be used when one or several products are introduced in the factory or in case of significant change in available machining capacity.

Multi-objective optimization methods
Most real world industrial applications have several, often conflicting, objectives.This has resulted in a surge in utilization of multi-objective optimization methods to obtain so-called efficient frontiers (see Ehrgott, 2006).The solutions on the efficient frontier help the decision maker understand the true trade-off between two or more objectives.Algorithms for bi-objective integer programming (BOIP) and bi-objective mixed-integer programming (BOMIP) are broadly classified into decision space search methods and criterion space search methods.Popular methods for decision space search include evolutionary multi-objective methods (e.g., NSGA-II; Deb et al., 2002), which has gained interest although it does not provide any measure of optimality or near-optimality.There are also some nonevolutionary metaheuristics, hybrid multi-objective metaheuristics, and parallel multiobjective optimization methods that have become popular in some applications (such nonstandard approaches are discussed in Talbi et al. (2012)).Such inexact methods are, however, extremely useful if the enumeration of all the nondominated points (NDPs) using exact methods is hindered by time or memory limitations.Some improvements in branch-and-bound methods for bi-objective mixed 0-1 linear optimization problems are presented in Vincent et al. (2013).
Our work focuses on criterion space search methods, motivated by the fact that they can exploit the power of integer programming solvers (Gurobi in our case), and that any future computational enhancement of the solver will be replicated for the end-user of our algorithm as well.Some popular methods for criterion space search are the weighted sum method (Aneja and Nair, 1979), the perpendicular search method (Chalmet et al., 1986), the augmented weighted Tchebycheff (AWT) method (Bowman, 1976;Steuer and Choo, 1983), and the -constraint method (Miettinen, 1988, p. 85).More recently proposed methods include the balanced box (BB) method (Boland et al., 2015a), which extends the traditional box method for BOIP (Hamacher et al., 2007), and the triangle splitting method for BOMIP (Boland et al., 2015b unsupported nondominated points (Ehrgott, 2005, Definition 8.7), the existence of which makes BOMIP (and BOIP) particularly hard to solve.Most algorithms for multi-objective integer programming (MOIP) possess one basic common operation, the so-called scalarization, the idea of which is to transform a MOIP into a sequence of single objective IPs (integer programs).Scalarization techniques for transforming a BOIP into a sequence of IPs include the augmentation method, which uses appropriate values for the objective weights (Mavrotas, 2009;Özlen and Azizoglu, 2009), and lexicographic minimization of the two objectives (Ehrgott, 2005, Section 5.1).Three important factors to consider when selecting an appropriate method are (i) the maximum number of scalarized problems to solve in order to identify all the nondominated points (NDPs), (ii) the quality of the feasible solutions (if) available for each scalarized problem, and (iii) the computational cost of a given scalarized problem.A combination of these factors along with information about the problem type and instance sizes help in identifying a functioning approach.For example, if a problem is computationally hard to solve unless it is initialized with a good feasible solution, then a decomposition method such as the balanced box (BB; Boland et al., 2015a) and the box method (Hamacher et al., 2007) may provide better (i.e., shorter) solution times as compared to the -constraint method.However, to identify a set G of NDPs, the BB method solves 3|G| scalarized problems (Boland et al., 2015a, Proposition 5), while the -constraint method solves only 2|G| + 1 scalarized problems (Chankong and Haimes, 1983).Hence, most decomposition methods solve more scalarized problems, which may result in large computing times for certain problem instances.Several attempts to combine different approaches to utilize their respective strengths and avoid weaknesses have been made.For one such example, as presented by Dai and Charkhgard (2018), the algorithm starts by employing the BB method and then switches to the -constraint method; this algorithm results in 25%-40% improvement in solution times over the pure BB method and the pure -constraint method on various popular benchmarking instances (Dai and Charkhgard, 2018, Table 1).In Leitner et al. (2016), switches between the perpendicular search and -constraint methods are done repeatedly during the course of the algorithm.The contrast between these two methods-apart from using different combinations-is that Dai and Charkhgard (2018) use a fixed switch rule while Leitner et al. (2016) do not.The idea of combining two or more approaches has, however, been around for a long time (to the best our knowledge, the first application appeared in Ulungu and Teghem (1994)).We conclude that different criterion space search methods have varying effects on computational performance, depending on the problem and instances.

Problem description and properties
The tactical resource allocation problem (TRAP) is defined as follows; notations are listed in Table 2.
Definition 1 (Tactical Resource Allocation Problem (TRAP)).Given a set J of job types (tasks) and a set K of machines, let p jk be the processing time (including set-up time) of job type j ∈ J when performed in a compatible machine k ∈ K j ⊆ K.Each machine k ∈ K has the capacity C kt (time units) in time period t ∈ T and a relative loading threshold ζ k ∈ [0, 1].The demand a jt of each job type j ∈ J in time period t ∈ T must be met.The number of machines allocated to the same job type in each time period may not exceed the value of the parameter τ ∈ Z + .For assignments ( j, k), © 2022 The Authors.

International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies
Table 2 Notation for the tactical resource allocation model Sets Description J = {1, . . ., J} set of job types to be performed on the products K = {1, . . ., K} set of machines K j ⊆ K set of machines feasible for job type j ∈ J N j ⊆ K j set of machines feasible, but not qualified for job type j ∈ J T = {1, . . ., T } set of time-periods Variables Description bold notations representing vectors of the corresponding indexed variables Parameters Description cost associated with qualifying machine k ∈ N j for job type j ∈ J γ ∈ Z + upper limit on the number of qualifications in a single time period τ ∈ Z + maximum number of alternative machines for each job type in a single time period loading threshold for machine k ∈ K such that k ∈ N j and j ∈ J , so-called qualifications are required that generate additional one-time costs.It holds that N j ⊆ K j for all j ∈ J ; for the case of a new job type (associated with a new product) j, K j = N j holds.For a job type j ∈ J , the machines in the set K j \ N j do not require any qualifications (pre-qualified machines).The total number of qualifications performed per time period t may not exceed the value of the parameter γ ∈ Z + .The objectives considered are to minimize the sum (over time periods) of maximum excess resource loading above given thresholds and to minimize the qualification costs (see Section 3.1 for more details).
Let us consider the example illustrated in Table 1, with two job-types, T and M. Using the notations (see Table 2 Let us denote the solution (the variable x in Tab. 2) illustrated in Table 1 for the job type milling (M) as and for turning (T) as (with the same row (t) and column (k) labels).The excess resource loading (loading level above the threshold ζ k ) for t = 1 and k = 1 is defined as (column) then is 0.3 0 0.05 0.3 0 0 0 0 0 , which implies that for t = 1 the maximum resource loading over all machines is 0.3.Analogously computed, for t = 2 the maximum is 0.3 and for t = 3 the maximum is 0. Hence, the sum over the time periods of the maximum excess resource loading equals 0.6 and the qualification cost equals 2 (machine 2 is used for M and machine 3 is used for T, but N M = {2} and N T = {3}).In another solution, if we do not consider qualifications, the corresponding solutions , which yields the excess resource loading matrix 0.3 0 0.05 0.3 0 0 0.05 0 0 .
Hence, the sum of the maximum excess resource loading equals 0.65 and the qualification cost is 0 (machine 2 is not used for M and machine 3 is not used for T).

A bi-objective mixed integer programming model of the TRAP
The constraints defining feasible resource allocations.The production should equal the demand for each job type j in each time period t, as expressed in (1a).The constraints (1b) ensure that no job is performed in any machine and time period to which it is not allocated.The constraints (1c) limit, for each job type and time period, the number of alternative allocated machines to τ , the value of which is given as input by the user; a too small value of τ may result in an empty set of feasible solutions; a too large value may lead to an increased complexity of the product routes.The constraints (1d) ensure that the machining hours do not exceed the respective machine capacities in any time period.The constraints (1e) imply that if a job of type j is assigned to machine k ∈ N j in time period l ∈ T , then machine k must be qualified for that job type at the latest in time period l.The constraints (1f) limit the number of scheduled qualifications in each time period to γ (due to limited number of skilled professionals for completing new qualifications).The constraints (1g)-(1j) define the allowed values of the variables x jkt , s jkt , z jkt , and n t .The feasible solutions to the TRAP are thus defined as the variable vectors (x, s, n, z) fulfilling the constraints Definition 2 (Set of feasible solutions and its image in the criterion space).Defining2 y := (x, s, n, z), for any values of τ, γ ∈ Z + , the set of feasible solutions to the model ( 1) is denoted as The two objectives considered represent the preferences of the planners.The first objective is to minimize the excess resource loading, expressed as to minimize x,s,n,z where n t is limited by the constraints (1d) and (1j).It considers the sum over the time periods t ∈ T of the excess resource loading of the machines (i.e., n t ≥ 0), which is defined as the maximum (over the machines) ratio between the allocated machining hours and the available hours (i.e., 1 C kt j∈J p jk x jkt ) minus the loading threshold ζ k ∈ [0, 1] for the machine.Therefore, in a solution (x, s, n, z) that minimizes the objective g 1 , the equality n t = max{0; max k∈K { 1 C kt j∈J p jk x jkt − ζ k }} will hold for every t ∈ T .In the context of a bi-objective mixed integer programming (BOMIP), this objective is thus defined by (2a), (1d), (1g), and (1j).The constraints (1d) prevent the loading level of each machine in each time period from exceeding one, hence, imposing a capacity limit on each machine.Minimizing the excess resource loading will result in an increased capacity of buffers, which in turn enables reduced lead times for the products.
An alternative to the min-sum objective in (2a), which amounts to solving min y∈Y (τ,γ ) { t∈T n t }, would be to choose a min-max objective, that is, to solve min y∈Y (τ,γ ) {max t∈T n t }.The latter may result in more reasonable solutions, but is not chosen due to two reasons.(i) First, since there is no clear priority among the excess resource loading n t over the time periods t ∈ T , let us hypothetically consider T objective functions (n t ) defining a T -dimensional criterion space.While the min-sum objective is guaranteed to always yield an efficient solution in this T -criterion space, the min-max objective (also known as max-ordering) is not (Ehrgott, 2005, Proposition 9).The min-max objective can, however, be utilized to calculate an efficient solution by performing a lexicographic max-ordering optimization (Ehrgott, 2005, Chapter 5.3).But since the order in which the objectives are minimized is not known a priori (generally T ≥ 16 holds for all our problem instances) identifying a nonincreasing ordered sequence of the objective functions' values will significantly increase the computing time.(ii) From a practical standpoint, in some time periods the demand for products might be consistently higher, leading to higher loading levels of the machining resources.Hence, a min-max objective-focusing on high-demand time periods and ignoring the others-can be counterintuitive for production planners.
The second objective is to minimize the total qualification cost, defined as the sum of the one-time costs incurred by qualifying machines for job types, over the time periods.It is expressed as to minimize x,s,n,z Note that an increase in the number of qualifications (g 2 ) will enable a reduction of the excess loading of the machines (g 1 ).
The vector-valued function of objectives is denoted as g(y) := (g 1 (y), g 2 (y)), and the image of the set Y (τ, γ ) in the criterion space is defined as g(Y (τ, γ ) such that the relations g(y) ≤ g(ȳ) and g(y) = g(ȳ) hold.If ȳ is efficient in the decision space, then g(ȳ) is a nondominated point in the criterion space.The set of efficient solutions to the bi-objective optimization problem ( 1)-( 2) is denoted as Y eff (τ, γ ).

Computational complexity
We show that the TRAP (Definition 1) is NP-hard, by comparing with the makespan minimization problem.
Definition 4 (Makespan minimization).Given n jobs and m machines, let d jk be the processing time of job j ∈ {1, . . ., n} if it is assigned to machine k ∈ {1, . . ., m}.No machine can process two jobs at the same time.After a job has started getting processed, it must be completed on that machine without any interruption.The makespan minimization of unrelated parallel machine is the problem of assigning each job to a machine such that the maximum completion time (makespan) C max over the machines is minimized.This problem is referred to as R||C max ; see (Lawler et al., 1993, Chapter 9).The problem is NP-hard in the strong sense, since its version where d jk = d, that is, the problem denoted P||C max , is also NP-hard; Garey and Johnson (1979).
Proposition 1 (Problem reduction).The makespan minimization of the unrelated parallel machine scheduling problem, that is, R||C max , is polynomially reducible to the TRAP of Definition 1.
Proof.Letting, for j ∈ {1, . . ., n} and k ∈ {1, . . ., m}, jk = 1 if job j is assigned to machine k, and (3) Consider then an instance of TRAP (see Definition 1) given by the following: The corresponding special case of our model ( 1)-( 2), is then expressed as the single-objective MIP to minimize The instance of TRAP, expressed in (4), is an R||C max , as expressed in ( 3).An optimal solution (x * , n * 1 ) to ( 4) is equivalent to an optimal solution to this instance of the TRAP (1)-( 2) and given by (x, s, Since the R||C max decision problem is NP-complete, it follows that the TRAP is NP-hard.

Considering a third objective: routing efficient solutions
The constraints (1c) limit the number of alternative machines allocated to a job type in each time period to τ , the practical effect of which is reducing the number of possible routings (i.e., sequences in which machines are visited) used by a product.A decision maker prefers efficient solutions to the bi-objective model ( 1)-( 2) that make use of low numbers of alternative machines.Hence, we define the set of routing efficient solutions, Y reff (τ, γ ) ⊆ Y eff (τ, γ ) (see Definition 3), by means of the function g 3 (y) := t∈T j∈J k∈K j s jkt , representing the total number of alternative machines used.Assigning g 3 a lower priority than g 1 and g 2 , its inclusion will not change the nondominated points (NDPs) in the two-dimensional criterion space to the TRAP.The less importance of g 3 as compared to g 1 and g 2 is because g 3 is associated with internal transportation cost that has less priority.

Integrality of the variables z
For fixed values of the binary variables s, we define the following polyhedron in the z-variable space: Proposition 2 (Integrality of the variables z).For any s jkt ∈ {0, 1}, k ∈ N j , j ∈ J , t ∈ T , all extreme points to the polyhedron Z(s, γ ), defined in ( 6), are integral.

Property of the efficient frontier for the TRAP model
The efficient frontier of a general BOMIP has isolated points as well as closed, half-open, and open line segments; (Boland et al., 2015a, Theorem 3).For the TRAP model, however, the following proposition holds.Proof.It is sufficient to show that no line segment between any two NDPs in the criterion space exists.First, in any efficient solution, the variables (x, s, z) will take finite integer values, say (x , s , z ).Hence, from (1j) and (1d) follow that each variable n t will either take the value n t = 0 or n t = max k∈K { 1 C kt j∈J p jk x jkt − ζ k }, that is, n t takes only discrete values.Hence, the efficient solutions constitute a discrete set, which maps to a discrete set of (finite) NDPs in the criterion space.The proposition follows.

The starting heuristic
For solving the BOMIP (1)-( 2) , we use a MILP solver, the efficiency of which partly relies on an early pruning of branches owing to high-quality initial feasible solutions.We propose a simple heuristic that either quickly finds a feasible solution, returns no solution (but does not conclude infeasibility), or concludes that the problem is infeasible.A feasible solution (if found) is provided to the MILP solver while looking for a nondominated point g TOP (see Section 3.3).The heuristic relies on a decomposition of the model ( 1)-( 2) with respect to the time periods resulting in one (smaller) MILP per time period. 3Since we aim at a solution with an objective vector value as close as possible to g TOP , the heuristic employs an objective function which prioritizes minimizing the objective function g 1 , defined in (2a).The heuristic is described below and summarized in Algorithm 1.Let the vector y t := (x t , s t , n t , z t ) represent the corresponding variables for t ∈ T .The constraints (1e) connect the time periods-thus complicating the model-and are replaced by the time-disconnected constraints where the constants ξ t−1 jk ∈ {0, 1} are computed according to (12).Sequentially for t = 1, . . ., T , the time-disconnected feasible sets are then defined as 4 3 Another approach to finding a good starting feasible solution would be to use a Lagrangian heuristic, where, for example, the constraints (1e) are relaxed, resulting in the subproblem separating over the time indices t ∈ T as well as between the variables (x, s, n) and z.Each of the 2T subproblems would still be NP-hard, but considerably smaller than the TRAP.A Lagrangian heuristic does, however, require several dual iterations to find good enough dual variable values and, finally, a feasibility heuristic. 4For brevity, we define the notation (constraint reference) t , referring solely to the time index t of the corresponding constraints.

International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies
If, for t = 1, a feasible solution to the model (11) t is found, the index t is increased by one and the model (11) t is solved with ξ t−1 as input and ȳt as output.A solution obtained for a certain value of t may not be changed; hence, for any t > 1, the model (11) t may be infeasible even if the original feasible set-as defined in (1)-is nonempty.For an infeasible model (11) t , feasibility is achieved by an increase of the value of γ (temporarily, for this value of t).The feasibility loop terminates when either a feasible solution to the model (11) t is found or the constraint (1f) t is not in an irreducible inconsistent subsystem5 (IIS); see Parker and Ryan (1996, Section 2).In the latter case, the original model ( 1) is infeasible, since neither the constraints (1e) t nor (1f) t can be in IIS, which implies that one of the other time-disconnected constraints (1a) where the set Z([s t ] t∈T , γ ) is defined analogously as in ( 6).The LP ( 13) can be solved efficiently by any commercial solver; if infeasible, however, the heuristic terminates without any conclusion about the feasibility of the model (1).

Criterion space search method
According to Proposition 3, the TRAP has only isolated NDPs.With some care for the respective algorithmic properties, we may now use ideas from BOIP algorithms to find (almost) all NDPs of the TRAP (1)-(2).

Search for nondominated points by the -constraint and AWT methods: strengths and weaknesses
We briefly describe the -constraint and the AWT methods along with their main properties; below, h 1 and h 2 denote the objective functions, h(x) = (h 1 (x), h 2 (x)) , h ref is a reference point, 1 = (1, 1) , and X is an integer set.
The -constraint method yields the scalarization for a BOIP as min An optimal solution to a given scalarization is guaranteed to be weakly efficient; Chankong and Haimes (1983).Further, all efficient solutions can be found by setting appropriate values for .Since no feasible solution is available when solving a given scalarized problem, computing times may be longer as compared to the case of an available feasible solution.This issue can be resolved by a bi-directional (BD) variant of the -constraint method, denoted the BDmethod; see (Boland et al., 2015a, Section 5.1) and Algorithm 2.Then, the -constraint method is applied alternately to the two objectives, such that a previously identified efficient solution is Algorithm 2. Modified bi-directional -constraint available when solving the scalarized problem.The bi-directional approach does not require solving an infeasible problem as the last scalarized problem, cf. the original -constraint method.
The AWT method yields the scalarization for a BOIP as min The method recursively searches for a yet-unknown nondominated point by minimizing the maximum weighted (by The second term (the l 1 -norm) is added to avoid generating too many weakly nondominated points (NDPs), using a suitable value of λ ≥ 0. The AWT method can identify all the NDPs and is used in many interactive multi-objective optimization approaches (Steuer and Choo, 1983).Dächert et al. (2012) establish that too small values of λ cause numerical difficulties, while too large values may result in oversight of some NDPs; they also derive problem dependent (adaptive) formulae for calculating α 1 , α 2 , and λ; these are, however, applicable only for pure BOIPs with integer valued objective functions.Dächert et al. (2012) has shown computational superiority of their adaptive formulae (for the algorithm presented in Ralphs et al. ( 2006)) over using fixed parameters instead.Steuer and Choo (1983) present a lexicographic weighted Tchebycheff method that guarantees that all NDPs are found, without requiring input parameters α 1 , α 2 , and λ; it requires, however, the solution of two optimization problems and is computationally expensive, as empirically verified by Miettinen et al. (2006); to linearize the max-function additional variables and constraints are included, which may increase the solution times as well.
As discussed in Ehrgott (2006, Section 4.4), a drawback of the -constraint method for solving a BOIP is that the -constraint is a knapsack constraint; this may reduce the computational tractability of the scalarized problem (Ehrgott, 2006, Sections 2.2 & 3.2).In a branch-and-bound © 2022 The Authors.

International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies
tree corresponding to a given -constraint scalarization, the contradictory nature of the two objectives implies that the -constraint will be active in an optimal solution at several nodes if possible; hence, many branch-and-bound nodes need be explored before reaching a node with inactiveconstraints (see Appendix A.1). Ehrgott (2006, Section 4) suggests replacing the -constraint by a so-called elastic constraint.The idea is to allow violation of the -constraint at the expense of a weighted penalty variable added to the objective function.The approach has been successful in several applications (see, e.g., Ehrgott and Ryan, 2002 for an implementation for a crew-scheduling problem).Finding appropriate values for the coefficient of the penalty variable is, however, tricky and the computational efficiency depends hugely on the type of problem and instance.
We propose and investigate a modified version of the bi-directional -constraint method, in which the second stage uses the reference-point based AWT method, which also avoids the -constraints (an example showing one of the differences between the scalarization of the two methods is highlighted in Appendix A.1).

Implementation of the modified bi-directional -constraint method
Our proposed approach, the modified bi-directional -constraint method, combines the BDmethod with the AWT method; it is described below and summarized in Algorithm 2.
Description of the algorithm.The input to Algorithm 2 are the two initial nondominated points g TOP and g BOT , defined in (5), and which constitute the initial points in the set NDP of nondominated points.The two points g TOP and g BOT are used to initialize g 1 and g 2 , respectively.The parameter φ ∈ (0, 1) determines when to switch from the BD-method to the AWT, while the parameter ψ 1 (ψ 2 ) denotes the value by which g 2 1 (g 1 2 ) is reduced while applying the -constraint on g 1 (y) (g 2 (y)).Since the qualification cost parameter β jk is integral, we set ψ 2 := 1.We propose switching to the AWT method if the area of the rectangle B(g 1 , g 2 ) to be explored for yet-unknown nondominated points fulfills A(g 1 , g 2 ) ≤ φ • A(g TOP , g BOT ) (see Fig. 1).We made computational tests for φ ∈ {0.25, 0.35} and ψ 2 = 0.01, which are reasonable for our problem instances.We set a time limit, t lim of 5000 seconds for solving each scalarized problem, after which the solver terminates the computations (see Section 6, §1 for details about the termination criteria used).
Algorithm 2 starts by the BD-method (i.e., if , which is expressed by the scalarizations referred to as TOP and BOT; TOP and BOT apply an -constraint on the qualification cost and the excess resource loading, respectively.The model TOP is solved and the value of g 1 is updated; if the solver does not verify optimality in a user-defined time limit it terminates (see Sections 6 and 1 for criteria to verify optimality).The model BOT is solved and the value of g 2 is updated.(In Fig. 1, g 2 1 = g BOT 1 = 2.48 and ψ 1 = 0.01, while g 1 2 = g TOP 2 = 24 and ψ 2 = 1.)If the area of the remaining rectangle fulfills A(g 1 , g 2 ) ≤ φ • A(g TOP , g BOT ), Algorithm 2 switches to the AWT method, which is provided two nondominated points, g 1 and g 2 .The AWT method searches for yetunknown nondominated points in the interior of B(g 1 , g 2 ); cf.Section 5.1 The adaptive formulae used to calculate α and λ are found in (Dächert et al., 2012, Table 2), and are guaranteed to yield only efficient solutions when both objective functions take only integer values (for this purpose, one could alter n t such that it represents exact time units above a given, integer valued capacity threshold).The method AWT (Section 5.1) identifies and includes all nondominated points in the set NDP (for details, see Steuer and Choo, 1983) provided that the parameters α and λ are set appropriately.
The while-loop ends when at least one side of the rectangle B(g 1 , g 2 ) is below a certain limit (i.e., Harvesting efficient solutions.Algorithm 2 only searches for yet-unknown nondominated points in the interior of B(g TOP , g BOT ).To initialize the algorithm, the solution ( x, s, n, ẑ) from the heuristic (Algorithm 1) is provided as input to computing g TOP = g(y TOP ), as defined in (5a).The solution y TOP is then provided as input to computing g BOT = g(y BOT ), as defined in (5b).Setting y 2 := y BOT , the solution y 2 is used as a starting feasible solution for solving the model TOP and computing g 1 = g(y 1 ); then the solution y 1 is used as a starting feasible solution for solving the model BOT and computing g 2 = g(y 2 ).This procedure is repeated, providing updates of the nondominated points g 1 and g 2 and the corresponding solutions y 1 and y 2 ; the former are included in the set NDP of efficient solutions, while the latter are provided as initial feasible solution to the solver.

Computational details, tests, and results
We generate 60 instances expected to sufficiently capture most of the realizations of actual data that the model might encounter at GKN Aerospace.All the computations are performed in Python 3.7 using Gurobi 9 on a system with 1.70 GHz processor, 4 cores, and 16 GB RAM.Algorithm 2 terminates when at least one of the following criteria is met: (i) the running time for a scalarized problem exceeds 5000 seconds; (ii) the (relative) MIP duality gap is < 0.0005; (iii) MIP absolute gap limit is ≤ 0.01 and there has been no improvements in the duality gap for the previous 1000 nodes of the branch-and-bound tree. 6Note that optimality is verified (see Algorithm 2) only if the solution satisfies (ii) or (iii).

Generating the industrial test instances
We have the data for most of the parameters and sets mentioned in Table 2, however, we do not have processing times p jk of jobs j ∈ J in machines k ∈ N j (qualification required), and the qualification cost parameters β jk , j ∈ J , k ∈ N j .In order to generate instances that represent possible realizations of the actual data, we introduce the following distributions that are based on knowledge of the managers.
Skewness of processing times.The skew normal distribution is a generalized normal distribution allowing for nonzero skewness (Weisstein, 2021). 7We generate processing times p jk , k ∈ N j , j ∈ J , for newly qualified machines from three differently skewed normal distributions with mean μ and skewness/shape parameter α: positive skew (α = 1 > 0), negative skew (α = −1 < 0), and zero skew (α = 0).A location parameter/mean μ is based on the expected processing time of a given job type j on an already qualified machine k ∈ K j \ N j and which is similar to that of the machine being qualified.For all these distributions, we set the scale parameter σ := 0.1 • μ; according to the internal statistical process control data (and managerial experience), processing times of newly qualified allocations have a standard deviation of 10% of the expected value.
Qualification cost.The exact cost for qualifying a machine for a job type is not known a priori and accurate predictions require detailed simulation work by the engineering team.Thus, the input received are so-called cost levels, assigned to each qualification.For testing our model and proposed modifications, we define 20 cost levels, H = {1, . . ., 20}, and select the qualification costs from different discrete distributions over the discrete domain H. Letting π h be the frequency of cost level h ∈ H, its relative frequency is πh := ( i∈H π i ) −1 π h ; we also define π0 = 0. To determine a cost β jk , a sample α is drawn from the interval [0,1].Then, The frequency distributions are defined as follows.For each h ∈ H, 6 Criterion (iii) is particularly useful while minimizing g 1 as the corresponding lower bounds are often close to 0.

Constant data across 60 instances
The demand is from quarterly forecasts made at GKN Aerospace in January 2015 (denoted by J 15 ) and December 2015 (denoted by D 15 ) for the period 2016-2019.The minimum, maximum, and median values of the demand of job types (# of orders) are 1, 172, and 11, respectively.For processing times over the job types and machines, the corresponding values are 0.1, 89.7, and 5.63 hours, respectively.Each machine has a yearly capacity of 5000 hours; the available hours per machine and quarter are thus 1250.The planning period of four years with quarterly time buckets yields T = 16 time periods.There are K = 125 machines, and J = 517 unique job types, each having integral demand per time period.The number of possible assignments of job types to machines during the entire planning period thus amounts to ∼ 10 6 .We employ the parameter values τ = 3, γ = 4, and ζ k = 0.7, k ∈ K.The instances vary based on (i) distributions used to draw samples of qualification cost parameters (denoted by β ∈ {Left,Right,Symmetric,Uniform,Bimodal}), and processing times (denoted by p ∈ {skew+,skew-,skew0}), (ii) demand forecast (denoted by ā ∈ {J 15 , D 15 }), and (iii) whether a new product is included (with new) or not (without new).Consequently, we have 60 instances.Figure 2 shows box plots for the effect of the above mentioned variations8 on the range of g 2 (•), that is, g TOP 2 − g BOT 2 , which also sets the upper limit g TOP 2 − g BOT 2 + 1 on the number of nondominated points.

Results and insights
We compare various state-of-the-art solution approaches with our proposed approach.For each approach, the objectives are combined using augmentation (Aug) or lexicographically (Lex).The criterion space search method used is either the bi-directional -constraint (BD-) or the balanced box (BB) method.A switch to the AWT method is denoted by AWT; an ∅ means that there is no such switch.The solution approaches compared are defined by the 3-tuples (Aug,BD-,AWT) (for the values of φ ∈ {0.25, 0.35}), (Aug,BD-, ∅), (Aug,BB,∅), and (Lex,BD-,AWT), in total five variants.
Performance.Table 3 compares solution times and provides numbers of nondominated points9 |NDP| for the methods (Aug,BD-,AWT), labeled φ = 0.25, and (Aug,BD-, ∅), labeled φ = 0.Each row consists of two instances having equal distributions and forecasts, one with new product and one without.Solution times refer to total wall clock time used to find the nondominated points-excluding the computation of g TOP and g BOT . 10Comparisons for solution times should be made between columns labeled as φ = 0.25 and φ = 0, and for the same instance.A "−" indicates that the solution times for (Aug,BD-,AWT) and (Aug,BD-, ∅) are equal, since the switch to the AWT method did not occur.There are 23 instances for which the switch to the AWT method did occur.Six of these 23 instances have either moderately shorter or equal solution times for (Aug,BD-, ∅) as compared to (Aug,BD-,AWT); for the remaining instances, the (Aug,BD-,AWT) method is computationally superior.
Figure 3 compares the five variants of solution methods mentioned above, as applied to the 60 instances by means of performance profiles (see Dolan and Moré, 2002), with the performance ratio defined as r ps := t ps min r∈S {t pr } , and t ps denoting the time used by the method s ∈ S for solving instance p ∈ P. The probability that method s solves each of the instances within the ratio ω of the time spent by the fastest method for that instance (i.e., the cumulative distribution function for the performance ratio of method s) is defined as The method (Aug,BB,∅) performs quite badly for our class of problems.One explanation is that BB solves at most 3|Y eff | scalarized problems (Boland et al., 2015a, Proposition 5), while BD-solves at most 2|Y eff | + 1 scalarized problems (see Chankong and Haimes, 1983), where Y eff is the set of nondominated points.Figure 4 presents a performance profile to highlight the effect of using the starting heuristic (Algorithm 1) within the computation of g TOP , which along with the corresponding solution y TOP is the input to Algorithm 2 (i.e., (Aug,BD-,AWT) with φ = 0.25).Additionally, we compare the performance of our algorithm with NSGA-II, implemented using the python library pymoo (described in Blank and Deb, 2020).For all the instances, our algorithm resulted in larger hypervolume (as defined in Zitzler et al., 2003) than NSGA-II that implies that our algorithm yields a better approximation of the efficient frontier.The ratio of hypervolumes is presented in Fig. 5 for all the instances.Since all the nondominated points were not obtained by NSGA-II, we have not compared it in Fig. 3.For details on parameters used for NSGA-II, we provide access to a jupyter notebook in https://bit.ly/3JKE4DAand corresponding data in our github repository  x, s, n, ẑ) is provided (H) and when no such solution is provided (WH), while searching for g TOP .
https://bit.ly/3iAV4Ao.We have used the feasible solutions obtained from the starting heuristic Algorithm 1 in the initial population and terminated NSGA-II after 3600 seconds.As a quality check for the starting heuristic, we compare the resulting values for the objectives g 1 and g 2 with g TOP , which is the nondominated point possessing the minimum possible value for g 1 (see definition in Section 3.3; cf. also the objective in ( 11)). Figure 6 shows normalized (with respect to g TOP − g BOT 1 ) distances between g 1 (•) and g TOP , where (x, s, n, ẑ) denotes the solution obtained from Algorithm 1.Note that, out of the 60 instances, the normalized distance equals 0 for six instances of the objective g 1 , and for 36 instances of g 2 .3 for which the switch to AWT did occur (φ = 0.25).The number of instances in each category is specified above the corresponding box.
Variations among the instances.has an impact on the magnitude of computational superiority of switching to the AWT.However, for all the 23 instances where AWT is used, we see the ratio is always either close to 1 or significantly greater than 1.It is to be noted that in all tests performed, most of the computation time is spent on improving the lower bounds (verifying optimality).In Fig. A4, we illustrate the efficient frontier for instance 5y.

Conclusion
A new resource allocation model for a large tier-1 supplier of aerospace engine systems is presented and analyzed.We have proved various mathematical properties of the tactical resource allocation problem (TRAP) model, one of them resulting in a relaxation of the integrality constraints on some of the variables that is utilized in a specialized heuristic for the TRAP (see Algorithm 1).The solution from the heuristic is used as a starting feasible solution that resulted in significant reductions of solution times (Fig. 4).We tested a tailored modified bi-directional -constraint method, which reduces the solution time for a majority of the instances.The algorithm attempts to use the strengths of the -constraint method while avoiding some of its drawbacks by switching to the augmented weighted Tchebycheff AWT in later stages when the -constraint is not likely to be tractable.Further tests should be conducted to check the performance of our algorithm as applied to other problem classes.
Future work may involve further polyhedral analysis of the problem aiming at a tighter formulation, and investigating exact decision space search methods, such as branch-and-bound for BOIP.The model ( 1)-( 2) should be extended to a robust formulation, incorporating uncertainty in processing times and qualification costs of newly qualified jobs.To incorporate more preferences of the decision makers, new objectives and/or alternative combinations of the current objectives can be explored as well.Further, we intend to investigate how the inclusion of inventory modeling at different stages of the production planning affects the trade-off between objectives.The points g TOP = (1, 4) and g BOT = (4, 0) are known and the ideal (reference) point is g ref = (1, 0) ; see Fig. A1.For comparison of the (bi-directional) -constraint (BD-) method with the AWT method, we examine the branch-and-bound trees corresponding to their respective scalarized problems.In the BD-method, the first step is to add the -constraint x 2 ≤ 3 and solve the scalarized problem where X is the set defined by the constraints in (A1) and the coefficient 1 5 is small enough to yield at least a weakly efficient solution.In the AWT method, the first scalarized problem to be solved is defined by the reference point g ref = (1, 0) , as min where the parameters λ, α 1 , and α 2 are derived as in (Dächert et al., 2012, Table 2).The branchand-bound trees for the two scalarized problems in (A2) are shown Figs.A2 and A3. 11At the root node to (A2), x * 2 = 3, such that the -constraint x 2 ≤ 3 is active, whereas at the root node to (A3), x * 2 = 2.05.The conflicting nature of the two objectives typically keeps the -constraint active for several nodes (in our small example its just once) while solving a scalarization in the BD-method, while this is not the case for the AWT method.

A.2. The resulting efficient frontier for instance 5y
In Fig. A4, we show nondominated points identified for instance 5y.The two rectangles marked in magenta have area less than φ • A(g TOP , g BOT ), where φ = 0.25, g BOT = (2.25,17) and g TOP = (0, 28) .Hence, AWT is applied in this instance.

©
Fig.1.BD-with switch to the AWT.The green points are nondominated, the red-dashed lines represent -constraints, and the purple-dashed rectangles represent the switch to the AWT method for φ = 0.35 and 0.25, respectively.The respective latest updated values of g 1 and g 2 are defined as functions of φ.

©
2022 The Authors.International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies

7
Fig.2.Distribution of the range for 2 , which affects the maximum number of nondominated points, for different instances.

Fig. 3 .
Fig. 3. Performance profiles of the different solution methods while exploring the strict interior of B(g TOP , g BOT ).

10
Fig.4.Performance profiles when ( x, s, n, ẑ) is provided (H) and when no such solution is provided (WH), while searching for g TOP .

Fig. 6 .Fig. 7 .
Fig.6.Normalized distance of the objective values provided by the heuristic, with respect to g TOP for 60 instances.The average value the normalized distances equal 0.1 and 0.27 for g 1 and g 2 , respectively.See Section 6.1 for details on instances.
Figure 7 shows box plots for the ratio of solution times for (Aug,BD-, ∅) versus (Aug,BD-,AWT) (for φ = 0.25).The top-left and top-right plots divide the 23 instances with respect to the skewness of processing times and the distributions of qualification costs, respectively.Furthermore, the bottom-left and bottom-right plots indicate observed variations for new product's inclusion and demand forecast, respectively.It is clear that the instance type © 2022 The Authors.International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies
). Methods such as weighted sum are not able to find .1), as © 2022 The Authors.International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies y TOP := arg lexmin y∈Y(τ,γ ) (Efficient frontier for the TRAP model).The efficient frontier of the model (1)-(2) contains only isolated nondominated points(NDPs), and no (closed, half open, or open)line segments, irrespective of the values of the parameters β jk , k ∈ N j , j ∈ J .
© 2022 The Authors.International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies Proposition 3 t -(1d) t or variable bounds (1g) t -(1j) t are in IIS, which cannot be resolved.Hence, the model is infeasible.After exiting the outer loop the values [s t ] t∈T are used for resolving any infeasibility in [z t ] t∈T through the time-connected linear program (LP) (see Proposition 2) defined as Ẑ([s t ] t∈T , γ ) := arg min z∈Z([s t] t∈T ,γ )