Robust optimization of a bi-objective tactical resource allocation problem with uncertain qualification costs

In the presence of uncertainties in the parameters of a mathematical model, optimal solutions using nominal or expected parameter values can be misleading. In practice, robust solutions to an optimization problem are desired. Although robustness is a key research topic within single-objective optimization, little attention is received within multi-objective optimization, i.e. robust multi-objective optimization.This work builds on recent work within robust multi-objective optimization and presents a new robust efficiency concept for bi-objective optimization problems with one uncertain objective. Our proposed concept and algorithmic contribution are tested on a real-world multi-item capacitated resource planning problem, appearing at a large aerospace company manufacturing high precision engine parts. Our algorithm finds all the robust efficient solutions required by the decision-makers in significantly less time than the approach of Kuhn et al. (Eur J Oper Res 252(2):418–431, 2016) on 28 of the 30 industrial instances.


Introduction
Two common challenges while dealing with complex real-world decision-making problems are (i) uncertain input parameters and (ii) multiple conflicting goals. The former may be a consequence of measurement errors, imprecise forecasts, unwanted disturbances, and a lack of historical data. The latter can be due to the involvement of multiple stakeholders/agents with competing goals. This work deals with bi-objective discrete optimization problems where one of the objective functions contains uncertain parameters. Few applications of uncertain bi-objective optimization are considered in the literature. Kuhn et al. [21,Sect. 8.1], consider a flight route planning problem (first presented in [20]), the two objectives being route efficiency-in terms of minimizing travel time-and risk-in terms of minimizing exposure to weather conditions; the latter (weather condition) is considered uncertain. A hazardous materials transportation problem is studied by Kuhn et al. [21,Sect. 8.2], with the objectives to minimize travel times and the number of people exposed in case of an accident, the travel times being considered uncertain, especially in urban areas.
One such application is presented in this work: it is focused on the planning of machining capacity for a large aerospace tier-1 [5, p. 64] engine system manufacturer with uncertain qualification costs, i.e. non-recurring fixed costs for verifying or assisting a machine in performing a specific task. The problem is the so-called tactical resource allocation problem (TRAP), a multi-item capacitated resource planning problem with a medium-to-long term planning horizon. A corresponding model is discussed in detail in [14,Chap. 2] and [13]. A solution to the TRAP is the allocation of tasks to machines and a schedule for the qualifications to be performed. The two objectives considered are minimizing the maximum excess resource loading above a given threshold and minimizing the total qualification costs incurred due to man-hours needed to verify machines for tasks or buying new fixtures/tools to adjust the machine for a given task. The uncertainty of the coefficients of latter objective function is not considered in [13]. In this work, we assess the impact of uncertainty in the qualification costs on the efficient frontier for 30 industrial instances.
Various methods have been proposed for solving uncertain multi-objective optimization problems (MOOPs) using stochastic programming and robust optimization. Stochastic programming for MOOP (see [16]) is used when enough data is available. A drawback of the stochastic approach is that for some problems so-called long-run optimality is not relevant, as the repeatability element of the decisions is missing; the decision-maker has to live with the consequences of the decisions made once. Since in the TRAP, qualification costs are incurred only once, it is evident that combining robust optimization and multi-objective optimization has certain benefits over other approaches. Moreover, the absence of historical data or expert opinions as bases for assumptions about the underlying probability distributions does not encourage us to use stochastic programming. Although a lot of work has been done in the field of single-objective robust optimization, its generalization to MOOP is relatively new. [1,3] present some basic concepts for single-objective robust optimization which have been used for MOOPs (e.g. [15,29,30]). Each of these relies, however, on some form of a priori scalarization and does not assess the impact of uncertainty on the efficient frontier in different scenarios. Consequently, so-called robust efficiency concepts for generalizing efficiency in MOOPs to uncertain MOOPs have been developed (see the discussion in Sect. 2). Most of the main concepts are described in the survey [18]. Our focus in this work is on bi-objective optimization problems with one uncertain objective function and a deterministic feasible set. Specifically, we build on the contribution of Page 3 of 31 36 [21] by suggesting a new robust efficiency concept and an algorithm that shows improved results on the industrial instances considered.

Relevance to multi-objective multi-agent decision problems
A multi-agent system (MAS) contains multiple agents deployed in a shared environment. Even though several multi-agent systems have multiple objectives, they are commonly (often erroneously) modeled as a single-objective decision problem. There has been some interest in using multi-objective approaches for MAS. [26] present several multi-objective approaches used for decision problems in MAS. To establish a relation between our contribution-on a bi-objective optimization problem with an uncertain objective function-to multi-objective MAS, we present an example of the Multi-Objective Normal Form Game (MONFG). 1 Suppose two agents wish to commute from an origin to a final destination. There are two modes of transportation: taxi and train. When both agents take a taxi, they split the cost. However, if both choose to take a train, they pay for their individual tickets. The exact costs of the train tickets and the taxi fare are available, but the travel time is uncertain for both modes of transportation. We assume two scenarios: expected/nominal (nom) and a worst-case (wc). The costs and travel times are given by cost time(nom) time(wc) taxi 20 10 15 train 5 30 32 and a corresponding individual payoff matrix for the commuting MONFG is given in Table 1 (the negative values representing costs). 2 The actions of the two agents result in a reward vector, representing cost and travel time. For instance, if agent X takes a taxi and agent Y takes a train, the reward received by the two agents in the nominal scenario is (−20, −10) and (−5, −30) , respectively. In the worst-case scenario, the two agents receive the reward (−20, −15) and (−5, −32) , respectively. Assuming that the two agents are cooperating (considering a team reward) and that there is no predefined utility function, it would make sense to identify actions that are equally good. However, if the two agents choose different modes of transportation the individual reward vectors of both agents have lower values than if both agents take a taxi, in the nominal as well as the worst-case scenarios. Hence, the actions of both agents choosing the same mode of transportation are of interest to the system planner. This type of problem can be modeled as bi-objective optimization problems with one uncertain objective (here, the team reward for travel time is uncertain), which is the topic of our work. Since we assume that the two agents are cooperating, the rewards for the two agents should be combined, e.g. summing costs and travel times, respectively. For instance, when both agents take a taxi, the total reward is (−20, −20) in the   In real-world problems there are  often several possible actions (or feasible solutions); hence, an explicit payoff matrix cannot be established, which calls for a combination of multi-objective and robust optimization.

Contributions and outline
The contribution of this work is in the field of bi-objective optimization with one uncertain objective function with well-defined nominal (most likely) and worst-case scenarios. Our first contribution is the new robust efficient (RE) concept of positive robustness -representative lightly RE solutions and a measure to evaluate the net gain by the inclusion of a worst-case scenario to the optimization problem. Our second contribution is the suggestion of a new approach to identify relevant RE solutions to the so-called tactical resource allocation problem (TRAP), leading to computational gains over existing approaches such as the algorithm presented in [21]. In Sect. 1.3, we present a mathematical model for the TRAP. In Sect. 2 we present some of the existing robust efficiency concepts from the literature. In Sect. 3, we present a new robust efficiency concept for bi-objective robust optimization problems with one uncertain objective function. Further, we present a measure for assessing the net gain by the inclusion of a worstcase scenario. In Sect. 4, we present a so-called 3-stage method and perform numerical tests on 30 industrial instances of the TRAP.

Problem description
Formally, the TRAP is defined as follows (see Table 2 for notations). Definition 1 (TRAP) Given a set J of job types (tasks) and a set K of machines, let p jk be the average 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 must not exceed the value of the parameter ∈ ℤ + . For assignments (j, k), such that k ∈ N j ⊆ K j and j ∈ J , so-called qualifications are required which generate additional one-time costs ( q jk , where q ∈ Q is the index of a scenario). For a job type j ∈ J , the machines in the set K j ⧵ N j do not require any qualifications. The total number of qualifications performed per time period t must not exceed the value of the parameter ∈ ℤ + .

Model description
The two objectives considered and the constraints defining the feasible set for the TRAP is expressed as (the index q ∈ Q denotes a scenario) (1a) minimize , , , subject to (1b) minimize , , , ∑ t∈T∶t≤l z jkt ≥ s jkl , k ∈ N j , j ∈ J, l ∈ T, Variables Description n t ∈ ℝ + Maximum resource loading above thresholds k , k ∈ K , in time period t ∈ T = ( , , , ) Bold notations representing vectors of the corresponding indexed variables Parameters Description a jt ∈ ℤ + Demand of orders of job type j ∈ J in time period t ∈ T p jk ∈ ℚ + Average machining time in machine k ∈ K j for job type j ∈ J C kt ∈ ℤ + Capacity (hours) available in machine k ∈ K in time period t ∈ T q jk ∈ ℤ + Qualification cost in scenario q ∈ Q , for machine k ∈ N j for job type j ∈ J ∈ ℤ + Upper limit on the number of qualifications in a single time period ∈ ℤ + Upper limit on number of alternative machines for each job type in a single time period We denote the set of feasible solutions as The set of job types ( J ) contains unique tasks (operations such as milling, turning, and grinding) to be performed on parts/products. Each job type j ∈ J has the demand a jt in time period t ∈ T , and the production should equal this demand, as expressed in (2a). The constraints (2b) ensure that no job is performed in any machine and time period to which it is not allocated. The constraints (2c) limit, for each job type and time period, the number of alternative machines to , the value of which is given as an input by the user; a too small value of may result in an empty set of feasible solutions and a too high value may be inconvenient for the production planners due to a resulting increased complexity of the product route. The constraints (2d) and (2j) make sure that the allocated machining hours do not exceed the capacity of the machines at any time period. The constraints (2d) are also used to define the objective function (1a), which is to minimize 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. The constraints (2e) implies that if a job type j is scheduled in machine k during time period t, then qualification of the machine for this job type must be done once during the time periods {1, … , t} . The constraints (2f) limit the number of qualifications allowed to be scheduled in each time period to ; this is due to a limited number of skilled professionals for completing new qualifications. The constraints (2g)-(2j) define the allowed values of the variables.

Qualification cost
The qualification cost is the one-time cost incurred by the company to qualify a machine for a job type. Once qualified, the new allocation can be used in all the subsequent time periods. For each allocation (j, k), j ∈ J , k ∈ N j , in scenario q, the qualification cost parameter q jk is represented by a natural number; hence, the uncertainty set is finite. It is common in the robust optimization literature to assume a nominal or most-likely scenario (see [3]). We denote the uncertainty set by Q ∶= {q,q} , where the indices q and q represent the nominal and worst-case scenario, respectively; for each j ∈ J and k ∈ N j it holds that q jk ≤̃q jk . Since the qualification costs of any two allocations (j � , k � ) and (j, k) are independent, the nominal and worst-case values, q and q , respectively, are well-defined. For uncertainty sets with more than two scenarios, most of the robust efficiency concepts to be defined will be retained, unless otherwise mentioned.

Preliminaries and an example to highlight robust efficiency concepts
Both robust and multi-objective optimization have, respectively, rich sets of literature. We introduce most of the key concepts in robust multi-objective optimization, aided by an instance of the TRAP, starting with some necessary notations ( [10]). For any two vectors , ∈ ℝ b it holds that In a robust counterpart to a deterministic version of the TRAP, the objective function is defined as ∶ Y × Q ↦ ℝ 2 + , i.e. the scenarios corresponding to the indices in Q affect the objective values. An uncertain bi-objective TRAP is defined as a set of parametrized problems, as where P(q) denotes the instance (corresponding to scenario q) of a bi-objective optimization problem, and the functions g 1 ∶ Y ↦ ℝ + and g 2 ∶ Y × Q ↦ ℤ + are defined by Note that the feasible set Y is independent of the scenario q. As compared to the cases of single-objective robust optimization 3 (SO-RO) and single-scenario (hence, deterministic) multi-objective optimization 4 (SS-MOOP), it is not evident what properties characterize good solutions. Consequently, we evaluate some of the key concepts mentioned in [18,21], and [11], and relate these to our model. However, most of the concepts apply to general cases as well i.e. for robust multi-objective optimization problems (robust MOOPs) with more than two objective functions. All set notations defined are listed in Table 4.
3 A field of optimization theory that requires a certain measure of solution robustness against uncertainty in parameters with a single objective function 4 A multi-objective optimization problem with no uncertain parameters is a non-dominated point in the criterion space corresponding to the objective functions g 1 (⋅) and g 2 (⋅, q) . The set of efficient solutions to the bi-objective optimization problem P(q) is denoted as Y ef f (q) , for the objective functions (⋅, q) = (g 1 (⋅), g 2 (⋅, q)).
As an illustrative example, we assume eight feasible solutions to a given instance of the TRAP denoted as i = ( i , i , i , i ) ∈ Y , i ∈ {1, … , 8} , and two scenarios Q ∶= {q,q} . The corresponding objective vectors are listed in Table 3 and visualizedin the criterion space-in Fig. 1.
Since in a bi-objective problem with one uncertain objective function with two scenarios each solution maps to two points in the criterion space, the concept of efficiency is not well-defined (see Fig. 1). We next present the main concepts of robust efficiency, the counterparts of efficiency for SS-MOOP.

Definition 3 (Flimsily Robust Efficient (FRE); [4]) A point ∈ is called flimsily RE
(FRE) if it is an efficient solution to the deterministic MOOP P(q) , for at least one scenario q ∈ Q . The set of FRE solutions is defined as where Y ef f (q) is the set of efficient solutions for the SS-MOOP P(q) (see (5b)). In our recurring example, Y FRE = { 1 , 2 , 3 , 4 , 7 , 8 } (cf. Fig. 1).

Definition 4 (Highly Robust Efficient (HRE); [18]) A point ∈ is called highly RE
(HRE) if it is an efficient solution to each of the deterministic MOOPs P(q) , q ∈ Q . The set of HRE solutions is defined as Unless some necessary conditions are satisfied [18,Le. 9], it is not guaranteed that every instance of a robust MOOP possesses HRE solutions.
In our recurring example, Y HRE = { 1 , 3 } (cf. Fig. 1). We next present a concept that mitigates the conservativeness of choosing the best solutions for the worst-case scenario. It is based on the light robustness concept for SO-RO problems ( [12]). It characterizes solutions that, w.r.t. an -neighbourhood (in the criterion space) of an efficient solution in the nominal scenario, are efficient in the worst-case scenario.
Definition 5 ( -lightly RE; [27]) For a parameterized uncertain MOOP P(Q) with predefined nominal scenario ( q ), and ∈ ℝ 2 + , the set Y light (̂ , ) of -lightly robust efficient solutions w.r.t. an efficient solution ̂ ∈ Y ef f (q) to the deterministic MOOP P(q) is defined as the set of efficient solutions to where the equivalence holds when the worst-case scenario q is well-defined, and defines a neighborhood of ̂ .
The set of all -lightly RE solutions is defined as For each efficient solution in the nominal scenario, there may exist multiple -lightly RE solutions. Hence, the number of solutions to be evaluated by the DM may be quite large if the entire set Y light (Y ef f (q), ) is considered. To reduce this number, the following definition identifies a representative set of -lightly robust efficient solutions. Definition 6 ( -representative lightly RE; [21]) For a parameterized uncertain MOOP P(Q) with pre-defined nominal ( q ) and worst-case ( q ) scenarios, and ∈ ℝ 2 + , the set of -representative lightly robust efficient solutions w.r.t. an efficient solution ̂ ∈ Y ef f (q) to the deterministic MOOP P(q) is defined as The set of all the -representative lightly RE solutions is defined as In our recurring example (Fig. 1 [17] prove that the solutions obtained using only the traditional concept of strict robustness for SO-RO, 5 which focuses solely on the worst-case scenario can be dominated; instead, they introduce the concept of Pareto Robust Optimal (PRO) solutions for SO-RO problems. [21] suggest a generalization of PRO solutions from SO-RO to robust MOOPs, called PRO Robust Efficient (PRO-RE) solutions.
Defining a vector valued function  Table 3), it holds that ( 7 ) ⪯ ( 8 ) . Hence, 8 is not a PRO-RE solution. HRE solutions are, however, PRO-RE, since they are efficient in all scenarios.
The RE concepts presented in [17,21,1,11,27,7], and [18] provide a sufficient basis for our contributions. Furthermore, as per [21, Def. 9]-for bi-objective-and [7]-for general multi-objective-robust optimization problems, flimsily, highly, -lightly, and -representative lightly RE solutions must be PRO-RE, (cf. Def. 7). The sets corresponding to each type of RE while satisfying the constraints (12) are defined as Similarly, the corresponding set of PRO solutions that are efficient in scenario q is defined as

Significance of uncertainty
As the numbers of uncertain objective functions or scenarios increase, the number of objective functions in the corresponding deterministic MOOP increases [17], hence increasing the computational effort needed to find all PRO-RE solutions. E.g., the bi-objective TRAP with one uncertain objective and two scenarios results in solving a tri-objective deterministic MOOP to find all PRO-RE solutions (see Def. 7). [8] established theoretically that the number of non-dominated points grows exponentially with the number of objective functions. Hence, restricting the number of uncertain objective functions or scenarios make the problem more tractable. However, ignoring uncertainties may also result in sub-optimal or too sensitive solutions. Consequently, it is important to understand whether the inclusion of scenarios reveals any deficiency caused by only using the nominal scenario for instances of a given problem class.

Motivation
To assess the gain from the inclusion of a worst-case scenario one should perform numerical tests (over a fairly large number of numerical instances representing various possible realizations of the TRAP). We consider two instances of the TRAP, illustrated 6 in Fig. 2. In Fig. 2 (left), all the efficient solutions ∈ Y PRO ef f (q) are highly PRO-RE (i.e. Y PRO ef f (q) ⊆ Y PRO HRE ). Thus, nothing is gained by including the worst-case scenario for this particular instance. However, for the instance #18 (see Sect. 5.1 and 7.2 for details on instances) in Fig.2 (right), none of the efficient solutions (i.e. Y PRO ef f (q) ) are in Y PRO HRE . Figure 2 considers ̂ ∈ Y PRO ef f (q) and ̃ ∈ Y PRO ef f (q) , corresponding to the vectors (̂ ) = (0, 14, 54) ⊤ and (̃ ) = (0, 17, 43) ⊤ , respectively (see Def. 7). Hence, the solution ̂ maps to the two points (̂ ,q) = (0, 14) ⊤ (left-most green-cross) and (̂ ,q) = (0, 54) ⊤ (left-most green-square), while the solution ̃ yields (̃ ,q) = (0, 17) ⊤ (left-most red cross) and (̃ ,q) = (0, 43) ⊤ (left-most red square). It follows that three units of qualification cost are lost in the nominal scenario, while in the worst-case, the qualification cost is reduced by eleven units. Hence, a total net reduction of eight units of qualification cost is obtained by this swap. A DM relying on its tolerance for loss of nominal quality may find this swap useful. We formalize this idea in the next section.

Positive robustness
We next present a new RE concept for the TRAP to identify solutions having a positive effect on mitigating risks due to the worst-case scenario.
Definition 8 (positive robustness -representative lightly RE solution) For a parameterized uncertain bi-objective MOOP P(Q) with one uncertain objective function g 2 , predefined nominal ( q ) and worst-case ( q ) scenarios, ∈ ℝ 2 + , and 0 < ≪ 1 , the set of positive robustness -representative lightly robust efficient solutions for the solution ̂ ∈ Y ef f (q) is defined as The set of all positive robustness -representative lightly RE solutions is defined as Positive robustness refers to the expression v( ,̂ ) > 0 , measuring the net reduction of qualification costs from swapping ̂ ∈ Y ef f (q) and ∈ Y nb (̂ , ) . The constraint ∈ Y nb (̂ , ) limits the loss of nominal value while also yielding a representative set for each efficient solution in the nominal scenario, which holds also for -representative lightly RE solutions (cf. Def. 6).
A corresponding PRO-RE set is denoted as

Remark 1
The set Y PRO r−light (Y ef f (q), , ) can be considered by the DM instead of the -representative lightly PRO-RE set Y PRO r−light (Y ef f (q), ) , as solutions with positive robustness are not guaranteed to exist in the latter. In Fig. 3, three solutions are in the set Y PRO ef f (q) ⧵ Y PRO ef f (q) (corresponding points in the criterion space are distinct red-cross and -square marks) two of which belong to Y r−light (Y PRO ef f (q), ) ; marked by black-dashed rectangles. None of these two solutions, however, possess the so-called positive robustness as compared with the respective closest nominal PRO-RE solutions (marked by a green cross at the lower-left vertex of the respective rectangle). Hence, solutions without positive robustness have a limited utility for the DM.

Inclusion of a worst-case scenario
The topic of assessing the importance of including a worst-case scenario for an uncertain MOOP, over a fairly large number of numerical instances has, to the best of our knowledge, not been discussed in the existing literature. While the positive robustness -representative lightly RE solutions rely on the values of , we next present a measure that does not require specifying .
A worst-case scenario and the corresponding efficient solutions in Y PRO ef f (q) add a significant value to a robust bi-objective optimization problem if the set ∪ ∈ P ef f (̃) { ( ,q)} is a good approximation 7 of the efficient frontier ∪ ∈ P ef f (̂) { ( ,q)} . Further, ∪ ∈Y PRO ef f (q) { ( ,q)} is not a good approximation of the efficient frontier ∪ ∈Y PRO ef f (q) { ( ,q)} for the worst-case scenario.
Definition 9 (Scenario approximation) For each r, q ∈ {q,q} , the region in the criterion space of (g 1 (⋅), g 2 (⋅, r)) , defined by the set where is dominated by the worst-case ( r =q ) or nominal ( r =q ) objective vectors corresponding to PRO-RE solutions that are efficient in the worst-case ( q =q ) and nominal ( q =q ) scenario, respectively.
The nominal ( q ) and worst-case ( q ) scenario approximation area is defined as (A(q,q) ⧵ A(q,q)) and (A(q,q) ⧵ A(q,q)) , respectively, where (B) measures the area covered by a set B ⊂ ℝ 2 .
In Fig. 4 the grey shaded area illustrates the set A(q,q) ⧵ A(q,q) , i.e. the region in the criterion space that is not dominated by any non-dominated point in the worst-case scenario (the set ∪ ∈ P ef f (̃) { ( ,q)} ), but is dominated by non-dominated points in the nominal scenario (the set ∪ ∈ P ef f (̂) { ( ,q)} ). Analogously, the red shaded area illustrates the set A(q,q) ⧵ A(q,q) , i.e. the region that is not dominated by any non-dominated point in the nominal scenario (the set ∪ ∈ P ef f (̂) { ( ,q)} ), but is dominated by non-dominated points in the worst-case scenario (the set ∪ ∈ P ef f (̃) { ( ,q)}).

Definition 10
The difference between the worst-case and nominal scenario approximations indicates the utility of including the worst-case scenario. We suggest a measure defined as the ratio of the difference in the areas of the worst-case and nominal scenario approximations and the maximum of these two values, as 8 If Δ dif f ≤ 0 , the worst-case scenario is not expected to add significant gains, whereas if Δ dif f > 0 , it is expected to result in solutions that provide significant gain. We validate this claim empirically in Sect. 5.3.

Solution approaches for bi-objective robust optimization with one uncertain objective function
From the definitions in (13), RE solutions are required to be PRO-RE. To the best of our knowledge, the only fully autonomous algorithm suggested for robust bi-objective optimization problems with one uncertain objective function is the one presented in [21,Sect. 6]. Some interactive methods-relying on users providing their preferences-are, however, suggested in [23,28], and [24, pp. 27-57]. The algorithm by [21] first finds a minimal set 9 of efficient solutions in Y ef f to a deterministic tri-objective (for |Q| = 2 ) optimization problem (cf. Def. 7). Then, it filters the resulting minimal subset of Y ef f for flimsily/highly/ -representative lightly PRO-RE solutions.
Note that it requires computationally expensive tri-objective optimization to identify all the PRO-RE solutions. Further, the set Y ef f may contain more solutions than required by the DMs (typically the ones in (13)).

A new approach to find required PRO-RE solutions
We propose a 3-stage method (Alg. 1): (i) compute Y PRO ef f (q) and Y PRO ef f (q) ; (ii) identify Y PRO HRE and Y PRO FRE from solutions identified in (i); (iii) for a given , compute Y r−light (Y PRO ef f (q), , ) . The set Y r−light (Y PRO ef f (q), , ) obtained is not guaranteed to be PRO-RE. 10 For instance, in Fig. 5 the second worst-case efficient solution from the left (marked red-cross and redsquare) dominates the solution corresponding to the second positive robustness -representative lightly RE solution from the left (marked blue-cross and blue-square). 11 In Sect. 4.3 we suggest a procedure to prevent this issue.

Finding nominal and worst-case PRO-RE solutions: Lex-BiObjective
In the first stage of Alg. 1 two bi-objective optimization problems are solved, to find Y PRO ef f (q) and Y PRO ef f (q) , respectively. The corresponding bi-objective optimization problem is denoted as P PRO (q) and is defined as The two minimization problems in (19) corresponding to q and q are bi-objective mixedinteger programming problems; 12 any efficient solution to the two is PRO-RE, as per Def. 7. Thus, a non-dominated point must possess a minimal value w.r.t. the other objective (i.e. g 2 (.,q) for P PRO (q) and g 2 (.,q) for P PRO (q)).
A scalarized problem within an -constraint method to solve an instance of P PRO (r) , for r ∈ Q , is given by where u ∈ Q ⧵ r , the constant is an upper bound used to explore the criterion space during each iteration of the -constraint method (updated in each iteration), and w(r) > 0 is appropriately small to ensure at least a weakly efficient solution w.r.t. g 1 and g 2 ( , r) (see [22]). Lexicographic minimization (lexmin) is performed to obtain a PRO-RE solution (for details on lexicographic minimization, see [10,Sect. 5.1]). (19) min 12 Since , , and are integral, in an efficient solution, the variables will take a finite number of discrete values. Hence, the efficient frontier contains a finite number of non-dominated points and no continuous line segments  If the solution to (23) from the first iteration of constraint generation is such that 1 = 0 3 , then the solution ̄ is PRO-RE; otherwise, the solution 1 (for f = 1 ) fulfills ( 1 ) ⪯ (̄ ) . In the latter case, the variables m if ∈ {0, 1} and the constraints (24) are added to the models (21) and (22), for the k th iteration of the constraint generation:

Tests and results
To investigate the computational efficiency of our proposed approach, we generated 30 instances, which are expected to represent some of the possible realizations of actual data that the model might encounter. Hence, allowing for confident conclusions about the benefits of including worst-case scenarios and computational performance of our proposed approach.
All the computations are done in Python 3.7 using Gurobi 9 on a system with 1.70GHz processor, 16 GB RAM, and 4 cores. We set a time limit of 5000 seconds for each scalarized optimization problem, and a global time limit on the algorithm is set to 20000 seconds for the entire process. We also terminate if the MIP duality gap is less than 0.05% . Another termination criterion records and updates the best incumbent solution after every 200 node explorations; if there is no improvement in the optimality gap, the solver terminates the process. However, this termination criterion is only activated after 1500 seconds. For a fair comparison these termination criteria are implemented for all the algorithms used to benchmark against our proposed approach. Moreover, all computed sets reported are minimal sets.

Industrial instances
We have used real industrial data for most of the parameters and sets indicated in Table 2. However, processing times p jk of job type j ∈ J in machines k ∈ N j (qualification required), and the qualification cost parameters q jk , j ∈ J , k ∈ N j , q ∈ Q , were not available. In order to generate instances which represent possible realizations of the actual data we introduce the following distributions, which are based on knowledge of the managers.
Skewness of processing times. The skew normal distribution is a generalized normal distribution allowing for non-zero skewness; [31]. We 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 ( > 0 ), negative skew ( < 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 ∶= 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.
Range of values for the qualification cost: The exact cost for qualifying a machine for a job type is not known a priori, and for prediction the engineering team has to spend time on simulations. Hence, the input received is the so-called cost levels, assigned to each qualification. For testing our model and proposed modifications, we use two sets of cost levels: H = {1, … , max } , with max ∈ {20, 50} . The qualification costs are selected from different discrete distributions over the discrete domain H.
Nominal qualification cost. 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

Constant data for the 30 instances
The demand is from quarterly forecasts made at GKN Aerospace in January 2015 ( J 15 ) for the period 2016-2017. The minimum, maximum, and median values of the demand are 1, 172, and 11, respectively. For processing times p jk , j ∈ J, k ∈ K j , the minimum, maximum, and median are 0.1, 89.7, and 5.63 hours, respectively. Each machine has a yearly capacity of 5000 hours which can be equally divided among four quarters in a year. The planning period of two years with quarterly time buckets yields T = 8 time periods. There are K = 125 machines, and about J = 510 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 5 . The parameter values = 3 , = 4 , and k = 0.7 , k ∈ K , are used for all 30 instances. The instances are denoted as ( ̄ , p , max ), where ̄∈ {Lef t, Right, Symmetric, Uniform, Bimodal} (representing distributions for the qualification costs), p ∈ {skew+, skew−, skew0} (representing distributions for processing times p jk , k ∈ N j ), and lastly, max ∈ {20, 50} . The details of each instance are indicated in the captions of the subfigures of Fig. 11. All 30 public instances are available at https:// www. short url. at/ uGVX3.

Results
We present two main results, the first regarding the difference of the areas of the worst-case and nominal scenario approximations (Def. 9). The second result concerns the benefits of using the 3-stage method as opposed to the approach of [21] for identifying relevant PRO-RE solutions. when Δ dif f ≤ 0 , which is the case for remaining twelve instances (the twelve left-most in Fig. 7), only four instances (16, 11, 1, and 7) have solutions with positive robustness values ( |Y PRO r−light (Y PRO ef f (q), , )| > 0 ). A possible explanation is that Δ dif f is an aggregate measure, hence these instances have positive robustness values even though Δ dif f ≤ 0 (see Fig. 11). Hence, Fig. 7 indicates that the worst-case scenario does have an effect for a significant proportion of the instances of the TRAP. Note that the values of influence the value of |Y PRO r−light (Y ef f (q), , )| . For our problem it is based on the range of values for the two objective functions ( 1 = 0.12 ⋅ (g BOT 1 (q) − g TOP 1 (q)) and 2 = 0.12 ⋅ (g TOP 2 (q) − g BOT 2 (q))).

Solution time comparisons
We compare the solution times obtained from using 3-stage method (Alg. 1) with the results obtained if algorithm from [21] is used for the 30 instances. As mentioned in Sect. 4 §2 the first step in [21] requires a tri-objective optimization. For this purpose a tri-objective criterion space search method called Quadrant Shrinking Method (QSM) [6] is used. We implemented the QSM [6, Alg. 1, p. 877] in Python 3.7 using the data structures suggested in [6]. 14   instances our 3-stage method is computationally superior. This validates (empirically) that solving two bi-objective optimization problems is faster than solving one tri-objective optimization problem for the given instances of the TRAP. (b) For the remaining 13 instances the two methods found different numbers of PRO-RE solutions; only for instance 11 the QSM (tri-objective) method found fewer PRO-RE solutions than the 3-stage method. 15 Our approach does not find all the PRO-RE solutions, but manages to find all those that are required by the DMs in significantly less time than the solution approach mentioned in [21] (using the QSM in [6] for tri-objective optimization).

Computational complexity
We compare the computational complexity of the algorithm in [21,Sect. 6] with that of our 3-stage method Alg. 1; 16 see Table 5 for a summary.

RE solutions
Y ef f set of efficient solutions with objective functions Def. 7 ∶= (g 1 (⋅), g 2 (⋅,q), g 2 (⋅,q)) ; also called PRO-RE solutions 15 The solution time limit of 20000 seconds was reached for instance 11 16 For the set notations, see Table 4 17 Using the QSM method Alg. 1 involves three stages: In stage 1 the sets Y PRO ef f (q) , q ∈ Q , are computed, i.e. |Q| bi-objective optimization problems are solved; the use of an -constraint method then leads to solving O(2 ∑ q∈Q �Y PRO ef f (q)� + �Q�) single-objective MILPs. 18 (see [9]). Stage 2 filters the sets Y PRO FRE and Y PRO HRE , by identifying the intersection and union, respectively, of the sets Y PRO ef f (q) , q ∈ Q ; since both procedures utilize a sorting step, the total number of elementary operations is O( ∑ q∈Q (�Y PRO ef f (q)�⋅log �Y PRO ef f (q)�)) . Stage 3 makes, say b, calls to Alg. 2, each of which comprising an iterative algorithm that requires solving ∈ [1, |Y nb (̂ , )|] MILPs (if ≠ else = 0 ), where ̂ ∈ Y ef f (q) ; for fairly small values of , the upper bound on is not expected to be large (see (9)).
Obtaining estimates of b and (number of calls to Alg. 2 and iterations, respectively) is quite tricky. These measures are illustrated for our instances in Fig. 9, with the values of b on the left axis, and the share of the total solution time of Alg. 1 spent on stages 1 and 2 on the right axis. For 15 out of the 30 instances, there are no calls to Alg. 2 (hence, the share equals 1). For 25 instances the share is greater than 0.9. However, for a few instances (6, 11, 18, 23, 28, and 30) the calls to Alg. 2 (i.e. stage 3) constitute a significant portion of the solution time. E.g. for instance #11, stages 1 and 2 use only 56.55% of the total time, the remaining 43.45% being spent on calls to Alg. 2 (stage 3).

Conclusion
We propose an approach for solving robust bi-objective optimization problems, with one uncertain objective function including two well-defined scenarios; a nominal (most likely) and a worst-case scenario. We present a new measure to assess the net gain from the inclusion of the worst-case scenario.
For understanding the general applicability, we discuss the two ideas separately. We first consider the positive robustness -representative lightly robust efficient (RE) concept, which is similar to that of -representative lightly RE, except that the former satisfies an additional constraint (cf. (15)). As -representative lightly RE solutions [21], our proposed RE solutions apply to MOOPs with one uncertain objective function and any number of scenarios as long as nominal and worst-case scenarios are well-defined. Second, the difference between worst-case and nominal scenario approximations can be used to filter out instances that are expected not to gain significantly from the addition of a worst-case scenario, and thus avoiding the need to find solutions in the set Y PRO r−light (Y PRO ef f (q), , ) (Def. 8). Our proposed measure Δ dif f (Def. 10) is applicable to bi-objective optimization problems, with one uncertain objective function having any number of scenarios, provided that the worst-case and nominal scenarios are well-defined.
Extensions and future research. For most of our 30 instances, the computation time is too long (see details at https:// www. short url. at/ nrIJQ) and may discourage decisionmakers from performing multiple what-if/sensitivity analyses. Hence, in-exact computing approaches could be investigated, especially if more scenarios are included. There is also a potential for generalizing our concepts to general robust multi-objective problems with several uncertain objectives. 18 Note that the inequality      (Y PRO ef f (q)) wc policy P-policy (Y PRO ef f (q)) nom policy P-policy (i) Effect of policies in the nominal scenario for g 1 (∆g 1 ) (ii) Effect of policies in the nominal scenario for g 2 (∆g 2 ) (iii) Effect of policies in the worst-case scenario for g 1 (∆g 1 ) (iv) Effect of policies in the worst-case scenario for g 2 (∆g 2 ) 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://creativecommons.org/licenses/by/4.0/.