Scholarly article on topic 'Human – Industrial Robot Collaboration: Application of Simulation Software for Workstation Optimisation'

Human – Industrial Robot Collaboration: Application of Simulation Software for Workstation Optimisation Academic research paper on "Economics and business"

CC BY-NC-ND
0
0
Share paper
Academic journal
Procedia CIRP
OECD Field of science
Keywords
{"Human Robot Collaboration" / HRC / "workstation design" / assembly / "man machine" / "workstation optimisation" / "handover position"}

Abstract of research paper on Economics and business, author of scientific article — Fredrik Ore, Bhanoday Reddy Vemula, Lars Hanson, Magnus Wiktorsson

Abstract The simulation possibilities of Human – Industrial Robot Collaboration (HIRC) are limited in commercial software and published research. In order to meet this a demonstrator software has been developed. This paper presents the combination of the quantitative output from the software (measuring operation time and biomechanical load) together with existing optimisation techniques used to design the optimal HIRC workstation. An industrial case is used as an example where the optimal geometric handover position between robot and human is found. From the simulation software metamodels were created in order to represent the investigated workstation. The model was used in a multi-objective optimisation problem and resulted in a trade-off chart between operation time and biomechanical load. The result shows one example of the possibilities to combine the quantitative results from the simulation with optimisation in order to get the best solution to a HIRC design problem.

Academic research paper on topic "Human – Industrial Robot Collaboration: Application of Simulation Software for Workstation Optimisation"

Available online at www.sciencedirect.com

ELSEVIER

6th CIRP Conference on Assembly Technologies and Systems (CATS)

Human - Industrial Robot Collaboration: Application of simulation software for workstation optimisation

Fredrik Orea,b*, Bhanoday Reddy Vemulaa, Lars Hansonb'c,d, Magnus Wiktorssona

a Mälardalen University, School of Innovation, Design and Engineering, Eskilstuna 631 05, Sweden b Scania CV AB, Global Industrial Development, Södertälje 151 87, Sweden c University of Skövde, School of Engineering Science, Skövde 541 28, Sweden d Chalmers University of Technology, Product and Production Development, Göteborg 412 96, Sweden

* Corresponding author. Tel.: +46-8-553-51237; E-mail address: fredrik.ore@scania.com

Abstract

The simulation possibilities of Human - Industrial Robot Collaboration (HIRC) are limited in commercial software and published research. In order to meet this a demonstrator software has been developed. This paper presents the combination of the quantitative output from the software (measuring operation time and biomechanical load) together with existing optimisation techniques used to design the optimal HIRC workstation. An industrial case is used as an example where the optimal geometric handover position between robot and human is found. From the simulation software metamodels were created in order to represent the investigated workstation. The model was used in a multi-objective optimisation problem and resulted in a trade-off chart between operation time and biomechanical load. The result shows one example of the possibilities to combine the quantitative results from the simulation with optimisation in order to get the best solution to a HIRC design problem. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of the organizing committee of the 6th CIRP Conference on Assembly Technologies and Systems (CATS)

Keywords: Human Robot Collaboration; HRC; workstation design; assembly; man machine; workstation optimisation; handover position

ScienceDirect

Procedia CIRP 44 (2016) 181 - 186

www.elsevier.com/looate/procedia

1. Introduction

In order to stay competitive in the globalised world there is a need for manufacturing companies to steadily increase their productivity. Industrial robots have been introduced in industry in the last 50 years and have assisted in huge productivity improvements since then. One way to future improvements is to let humans and robots work closer together. By removing the fences surrounding the industrial robots of today, production systems can be designed that combine the robot's strength and repeatability with the human's tactile sense and flexibility [1]. The best utilisation of these systems can be achieved by introducing robots in today's manual workstations [2]. Compared with manual stations, this should not only increase productivity [3] but also improve human ergonomics [4]. These systems imply that it is possible to work in manufacturing industries for a long time without injuries, thus meeting the

demographic change that is a reality in developed countries [5, 6].

These fences surrounding industry robots may be physical or in the form of certified optical systems. They exists to ensure that no human gets in the way of a moving robot. Major current research efforts in Human - Industrial Robot Collaboration (HIRC) focus on developing systems that prevent the human from getting injured by the robot in a fenceless environment. Less research is found on design and simulation of HIRC systems. In order to meet this gap a demonstrator software has been developed [7, 8]. The software provides the possibility to analyse operation time and human biomechanical load in HIRC assembly workstation design tasks. These evaluations are connected to the proposed benefits of higher productivity and lower ergonomic load. The quantitative results can, together with the exact geometrics of both the product and the workstation, open up possibilities to use mathematical

2212-8271 © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of the organizing committee of the 6th CIRP Conference on Assembly Technologies and Systems (CATS) doi:10.1016/j.procir.2016.02.002

optimisation techniques to design the most favourable workstation [9].

Optimisation techniques have been used in workstation design research for many years. The general objective has for any company always been to optimise profit, and in workstation design that has resulted in operation time optimisation, such as in the design system EMMA [10]. In that system, ergonomic considerations are included, in addition to the time, in the workstation design problem. Manual evaluation of the solutions are however demanded.

Ben-Gal et al. [11] present a method to find the optimal layout considering economic and ergonomic goals by including the statistical technique Response Surface Methodology. In total four objectives are considered and a weight is put on each of them in order to find one optimised result. The same kind of weights on the optimisation objective are set in [12], but also include the digital human modelling tool.

The aim of this research paper is to present a method for applying existing optimisation techniques together with the HIRC demonstrator software to create productive and healthy workstations. These HIRC workstation design problems can be formulated as different optimisation problems. The software analyses both operation time and biomechanical load and the optimal solution must consider both. Thus the task can be describes as an multi-objective optimisation problem. As a proof of concept, this method is applied to optimise a real industrial workstation design case.

2. Method

The optimisation method use the developed HIRC simulation demonstrator software to create virtual simulation of a human and industrial robot collaborative workstation. This quantitative input and output data from the simulation are used to find the input which can lead to optimal workstation design, by using metamodel based multi-objective optimization approach. This process is described in detail by applying it on an industrial case in the following chapters.

2.1 HIRC .simulation .software used

The demonstrator software [7, 8] was used to generate the HIRC simulation used in this study. It combines a digital human modelling tool that evaluates biomechanical load with an industrial robot simulation tool. In this demonstrator software a HIRC work system can be designed, visualised, and evaluated. The evaluation is performed to analyse operation time and biomechanical load in these systems.

Operation time of the robot task is calculated by the simulation software, using robot and environmental data. Human operation time is also calculated by the software using a predetermined time system called SAM [13], based on the MTM (Methods-Time Measurement) system [14]. A predetermined time standard is used to analyse a standard time to perform a manual task. The MTM system was developed in the US in mid 1940's and is one of the most widespread standards [13].

Biomechanical load is analysed with RULA (Rapid Upper Limb Assessment). It is an observational based posture, force,

and muscular assessment tool used to evaluate the risk factors connected with work tasks [15]. It is developed in beginning of the 1990's and consists of an assessment worksheet that is a method used to analyse individual human poses within a work task. It focuses on the positions of individual limbs of the human body and grade them according to the assessment worksheets. The tool quantifies the risk of musculoskeletal injuries on a human posture on a scale from one to seven, where a high score represents a high risk of future injuries.

2.2 Industrial case description

In this study an assembly workstation at a heavy vehicle manufacturer was optimised. It is a flywheel cover assembly station in the engine assembly plant. In this station the flywheel cover, with a weight of up to 60 kilograms and a diameter of 0.6 metres, is assembled on the engine block. The assembly is currently done by a human using an overhead conveyor system and a pneumatic lifting device.

Five operations are performed at the station. Figure 1 shows the workstation layout and the five operations. First, the assembly process starts as the flywheel cover is fetched from a rack of incoming material. Second, the cover is placed in the silicone applying machine, where a small silicone string is applied. Third, the cover is placed at the assembly position on the engine block. In the fourth operation, the flywheel cover is secured on the engine block by 12 screws, using an electric nut runner. The assembly process is then finished as the robot returns to its home position and a new cycle begins with a new product being fetched from the racks.

2.3 HIRC optimisation problems

Multi-objective optimisation aims at seeking a set of optimal trade-off solutions with respect to multiple conflicting objectives [16]. The input data (factors) in the model are assigned to be either parameters (that are chosen to be stable) or variables (that are to be optimised), in order to present the output of the experiment (the response) [17].

The multi-objective optimisation of the HIRC simulation presented is based on two objectives: to minimise the operation time and to minimise the biomechanical load. This HIRC

Fig. 1. Layout of the HIRC assembly station of flywheel cover, with five operations of interest.

system optimisation problem can be defined in numerous ways. Three of them are presented below:

i) Optimise the layout

The positions of all workstation equipment and the industrial robot in three dimensions need to be optimal. The layout should result in short operation times and good biomechanical solutions.

All operations at the flywheel cover assembly station are today performed manually using an overhead conveyor system and a pneumatic lifting device. In [7, 8] an industrial robot (ABB IRB 6620) was placed in the current assembly station so that it could reach all operation positions, visible in figure 1. The manikin used in the simulations was a single male manikin with stature as the key anthropometric variable. A 50th percentile male from a Swedish anthropometric database [18] was used with a weight of 77 kilograms and stature of 179.1 centimeters.

The position of the robot and equipment in the work cell could be optimised further, but the station layout presented in figure 1 will be used in the industrial case in this paper.

ii) Optimise task allocation and sequence between human and industrial robot

To get full use of the HIRC system, the assembly task division between human and industrial robot must be at an optimum.

Table 1 presents the options of performing the assembly operations manually by a human or fully automatically by an industrial robot. The table is used to define the best task allocation and sequence in this example. Since robot movements are faster than human, they should be preferred from an operation time point of view. However, the "move product to assembly" operation does include interaction with two other operators that simultaneously work on the engine. This makes it impractical from a safety point of view to automate the movement to assembly since that would imply that the robot moves in a human work area without human control. The "assemble product with nut runner" demands manual interaction since the robot cannot handle the nut runner and the product simultaneously. This discussion results in a HIRC system where the robot performs operations 1, 2 and 5 and the human in collaboration with the robot operations 3 and 4

This optimal division of task allocation in an operation could also be achieved through optimisation techniques.

iii) Optimise the handover position

The geometrical position where the human takes control of the assembly operation from the industrial robot also needs to be optimal.

This handover position optimisation problem was chosen to be used as a case in this paper. The handover from robotic control to human control is a position between operations 2 and 3. The return of control to the robot is between operations 4 and 5 and has the same geometric position as the first handover. The rest of this paper considers this handover position optimisation problem.

Table 1. Operations to be performed and the resources capable of executing them.

Operations in HIRC system

2. Move

1 g t product to/ 3. Move 4. Assemble

, ^ from product to with nut- 5. Go home

product

silicone assembly runner _machine_

Manual X X X X X

Robotic XX X

2.4 Optimisation process

The optimisation process used in this study is visualised in figure 2. The following sections describe each of the three processes in the boxes in the figure.

2.4.1 Data collection/simulation

The HIRC demonstrator software was used for data collection for the optimisation of the metamodel. The software performed the simulation and transformed factors into responses.

Constraints were put on the available solution space where the optimal handover position was to be found. The handover position was defined as polar coordinates in the x-y plane and a height in the z direction. In order to facilitate other human operations around the engine, the handover position must be in a safe distance from the engine, in this case 1.2 metres. Since the robotic speed is higher than the human, the most beneficial operation time will be achieved when the human motion compared with robotic motion is as short as possible. Hence the length of the polar coordinate is considered to be a parameter of 1.2 metres. The handover position is at this distance in a 90° arc from the engine, leaving the angle as a variable a (0° < a < 90°). The height of the possible handover positions is set to 1.15 +/-0.2 metres, as 1.15 metres is the elbow height of the 50th percentile male manikin that is used in this case. This gives the second variable height from floor, Z (0.95 m < Z < 1.35 m). Figure 3 shows these constraints and variables and presents the geometrical solution space.

An optimal Latin hypercube sample was used to define ten combinations of Z and a that were used in creating the simulated experiments. The four boundary condition factors (combination of max and min Z and a) were also simulated resulting in a total of 14 simulations.

Fig. 3. Graphical representation of parameters and variables that create the solution space (visualised in the red arc).

The simulations performed in the demonstrator software include robotic movements from silicone application machine to one of the 14 simulated handover positions (factors). It also include the manual movement from the engine to the fly wheel cover in order to grasp it and assemble it while the robot carries the load. The final manual return of the robot to the same handover position is also included. This simulation was performed for each of the 14 factors and the resulting responses considering biomechanical load (in terms of a RULA score) and operation time was collected.

2.4.2 Metamodelling

Metamodels (also known as surrogate models) use a limited number of experiments to create a model that represents the system's behaviour. It is described in the form of mathematical expressions that describe the relation between inputs and outputs of the modelled system [19]. These models, if verified in terms of accuracy, can replace the time-consuming simulation process in performing highly iterative evaluation and analysis processes such as an optimisation process [20, 21]. This is needed since simulation of one of the factors in the demonstrator software demands approximately one hour of work. The large majority of this time are depends on manual input and adjustments needed in the CAD data (e.g. to change positions of geometries, and move data between software).

In this study, a set of 14 factors (handover positions) and their corresponding responses from the simulation process were used to develop two metamodels in order to predict time and RULA in the full solutions space. These metamodels are traditionally generated through different forms of polynomial regression. In this paper are the modelling technique kriging used. Kriging adds a stochastic error element as a function of the factors on top of the polynomial function [19]. These models were developed using the Dace Matlab kriging toolbox [22]. The validation of the metamodels was performed through mean Absolute Relative Error (ARE) [23]:

mean ARE = ■

ISO, - MOj I

I So; I

min F(x) = Xftime(x)) + (1 - X)(fRUIA(x)) (2)

x = (Z,a) 0 < a < 90 0.95 < Z < 1.35 0 < X > 1, X = 0:1/N: 1

3. Result

The 14 simulated cases and their resulting time and biomechanical RULA score is presented as black nodes in the metamodels visualised in figure 4.

Two metamodels were created using the Dace Matlab kriging toolbox [22] to predict the values of time and RULA in the whole solution space. The correlation of these values across the grid of factors within the design space is presented in figure

4. In figure 4 are also the optimal values for each of the models presented as highlighted values.

The accuracy of the metamodels was evaluated by comparing the metamodel responses with the HIRC simulation responses of two randomly chosen factors. The results are presented in Table 2. These data in eqn. (1) gives the mean ARE of 1.3% and 3.0% for time and RULA metamodels respectively.

The optimal factors and objectives identified from the metamodels are also presented in Table 2 together with their simulation responses.

Even though the optimal time factors give a close to optimal RULA score, a multi-objective optimisation was performed in order to present the optimisation process. In the weighted sum approach from eqn. (2), 200 different trade-off situations (N= 200) were analysed by varying the weighting factor (X). This resulted in a set of Pareto-optimal solutions presented in figure

5, where each dot represent one solution to the multi-objective optimisation problem in eqn. (2). Each of these dots represent a position on the red arc in figure 3. From figure 5 is it possible to choose the most appropriate trade-off considering the strategies and priorities from the company.

where So and Mo are the simulation and metamodel outputs. 2.4.2 Multi-objective optimisation

Time reduction and biomechanical load reduction have been identified as two objectives that an optimal HIRC layout design has to fulfil. These objectives are conflicting as there don't exist a common minimal handover position for both. Hence this layout design problem is described as a multi-objective optimisation problem [24]. In order to solve this the concept of Pareto optimality [25] is used, where a set of Pareto-optimal solutions containing potential solutions with respect to all possible trade-off situations is generated. This facilitates choosing the most preferable layout design by exploring the trade-off relations between time and biomechanical load in the HIRC system.

The trade-off relations for the optimal layout design are accomplished by using a weighted sum objective function as given in eqn. (2). Multiple optimisation runs of this eqn. are conducted in a complex optimisation algorithm [26] using different weighting factors (A) in order to create Pareto-optimal solutions for (N) possible trade-off situations.

4. Discussion

4.1 The industrial case results

Good accuracy of the metamodels is one prerequisite in applying them in optimisation. The error in the presented HIRC

Table 2. Random values used to calculation accuracy, and values of optimal positions in the developed metamodels.

Factors

Metamodel response

Simulation response

Angle (°) Height (m) Time (s) RULA Time (s) RULA

Random1 76.15 1.206 14.09 3.459 13.8232 3.35

Random2 27.69 1.073 14.22 3.444 14.1362 3.355

Optimal time 34.62 1.145 13.56 3.265 13.5665 3.270

Optimal RULA 23.08 1.165 13.84 3.255 14.5664 3.260

Fig.4: Metamodels for time and RULA. The surfaces represent the functions ftime and fRULA used in the optimisation. The black dots represent the simulated data and the highlighted values are the optimal factors and responses.

workstation design was 1.3 and 3.0%. These numbers shows that the metamodels describe the simulation to a good extent and that the models can be used in the optimisation processes.

The trade-off curve for the presented HIRC workstation design problem are presented in figure 5. The decision can for example be made on constrains in one of the factor, e.g. a maximum allowed RULA score. However, the small difference between minimum and maximum RULA score in this example (3.254- 3.268) have a slight difference in practice. Thus position A marked as a red cross in figure 5 that minimise time is suggested as the optimal position for this problem.

Using optimisation in the example presented results in small savings in operation time or biomechanical load compared with any random selection. From figure 4 can it be seen that the interval for time response was 13.6 to 14.7 s and for RULA 3.2 to 3.9. It can be questioned whether this is a valid way to use the

Fig. 5 : Pareto-optimal solutions presenting the trade-off selection between time and RULA. A presents an suggested optimal trade-off.

workstation design resources in this process. It might be more resource-efficient to focus on the ergonomie height only and pick a handover position somewhere in the middle of the arc. It is however difficult to know the relationships between the factors and the responses before the metamodels are developed, so the simulation and metamodel process in figure 2 have to be performed. The aim with this study was to present the full multi-objective optimisation method and that is why the last multi-objective optimisation process 2 also is included in the text.

4.2 HIRC optimisation method

The study shows that the optimisation method is possible to use when designing HIRC workstations with the demonstration software as a simulation tool. Few workstation design problems are similar to each other. This results in infinite potential problems to address with this method. Thus the metamodel optimisation results presented in this paper can not be generalised and must be developed from the scratch for every problem, thus requiring substantial resources for manual simulation and optimisation tasks.

Three of the potential optimisation problems in HIRC systems are described in Section 2.3. All these can be analysed and optimised with the proposed method. The challenging task is to formulate a question and identify which factors to compare.

A difference between the reviewed multi-objective optimisation work presented in the introduction and the method presented in this paper are when the weights between the objective are chosen. In the earlier work must the weights be selected before the multi-objective optimisation are performed, while this selection is done at the end in the presented method. This enable the decision maker to more easily understand the impact of their choice and base it on their priorities and the resulting trade-off curve.

One way to calculate a general optimum HIRC workstation design is to translate operation time, and investments into one general economic cost parameter. The parameters can include both investment (e.g., industrial robot, control system, materials handling equipment, floor utilisation) and running production (e.g., labour, maintenance) costs. This general optimum production cost solution is not included in this work but is a possible future research area.

4.3 Use of method for wider workstation design tasks

The calculation of optimal HIRC workstations through multi-objective optimisation can be applied in any optimisation problem with conflicting constraints. Multi-objective optimisation problems exist in numerous practical situations in industry, and modern simulation tools offer the possibilities to get a wide variety of quantitative data early in the design process. These can then be used to design the optimum layout given the constraints and model available.

5. Conclusion

This paper presents a method to apply existing optimisation techniques with virtual simulation of human industrial robotic collaboration (HIRC) workstation design. An industry workstation design case is used to demonstrate the optimisation method that includes metamodelling and multi-objective

optimisation. The final result of the design case is a trade-off graph between the objectives of minimising the operation time and biomechanical load that can be used in the decision making process of the workstation. The optimisation possibilities in HIRC systems together with the simulation software will be further developed in order to drastically reduce the simulation time, thus making it possible to include a higher number of factors in the metamodels. This will in turn increase the quality of the models and the optimisation.

Acknowledgements

The research work was funded by the Swedish Knowledge Foundation (for the INNOFACTURE Research School), Scania, and Malardalen University. The research is also supported by the research project HIRC (Virtual Verification of Human Robot Collaboration), funded by the Swedish Governmental Agency for Innovation Systems (VINNOVA). The research was conducted in the context of the XPRES research and education environment at Malardalen University.

Reference

[1 ] Schraft RD, Meyer C, Parlitz C, and Helms E. PowerMate - A Safe and Intuitive Robot Assistant for Handling and Assembly Tasks. in Proceedings of 2005 IEEE International Conference on Robotics and Automation: 4074-4079.

[2] Krüger J, Nickolay B, Heyer P, and Seliger G. Image based 3D Surveillance for flexible Man-Robot-Co operation. CIRP Annals -Manufacturing Technology 2005; 54(1): 19-22.

[3] Krüger J, Lien TK, and Verl A. Cooperation of human and machines in assembly lines. CIRP Annals - Manufacturing Technology 2009; 58(2): 628-646.

[4] Oberer-Treitz S, Dietz T, and Verl A. Safety in industrial applications: From fixed fences to direct interaction. in 2013 44th International Symposium on Robotics (ISR), 1-4.

[5] United Nations, Department of Economic and Social Affairs, Population Division. World Population Prospects: The 2012 Revision. Highlights and Advance Tables. Working Paper No. ESA/P/WP.228, New York, NY, USA: United Nations, 2013.

[6] Zaeh MF and Prasch M. Systematic workplace and assembly redesign for aging workforces. Production Engineering 2007; 1(1): 57-64.

[7] Ore F, Hanson L, Delfs N, and Wiktorsson M. Virtual evaluation and Optimisation of industrial Human-Robot Cooperation: An Automotive Case Study, presented at the 3rd Digital Human Modeling Symposium (DHM2014). Tokyo, Japan, 2014.

[8] Ore F, Hanson L, Delfs N, and Wiktorsson M. Human-Industrial Robot Collaboration - development and application of simulation software. International Journal of Human Factors Modelling and Simulation. In press.

[9] Michalos G, Makris S, and Mourtzis D. An intelligent search algorithm-based method to derive assembly line design alternatives. International Journal of Computer Integrated Manufacturing 2012; 25(3): 211-229.

[10] Braun WJ, Rebollar R, and Schiller EF. Computer aided planning and design of manual assembly systems. International Journal of Production Research 1996; 34(8): 2317-2333.

[11] Ben-Gal I and Bukchin J. The ergonomic design of workstations using virtual manufacturing and response surface methodology. IIE Transactions 2002; 34(4): 375-391.

[12] del Rio Vilas D, Longo F, and Monteil NR. A general framework for the manufacturing workstation design optimization: a combined ergonomic and operational approach. Simulation 2013; 89(3): 306-329.

[13] Laring J, Forsman M, Kadefors R, and Ortengren R. MTM-based ergonomic workload analysis. International Journal of Industrial Ergonomics 2002; 30(3): 135-148.

[14] Maynard HB, Stegemerten GJ, and Schwab JL. Methods Time Measurement. New York: McGraw-Hill; 1948.

[15] McAtamney L and Corlett EN. RULA: a survey method for the investigation of work-related upper limb disorders. Applied Ergonomics 1993; 24(2): 91-99.

[16] Deb K. Multi-Objective Optimization Using Evolutionary Algorithms. Chichester, UK; John Wiley & Sons: 2001.

[17] Kleijnen JP. Regression models and experimental designs: A tutorial for simulation analysts. in Proceedings of the 39th conference on Winter simulation: 40 years! The best is yet to come 2007; 183-194.

[18] Hanson L, Sperling L, Gard G, Ipsen S, and Olivares Vergara C. Swedish anthropometrics for product and workplace design. Applied Ergonomics 2009; 40(4): 797-806.

[19] Simpson TW, Poplinski JD, Koch PN, and Allen JK. Metamodels for Computer-Based Engineering Design: Survey and recommendations. Engineering with Computers 2001; 17(2): 129-150.

[20] Tarkian M, Persson J, Olvander J, and Feng X. Multidisciplinary design optimization of modular industrial robots. in ASME 2011 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference 2011; 867-876.

[21] Jin R, Chen W, and Simpson TW. Comparative Studies of Metamodeling Techniques Under Multiple Modeling Criteria. Structural and Multidisciplinary Optimization 2001; 23(1): 1-12.

[22] Lophaven SN, Bruun Nielsen H, and Sondergaard J. DACE:—a Matlab Kriging Toolbox: Version 2.0, Technical Report IMM-TR-2002-12. Kgs. Lyngby, Denmark; Technical University of Denmark: 2002.

[23] Kleijnen JP and Sargent RG. A methodology for fitting and validating metamodels in simulation. European Journal of Operational Research 2000; 120(1): 14-29.

[24] Hwang C-L, Paidy SR, Yoon K, and Masud ASM. Mathematical programming with multiple objectives: A tutorial. Computers & Operations Research 1980; 7(1): 5-31.

[25] Kasprzak EM and Lewis KE. An Approach to Facilitate Decision Tradeoffs in Pareto Solution Sets. Journal of Engineering Valuation and Cost Analysis 2000; 3(1): 173-187.

[26] Box MJ. A New Method of Constrained Optimization and a Comparison With Other Methods. The Computer Journal 1965; 8(1): 42-52.