# Simulation-Based Multiobjective Optimization of Open-Pit Mine Haulage System: A Modified-NBI Method and Meta Modeling Approach

Authors: Milad Abolghasemian; Armin Ghane Kanafi; Maryam Daneshmand-Mehr
Journal: Complexity
 (2022)

---

## Abstract

A large number of engineering problems involve several conflicting objectives, which today are often solved through expensive simulation calculations. Methods based on meta-models are one of the approaches to solving this group of problems. In this paper, multiobjective optimization in the extraction system of a copper open-pit mine complex is presented by the modified-NBI optimization method and regression meta-model. For this purpose, two objective functions of maximizing the amount of total extraction, which is the sum of the extraction of sulfide, oxide, low-grade ores, and waste in this mine, and minimizing the transport time of haulage according to the limitation of its storage capacity, transport equipment, and budget are considered. The Central Composite Design (CCD) method is used to build the Design of Experiments (DOE) for the design variables. The considered design variables are the number of trucks of 120 tons, 240 tons, 35 tons, and 100 tons. The number of targets considered in each design combination is considered the response surface. The suitable meta-model to maximize the total extraction rate and minimize the transport time of the haulage, two modified functions of nonlinear regression have been determined. The accuracy of the models for selection has been done using PRESS andR2 statistics. The most common PRESS error has also been used to validate the meta-models. Then the multiobjective optimization problem was solved using the modified-NBI method. Finally, Pareto and optimal solutions using the proposed approach were presented and discussed.


---

## Body

## 1. Introduction

Mines have been considered one of the costliest and complicated industries for many years, and various studies have been done on different parts of this industry such as geology, drilling planning, and operational processes [1]. Without a doubt, the proper exploitation of the country’s mines is considered an important and positive factor in economic growth and development [2]. Because each ton of copper ore loaded in trucks is worth nearly 100 thousand dollars [3, 4]. Mines contain several uncertain parameters that make their modeling by traditional techniques very complicated. Simulation models are a powerful tool for solving estimation problems that can create flexible models for systems without considering many assumptions [5, 6]. Considering today’s competitive world, companies emphasize finding ways to produce products faster, cheaper, and more effectively. Therefore, the use of simulation techniques is increasing to investigate system behavior and design effects on system performance [7]. In the real world of engineering design, optimization processes are often performed with more than one objective, which is called multiobjective optimization [8, 9]. Multiobjective optimization for an engineering problem involves several extensive evaluations of each objective in the design space, which leads to a large number of simulation runs, each run requiring hours of computation to find the optimal solution through optimization-based simulation. Despite the high execution time required for simulation calculations, it can be accepted that all computational time costs should be spent on simulation based on multiobjective optimization to find a set of optimal solutions. This is while all the results obtained from an optimization implementation may not be acceptable, in which case it is necessary to make corrections such as changing the formulation, parameters, and constraint [10, 11].One of the strategies for quick and accurate estimation of complex and expensive models is the simulation-based optimization meta-modeling approach [2]. The meta-models develop a relationship between the input variables and the response level to predict the simulation calculation model [12]. Meta-models are mathematical estimation models for simulation models [13, 14]. Various meta-model methods have been developed to solve optimization problems based on meta-models. Response surface method (RSM) or polynomial regression [15], Kriging [16], and artificial neural networks [17] are several known functions. Various studies have compared different meta-models in terms of accuracy, efficiency, stability, and effectiveness. But by reviewing the literature, we can conclude that there is no specific method superior to other techniques in terms of performance, and the choice of meta-model is chosen arbitrarily [18, 19]. Usually, in related studies, low-order polynomials such as quadratic polynomials are used, in which the unknown coefficients are obtained by minimizing the error of the residuals between the fitted values and the value of the objective function [20]. The response surface method is a set of statistical and mathematical methods that can optimize probabilistic functions such as simulation models. Recently, the response surface method has been widely used in the field of engineering to design a new product or redesign a product or develop a new product [7]. For example, Dengiz et al. [5], by presenting a response surface meta-model based on analytical modeling to increase productivity in the automotive industry in Turkey, were able to increase the daily production rate by 15%. Amouzgar et al. [12] considered a potential advantage for Meta model-based multiobjective optimization in machining operations. The multiobjective optimization method based on the meta model has been useful in reducing the calculation time in this study. In addition, it can find more infinite points as a solution than other existing methods.Therefore, the tendency to use meta-models in multiobjective optimization is very important. Because, in engineering problems, generally more than one goal is considered, and considering that the goal functions conflict with each other, there is no optimal solution for them, but instead, a set called Pareto solutions. It seems that multiobjective optimization based on a meta-model is an effective approach both in multiobjective optimization and in the design of complex products, whose main goal is to determine a suitable functional relationship between input and output in the system. Therefore, in this paper, for the simulation-based optimization framework, an estimated function is substituted for the complex simulation model. Therefore, the main contribution of the paper is as follows:(i)
Presenting a comprehensive framework for multiobjective simulation optimization based meta-modeling in an open-pit mine,(ii)
Determining an optimal production plan through effective haulage equipment control,(iii)
Controlling the duration of using effective haulage equipment in open-pit mine.The paper is organized as follows: Section2 presents literature review. Section 3 presents research methodology. For this purpose, the definitions, concepts, and details of the modified-NBI, multiobjective optimization method, formulated problem structure, and the meta-modeling method are explained. Section 4 presents statistical analysis, optimization of the mathematical model, and sensitivity analysis. Section 5 presents managerial insight. Finally, Section 6 contains conclusions and some suggestions for future research studies.


## 2. Literature Review

Multiobjective optimization is one of the attractive research fields in the branch of optimization methods, especially the use of interactive methods. However, a small number of researches in the literature on multiobjective optimization based on the simulation of interactive algorithms has been done. This is although using a lot of evolutionary algorithms has been seen. For example, Syberfeldt et al. [21] presented an evolutionary algorithm-based method for simulation-based multiobjective optimization in a manufacturing problem to improve cell manufacturing in VOLVO in Sweden. The results have shown that by using simulation and evolutionary algorithms, it is possible to increase the amount of cell usage and also reduce the delay components. In another study, Syberfeldt et al. [22] presented a simulation-based multiobjective optimization using evolutionary algorithms for the personnel planning system of the post office in Sweden. The purpose of this study is to determine the best work schedule for personnel to reduce working time and administrative work pressure. The NSGA-II algorithm is used for multi-objective optimization. The results of this research show that the algorithm can be easily implemented in optimization. Moussavi et al. [23] presented an integer multiobjective programming problem to implement an ergonomic work cycle in a truck assembly production system. The main goal of this study is to balance the workload of workers and reduce the production cycle time in the study. This model is programmed using Goal Programming and solved using the Gurobi algorithm. The results show that the proposed model can optimize for both purposes. Amouzgar et al. [12] provide an effective framework for multiobjective optimization for metal cutting machining processes. The goal of multiobjective optimization is to minimize the tool-chip temperature and wear depth while maximizing the removal rate. In this study, by performing a knowledge discovery and data weighting style, the nondominated solutions are analyzed using data mining techniques to gain a deep understanding of the metal cutting process. Das and Pratihar [24] presented an approach to increase the accuracy of the solutions of multiobjective optimization evaluation algorithms. In this study, after obtaining a set of Pareto points using a weighted multiobjective evaluation algorithm, it is used in a neural system. Then, using this neural system, modified Pareto solutions are obtained. The presented algorithm provides analysts with valuable information for analyzing engineering problems. Karmellos et al. [25] compared multiobjective optimization frameworks for the design of energy distribution systems under uncertainty. For this purpose, they presented two multiobjective models for the design of energy system distribution to identify places in need of heating, cooling, and electricity, taking into account uncertainty parameters such as energy price, interest rate, solar radiation, wind speed, and energy demand. The results of the research show that by using this method, the decision maker can make an informed choice to determine the energy distribution system under conditions of uncertainty. Russell and Taghipour [26]; a new solution method using multiobjective optimization is presented to solve the complex scheduling problem in low-volume production systems. For this purpose, using integer multiobjective linear mathematical programming models, the scheduling problem in low-volume production systems has been modeled. The models presented in this paper have been used for compatibility in the real world in a case study in the aerospace industry, through which the reliability of the models is confirmed. Zhang et al. [27] have used multiobjective optimization to determine concrete mix ratios with several objectives and under nonlinear constraints. In this study, an optimization method based on machine learning using metaheuristic algorithms is presented. The results show that the multiobjective optimization model can help as a design guide to facilitate decision-making before the construction stage.As mentioned above simulation-based optimization is an efficient method. The idea of optimization based on simulation is presented to find optimal or near-optimal solutions. Choosing a suitable approach for optimization depends on the characteristics of the problem [28]. In terms of choosing the approach, we classify optimization based on simulation into two types. The first type is a common optimization that generally considers one or more objective functions based on several constraints that can be linear or non-linear. For example, Dengiz and Belgin [7] presented a simulation-based optimization for a painting production line in the automotive industry using response surface methodology. The Meta-model estimated in this study can reduce the deviation in results and current costs in the system. Shishvan et al. [28] have presented a new approach for simulation optimization to solve the problems of transportation and job-shop scheduling. Based on the obtained results, the quality of the obtained solutions increases compared to other considered algorithms. Burak and Kumral [29] presented a simulation-based optimization for a truck-shovel system in an open pit mine. For this purpose, aimed to maximize the use of the truck-shovel system. This approach has a good ability to increase the productivity of the truck-shovel system. Based on this approach, the movement of materials in the system increases by 6 k tons. Jahangiri et al. [13] presented a simulation-based optimization approach to evaluate the emergency department of a public hospital in Iran during the pandemic of COVID-19. By using variables influencing the flow of patients’ admission, they determined the optimal combination of resources to obtain the minimizing waiting time for patients. Moniri-Morad et al. [30] developed a simulation-based optimization algorithm to determine the most optimal handling equipment by considering influencing factors such as availability and maintenance analysis, production scheduling, material flow rates, and random environmental and operational phenomena. The proposed method is used to size the transportation fleet in one step by developing a parallel combination of mixed integer programming and discrete event simulation. Finally, the proposed approach has been implemented in the Sungun copper mine complex in Iran. The second type is hierarchical optimization based on simulation. Generally, these types of issues are divided into two levels. At the high level, the main goal is considered, and at the low level, the number of influencing variables to achieve the high-level goal is considered. For example, Nageshwaraniyer et al. [3] presented a simulation-based two-level hierarchical optimization framework for time scheduling in a coal mine. At the top level, the direct flow of coal from the pit to the trains is considered, and the problem of scheduling the machines is solved at this stage. At the low level, using OptQuest®, an optimization problem has been solved to determine haulage variables such as trains, trucks, and conveyors. Based on the obtained results, the travel and loading time of trucks has decreased and the rate of using machines has increased. In Table 1, above mentioned literature categorized.Table 1
Literature review.
AuthorOptimizationDesign of experimentMeta-modelObjective functionMeta-heuristicExactLHS2kCCDANNKrigingRegressionMultiSingleSyberfeldt et al. [21]∗∗Dengiz and Belgin [7]∗∗∗∗Syberfeldt et al. [22]∗∗Moussavi et al. [23]∗∗Amouzgar et al. [12]∗∗∗Shishvan et al. [28]∗∗∗Pratihar [24]∗∗∗Karmelos and Mavrotas [25]∗∗Rassel and Taghipour [26]∗∗Zhang et al. [27]∗∗Jahangiri et al. [13]∗∗∗Moniri-Morad et al. [30]∗This research∗∗∗∗


## 3. Research Methodology

### 3.1. Problem Statement

Sarcheshmeh open-pit copper mine complex is located in the Kerman Province, southeast of Iran. Sarcheshmeh is a large open copper mine, considered to be the second largest copper deposit worldwide. It is located at 65 kilometers off the southwest of Kerman city and 50 kilometers from south of Rafsanjan. The average altitude of the region is about 2600 m, and the highest spot is approximately 3000 m. The extracted deposits can be categorized into four groups:(1)
Sulfide ore (grade of copper greater than 0.7%)(2)
Oxide ore (grade of copper between 0.25% and 0.7%)(3)
Low grade ore (grade of copper between 0.15% and 0.25%)(4)
Waste (grade of copper less than 0.15%)The proportion of the amount of these rocks to the whole amount is 45%, 5%, 44%, and 6%, respectively. Based on the different kinds of ore, a transportation strategy and the way of storing ores would be chosen. The first kind of mineral substance, sulfide ore, is transferred to a crusher station. There is a crusher machine with a capacity of 60,000 tons per day. Then, the substance is moved to harp copper storage with the capacity of 150,000 tons. After harping, the substance is stored in a soft copper storage. Oxide ore, low grade ore, and waste are transferred to their respective dumping station. The conceptual model of the haulage system at Sarcheshmeh copper complex is demonstrated in Figure1.Figure 1
The conceptual model of the haulage system at sarcheshmeh copper mine.In order to transport ores, a number of trucks are assigned to the loading station. The mineral substance is loaded onto a truck with a shovel, and when the truck is filled, it is led to a dump. The main purpose of this research is determining the optimal number of haulage system equipment especially number of trucks in order to maximize sulfide ore and maximize loaded ores in the trucks according to equipment and storage capacity and budget. The key resources in the Sarcheshme copper mine are as follows:(1)
Truck 120 TonsX1(2)
Truck 240 TonsX2(3)
Truck 35 TonsX3(4)
Truck 100 TonsX4Currently, Sarcheshmeh open copper mine has nine trucks of 35 tons, 36 trucks of 100 tons, 20 trucks of 120 tons, and two trucks of 240 tons. Trucks which are used to transfer oxide ore and wastes are varied between 35 tons and 100 tons. Moreover, 120 tons to 240 tons’ trucks are used to move sulfide and low-grade ore. It is possible to assign each shovel to every kind of ores. The hourly operating costs of trucks are shown in Table2.Table 2
Hourly operating cost of trucks (cost unit) [1].
Truck (tons)DepreciationOverheadRepairMaintenanceGasLubricantsOtherTotalSpare partsSalarySpare partsSalary35110102103422100211215115945120321417118156915030141711816672404325110221143118In this paper, the application of the simulation modeling approach using the Arena software® to model the haulage system in the Sarchesmeh copper open-pit mine, which was developed by Eskandari et al. [1] is considered. In the developed model for all cases, there is no unacceptable difference between the results at the 95% confidence level. So we conclude that the model was built correctly. To run the model, it is first necessary to determine the simulation parameters such as length and number of repetitions. The working schedule of the haulage system in the Sarcheshmeh copper mine is 24 days a month. Therefore, the length of each repetition is 1 month. The number of repetitions is determined by considering the half-width of trucks as the main value for measuring system performance. The results show that the number of repetitions should be considered 10. To achieve stability in the system, a warm-up period is considered. Experiments show that after four days of the warm-up period, the performance of the system shows a steady state. In Table 3, we show the lower bound, the upper bound, and the current bound of trucks. The mentioned settings are included in the model and implemented on a personal computer with Intel Core i3 1.8 GHz CPU specifications and 4 GB RAM.Table 3
Number of trucks combination at each type.
TypeCombination (lower bound, current, upper bound)Truck 120 tons(12, 20, 28)Truck 240 tons(3, 4, 5)Truck 35 tons(9, 15, 25)Truck 100 tons(25, 36, 45)Based on the previous information, the managers of the Sarcheshmeh copper mine complex tend to optimize the combination of key haulage resources based on two objective functions: maximizing the total extraction amount, which is the sum of the extraction amount of sulfide, oxide, low-grade ores, and waste in this mine and minimizing the travel time of the haulage. Find the relocation according to the storage capacity limit and the haulage and budget considered. The optimization problem is mathematically formulated as follows:(1)maxf1X1;X2;X3;X4,minf2X1;X2;X3;X4,subject to,∑i=14cixi≤B,∑i=14ci′xi≤C,Li≤xi≤Uifori=1;2;3;4,xiinteger.Equation (1) is an integer multiobjective optimization problem. The functions of this problem are unknown and we do not have an analytical mode. They must be evaluated through simulation according to the proposed framework. ci is the cost of each truck. B is the total available budget. ci′ is the capacity of each key resource. C is the total storage capacity in the system. Li and Ui are respectively the lower and upper bounds of resources in the mine complex.


### 3.2. Preliminary Definitions

In this subsection of the research, the preliminary definition of the concepts that will be used in the following is considered.Definition 1.
Multiobjective optimization problem (MOP)
A multiobjective optimization problem (MOP) is shown in the following equation:(2)minf1x,…,fpx,p≥2,s.t,gx≤0,hx=0,xil≤xi≤xiui=1,…,N,where, F:RN⟶Rp;h:RN⟶Reandg:Rn⟶Ri are twice continuously differentiable mappings and xil∈R∪−∞N, xiu∈R∪∞N, N being the number of variables, p the number of objectives, and e and i the number of equality and inequality constraint.Definition 2.
Pareto Set
If none of the objective functions can be improved by a feasible solution without worsening at least one of the other objectives, then a nondominated or Pareto optimal solution is feasible in the design space. The set of feasible solutions that are nondominated is also called the Pareto optimal or nondominated set. If there is a solution that does not belong to this set, it is called a dominated solution.Definition 3.
Convex hull of individual minima (CHIM)
Letxi∗ be the respective global minimizers of fix,i=1,…,n over x∈C. Let Fi∗=Fxi∗,i=1,…,n. Let ϕ be the n×n matrix whose ith column is Fi∗−F∗ sometimes known as the pay-off matrix. Then the set of points in Rn that are convex combinations of Fi∗−F∗, ϕβ:β∈Rn;∑i=1nβi=1,βi≥0 is referred to as the convex hull of individual minima [31].Definition 4.
Meta-Model
A meta-model or surrogate model is a mathematical approximation of a simulation model. Therefore, meta-model is an abstract model for simulation model.Definition 5.
Validation of meta-model
Validation means whether the model is designed correctly or not. There are various methods to provide the validity of the meta-models. In this paper, we use the most common PRESS error is the root mean square PRESS denoted asRMSEPRESS Calculated by PRESS/n where, n is number of test points selected to evaluate the model. It is obvious that a value of zero for RMSE is the optimal desired value.


### 3.3. Solution Approach

#### 3.3.1. Meta-Modeling Approach

The main elements of the proposed framework are meta model-based optimization [32] in the form of identifying the shape of the meta model, designing experiments to adapt the Meta model, performing simulation experiments, fitting the Meta model, and verifying its accuracy, and optimization considering the meta model in the problem. The algorithm for finding a suitable meta model is shown below. We have used this algorithm for multi-objective optimization based on simulation, the steps of which are fully presented in the next sections. (Algorithm1)Algorithm 1: Meta modelling algorithm.
Input: Objectives, Variables and ParametersOutput: Meta Model and a non-dominated set//Start//Step1: Simulation Model DevelopmentStep 2: Design of Experiment (DOE)Step 3: Run Simulation Model ExperimentStep 4: Fit the Meta ModelStep 5: Meta Model ValidationStep 6: Apply optimization algorithmStep 7: Quality of Results//Finish//In the first step of the algorithm, a discrete event simulation model is developed. If the built model has the necessary validity, we go to the next step; otherwise, the model is accompanied by modifications to obtain the necessary validity. In the second step, a suitable design for an experiment is developed. In the third step, the scenarios designed in step 2 are implemented in the simulation model developed in step 1 to determine the dependent variable of the model. In the fourth step, the best meta model is selected and its unknown coefficients are determined by performing statistical analysis. In this step, we should run the simulation model to determine the response surface for fitting the meta-model. Then, from the data, we obtain an approximate value for the parameter value of the meta-model. Finally, we evaluate these estimates using mathematical and statistical criteria. The fifth step is the answer to the question of whether the meta model built in the fourth step can sufficiently predict the performance of the system or not. In the case of lack of validity, we go to step 2 and change the design of the experiment or go back to step 4 and change the type of meta model. In the sixth step, by applying management constraints, a set of nondominant solutions for the multiobjective optimization problem is obtained using the modified-NBI method. In the seventh step, we will compare the obtained results with the existing situation to identify the improved level.


#### 3.3.2. Modified Normal-Boundary Intersection (Modified-NBI) Method

The first goal of the modified-NBI method is to determine the Pareto frontier by solving an optimization problem for multiobjective optimization problems with a continuous or piecewise Pareto frontier. This goal is accomplished in two stages in this method. The first step is to use the modified CHIM in the optimization problem compared to the original NBI algorithm. The second step is to control the iterations of the optimization problem, which are solved by the modified CHIM. According to the modified-NBI algorithm, it is necessary to normalize the objective functions. Therefore, all objective functions have a minimum value of zero and a maximum value of one. If the objective function is unlimited or its maximum value cannot be determined, the user can impose an upper limit value on the algorithm. The first step in the modified-NBI is to use the modified CHIM, which is described below [33]. The second step of the modified-NBI method is to control the iterations while solving the optimization problem to obtain the Pareto frontier. One of the methods used to solve the optimization problem is the quasi-Newton method. Using this method, a relative minimum is obtained for a multi-objective optimization problem. This method has differences compared to meta-heuristic method such as genetic algorithm (GA) that are based on primary populations such as genetic algorithm (GA). For example, a GA requires the selection of parameters such as population size, type of crossover and probability, mutation probability, and number of generations. Also, due to the random nature of population-based methods, even with the same settings, similar answers are often not produced. Hence, these algorithms are often run many times to obtain a reliable set of solutions. The modified-NBI optimization method involves choosing only one parameter to influence the number of generated Pareto points. In addition, similar results are produced each time. Therefore, the algorithm does not need a large number of execution times.Let us show through mathematical formulation how any such boundary point can be found by solving an optimization problem. Givenβ, that ϕβ shows a point in the CHIM. Suppose that n^ represent the normalization unit to the CHIM simplex pointing towards the origin. Then, ϕβ+tn^,t∈R shows the set of points on that normal. Therefore, normal points and boundary of F closet to the origin are the global solution according to the following problem:(3)maxt,s.t,ϕβ+tn^=Fx,gx≤0,hx=0,xil≤xi≤xiui=1,…,N.The vector constrainϕβ+tn^=Fx ensures that the point x is actually mapped by F to a point on the normal. While the remaining constraints ensure feasibility of x. Instead of ϕβ+tn^=Fx we can use Fx≤u+tv. Where, v is normal vector and u is point of origin of the normal which are user-defined values. The general algorithm of the modified NBI optimization mathematical method is written below. (Algorithm2)Algorithm 2: Modified-NBI algorithm.
Input: Number of significant digitsNumber of initial pointsOutput: non-dominated set//Start//Step 1: Normalize objective FunctionsStep 2: Generating Pareto PointStep 3: Initiate objective function optimizationStep 4: Accurate Pareto setStep 5: Fix objective function ValuesStep 6: Remove Pareto dominated PointsStep 7: Optimal Stage//Finish//In the first stage, the objective functions are placed between the minimum value of zero and the maximum value of one. Therefore, we determine the objective functions in the interval0;1. This causes the Pareto frontier to be placed inside a bound box. In the second step, a V is selected for each pair in this step to generate the space of Pareto points. In the third step, the first optimization problem with the starting minimum point for f1;t=1;β=0 starts. If at the end of optimization t=0;β=1, then go to step 4. Because in this case, a Pareto point has been estimated. Otherwise, use t+V and β+V as starting points for the next optimization, and this step is repeated. In the fourth step, if the estimated Pareto set needs more accuracy, a smaller Vm is used to generate more Pareto points and we return to step 3. Otherwise, we go to step 5. In the fifth step, for N > 2, we use different values in [0, 1]in order to construct the values of objective functions to determine the results of multiple objectives. Then we return to step 1 and determine the objective functions in the interval 0;1. If all combinations are determined, we stop. Otherwise, go to step 6. In the sixth step, based on the described filter [34], the set of inferior Pareto points is removed from the set of generated Pareto points. In the seventh step, choose the best solution from the Pareto set as the optimal solution. For this, choose a range for V for more precision of the generated Pareto set and repeat steps 2 to 6 for all V values in the range. In Figure 2, research solution approach framework are shown.Figure 2
Research solution approach framework.


## 3.1. Problem Statement

Sarcheshmeh open-pit copper mine complex is located in the Kerman Province, southeast of Iran. Sarcheshmeh is a large open copper mine, considered to be the second largest copper deposit worldwide. It is located at 65 kilometers off the southwest of Kerman city and 50 kilometers from south of Rafsanjan. The average altitude of the region is about 2600 m, and the highest spot is approximately 3000 m. The extracted deposits can be categorized into four groups:(1)
Sulfide ore (grade of copper greater than 0.7%)(2)
Oxide ore (grade of copper between 0.25% and 0.7%)(3)
Low grade ore (grade of copper between 0.15% and 0.25%)(4)
Waste (grade of copper less than 0.15%)The proportion of the amount of these rocks to the whole amount is 45%, 5%, 44%, and 6%, respectively. Based on the different kinds of ore, a transportation strategy and the way of storing ores would be chosen. The first kind of mineral substance, sulfide ore, is transferred to a crusher station. There is a crusher machine with a capacity of 60,000 tons per day. Then, the substance is moved to harp copper storage with the capacity of 150,000 tons. After harping, the substance is stored in a soft copper storage. Oxide ore, low grade ore, and waste are transferred to their respective dumping station. The conceptual model of the haulage system at Sarcheshmeh copper complex is demonstrated in Figure1.Figure 1
The conceptual model of the haulage system at sarcheshmeh copper mine.In order to transport ores, a number of trucks are assigned to the loading station. The mineral substance is loaded onto a truck with a shovel, and when the truck is filled, it is led to a dump. The main purpose of this research is determining the optimal number of haulage system equipment especially number of trucks in order to maximize sulfide ore and maximize loaded ores in the trucks according to equipment and storage capacity and budget. The key resources in the Sarcheshme copper mine are as follows:(1)
Truck 120 TonsX1(2)
Truck 240 TonsX2(3)
Truck 35 TonsX3(4)
Truck 100 TonsX4Currently, Sarcheshmeh open copper mine has nine trucks of 35 tons, 36 trucks of 100 tons, 20 trucks of 120 tons, and two trucks of 240 tons. Trucks which are used to transfer oxide ore and wastes are varied between 35 tons and 100 tons. Moreover, 120 tons to 240 tons’ trucks are used to move sulfide and low-grade ore. It is possible to assign each shovel to every kind of ores. The hourly operating costs of trucks are shown in Table2.Table 2
Hourly operating cost of trucks (cost unit) [1].
Truck (tons)DepreciationOverheadRepairMaintenanceGasLubricantsOtherTotalSpare partsSalarySpare partsSalary35110102103422100211215115945120321417118156915030141711816672404325110221143118In this paper, the application of the simulation modeling approach using the Arena software® to model the haulage system in the Sarchesmeh copper open-pit mine, which was developed by Eskandari et al. [1] is considered. In the developed model for all cases, there is no unacceptable difference between the results at the 95% confidence level. So we conclude that the model was built correctly. To run the model, it is first necessary to determine the simulation parameters such as length and number of repetitions. The working schedule of the haulage system in the Sarcheshmeh copper mine is 24 days a month. Therefore, the length of each repetition is 1 month. The number of repetitions is determined by considering the half-width of trucks as the main value for measuring system performance. The results show that the number of repetitions should be considered 10. To achieve stability in the system, a warm-up period is considered. Experiments show that after four days of the warm-up period, the performance of the system shows a steady state. In Table 3, we show the lower bound, the upper bound, and the current bound of trucks. The mentioned settings are included in the model and implemented on a personal computer with Intel Core i3 1.8 GHz CPU specifications and 4 GB RAM.Table 3
Number of trucks combination at each type.
TypeCombination (lower bound, current, upper bound)Truck 120 tons(12, 20, 28)Truck 240 tons(3, 4, 5)Truck 35 tons(9, 15, 25)Truck 100 tons(25, 36, 45)Based on the previous information, the managers of the Sarcheshmeh copper mine complex tend to optimize the combination of key haulage resources based on two objective functions: maximizing the total extraction amount, which is the sum of the extraction amount of sulfide, oxide, low-grade ores, and waste in this mine and minimizing the travel time of the haulage. Find the relocation according to the storage capacity limit and the haulage and budget considered. The optimization problem is mathematically formulated as follows:(1)maxf1X1;X2;X3;X4,minf2X1;X2;X3;X4,subject to,∑i=14cixi≤B,∑i=14ci′xi≤C,Li≤xi≤Uifori=1;2;3;4,xiinteger.Equation (1) is an integer multiobjective optimization problem. The functions of this problem are unknown and we do not have an analytical mode. They must be evaluated through simulation according to the proposed framework. ci is the cost of each truck. B is the total available budget. ci′ is the capacity of each key resource. C is the total storage capacity in the system. Li and Ui are respectively the lower and upper bounds of resources in the mine complex.


## 3.2. Preliminary Definitions

In this subsection of the research, the preliminary definition of the concepts that will be used in the following is considered.Definition 1.
Multiobjective optimization problem (MOP)
A multiobjective optimization problem (MOP) is shown in the following equation:(2)minf1x,…,fpx,p≥2,s.t,gx≤0,hx=0,xil≤xi≤xiui=1,…,N,where, F:RN⟶Rp;h:RN⟶Reandg:Rn⟶Ri are twice continuously differentiable mappings and xil∈R∪−∞N, xiu∈R∪∞N, N being the number of variables, p the number of objectives, and e and i the number of equality and inequality constraint.Definition 2.
Pareto Set
If none of the objective functions can be improved by a feasible solution without worsening at least one of the other objectives, then a nondominated or Pareto optimal solution is feasible in the design space. The set of feasible solutions that are nondominated is also called the Pareto optimal or nondominated set. If there is a solution that does not belong to this set, it is called a dominated solution.Definition 3.
Convex hull of individual minima (CHIM)
Letxi∗ be the respective global minimizers of fix,i=1,…,n over x∈C. Let Fi∗=Fxi∗,i=1,…,n. Let ϕ be the n×n matrix whose ith column is Fi∗−F∗ sometimes known as the pay-off matrix. Then the set of points in Rn that are convex combinations of Fi∗−F∗, ϕβ:β∈Rn;∑i=1nβi=1,βi≥0 is referred to as the convex hull of individual minima [31].Definition 4.
Meta-Model
A meta-model or surrogate model is a mathematical approximation of a simulation model. Therefore, meta-model is an abstract model for simulation model.Definition 5.
Validation of meta-model
Validation means whether the model is designed correctly or not. There are various methods to provide the validity of the meta-models. In this paper, we use the most common PRESS error is the root mean square PRESS denoted asRMSEPRESS Calculated by PRESS/n where, n is number of test points selected to evaluate the model. It is obvious that a value of zero for RMSE is the optimal desired value.


## 3.3. Solution Approach

### 3.3.1. Meta-Modeling Approach

The main elements of the proposed framework are meta model-based optimization [32] in the form of identifying the shape of the meta model, designing experiments to adapt the Meta model, performing simulation experiments, fitting the Meta model, and verifying its accuracy, and optimization considering the meta model in the problem. The algorithm for finding a suitable meta model is shown below. We have used this algorithm for multi-objective optimization based on simulation, the steps of which are fully presented in the next sections. (Algorithm1)Algorithm 1: Meta modelling algorithm.
Input: Objectives, Variables and ParametersOutput: Meta Model and a non-dominated set//Start//Step1: Simulation Model DevelopmentStep 2: Design of Experiment (DOE)Step 3: Run Simulation Model ExperimentStep 4: Fit the Meta ModelStep 5: Meta Model ValidationStep 6: Apply optimization algorithmStep 7: Quality of Results//Finish//In the first step of the algorithm, a discrete event simulation model is developed. If the built model has the necessary validity, we go to the next step; otherwise, the model is accompanied by modifications to obtain the necessary validity. In the second step, a suitable design for an experiment is developed. In the third step, the scenarios designed in step 2 are implemented in the simulation model developed in step 1 to determine the dependent variable of the model. In the fourth step, the best meta model is selected and its unknown coefficients are determined by performing statistical analysis. In this step, we should run the simulation model to determine the response surface for fitting the meta-model. Then, from the data, we obtain an approximate value for the parameter value of the meta-model. Finally, we evaluate these estimates using mathematical and statistical criteria. The fifth step is the answer to the question of whether the meta model built in the fourth step can sufficiently predict the performance of the system or not. In the case of lack of validity, we go to step 2 and change the design of the experiment or go back to step 4 and change the type of meta model. In the sixth step, by applying management constraints, a set of nondominant solutions for the multiobjective optimization problem is obtained using the modified-NBI method. In the seventh step, we will compare the obtained results with the existing situation to identify the improved level.


### 3.3.2. Modified Normal-Boundary Intersection (Modified-NBI) Method

The first goal of the modified-NBI method is to determine the Pareto frontier by solving an optimization problem for multiobjective optimization problems with a continuous or piecewise Pareto frontier. This goal is accomplished in two stages in this method. The first step is to use the modified CHIM in the optimization problem compared to the original NBI algorithm. The second step is to control the iterations of the optimization problem, which are solved by the modified CHIM. According to the modified-NBI algorithm, it is necessary to normalize the objective functions. Therefore, all objective functions have a minimum value of zero and a maximum value of one. If the objective function is unlimited or its maximum value cannot be determined, the user can impose an upper limit value on the algorithm. The first step in the modified-NBI is to use the modified CHIM, which is described below [33]. The second step of the modified-NBI method is to control the iterations while solving the optimization problem to obtain the Pareto frontier. One of the methods used to solve the optimization problem is the quasi-Newton method. Using this method, a relative minimum is obtained for a multi-objective optimization problem. This method has differences compared to meta-heuristic method such as genetic algorithm (GA) that are based on primary populations such as genetic algorithm (GA). For example, a GA requires the selection of parameters such as population size, type of crossover and probability, mutation probability, and number of generations. Also, due to the random nature of population-based methods, even with the same settings, similar answers are often not produced. Hence, these algorithms are often run many times to obtain a reliable set of solutions. The modified-NBI optimization method involves choosing only one parameter to influence the number of generated Pareto points. In addition, similar results are produced each time. Therefore, the algorithm does not need a large number of execution times.Let us show through mathematical formulation how any such boundary point can be found by solving an optimization problem. Givenβ, that ϕβ shows a point in the CHIM. Suppose that n^ represent the normalization unit to the CHIM simplex pointing towards the origin. Then, ϕβ+tn^,t∈R shows the set of points on that normal. Therefore, normal points and boundary of F closet to the origin are the global solution according to the following problem:(3)maxt,s.t,ϕβ+tn^=Fx,gx≤0,hx=0,xil≤xi≤xiui=1,…,N.The vector constrainϕβ+tn^=Fx ensures that the point x is actually mapped by F to a point on the normal. While the remaining constraints ensure feasibility of x. Instead of ϕβ+tn^=Fx we can use Fx≤u+tv. Where, v is normal vector and u is point of origin of the normal which are user-defined values. The general algorithm of the modified NBI optimization mathematical method is written below. (Algorithm2)Algorithm 2: Modified-NBI algorithm.
Input: Number of significant digitsNumber of initial pointsOutput: non-dominated set//Start//Step 1: Normalize objective FunctionsStep 2: Generating Pareto PointStep 3: Initiate objective function optimizationStep 4: Accurate Pareto setStep 5: Fix objective function ValuesStep 6: Remove Pareto dominated PointsStep 7: Optimal Stage//Finish//In the first stage, the objective functions are placed between the minimum value of zero and the maximum value of one. Therefore, we determine the objective functions in the interval0;1. This causes the Pareto frontier to be placed inside a bound box. In the second step, a V is selected for each pair in this step to generate the space of Pareto points. In the third step, the first optimization problem with the starting minimum point for f1;t=1;β=0 starts. If at the end of optimization t=0;β=1, then go to step 4. Because in this case, a Pareto point has been estimated. Otherwise, use t+V and β+V as starting points for the next optimization, and this step is repeated. In the fourth step, if the estimated Pareto set needs more accuracy, a smaller Vm is used to generate more Pareto points and we return to step 3. Otherwise, we go to step 5. In the fifth step, for N > 2, we use different values in [0, 1]in order to construct the values of objective functions to determine the results of multiple objectives. Then we return to step 1 and determine the objective functions in the interval 0;1. If all combinations are determined, we stop. Otherwise, go to step 6. In the sixth step, based on the described filter [34], the set of inferior Pareto points is removed from the set of generated Pareto points. In the seventh step, choose the best solution from the Pareto set as the optimal solution. For this, choose a range for V for more precision of the generated Pareto set and repeat steps 2 to 6 for all V values in the range. In Figure 2, research solution approach framework are shown.Figure 2
Research solution approach framework.


## 3.3.1. Meta-Modeling Approach

The main elements of the proposed framework are meta model-based optimization [32] in the form of identifying the shape of the meta model, designing experiments to adapt the Meta model, performing simulation experiments, fitting the Meta model, and verifying its accuracy, and optimization considering the meta model in the problem. The algorithm for finding a suitable meta model is shown below. We have used this algorithm for multi-objective optimization based on simulation, the steps of which are fully presented in the next sections. (Algorithm1)Algorithm 1: Meta modelling algorithm.
Input: Objectives, Variables and ParametersOutput: Meta Model and a non-dominated set//Start//Step1: Simulation Model DevelopmentStep 2: Design of Experiment (DOE)Step 3: Run Simulation Model ExperimentStep 4: Fit the Meta ModelStep 5: Meta Model ValidationStep 6: Apply optimization algorithmStep 7: Quality of Results//Finish//In the first step of the algorithm, a discrete event simulation model is developed. If the built model has the necessary validity, we go to the next step; otherwise, the model is accompanied by modifications to obtain the necessary validity. In the second step, a suitable design for an experiment is developed. In the third step, the scenarios designed in step 2 are implemented in the simulation model developed in step 1 to determine the dependent variable of the model. In the fourth step, the best meta model is selected and its unknown coefficients are determined by performing statistical analysis. In this step, we should run the simulation model to determine the response surface for fitting the meta-model. Then, from the data, we obtain an approximate value for the parameter value of the meta-model. Finally, we evaluate these estimates using mathematical and statistical criteria. The fifth step is the answer to the question of whether the meta model built in the fourth step can sufficiently predict the performance of the system or not. In the case of lack of validity, we go to step 2 and change the design of the experiment or go back to step 4 and change the type of meta model. In the sixth step, by applying management constraints, a set of nondominant solutions for the multiobjective optimization problem is obtained using the modified-NBI method. In the seventh step, we will compare the obtained results with the existing situation to identify the improved level.


## 3.3.2. Modified Normal-Boundary Intersection (Modified-NBI) Method

The first goal of the modified-NBI method is to determine the Pareto frontier by solving an optimization problem for multiobjective optimization problems with a continuous or piecewise Pareto frontier. This goal is accomplished in two stages in this method. The first step is to use the modified CHIM in the optimization problem compared to the original NBI algorithm. The second step is to control the iterations of the optimization problem, which are solved by the modified CHIM. According to the modified-NBI algorithm, it is necessary to normalize the objective functions. Therefore, all objective functions have a minimum value of zero and a maximum value of one. If the objective function is unlimited or its maximum value cannot be determined, the user can impose an upper limit value on the algorithm. The first step in the modified-NBI is to use the modified CHIM, which is described below [33]. The second step of the modified-NBI method is to control the iterations while solving the optimization problem to obtain the Pareto frontier. One of the methods used to solve the optimization problem is the quasi-Newton method. Using this method, a relative minimum is obtained for a multi-objective optimization problem. This method has differences compared to meta-heuristic method such as genetic algorithm (GA) that are based on primary populations such as genetic algorithm (GA). For example, a GA requires the selection of parameters such as population size, type of crossover and probability, mutation probability, and number of generations. Also, due to the random nature of population-based methods, even with the same settings, similar answers are often not produced. Hence, these algorithms are often run many times to obtain a reliable set of solutions. The modified-NBI optimization method involves choosing only one parameter to influence the number of generated Pareto points. In addition, similar results are produced each time. Therefore, the algorithm does not need a large number of execution times.Let us show through mathematical formulation how any such boundary point can be found by solving an optimization problem. Givenβ, that ϕβ shows a point in the CHIM. Suppose that n^ represent the normalization unit to the CHIM simplex pointing towards the origin. Then, ϕβ+tn^,t∈R shows the set of points on that normal. Therefore, normal points and boundary of F closet to the origin are the global solution according to the following problem:(3)maxt,s.t,ϕβ+tn^=Fx,gx≤0,hx=0,xil≤xi≤xiui=1,…,N.The vector constrainϕβ+tn^=Fx ensures that the point x is actually mapped by F to a point on the normal. While the remaining constraints ensure feasibility of x. Instead of ϕβ+tn^=Fx we can use Fx≤u+tv. Where, v is normal vector and u is point of origin of the normal which are user-defined values. The general algorithm of the modified NBI optimization mathematical method is written below. (Algorithm2)Algorithm 2: Modified-NBI algorithm.
Input: Number of significant digitsNumber of initial pointsOutput: non-dominated set//Start//Step 1: Normalize objective FunctionsStep 2: Generating Pareto PointStep 3: Initiate objective function optimizationStep 4: Accurate Pareto setStep 5: Fix objective function ValuesStep 6: Remove Pareto dominated PointsStep 7: Optimal Stage//Finish//In the first stage, the objective functions are placed between the minimum value of zero and the maximum value of one. Therefore, we determine the objective functions in the interval0;1. This causes the Pareto frontier to be placed inside a bound box. In the second step, a V is selected for each pair in this step to generate the space of Pareto points. In the third step, the first optimization problem with the starting minimum point for f1;t=1;β=0 starts. If at the end of optimization t=0;β=1, then go to step 4. Because in this case, a Pareto point has been estimated. Otherwise, use t+V and β+V as starting points for the next optimization, and this step is repeated. In the fourth step, if the estimated Pareto set needs more accuracy, a smaller Vm is used to generate more Pareto points and we return to step 3. Otherwise, we go to step 5. In the fifth step, for N > 2, we use different values in [0, 1]in order to construct the values of objective functions to determine the results of multiple objectives. Then we return to step 1 and determine the objective functions in the interval 0;1. If all combinations are determined, we stop. Otherwise, go to step 6. In the sixth step, based on the described filter [34], the set of inferior Pareto points is removed from the set of generated Pareto points. In the seventh step, choose the best solution from the Pareto set as the optimal solution. For this, choose a range for V for more precision of the generated Pareto set and repeat steps 2 to 6 for all V values in the range. In Figure 2, research solution approach framework are shown.Figure 2
Research solution approach framework.


## 4. Computational Results

In this paper, using the proposed framework, the multi-objective optimization of the extraction rate and travel time of moving equipment in a copper mine based on a meta modeling approach has been done. The central composite design is used for sampling and determining the objective values. In this paper, maximizing the amount of total extraction, which is the sum of the extraction of sulfide, oxide, low-grade rocks, and waste in this mine, and minimizing the travel time of hauling according to the limitations of its storage capacity, haulage, and budget are our goals. In addition, we have considered only several important sources such as 120-ton trucks, 240-ton trucks, 35-ton trucks, and 100-ton trucks in the mining complex and used other sources such as shovels is ignored. Because these resources do not directly affect our goals or have little effect on the existing process steps in the copper mining complex. Therefore, by not considering these sources, the size of the design space may decrease. Our approach starts with the design of the operational process of the copper mining complex through the discrete event simulation model. In multi-objective optimization, meta model-based simulation using design of experiment (DOE) is used to analyze scenarios. The design of experiments is used as a valuable set of mathematical techniques for statistical modeling and systematic analysis of a problem with the desired answer to optimize variables [35]. The first step for creating a meta model in the DOE section is the selection of input variables and their considered levels in the system limitation. These variables and their levels are shown in Table 4.Table 4
Predetermined parameters of experimental design.
Variable IDVariable nameMinimumMaximum−α+αX1Trucks 120 tons1228436X2Trucks 240 tons3526X3Trucks 35 tons925133X4Trucks 100 tons25451555These variablesX1;X2;X3;X4 are independent variables that are used as the input value of the simulation model to make the dependent variables of the extraction rate of minerals and the duration of moving trucks. Distances X1;X2;X3;X4 is 17, 5, 21 and 21, respectively, which is considered for each combination. By using a regression meta-model, instead of 17×5×21×21=37485 combinations for only one objective and 2×37485 for both objectives, all combinations of input variables can be shown. A central composite design (CCD) with 25 experiments for both objectives is employed for this purpose. CCD is the most famous design of the response surface method. A CCD consists of a two-stage fractional or full factorial design with central points to which several points called noncenter points have been added. If the distance of the center of the design to the factorial points is considered to be ±1 for each variable, the distance of the center of the design to the noncentered points will be ±α where α>1. The reason for using this design is the proper estimation of curvature in the system model. Then, each combination in this plan is repeated 10 times, and the average of each performance is determined as the dependent variable. Then, the best and most qualitative meta model is selected through statistical analysis.The problem consists of two objectives: maximizing the extraction rate of mineral stones and minimizing the time of moving haulage in the mine. Both goals are calculated using simulation results. Before fitting, we must determine the accuracy of the functions for each objective. UsingR2 and P value statistics for candidate meta models, the best prediction function is selected for each of the objectives. The R2 statistic indicates the difference between the experimental and predicted values. The higher the value, the more significant it means that there is no significant difference between these two values. In Table 5, the validation of the candidate models for each objective has been examined.Table 5
Accuracy model of responses for three candidate regression functions.
ResponseFunctionP valueR2StatusTotal ores productionLinear model0.00040.62SignificantTwo factor interactions model0.02490.69SignificantModified model0.00010.93Significant > selectedTotal haulage transport timeLinear model0.00030.63SignificantTwo factor interactions model0.03350.67SignificantModified model0.00010.94Significant > selectedAccording to the results obtained from evaluating the accuracy of the model, which are shown in Table4, the modified model has sufficient accuracy to predict performance on both response surfaces. Therefore, it is necessary to estimate the coefficients of the significant effects of the model in both objectives to fit the model. Statistical analysis, effect identification, and estimated coefficients for the total production and the trucks travel time in the copper mine complex are presented in Table 6 and 7. The total ore production model F-Value of 3.24 implies significant. Therefore, there is sufficient agreement between experimental and predicted values. Also, the trucks travel time model F value of 19.62 implies the model is significant too. Thus, there is only a 0.01% chance that an F Value this large cloud occurs due to noise. Values of Prob>F less than 0.05 indicate model terms are significant. Also, values of t>1.96 the model terms are significant. Values greater than 0.1 indicate the model terms are not significant. Although X2 in the first response and X4 in second response are not meaningful, but used them for analysis in the model, because these variables are the system decision variables and we intend to calculate their optimal value in the future.Table 6
Estimated effects and coefficients for the total ores production.
TermCoefficientS.E coefficientP valuet valueIntercept0.470.0350.000113.4285X10.150.0330.00054.5454X20.0580.0330.10101.7575X30.200.0330.00016.0606X40.0530.0190.01402.7894X1X20.0690.0230.01053X120.0410.0190.05352.1578X220.0550.0190.01402.8947X1X2X3−0.0640.0230.0160−2.7826X12X20.140.0400.00383.15X12X3−0.180.0400.0007−4.5Table 7
Estimated effects and coefficients for the haulage transport time.
SourceCoefficientS.E coefficientP valuet valueIntercept0.710.0180.000139.4444X1−0.0820.0150.00015.4666X20.230.0260.00018.8461X3−0.0470.0150.0075−3.1333X40.0080.0150.56720.5333X1X40.0390.0180.05032.1666X1X2X30.0480.0180.0205−2.6666X1X2X4−0.0520.0180.0139−2.8888X12X2−0.120.0320.0021−6.6666X1X2X3X40.0730.0180.00154.0555X24−0.0150.0030.0008−5X44−0.0070.0030.0459−2.3333Based on the statistical analysis, the meta-model of the total amount of ore extraction and the total haulage transport time are formulated as follows:(4)Y1=0.47+0.15X1+0.058X2+0.2X3+0.053X4+0.069X1X2+0.041X12+0.055X22−0.064X1X2X3+0.14X12X2−0.18X12X3,(5)Y2=0.71−0.082X1+0.23X2−0.047X3+0.008X4+0.039X1X4−0.048X1X2X3−0.052X1X2X4−0.12X12X2+0.073X1X2X3X4−0.015X24−0.007X44.The copper mine complex can use the above Meta models to find the nondominated solutions subject to given constraint when all functions are validating. Simulation validity measures how well the model represents the real world system [5].


### 4.1. Meta-Model Validation Results

To provide the validity of the meta models built in our paper, we use the most common PRESS error, which is the root mean square PRESS denoted asRMSEPRESS Calculated by PRESS/n , where, n is number of test points selected to evaluate the model. It is obvious that a value of zero for RMSE is the optimal desired values [10]. In Table 8, the RMSEPRESS value obtained for modified models and other considered models for each objective is shown. Therefore, we conclude that the modified models can be used as an abstraction model of the simulation model.Table 8
Validation of the meta models.
Objective functionRMSEPRESSTotal ores extractionY1RMSEPRESSModified=0.16RMSEPRESSLinear=0.2RMSEPRESS2FI=0.26Total trucks travel timeY2RMSEPRESSModified=0.11RMSEPRESSLinear=0.17RMSEPRESS2FI=0.24


### 4.2. Mathematical Optimization

Mathematical problem considered as follows:(6)MaximizeY1;Textractionx,MinimizeY2;TTrucktransportationx,subject to,120x1+240x2+35x3+100x4≤60000,69x1+118x2+22x3+45x4≤3900,12≤x1≤28,1≤x2≤5,5≤x3≤25,25≤x4≤45,xiinteger fori=1;2;3;4.In the integer nonlinear multiobjective optimization (INMOO) problem is considered,Y1 is the function of the total amount of ores extraction in the open-pit mine, and its equation is specified by the symbol Textraction in the model. Y2 is the function of the total transportation time of haulages, which is represented by the symbol TTrucktransportation. x is a vector of design variables that has four components x1 number of trucks is 120 tons; x2 number of trucks is 240 tons; x3 number of trucks 35 tons; x4 The number of trucks is 100 tons. Both functions were obtained using the analysis described earlier. 120x1+240x2+35x3+100x4≤60000 is capacity inequality and 69x1+118x2+22x3+45x4≤3900 is cost inequality. The multiobjective optimization problem has been coded and solved through the modified-NBI method with two meta models with Maple software.Two-dimensional (2D) graphs have been used to show the Pareto frontier of both objectives, where each axis represents each objective. The Pareto frontier represents a surface covering all possible mass values. In engineering applications, including the case study in this paper, the relationship between objective functions and nondominated solutions in the relevant space is an essential issue. Investigating the difference between each nondominated solution and the effect of the difference in the objective functions can help to understand the multi-objective optimization problem. According to the presented algorithm of the modified-NBI method, in Figure3, the values of normalized points are shown. In Figure 3, the diagram of the normalized value points of both objective functions is shown in two dimensions. According to the normalization process, both objective functions are fixed at a value between zero and one.Figure 3
2D diagram of normalized value.Based on the optimization algorithm, 22 nondominated points are produced, and the overall set of Pareto points is created using these points. In Table9, 22 nondominated points and their objective values are shown. Also, in Figure 4 shows the diagram of the space created by these 22 points. According to Figure 4, the convergence of the obtained nondominant points is desirable.Table 9
Nondominated points obtained.
RowSolutionsObjective functionx1x2x3x4Objective 1Y1 (tons)Objective 2Y2 (hourly)1125204521800207821249452170032133155154522850165641252545214001399512325452017016006135254521690116771342545210801277814325452075011579141254519520133810153254521050957.311162254520720842.212171254520400738.813182254521300492.814201254521270278.215211254521270278.216221254521570173.71723125452186089.871824225452215032.491926325452425098.4120164254521960682.3211432545207501157221322545198401482Figure 4
Space of the generated nondominate points.Using the nondominant set created, the optimization phase begins. After all the sets of nondominated solutions are obtained, the best nondominated solutions are reported as the Pareto set of the problem in Table10. The final Pareto set, which includes 13 nondominated solutions out of 22 nondominated solutions, is obtained from the optimization of the multiobjective optimization problem. The nondominated final solutions that make up the Pareto frontier are shown in Figure 5.Table 10
Pareto points obtained.
RowSolutionsObjective functionx1x2x3x4Objective 1Y1 (tons)Objective 2Y2 (hourly)1125204521800207821252545214001399313525452169011674134254521080127751432545207501157614125451952013387153254521050957.38162254520720842.29171254520400738.810182254521300492.811221254521570173.712164254521960682.3131322545198401482Figure 5
Space of the generated pareto points.The accuracy of the solution obtained by the built models compared to the existing situation is shown in Figure5. Based on the obtained results, it is clear that by using the built model and the optimization method, appropriate solutions have been obtained compared to the existing situation. In Figure 5, the red-highlighted dot shows the values of the objective functions in the current state. According to this situation, the value of the objective functions in all Pareto points is better than in the existing situation. Therefore, a designer or engineer chooses the best design of variables that can fulfill the demands of the goals through the Pareto frontier and the set of nondominated solutions.


### 4.3. Optimal Solution

The mathematical programming problem of INMOO is shown in equation (7). This problem can be converted into a mixed integer nonlinear programming (MINLP) according to equation (8). Because, in INMOO mode, instead of calculating an optimal solution, we deal with a set of Pareto solutions. While we want to choose an optimal solution from the set of Pareto solutions. Therefore, it is necessary to convert the INMOO problem into a MINLP problem. For this purpose, first each problem is solved separately, and then the optimal value of each objective function is obtained. Then, by defining the new variable SS instead of the objective functions according to equation (8), we convert the model into a single objective programming problem.(7)MaximizeY1;Textractionx,MinimizeY2;TTrucktransportationx,subject to,hx=0,gx≤0,xil≤xi≤xiu,xiinteger fori=1;2;3;4,where, in above INMOO model Y1 and Y2 are objective functions, hx and gx are equality and inequality, xil and xiu are lower and upper bound of decision variables.(8)MaximizeS,subject to,Y1∗×S≥Y1,Y2∗×S≥Y2,hx=0,gx≤0,xil≤xi≤xiu,Sfree andxiinteger fori=1;2;3;4,where in the above MINLP model S is objective function and is free in sign, Y1∗ and Y2∗ are the optimal value of the objective function, hx and gx are equality and inequality, xil and xiu are lower and upper bound of decision variables. Equation (8) is a nonlinear model. The nonlinear feature of this model is inherent and cannot be fixed. This model is solved using Lingo software.Finally, an optimal solution according to the proposed algorithm is determined among the obtained Pareto solutions. In Table11, the optimal solution is compared with the existing situation.Table 11
Optimum and existing situation.
SituationVariablesObjectivesImprovement (%)x1x2x3x4Objective 1 (tons)Y1Objective 2 (hours)Y2Y1Y2Existing204153612152332679% increase60% decreaseOptimum1252045218002078According to the obtained results, with the change in the combination of the considered resources, there was a 79% increase in the amount of extraction and a 60% decrease in the duration of transport the haulage. Therefore, if the number of trucks with 120 tonsx1 changes from 20 to 12, the number of trucks with 240 tons x2 changes from 4 to 5, the number of trucks with 35 tons x3 changes from 15 to 20, and the number of trucks with 100 tons x4 changes from 36 to 45, the objective functions are improved.


### 4.4. Sensitivity Analysis

In this section, the changes in the important parameters of the problem are examined. For this purpose, two parameters of cost and capacity are considered. Table12 shows the changes of both parameters and their effects on the objective functions.Table 12
Sensitivity analysis.
ParametersShiftObjectivesObjective 1Objective 2Cost+1000229762235−1000207821977Capacity+2000228922175−2000181522789In Figures6 and 7 are shown the cost and capacity changes, respectively. According to the changes in cost and capacity, the optimality of the objective functions is disturbed. For example, if the capacity decreases, the objective function (1) (total extraction) is reduced. In addition, the objective function (2) (transportation time of haulage) is increased. The obtained state is opposite to the optimal situation. Also, if the cost increases, both of function 1and 2 increased, which is the opposite to the optimal situation. Because, function 2 is increased and it is an undesired condition.Figure 6
Cost changing effect on both functions.Figure 7
Capacity changing effects on both functions.


## 4.1. Meta-Model Validation Results

To provide the validity of the meta models built in our paper, we use the most common PRESS error, which is the root mean square PRESS denoted asRMSEPRESS Calculated by PRESS/n , where, n is number of test points selected to evaluate the model. It is obvious that a value of zero for RMSE is the optimal desired values [10]. In Table 8, the RMSEPRESS value obtained for modified models and other considered models for each objective is shown. Therefore, we conclude that the modified models can be used as an abstraction model of the simulation model.Table 8
Validation of the meta models.
Objective functionRMSEPRESSTotal ores extractionY1RMSEPRESSModified=0.16RMSEPRESSLinear=0.2RMSEPRESS2FI=0.26Total trucks travel timeY2RMSEPRESSModified=0.11RMSEPRESSLinear=0.17RMSEPRESS2FI=0.24


## 4.2. Mathematical Optimization

Mathematical problem considered as follows:(6)MaximizeY1;Textractionx,MinimizeY2;TTrucktransportationx,subject to,120x1+240x2+35x3+100x4≤60000,69x1+118x2+22x3+45x4≤3900,12≤x1≤28,1≤x2≤5,5≤x3≤25,25≤x4≤45,xiinteger fori=1;2;3;4.In the integer nonlinear multiobjective optimization (INMOO) problem is considered,Y1 is the function of the total amount of ores extraction in the open-pit mine, and its equation is specified by the symbol Textraction in the model. Y2 is the function of the total transportation time of haulages, which is represented by the symbol TTrucktransportation. x is a vector of design variables that has four components x1 number of trucks is 120 tons; x2 number of trucks is 240 tons; x3 number of trucks 35 tons; x4 The number of trucks is 100 tons. Both functions were obtained using the analysis described earlier. 120x1+240x2+35x3+100x4≤60000 is capacity inequality and 69x1+118x2+22x3+45x4≤3900 is cost inequality. The multiobjective optimization problem has been coded and solved through the modified-NBI method with two meta models with Maple software.Two-dimensional (2D) graphs have been used to show the Pareto frontier of both objectives, where each axis represents each objective. The Pareto frontier represents a surface covering all possible mass values. In engineering applications, including the case study in this paper, the relationship between objective functions and nondominated solutions in the relevant space is an essential issue. Investigating the difference between each nondominated solution and the effect of the difference in the objective functions can help to understand the multi-objective optimization problem. According to the presented algorithm of the modified-NBI method, in Figure3, the values of normalized points are shown. In Figure 3, the diagram of the normalized value points of both objective functions is shown in two dimensions. According to the normalization process, both objective functions are fixed at a value between zero and one.Figure 3
2D diagram of normalized value.Based on the optimization algorithm, 22 nondominated points are produced, and the overall set of Pareto points is created using these points. In Table9, 22 nondominated points and their objective values are shown. Also, in Figure 4 shows the diagram of the space created by these 22 points. According to Figure 4, the convergence of the obtained nondominant points is desirable.Table 9
Nondominated points obtained.
RowSolutionsObjective functionx1x2x3x4Objective 1Y1 (tons)Objective 2Y2 (hourly)1125204521800207821249452170032133155154522850165641252545214001399512325452017016006135254521690116771342545210801277814325452075011579141254519520133810153254521050957.311162254520720842.212171254520400738.813182254521300492.814201254521270278.215211254521270278.216221254521570173.71723125452186089.871824225452215032.491926325452425098.4120164254521960682.3211432545207501157221322545198401482Figure 4
Space of the generated nondominate points.Using the nondominant set created, the optimization phase begins. After all the sets of nondominated solutions are obtained, the best nondominated solutions are reported as the Pareto set of the problem in Table10. The final Pareto set, which includes 13 nondominated solutions out of 22 nondominated solutions, is obtained from the optimization of the multiobjective optimization problem. The nondominated final solutions that make up the Pareto frontier are shown in Figure 5.Table 10
Pareto points obtained.
RowSolutionsObjective functionx1x2x3x4Objective 1Y1 (tons)Objective 2Y2 (hourly)1125204521800207821252545214001399313525452169011674134254521080127751432545207501157614125451952013387153254521050957.38162254520720842.29171254520400738.810182254521300492.811221254521570173.712164254521960682.3131322545198401482Figure 5
Space of the generated pareto points.The accuracy of the solution obtained by the built models compared to the existing situation is shown in Figure5. Based on the obtained results, it is clear that by using the built model and the optimization method, appropriate solutions have been obtained compared to the existing situation. In Figure 5, the red-highlighted dot shows the values of the objective functions in the current state. According to this situation, the value of the objective functions in all Pareto points is better than in the existing situation. Therefore, a designer or engineer chooses the best design of variables that can fulfill the demands of the goals through the Pareto frontier and the set of nondominated solutions.


## 4.3. Optimal Solution

The mathematical programming problem of INMOO is shown in equation (7). This problem can be converted into a mixed integer nonlinear programming (MINLP) according to equation (8). Because, in INMOO mode, instead of calculating an optimal solution, we deal with a set of Pareto solutions. While we want to choose an optimal solution from the set of Pareto solutions. Therefore, it is necessary to convert the INMOO problem into a MINLP problem. For this purpose, first each problem is solved separately, and then the optimal value of each objective function is obtained. Then, by defining the new variable SS instead of the objective functions according to equation (8), we convert the model into a single objective programming problem.(7)MaximizeY1;Textractionx,MinimizeY2;TTrucktransportationx,subject to,hx=0,gx≤0,xil≤xi≤xiu,xiinteger fori=1;2;3;4,where, in above INMOO model Y1 and Y2 are objective functions, hx and gx are equality and inequality, xil and xiu are lower and upper bound of decision variables.(8)MaximizeS,subject to,Y1∗×S≥Y1,Y2∗×S≥Y2,hx=0,gx≤0,xil≤xi≤xiu,Sfree andxiinteger fori=1;2;3;4,where in the above MINLP model S is objective function and is free in sign, Y1∗ and Y2∗ are the optimal value of the objective function, hx and gx are equality and inequality, xil and xiu are lower and upper bound of decision variables. Equation (8) is a nonlinear model. The nonlinear feature of this model is inherent and cannot be fixed. This model is solved using Lingo software.Finally, an optimal solution according to the proposed algorithm is determined among the obtained Pareto solutions. In Table11, the optimal solution is compared with the existing situation.Table 11
Optimum and existing situation.
SituationVariablesObjectivesImprovement (%)x1x2x3x4Objective 1 (tons)Y1Objective 2 (hours)Y2Y1Y2Existing204153612152332679% increase60% decreaseOptimum1252045218002078According to the obtained results, with the change in the combination of the considered resources, there was a 79% increase in the amount of extraction and a 60% decrease in the duration of transport the haulage. Therefore, if the number of trucks with 120 tonsx1 changes from 20 to 12, the number of trucks with 240 tons x2 changes from 4 to 5, the number of trucks with 35 tons x3 changes from 15 to 20, and the number of trucks with 100 tons x4 changes from 36 to 45, the objective functions are improved.


## 4.4. Sensitivity Analysis

In this section, the changes in the important parameters of the problem are examined. For this purpose, two parameters of cost and capacity are considered. Table12 shows the changes of both parameters and their effects on the objective functions.Table 12
Sensitivity analysis.
ParametersShiftObjectivesObjective 1Objective 2Cost+1000229762235−1000207821977Capacity+2000228922175−2000181522789In Figures6 and 7 are shown the cost and capacity changes, respectively. According to the changes in cost and capacity, the optimality of the objective functions is disturbed. For example, if the capacity decreases, the objective function (1) (total extraction) is reduced. In addition, the objective function (2) (transportation time of haulage) is increased. The obtained state is opposite to the optimal situation. Also, if the cost increases, both of function 1and 2 increased, which is the opposite to the optimal situation. Because, function 2 is increased and it is an undesired condition.Figure 6
Cost changing effect on both functions.Figure 7
Capacity changing effects on both functions.


## 5. Managerial Insight

The establishment of the optimal combination in the Sarcheshmeh copper mine complex will improve the amount of total extraction by controlling the number and duration of haulage transportation. In this case, by increasing the amount of extraction, economic growth will be created in a country such as Iran.The presented framework can provide valuable knowledge to mine managers and be used in short-term planning, such as the mining activity in the next shift, and long-term planning, such as the entire life of the mine. The advantage of the proposed planning for the mine is that it facilitates good decisions for mine redesign, mining planning, production rate, and process method. In general, the proposed planning process in the mine is according to Figure8.Figure 8
Proposed mining planning structure.


## 6. Conclusion

In this paper, it is shown that the modified NBI method and regression models are powerful tools for performing multiobjective optimization of physical systems such as mining complexes. In this study, the control model of sulfide, oxide, low-grade and wastes ores extraction and haulage transportation time has been developed for the Sarcheshmeh copper mining complex in Iran. This paper explains how an engineer or designer can easily choose the best variable design that fully satisfies the desired objectives by determining the Pareto frontier and the set of nondominated solutions. By using multi-objective optimization, the response surface method (RSM), design of experiments (DOE), simulation modeling, and optimization based on the modified NBI mathematical method, the effect of input variables and the effect between them have been investigated. The response surface method provides several advantages such as a large amount of information from a small number of experiments that consume time. In addition, the effect of interaction between factors (input variables) on the response is easily revealed using RSM. According to the discrete event model used in this paper, the number of trucks of 120 tons, the number of trucks of 240 tons, the number of trucks of 35 tons, and the number of trucks of 100 tons are the factors that these systems are considered, and their permissible levels are also determined. Using the modified-NBI, the set of nondominating points has been found for the two objectives of maximizing the amount of ore extraction and minimizing the haulage transportation time. Compared to the existing situation, the nondominant points obtained are of very high accuracy in all the combinations obtained. The main results of the paper are as follows:(i)
Determining the best modified regression meta-model to estimate the objective functions and assessing validity to perform the optimization process.(ii)
Formulation of INMOO mathematical problem to discover nondominated solutions.(iii)
Determining 22 nondominant solutions for the problem using the modified-NBI optimization algorithm.(iv)
Discovering 13 Pareto solutions for the INMOO problem.(v)
Determining an optimal solution for the MINLP model using the proposed algorithmIn future studies, multiobjective optimization in the mining complex can be used to rank the final nondominated sets using multicriteria decision-making methods such as TOPSIS. Also, multi-period model can be considered as an attractive suggestion for further studies. [36–37].



---

*Source: 3540736-2022-10-11.xml*