Optimal Scheduling of Copper Concentrate Operations under Uncertainty

We propose a continuous-time scheduling model for the logistic and blending operations of copper concentrates with uncertain composition. The formulation, based on the MultiOperation Sequencing (MOS) model, gives rise to a large-scale nonconvex mixed-integer nonlinear programming (MINLP) model. We adopt a two-step MILP-NLP decomposition strategy and enhance the MILP relaxation to propose schedules that significantly reduce the optimality gaps. The bounded uncertainty in element composition of the concentrates is addressed by an extended robust MOS model, which combines robust optimization and flexibility analysis techniques. The effectiveness of the models and the solution strategy is validated with an illustrative example and an industrial case study.


Introduction
Copper is used extensively throughout domestic and industrial applications because of its excellent corrosion resistance, malleability, and thermal conductivity (Langner, 2011). For this reason, copper demand is constantly growing due to the increasing needs on information technology, energy supply, manufacturing, and other copper-intensive industries. Demand 5 growth has been pushing the copper industry to improve its production efficiency, creating opportunities for mathematical programming techniques to be applied in the copper industry to systematically improve the performance of different stages of the process. In this article, we optimize the logistic and blending operations of a copper smelting process, which is the first step of the refining process in the copper value-added chain shown 10 in Figure 1. Copper concentrates are solid materials obtained from copper ores that undergo a concentration process. The smelting process transforms the copper concentrates to copper mattes, increasing the copper content from 25-35% to 50-70%. One important challenge in this process is to find the best mixture of various raw concentrates and recycling materials, such that the gross margin of the process is maximized, the environmental standards are 15 met, and the quality of the products is guaranteed.
Scheduling problems are very prevalent in process industries such as chemicals, food, pharmaceuticals, and oil and gas (Harjunkoski et al., 2014). Generally, the objective of the scheduling problem is to maximize the profit, minimize the cost, or minimize the makespan.
The decisions in scheduling problems usually include (Maravelias, 2012): 20 • selection and lot-sizing of tasks • assigning tasks to resources (units) • timing of tasks • sequencing of tasks Due to the complexity of the scheduling problem, mathematical programming approaches 25 have been introduced to systematically improve the quality of the schedules. A number of review articles summarizing diverse mathematical models for scheduling have been pub-lished in the past (Méndez et al., 2006;Maravelias, 2012;Harjunkoski et al., 2014). They identify two distinctive types of models according to time representation: discrete-time and continuous-time formulations. Discrete-time formulations divide the time horizon into uni-30 form time intervals and often lead to large-size, tight formulations. Early works in process scheduling explored various discrete-time representations based on the State-Task Network model (STN) (Kondili et al., 1993) and the Resource-Task Network (RTN) model (Pantelides, 1994). Continuous-time formulations partition the time horizon using time-points, time intervals or event points, and the length of each time period is to be determined by 35 the optimization model. These formulations are often smaller in size because they need fewer time periods to represent the schedule, but they tend to have weaker relaxations.
Continuous-time models can be further divided into two types based on the coupling between task and unit events: single time grid (Castro et al., 2001) and multiple time grids (unit-specific event based) (Ierapetritou and Floudas, 1998a,b). 40 For the scheduling of copper concentrate operations, Song et al. (2018) used a discretetime formulation to model the process as an mixed-integer nonlinear programming (MINLP) problem; they compared the performance of a split fraction model and a process network model for the same problem. Lalpuria (2017)  erations. It has proved to be efficient for crude oil scheduling problems (Mouret et al., 2011).
An important challenge related to copper concentrate operations arises from the uncer-50 tainty in their element composition. After being mined, copper ores go through a complex process in which they are transformed into concentrates; the process include crushing, grinding, flotation, and drying. The complexity of these processes causes the element compositions of each concentrate to fluctuate within certain ranges around their nominal value. These fluctuations are critical for the preparation of the concentrate blend, since the copper refin-55 ing requires strict conditions on the composition of different elements. In order to comply with these conditions, it is important to take the variability in element composition into account by treating it as a bounded uncertainty within the optimization framework.
Research on optimization under uncertainty has been an active area since 1970s . Two major uncertainty modeling techniques are robust optimiza-60 tion (Bertsimas et al., 2011) and stochastic programming (Birge and Louveaux, 2011). Robust optimization follows the idea of optimizing for the worst realization of the uncertainty while guaranteeing the feasibility of the solution. On the other hand, stochastic programming optimizes the expectation of an objective function based on the probability distribution of the uncertain parameters. Results from robust optimization are more conservative than 65 those obtained from stochastic programming. Therefore, stochastic programming is suitable for long-term planning problems to maximize the expected profit, while robust optimization is better suited for process scheduling under uncertainty, as the feasibility of the solution is critical for short-term problems. Alonso-Ayuso et al. (2014) used stochastic programming to consider the effect of changing copper prices on the expected profit for a five-year 70 copper extraction planning problem. Li and Ierapetritou (2008) reviewed several methods for process scheduling under uncertainty, including both robust optimization and stochastic programming. In addition to these two methods, flexibility analysis  can also be a suitable approach for scheduling under uncertainty, since it focuses on feasibility by accounting for recourse actions. Flexibility analysis is an uncertainty modeling 75 technique that was originally developed for plant design under uncertainty, and it has been used in the process systems engineering community for more than thirty years. As shown in a recent article by Zhang et al. (Zhang et al., 2016), flexibility analysis can yield either identical or more rigorous results than robust optimization for linear problems, although it is computationally more expensive. 80 In this work, we propose a continuous-time scheduling formulation based on the MOS model. The goal is to schedule the logistic and blending operations necessary to feed the copper smelting process with suitable concentrate mixtures. The model includes robust optimization and flexibility analysis techniques to deal with the uncertainty in element composition of concentrates. A two-step MILP-NLP decomposition strategy is adopted and 85 enhanced to solve the scheduling models efficiently.
The remainder of this paper is organized as follows. The copper smelting process and the challenges derived from uncertainty in element composition are described in Section 2. The time representation and the mathematical formulation of the deterministic MOS model are presented in Section 3. In Section 4, the robust MOS model that considers uncertainty in 90 element composition is developed. Section 5 presents the solution strategy based on a twostep MILP-NLP decomposition strategy. Finally, an illustrative example and an industrial case study are described in Section 6; the corresponding results are shown in Section 7. The logistic and blending operations for a standard copper smelting process are represented in Figure 2. The concentrates available for smelting are located at storage facilities near the port since many of them arrive by maritime transportation. The concentrates are first unloaded from maritime vessels at the port and piled separately; then, they are transferred to a blending unit where an initial blend is prepared. The blends are then distributed 100 to several bins, from where they are simultaneously discharged into the smelter. In addition to concentrates, daily-arrival materials are also transferred to bins and processed at the smelter. These daily-arrivals are a special raw material consisting mainly of process leftovers and scraps; processing these materials improves resource utilization and contributes to energy saving of the smelting process.

Problem statement
to the operation rules, there are several restrictions on the element composition for the 125 blends feeding the smelter. These restrictions must be satisfied to guarantee the operating conditions necessary for the smelting and the successive downstream process.

Uncertainty in element composition
The natural variability in ore grades and small fluctuations in the concentration conditions are the main sources of uncertainty in the element composition of concentrates. Since 130 the smelting process is very sensitive to any changes in blend composition, this variability must be considered in the scheduling problem to avoid sub-optimal or even infeasible solutions yielding off-spec products.
For each concentrate, we consider a subset of elements whose composition fluctuates within a bound relative to its nominal value. Therefore, this variability can be formulated as a bounded symmetric uncertainty. By using f c,k to denote the mass fraction of element k in concentrate c, the uncertainty set can be expressed as follows: wheref c,k is the nominal mass fraction value of element k in concentrate c, δ is the relative bound of the fluctuation, and K U is the set of elements whose composition fluctuate. The

135
mass fraction (f c,k ) is related to various quality requirements of the process, which are formulated as a set of linear constraints in Section 3.

Multi-Operation Sequencing model
We use the MOS model to formulate the deterministic scheduling problem. The MOS model is a continuous-time, multiple-time-grid scheduling model that has been previously 140 applied to crude oil scheduling (Mouret et al., 2011); it has proved to be computationally efficient compared to other common time formulations. In this section, we first describe its time representation, then explain the mathematical formulation in detail.

Time representation
In the study by Mouret et al. (2011), the time representation of the MOS model is built 145 on the concept of priority slot, and the computational efficiency of the model is based on the application of cliques. However, the connection between the clique and the time representation in the model was not explicitly presented. Here, an extensive time representation for the MOS model is given to explain in detail the function of the cliques.
The time representation of the MOS model is based on assigning a priority slot i to 150 each execution of operation v. A graphic representation of the MOS model is presented in Figure 3, where the execution of each operation is represented by a segment labeled with a certain priority slot. The total number of priority slots |I| is specified a priori, which implies that each operation can be assigned to at most |I| priority slots. The time span of a certain slot for different operations are independent of each other. The priority slots specify the sequence in which the operations are executed. The priority mechanism works in two ways. For a single operation, the operation can be executed several times by assigning it to multiple priority slots. Each execution does not overlap with each other, and occurs in the order defined by the priorities. For multiple operations, the assignment of priority slots is based on the concept of cliques.

160
A clique represents a set of non-overlapping operations, i.e. operations that cannot be executed simultaneously, in the MOS model. In graph theory (Diestel, 2005), a clique refers to a subset of vertices in an undirected graph such that any two distinct vertices are connected to each other. In the MOS model, all operations can be mapped as vertices in a graph, and two vertices will be connected if the operations they represent are non-overlapping.

165
A small example of cliques is presented in Figure 4. Both Tank 1 and Tank 2 cannot be charged/discharged by more than one operation simultaneously, generating two sets of nonoverlapping operations ({v 1 , v 2 , v 3 } for Tank 1 and {v 3 , v 4 , v 5 } for Tank 2). Correspondingly, there are two cliques w 1 and w 2 in the process. Tank 1 Tank 2 For all the operations in a clique, at most |I| slots can be assigned to these operations, 170 and each slot can only be assigned to a single operation. These operations are executed in the priority order, and their execution time cannot overlap with each other. As an example, for clique w 1 in Figure 4, a potential schedule is presented in Figure 5. The three nonoverlapping operations are assigned to three different priority slots, and each slot is assigned to only one operation. This mechanism ensures that the executions of operations in the 175 same clique do not overlap with each other.  For instance, Figure 6 shows a feasible schedule that considers the interaction between both cliques (w 1 and w 2 ) presented in Figure 4. The schedule of operations in clique w 1 is covered by blue dots, and schedule of operations in clique w 2 is covered by a red block. In Figure 6, it can be observed that 3 slots are assigned to operations v 1 , v 2 and v 3 in clique w 1 , and 2 slots are assigned to operations v 3 and v 5 in clique w 2 . These assignments guarantee that 185 execution of operations within each clique do not overlap with each other. Since v 2 and v 5 belong to two different cliques, they can both be assigned to slot i 3 and their execution time overlaps with each other. In this schedule, operation v 3 connects cliques w 1 and w 2 . Notice that in this schedule, operation v 4 is not executed. Based on the cliques, the sequential relationship among non-overlapping operations is 190 established. The schedule for the whole process can be seen as a combination of the schedules of all cliques. In this way, each priority slot can be utilized by different cliques and the non-overlapping structure is captured in the time representation, leading to an efficient mathematical formulation.

Assignment constraints
Binary variable Z i,v denotes whether operation v is assigned to priority slot i. In clique W , at most one operation v can be assigned to slot i.

Cardinality constraints
The cardinality of operations in clique W should be within its lower and upper bounds.
For example, since vessels are required to unload all concentrates exactly once, the lower and upper bounds of cardinalities of unloading operations are both set to 1.

Time constraints
The end time E i,v of the operation v at slot i is the sum of its start time S i,v and its If operation v is assigned to slot i, then its start time, duration and end time should be within their lower and upper bounds respectively. Otherwise, the time variables should be zero.

Non-overlapping constraints
The non-overlapping constraints specify the time relationship between any two operations in a clique. In clique W , if slot i is assigned to an operation, then the gap between the end time (E j,v ) of operation at an earlier slot j and the start time (S i,v ) of operation at slot i should be at least the sum of durations corresponding to all slots between i and j. According to the assignment constraints (2) and time constraints (5)-(7), at most one operation in the clique W can have non-zero values for the time variables (E i,v , S i,v and D i,v ) at each slot.
Equation (8) is formulated as a big-M constraint, such that if slot i is not assigned to any operation, the time variables at slot i are all zero and the constraint is relaxed. Here, H is the length of the time horizon and serves as a big-M parameter.
In addition, it is necessary to ensure that executions of the same operation at different priority slots do not overlap with each other. This condition is enforced with Equation (9).

Mass balance constraints 200
Variable L c i,r,c denotes the accumulated level of concentrate c in unit r at slot i, and variable M c i,v,c denotes the amount of concentrate c in operation v at slot i. The mass balance constraints are based on these two variables at the concentrate level. L c i,r,c in unit r at slot i equals the sum of the initial level L c,initial r,c , the accumulation terms from inlet transfer operations, and the consumption terms from outlet transfer operations from all previous slots.
The level of concentrate c in unit r follows the priority sequence, i.e. the level of c in unit r is changing from L c,initial r,c to L c i1,r,c , L c i2,r,c , . . . , L c |I|,r,c . However, it is worth noting that only M c i,v,c follows the time representation for operation v at slot i, and there are no time variables calculated for L c i,r,c . In other words, there is no specific time information about the level variables, but only the relative sequence of them.

205
A disaggregated formulation of the mass balance constraints is presented in Equation (10 ). The level of c in unit r at slot i+1 equals the sum of the original level, the accumulation terms and the consumption terms M c i,v,c at slot i.
The final level L c,Final r,c equals the sum of the level, the accumulation terms and the consumption terms at the last slot.
Variable L t i,r denotes the total level of concentrates in unit r at slot i, and variable M t i,v denotes the total amount of concentrates in operation v at slot i. They equal the sum of corresponding variables at the concentrate level, respectively.
Similarly to Equation (10), L t i,r can also be represented as the sum of initial total level L t,initial r , accumulation terms and consumption terms from previous slots. This constraint is a linear combination of Equation (11)-(13) and is redundant to the model, but it can help to reduce the search space of the mixed-integer linear programming (MILP) subproblem.
Mass balances at the element level are implicitly included in the model. Since concentrates are solids and their element compositions remain constant through the process, the mass balance for the amount of each element can be formulated as a linear combination of mass balances at the concentrate level. This is different from the crude oil scheduling problem, because crude oils are liquids and blending will introduce bilinear terms in the 210 mass balance constraints.

Sequencing constraints
Different from other units in the process, bins and daily-arrival piles are allowed to be charged and discharged simultaneously, which makes the clique and Equations (2), (3), (8) non-applicable to inlet/outlet transfer operations of these units. Therefore, new constraints 215 need to be specified for the sequences of these operations.
The main rule ensures that outlet transfer operations of a particular blend are not executed before the preceding inlet operations. In order to simplify the model, outlet transfer operations are assigned to one slot after inlet transfer operations.
Furthermore, the start time of outlet transfer operations should be no earlier than the inlet transfer operations.
In addition, the inventory at these units must be non-negative throughout the time horizon. At slot i, the total outlet amount until the end time of the outlet operation is bounded by the maximum potential inlet amount, which can be calculated from the upper bounds of incoming flowrates (F v ) and the end time of the outlet operation (

Composition constraints
The proportions of different concentrates in a blend must remain consistent in the subsequent transfer operations. The fraction of a concentrate in the blend being transferred in ; similarly, the fraction of a concentrate in the blend being accumulated in unit r can be represented as . The consistency of the proportions can be achieved by setting these two terms equal to each other for the units and their corresponding outlet transfer operations, so that these proportions stay the same along blending units and the downstream operations.
To avoid singularity, the fractions can be reformulated as bilinear equations:

Capacity constraints
If operation v is assigned to slot i, then M t i,v should be within its lower and upper bounds.
If operation v is assigned to slot i, then the total flowrate of v should be within its lower and upper bounds.
This constraint can be reformulated as the linear inequality (20). If operation v is not The total level in unit r at slot i and at the end of time horizon must be within its lower and upper bounds.

Element concentration constraints
As stated in previous sections, there are several quality requirements in this process, 220 both for the main products and for the byproducts. These requirements are formulated as linear constraints that set restrictions on the ratios of different key elements in final blend feeding the smelter.
The proportion of key elements in the final blend of concentrates transferred into the smelter is restricted as follows: v∈V Final c∈C where V Final denotes the set of final operations that transfer concentrates from bins to the smelter, the numerator denotes the total amount of element k in the final mixture, andf k denotes the upper bounds of the fraction of element k in the final mixture. This constraint can be rearranged as a linear inequality:

Interdependency constraints
The ratio between two key elements k and k can be bounded as follows: This constraint can be rearranged as a linear inequality: Additional constraints are necessary to guarantee that the weighted proportion of a certain element k does not exceed its upper bound UE k , where KE k denotes the weighting coefficient for the elements. This constraint can also be rearranged as a linear inequality: v∈V Final c∈C

Operation rules for daily-arrival materials 225
The daily-arrival materials must be transferred to the daily-arrival piles continuously throughout the time horizon.
The total inventory of daily-arrival materials at the end of time horizon is restricted.
The sum of inventories of these materials in all the daily-arrival piles and bins should not exceed a given boundL daily .

Operation rules for bins
The concentrates in all bins must be transferred to the smelter simultaneously. Therefore, will not be assigned to slot i.
An additional constraint is necessary to ensure that the contribution of each bin feeding the smelter is maintained within a certain range.
This constraint can be rearranged as the following linear inequality:

Operation rules for smelter
The operation of the smelter is constrained by its processing capacity. In this context, the total flowrate of concentrates transferred to the smelter should be below its upper bound.
This constraint can be rearranged as a linear inequality:

Symmetry breaking constraints
A general symmetry breaking constraint can be applied to the MOS model (Mouret et al., 2011). The logic implies that an operation v should not be assigned to i if no other non-overlapping operations v are assigned to i − 1. This is because if the slot i − 1 is empty,

Objective function
The objective function is to maximize the gross margin of concentrates transferred to the smelter as presented in Equation (MOS), where G c denotes the unit price of concentrate c ($/ton).

230
The complete formulation for the deterministic MOS model is obtained by maximizing the objective function (MOS) subject to constraints (2)-(32). This model is a nonconvex MINLP problem due to the bilinear constraint (18).

Robust MOS model
As mentioned earlier, there is inherent uncertainty in element composition of concentrates 235 due to the upstream processes in which they are produced. This uncertainty is crucial to the quality requirements of the smelting process modeled with linear constraints (23)-(25).
If the MOS model is solved using the nominal values of element composition, there is a good chance that some of the quality requirements be violated, and the solution is no longer feasible. Therefore, it is necessary to extend the MOS model to gain robustness against the 240 impact of variability in the composition of concentrates.
In our process, the element composition of concentrates is unknown through the entire time horizon. Therefore, a static robust optimization technique, i.e. robust counterpart, is suitable to account for the uncertainty by optimizing under its worst realization (Bertsimas and Sim, 2003). In most robust optimization approaches, the value of the uncertainty 245 corresponding to the worst case is usually dualized and eliminated in explicit form. Such information is valuable for the robust model, as it can be used to verify the model validity based on our existing knowledge of the process. In this context, the flexibility test problem  can be included in the robust framework as a complement to reveal the values of element composition obtained in the robust counterpart.

250
Since there are multiple quality requirements corresponding to different products at different stages of the process, it becomes impossible to satisfy all of them simultaneously when the uncertainty range is too large, leading to violations in some of the quality requirements.
In these cases, the traditional robust model cannot obtain feasible solutions, or determine which requirements would be violated under this circumstance. To deal with this situation, 255 we build a bi-criterion robust MOS model that allows violations of quality requirements to keep the model feasible under high uncertainty, and to quantify the violations. The information of quantified violations should guide response actions in the real process.
In this section, we first introduce and derive the robust counterpart of our MOS model.
Then, we present the flexibility test problem. Finally, we discuss the robust MOS model 260 that allows constraint violations.

Robust counterpart
The basic idea of the robust counterpart is to guarantee that the solution is feasible for any realization in the uncertainty set, which is enforced with a deterministic equivalent of their worst case (Soyster, 1973). In order to avoid the solution to be overly conservative, a 265 budget parameter Γ, developed by Bertismas and Sim (Bertsimas and Sim, 2003), can be adopted in the formulation to constrain the uncertainty set.
The robust counterpart formulation is a constraint-wise method, which implies that every constraint containing uncertainty is treated independently. Therefore, the robust counterpart formulation for each quality requirements (Equations (23)-(25)) needs to be 270 derived separately. We use the element concentration constraint (23) to show the four steps involved in the derivation of the robust counterpart. The derivations for the other quality requirements (Equations (24)-(25)) are presented in Appendix A.

Defining uncertainty set
The uncertainty set is defined in Equation (1). The original set is a box uncertainty set, in which the worst case occurs when all uncertain parameters reach one of their extreme values at the same time. In the smelting process, this case corresponds to the situation in which the composition of every element in each concentrate deviates δ from its nominal value. This case is very pessimistic and it can cause the solution to be overly conservative.
To avoid the extreme case, a budget parameter Γ is introduced with the following idea: for certain element k, there can be at most Γ concentrates whose composition of k can reach the worst-case values simultaneously. By normalizing the deviation of f c,k as p c,k , the budget parameter can be implemented in the uncertainty set by requiring that the sum of normalized deviations be less than Γ for each k: The objective of the problem is to maximize the gross margin of the process, which is a linear function of the variables M c i,v,c that are restricted by the element concentration constraint (23). Therefore, the worst case is that the variables in the objective functions (M c i,v,c ) are most constrained by the uncertain parameters, resulting in a lower objective value.

285
The worst cases can be achieved by maximizing the positive terms on the left hand side of the inequality (23). The corresponding worst case is presented in Equation (34).

Building auxiliary optimization problem
The worst case in Equation (34) transforms the model into a bi-level optimization problem. To make the model tractable, the maximized term in the inequality above can be reformulated by building an auxiliary problem. Considering the maximization term in (34), the reformulation for every individual inequality (i ∈ I, k ∈ K U ) is given by (35):

Dualizing auxiliary problem
The auxiliary optimization problem can be dualized to transform the model back to a single-level problem. For the auxiliary problem (35), with subscripts i ∈ I, k ∈ K U , the corresponding maximization problem is equivalent to: where s i,c,k and q i,k are dual variables corresponding to the first and the second constraints of problem (35), respectively. By strong duality, the dual problem is bounded and the optimal objective of the dual problem equals the optimal objective value of the auxiliary problem (35).
Substituting the dualized problems into Equation (23), the model will automatically minimize s i,c,k , q i,k from the dualized problem. The robust counterpart formulation is then given by Equation (37).
Compared with the original constraint (23), the uncertainty is taken into consideration Similarly, the robust counterparts for Equations (24) with subscripts i ∈ I, (k, k ) ∈ KK , k, k ∈ K U are: The robust counterpart for Equation (25) with subscripts i ∈ I, k ∈ K U is:

Flexibility test problem
As presented above, the robust counterpart dualizes the uncertain parameters in the constraints to transform the bi-level optimization problems into single-level ones. In the final for-300 mulation, the uncertain parameter f c,k is replaced with dual variables q i,k , q i,k,k , s i,c,k , s i,c,k,k , and there are no direct ways to relate these new variables to values of the uncertain parameters. In order to reveal which parameters are selected to fluctuate in the robust solution, we introduce the flexibility test problem to complement the robust counterpart.
The flexibility test problem was developed to quantify the capability of a system to   For the given uncertainty range ±10% for k , the fluctuation of the ratio is so large that any value within [0.58, 0.64] will violate at least one of the bounds in the extreme cases.
In other words, the violation of constraint (24) will always occur for the given uncertainty range in the robust model. When such a case occurs, the model will be infeasible, since there is no solution that satisfies both bounds for the ratio simultaneously. To make the problem tractable, and to quantify the violations, the robust counterpart can be further modified.
We introduce a violation term, ε, to the robust counterpart (37)-(40) to allow violations of these terms: Although unavoidable, the total violation of quality requirements should be reduced as much as possible so that the quality of products are less affected. This can be formulated as a second objective function: where ε 1 stands for the 1 norm of the violation terms, The new objective along with the modified quality requirements can be used to formulate a bi-criterion robust MOS model: The bi-criterion model can be reformulated as a single-objective problem using the εconstraint method (Cohon, 1978): in which ε 1 is the maximal value of ε 1 that can be obtained by solving the bi-criterion model only with the objective of maximizing the gross margin, and w is the coefficient that determines the upper bound of ε 1 . It is obvious that the gross margin is positively 340 correlated to to the violation term ε 1 , as the smaller violation indicates that the solution is more restricted by the quality requirements.
By introducing ε, the bi-criterion robust model is capable of handling large range of uncertainty by relaxing the quality requirements. In addition, a Pareto-front between the gross margin and the violations can be obtained by adjusting the value of w, indicating the 345 trade-offs that can be achieved between these two criteria.

Remarks
1. The constraints to which the robust counterpart formulation is applied are all linear.
Although the deterministic MOS model is an MINLP model, the quality requirements where the robust counterpart is implemented are all linear constraints. 2. The robust counterpart formulation does not introduce any integer variables or nonlinearity. The increase of model size is moderate, since it only introduces dual variables, whose sizes are O(|I| · |C| · |K U | 2 + |I| · |C| · |KK |).
3. The 1 norm of the violation terms can be further modified to weigh each violation term to indicate the total economic cost of dealing with these violations.

Two-step MILP-NLP decomposition strategy
Solving the robust counterpart formulation directly using MINLP solvers is computationally expensive due to its size. In previous works (Mouret et al., 2009(Mouret et al., , 2011Song et al., 2018), a two-step MILP-NLP decomposition strategy was developed to tackle such problems.

360
Its scheme is outlined in Figure 9. Because the MILP subproblem is a relaxation of the original maximization problem, its optimal objective value yields an upper bound for the original problem. With the binary variables fixed, the locally optimal objective value of the NLP subproblem yields a lower bound for the original problem. The final solution is the solution of the NLP subproblem.
The method provides an optimality gap, which is the relative difference between objective values of MILP and NLP solutions.
It is worth noting that feasibility and optimality of the final solution are not guaranteed with this approach. By passing the fixed assignment variables, the sequences of operations 375 determined in the MILP subproblem is fixed in the NLP subproblem. These sequences may be infeasible to the NLP subproblem if the nonlinear composition constraints (18)  with small or no optimality gaps with short computational time (Mouret et al., 2009(Mouret et al., , 2011Song et al., 2018).

Enhancement of decomposition strategy
The original MILP-NLP decomposition strategy can be further enhanced to reduce the optimality gaps. By comparing the MILP subproblem and the NLP subproblem in Fig-385 ure 9, the presence of concentration constraints (18)  One way to enhance the decomposition strategy is to maintain the concentrate types before and after blending to eliminate improper MILP solutions like in Figure 10. We introduce a new binary variable Y i,v,c to indicate if concentrate c is transferred in operation v in slot i: in which ε is a sufficiently small parameter to enforce positive values.
The new binary variable has the following relationship with the original binary variable The concentrate types arriving and leaving the blending unit can be related by specifying Y i,v,c for these two parts. There are various ways to construct such relationships. Here, we assume that the transfer of blends to the smelter occurs exactly 2 slots after the transfer of concentrates to the blending unit, because "2 slots" is the minimum slots needed to transfer blends from the blending unit to the smelter (blending unit → bins → smelter). Then the concentrate types in the former operations should be identical to those in the latter operations. This relationship can be formulated in aggregated form as follows: We also assume that the amount of each concentrate, M c i,v,c , is the same between transfer of concentrates to the blending unit at slot i and the transfer of blends to the smelter at i + 2: After applying these new constraints, the previous example will have a new MILP solution For daily-arrival materials, a similar constraint can be formulated to maintain their amounts between inlet and outlet operations of bins. We assume that the transfer of dailyarrival materials to the smelter is exactly one slot after the transfer to the bins (minimum slots needed), then the amount of c transferred to the bins at i should be no less than its amount transferred to the smelter at i + 1: The new assumptions for constraints (49)-(51) together imply a specific transfer pattern for the blending unit: the unit can be filled with new concentrates only when it is empty.
As the same amount of concentrates are required to be discharged from the blending unit to the smelter, there will be no concentrate left in the blending unit after the outlet operations 415 and before the next inlet operations. This pattern eliminates some possible ways of transferring blends from the blending unit to the smelter. This may affect the objective value for the MILP subproblem, but will reduce optimality gap and the computational time (see Section 7).
The scheme of the enhanced solution strategy is outlined in Figure 13 where the new 420 linear constraints (47)-(51) are added to the MILP subproblem.

Case studies
We address an illustrative example and an industrial case study to evaluate the performance of the MOS model, the robust formulations and the decomposition strategies. The topology of the illustrative example is presented in Figure 14. This example considers 6 425 types of concentrates and 4 key elements within a 10-day horizon. The process consists of 2 maritime vessels, 5 piles at the port, 1 blending unit, 1 daily-arrival pile, 3 bins, and 1 smelter. To simplify the model, interdependency constraints (24) and (25)  The topology of the industrial case study is shown in Figure 15. The industrial case study considers 14 types of concentrates with 9 key elements and a time horizon of 15 days.
The process consists of 5 maritime vessels, 11 piles at the port, 1 blending unit, 3 dailyarrival piles, 6 bins, and 1 smelter. All daily-arrival materials in piles r 25 , r 26 and r 27 can 435 be discharged to bins r 21 , r 22 , and r 23 . The data for the process and the set of cliques are presented in Tables C.3   fraction values. For the illustrative example, the uncertainty elemental set is K U = {k 1 , k 2 } and the budget parameter Γ is set to 2. For the industrial case study, K U = {k 1 , k 2 , k 3 } and Γ is set to 5. The equations involved in each case and each model are shown in Table 1.
The two cases are solved on an Intel Core i7 (6 cores) 2.60 GHz processor with 16 GB   Table 2. The model using the enhanced solution strategy has a significantly larger number of discrete variables as a result of introducing the binary variable Y i,v,c . The number of equa-455 tions also increases, as the linear constraints (47)-(51) are added to reduce the optimality gap. However, larger model size does not necessary lead to longer computational time (see Table 3). The enhanced solution strategy yields the same objective values in the MILP and NLP subproblems. The original solution strategy obtains a higher objective value in the MILP  Table 4. The number of continuous variables and equations in the model is slightly higher than the ones in the deterministic model because of the addition of the robust counterparts.

475
The flexibility test problem is a small MILP problem compared with the main MOS model. The computational results for the robust model are presented in Table 5. The enhanced solution strategy manages to obtain objective values with zero-gaps for all cases with different values of δ. The original solution strategy achieves the objectives with optimality gaps around 7%. The enhanced solution strategy also yields smaller computa-480 tional time for all cases compared with the original solution strategy. The MILP subproblems still take most of the total computational time, whereas the flexibility test problem usually takes less than 0.1 s to solve. The objective profile is outlined in Figure 16. The enhanced solution strategy always yields a higher gross margin than the original solution strategy. In addition, the gross margin tends to decrease monotonically with the increase of the given 485 uncertainty bounds, because the higher uncertainty makes the feasible region of the model smaller. The values of element composition obtained from the flexibility test problem for the case with δ = 0.1 are given in Table 6

Industrial case
The robust MOS model is first solved for the industrial case using the enhanced solution 495 strategy. The model statistics are given in Table 7. The model is much larger than the   Table 8. Similar to the results of the robust model for the illustrative example, for δ = 0 ∼ 0.03, the objective value decreases monotonically with the increase of δ. Interestingly, the 500 enhanced solution strategy also yields zero-gaps for these cases. The uncertain case with δ = 0.01 has a significantly larger computational time compared with the deterministic case.
However, as δ increases, the computational time drops monotonically. This may be because the introduction of robust counterparts makes the solution process much harder at first; then, the larger δ requires the solution to be more conservative to be robust to larger ranges   Figure 17. The Pareto-front shows the trade-off between gross margin and 515 constraint violations. The Gantt chart for the case with δ = 0.03 is shown in Figure 18;

Conclusions
In this work, we have addressed the optimal scheduling of logistic and blending operations for copper smelting process so that the gross margin associated to the concentrates being processed in the smelter are maximized. This problem is very difficult to solve due to the complexity of the process operations rules and the quality requirements for products,

Appendix A. Derivation of robust counterpart
The robust counterpart formulation is constraint-wise, i.e. it is different for every constraint containing uncertainty based on the uncertain terms. The derivation of robust counterparts for constraints (24) and (25) are given here.

555
Appendix A.1. Defining uncertainty set This step is the same as the one shown in the paper. By introducing the budget parameter, the uncertainty set can be modified as:

. Building worst cases
The worst cases can be obtained by maximizing the positive terms and minimizing the negative terms on the left hand side of the constraints. The corresponding worst cases are: (i ∈ I, (k, k ) ∈ KK , k, k ∈ K U ), the corresponding maximization problem is equivalent to: It is worth noting that max For the minimized term min , the corresponding problem is equivalent to:

. Dualizing auxiliary problems
The auxiliary optimization problem can be dualized to transform the model back to single-level problem. The auxiliary problem (A.4) has the same dualized problem as (35).
For every problem with subindices (i ∈ I, k ∈ K U ), the corresponding maximization problem is equivalent to: For the auxiliary problem (A.5) with subindices (i ∈ I, k ∈ K U ), the corresponding maximization problem is equivalent to: where s i,c,k,k and q i,k,k are auxiliary variables of the dualized problem. By strong duality, the optimal objective value of the dual problem is equal to the optimal objective value of (A.5).