ISSN 1854-6250 APEM journal Advances in Production Engineering & Management Volume 14 | Number 4 | December 2019 University of Mari bor ■ Published by CPE apem-journal.org Advances in Production Engineering & Management Identification Statement APEM ISSN 1854-6250 | Abbreviated key title: Adv produc engineer manag | Start year: 2006 ISSN 1855-6531 (on-line) Published quarterly by Chair of Production Engineering (CPE), University of Maribor Smetanova ulica 17, SI - 2000 Maribor, Slovenia, European Union (EU) Phone: 00386 2 2207522,Fax: 00386 2 2207990 Language of text: English APEM homepage: apem-journal.org UniversityofMaribor University homePage: WWW.um.si APEM Editorial Editor-in-Chief Miran Brezocnik editor@apem-journal.org, info@apem-journal.org University of Maribor, Faculty of Mechanical Engineering Smetanova ulica 17, SI - 2000 Maribor, Slovenia, EU Desk Editor Martina Meh desk1@apem-journal.org Janez Gotlih desk2@apem-journal.org Website Technical Editor Lucija Brezocnik desk3@apem-journal.org Editorial Board Members Eberhard Abele, Technical University of Darmstadt, Germany Bojan Acko, University of Maribor, Slovenia Joze Balic, University of Maribor, Slovenia Agostino Bruzzone, University of Genoa, Italy Borut Buchmeister, University of Maribor, Slovenia Ludwig Cardon, Ghent University, Belgium Nirupam Chakraborti, Indian Institute of Technology, Kharagpur, India Edward Chlebus, Wroclaw University of Technology, Poland Franci Cus, University of Maribor, Slovenia Igor Drstvensek, University of Maribor, Slovenia Illes Dudas, University of Miskolc, Hungary Mirko Ficko, University of Maribor, Slovenia Vlatka Hlupic, University of Westminster, UK David Hui, University of New Orleans, USA Pramod K. Jain, Indian Institute of Technology Roorkee, India Isak Karabegovic, University of Bihac, Bosnia and Herzegovina Janez Kopac, University of Ljubljana, Slovenia Qingliang Meng, Jiangsu University of Science and Technology, China Lanndon A. Ocampo, Cebu Technological University, Philippines Iztok Palcic, University of Maribor, Slovenia Krsto Pandza, University of Leeds, UK Andrej Polajnar, University of Maribor, Slovenia Antonio Pouzada, University of Minho, Portugal Rajiv Kumar Sharma, National Institute of Technology, India Katica Simunovic, J. J. Strossmayer University of Osijek, Croatia Daizhong Su, Nottingham Trent University, UK Soemon Takakuwa, Nagoya University, Japan Nikos Tsourveloudis, Technical University of Crete, Greece Tomo Udiljak, University of Zagreb, Croatia Ivica Veza, University of Split, Croatia Limited Permission to Photocopy: Permission is granted to photocopy portions of this publication for personal use and for the use of clients and students as allowed by national copyright laws. This permission does not extend to other types of reproduction nor to copying for incorporation into commercial advertising or any other profit-making purpose. Subscription Rate: 120 EUR for 4 issues (worldwide postage included); 30 EUR for single copies (plus 10 EUR for postage); for details about payment please contact: info@apem-journal.org Cover and interior design: Miran Brezocnik Printed: Tiskarna Kostomaj, Celje, Slovenia Subsidizer: The journal is subsidized by Slovenian Research Agency Statements and opinions expressed in the articles and communications are those of the individual contributors and not necessarily those of the editors or the publisher. No responsibility is accepted for the accuracy of information contained in the text, illustrations or advertisements. Chair of Production Engineering assumes no responsibility or liability for any damage or injury to persons or property arising from the use of any materials, instructions, methods or ideas contained herein. Copyright © 2019 CPE, University of Maribor. All rights reserved. Advances in Production Engineering & Management is indexed and abstracted in the WEB OF SCIENCE (maintained by Clarivate Analytics): Science Citation Index Expanded, Journal Citation Reports - Science Edition, Current Contents - Engineering, Computing and Technology • Scopus (maintained by Elsevier) • Inspec • EBSCO: Academic Search Alumni Edition, Academic Search Complete, Academic Search Elite, Academic Search Premier, Engineering Source, Sales & Marketing Source, TOC Premier • ProQuest: CSA Engineering Research Database -Cambridge Scientific Abstracts, Materials Business File, Materials Research Database, Mechanical & Transportation Engineering Abstracts, ProQuest SciTech Collection • TEMA (DOMA) • The journal is listed in Ulrich's Periodicals Directory and Cabell's Directory journal University of Maribor Chair of Production Engineering (CPE) Advances in Production Engineering & Management Volume 14 | Number 4 | December 2019 | pp 403-510 Contents Scope and topics 406 Femtosecond laser helical drilling of nickel-base single-crystal super-alloy: Effect of 407 machining parameters on geometrical characteristics of micro-holes Yin, C.P.; Wu, Z.P.; Dong, Y.W.; You, Y.C.; Liao, T. Advanced risk assessment in reverse supply chain processes: A case study in 421 Republic of Serbia Aleksic, A.; Runic Ristic, M.; Komatina, N.; Tadic, D. A new method for mathematical and simulation modelling interactivity: A case study in 435 flexible job shop scheduling Ojstersek, R.; Lalic, D.; Buchmeister, B. Surface features of chromium alloyed carbon steel specimens after salt-spray tests in 449 NaCl solution Varga, G.; Torok, T.; Felho, C.; Orosz-Szirmai, G.; Rez, I. A new approach for product quality prediction of complex equipment by grey system theory: 461 A case study of cutting tools for CNC machine tool Pang, J.H.; Zhao, H.; Qin, F.F.; Xue, X.B.; Yuan, K.Y. Effect of the manufacturer quality inspection policy on the supply chain 472 decision-making and profits Hu, H.; Wu, Q.; Zhang, Z.; Han, S. Hybrid fuzzy multi-attribute decision making model for evaluation of advanced digital 483 technologies in manufacturing: Industry 4.0 perspective Medic, N.; Anisic, Z.; Lalic, B.; Marjanovic, U.; Brezocnik, M. High-performance end milling of aluminum alloy: Influence of different serrated cutting 494 edge tool shapes on the cutting force Burek, J.; Plodzien, M.; Zylka, L.; Sulkowicz, P. Calendar of events 507 Notes for contributors 509 Journal homepage: apem-journal.org ISSN 1854-6250 (print) ISSN 1855-6531 (on-line) ©2019 CPE, University of Maribor. All rights reserved. Scope and topics Advances in Production Engineering & Management (APEM journal) is an interdisciplinary refer-eed international academic journal published quarterly by the Chair of Production Engineering at the University of Maribor. The main goal of the APEM journal is to present original, high quality, theoretical and application-oriented research developments in all areas of production engineering and production management to a broad audience of academics and practitioners. In order to bridge the gap between theory and practice, applications based on advanced theory and case studies are particularly welcome. For theoretical papers, their originality and research contributions are the main factors in the evaluation process. General approaches, formalisms, algorithms or techniques should be illustrated with significant applications that demonstrate their applicability to real-world problems. Although the APEM journal main goal is to publish original research papers, review articles and professional papers are occasionally published. Fields of interest include, but are not limited to: Additive Manufacturing Processes Advanced Production Technologies Artificial Intelligence in Production Assembly Systems Automation Big Data in Production Computer-Integrated Manufacturing Cutting and Forming Processes Decision Support Systems Deep Learning in Manufacturing Discrete Systems and Methodology e-Manufacturing Evolutionary Computation in Production Fuzzy Systems Human Factor Engineering, Ergonomics Industrial Engineering Industrial Processes Industrial Robotics Intelligent Manufacturing Systems Joining Processes Knowledge Management Logistics in Production Machine Learning in Production Machine Tools Machining Systems Manufacturing Systems Materials Science, Multidisciplinary Mechanical Engineering Mechatronics Metrology in Production Modelling and Simulation Numerical Techniques Operations Research Operations Planning, Scheduling and Control Optimisation Techniques Project Management Quality Management Risk and Uncertainty Self-Organizing Systems Statistical Methods Supply Chain Management Virtual Reality in Production 406 Advances in Production Engineering & Management Volume 14 | Number 4 | December 2019 | pp 407-420 https://doi.Org/10.14743/apem2019.4.337 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper Femtosecond laser helical drilling of nickel-base single-crystal super-alloy: Effect of machining parameters on geometrical characteristics of micro-holes Yin, C.P.a, Wu, Z.P.a, Dong, Y.W.ab*, You, Y.C.a, Liao, T.a aSchool of Aerospace Engineering, Xiamen University, Xiamen, P.R. China bShenzhen Research Institute, Xiamen University, Shenzhen, P.R. China A B S T R A C T A R T I C L E I N F O Laser micro-hole processing has been widely used in industry. Many laser processing parameters can affect the processing results. The relationship between the geometrical shapes of micro-holes and the laser processing parameters has not been determined accurately. In this paper, experiments on the femtosecond laser drilling of the nickel-base single-crystal super-alloy (DD6) materials were conducted to determine the relationship between the parameters, such as the laser single-pulse energy, rotation rate, and downward focus rate, and the geometrical characteristics of the micro-holes, such as the diameter, and roundness. A group of orthogonal experiments were conducted to determine the effects of the comprehensive influencing factors on the geometrical characteristics of the micro-holes. After the experiments were conducted and analysed, the experimental results were modelled by a backpropagation neural network, and the mapping relationship between the laser parameters and the geometrical morphologies of the micro-holes was constructed. The model established by the backpropagation neural network could obtain accurate prediction results, and the predictions of the diameters of the micro-holes were better than those of the roundness. Keywords: Femtosecond laser; Micro-hole machining; Helical drilling; Nickel-base single-crystal super-alloy (DD6); Orthogonal experiment; Artificial neural networks (ANN) *Corresponding author: yiweidong@xmu.edu.cn (Dong, Y.W.) Article history: Received 26 November 2019 Revised 10 December 2019 Accepted 11 December 2019 © 2019 CPE, University of Maribor. All rights reserved. 1. Introduction Laser micromachining has many applications in material processing, and it has received considerable attention because it can be used in nearly all of the manufacturing sectors [1], such as in the fields of aerospace, the electronics industry [2], biomedical engineering [3], measurement instruments, and the automotive industry [4]. With the rapid development of all these manufacturing sectors, the accuracy requirements of the laser technology have become stricter. The short pulse or even ultarshort pulse laser (pulse duration < 10 ps) has become a reliable technology for manufacturing and production [5] and considered to be a better choice for micro-hole drilling due to its high efficiency, low loss, noncontact procession without any special fixtures for the workpiece [6]. Laser drilling is one of the earliest laser machining technologies used in industry. Since the end of the 1990s, femtosecond lasers have been powerful tools in solid materials machining [7]. Compared with traditional long pulse lasers, femtosecond lasers have advantages in ultra-fine processing, from the submicron scale to the nanoscale, due to their extremely high peak power and multiphoton process with the materials. This advantage lies in limiting the range of energy 407 Yin, Wu, Dong, You, Liao in a small interspace in the center of focus, and not the whole area of irradiation by adjusting the incident energy. Furthermore, another advantage is the versatility in the types of materials that can be processed by femtosecond laser micromachining, such as metals, glasses, polymers, plastics, and semiconductor [8]. The multiphoton absorption and ionization threshold only hinges upon the atomic properties of the materials, not the concentration of free electrons. Theoretically speaking, the range of the heat-affected zone in femtosecond laser micromachining processing is negligible [9], and thus, non-hot melt processing with high precision is available. Due to its ultrashort pulse duration acting on the materials in short periods of time with high power densities, the energy absorption can be controlled only in the machining region and cause materials to enter the plasma state while materials outside the machining region remains cold without heat diffusion. The laser drilling process is a preferred method due to its unique properties of no tool wear, flexibility, and high rates [10]. Laser processing parameters have important influences on quality, efficiency, and cost. The main research parameters include the wavelength, pules width, spot radius, pules energy, repeat frequency, and rotation rate [11, 12], which directly affects the dimensional and shape accuracy of the processed micro-holes. Furthermore, there are four common strategies in laser drilling: single pulse percussion drilling, multiple pulse percussion drilling, rotary trepanning and helical drilling [13] as shown in Fig. 1 [14]. In percussion drilling, laser pulses shaped as voxels of a particular size are focused on a desired location for the desired time to ablate materials with the same geometric dimensions as those of the voxels. For trepanning and helical drilling, the laser pulse moves in a circular route to the required diameter of the presupposed micro-hole [15]. Although the application of femtosecond laser in drilling has attracted considerable attention in academia and industry, there are few accurate quantitative models that can describe the relationship between the laser parameters and geometrical shape of the micro-holes. In practice, the selection of parameters always depends on experiences of previous experiments, and the geometrical parameters of the micro-holes are difficult to forecast before machining. Therefore, it is vital to establish certain mapping relations between the processing parameters and geometric features. Single pulse drilling Multi-pulsed drilling Trepanning drilling Helical drilling Fig. 1 Four micro-hole drilling strategies: From the left to right are single pulse percussion drilling, multiple pulse percussion drilling, rotary trepanning, and helical drilling 2. Materials and methods 2.1 Experimental equipment and methods In the first experiment, lasers with different single-pulse energies, rotation rates, and downward focus rates were used to perform helical drilling. In addition, a group of L25(54) orthogonal experiments are processed to determine the comprehensive influences of the single-pulse energy, rotation rate, and downward focus rate on the geometrical parameters of the micro-holes. The machined material was a nickel-base single crystal super-alloy (DD6) with a thickness of 1 mm. The basic parameters of the laser system are listed in Table 1. A KH7040A-1 5-axis numerical control machine with PHAROS high power and energy femtosecond lasers was used, as shown in Fig. 2. The schematic diagram of the femtosecond laser machining system is shown in Fig. 3. Although the polarization state of the light mainly affects the roundness and the surface accuracy of 408 Advances in Production Engineering & Management 14(4) 2019 Femtosecond laser helical drilling of nickel-base single-crystal super-alloy: Effect of machining parameters on geometrical ... the micro-holes, there was no significant difference between polarized light in the quality of the micro-holes from the helical drilling [16]. Therefore, the influence of the laser polarization state factor was not considered. In the experiment, the quarter-wave plate (B) was used to convert the laser light into circularly polarized light and then focused on the surface of the processed materials. A four-optical-wedge beam rotation apparatus (E) is used to deflect the laser beam to achieve a planar helical motion. The rotary cutting module was composed of the focusing lens (F), and the laser outlet nozzle (G) could move along the material surface perpendicularly. By controlling the vertical adjustment device (H), the downward focus rate could be controlled, and combined with the planar helical movement, helical drilling could be achieved. The projection of the machining path is shown in Fig. 4, with a 0.5 mm radius of the outermost ring. The features of the micro-holes were observed and measured using SUPRA55 field emission gun scanning electron microscope (SEM). After machining, the specimens were soaked in ethanol solution and cleaned ultrasonically to remove the liquid melt on the surfaces of the specimens before the SEM measurements. Table 1 Parameters of laser machining system Parameter Unit Range Average power W 0-15 Pulse width fs 250 Repetition rate kHz 60-600 Wavelength nm 1064 Focal distance mm 150 Air blowing pressure MPa 0-0.5 Rotation rate r/min 600-2400 Spot radius |m 15-25 Fig. 2 Photo of laser device and specimen t A Î B t C A. PHAROS femtosecond laser B. Quarter-wave plate C. Beam expander D. Mirror H E. Beam rotation apparatus F. Focusing lens G. Outlet nozzle H. Vertical adjustment device I. Workpiece Fig. 3 The schematic diagram of the femtosecond laser machining system Advances in Production Engineering & Management 14(4) 2019 409 Yin, Wu, Dong, You, Liao 0.5 0.4 0.3 0.2 £ o.i 0 -0.1 -0.2 -0.3 -0.4 -0.5 Fig. 4 The projection of the laser helical scanning path 2.2 Single factor experiments, helical micro-hole machining, trepanning machining The first experiment was to determine the influence of the single-pulse energy on the features of micro-holes. The parameters that describe the features include the diameters, the roundness of entrance and exit, and the taper of each hole. The roundness of a micro-hole Ar can be calculated by the difference of diameters of the outer contour Ct an d inner contour C2 [17]. Ar = C1-C2 (1) The outer contour is the minimum circumscribed circle of the hole, while the inner contour means the maximum inscribed circle of the hole. A smaller Ar means a more circular and less rough hole. The taper can be calculated by Eq. 2. Taper = 2 x tan'1 (2) where R1 and R2 denote the entrance and exit diameters, respectively. The focus of the laser was set on the platform of the specimen with a repetition rate of 100 kHz, rotation rate of 2400 r/min, and maximum processing time of 200 s. The single-pulse energy was set from 80 to 130 |iJ at intervals of 10 |iJ. The second experiment examined the relationship between the rotation rate and the features of the micro-holes. Similar to the first experiment, the focus was set on the platform of the specimen with a repetition rate of 100 kHz, single pulse energy of 100 J and maximum processing time of 200 s. The rotation rate was set from 600 to 2400 r/min at intervals of 300 r/min. The third experiment aims to explore the effect of the focus downward rate on helical laser micro-holes drilling. The downward focus rate was set from 0 to 0.03 mm/s at intervals of 0.005 mm/s. The repetition rate was set to 100 kHz, the single-pulse energy was set to 100 J and the rotation rate was set to 2400 r/min. The processing time was 100 s. For all of the experiments above, blowing pressure was set to 0.5 MPa. For each set of data, three drilling processes were conducted and the average was obtained. 2.3 Orthogonal experimental design (OED) To investigate the effect and significance of the single-pulse energy, the rotation rate and the downward focus rate on the geometric characteristics, such as the diameters and roundness of the entrance and exit of the micro-hole, the backpropagation (BP) neural network (BPNN) was used to obtain the mapping relationship between the geometric characteristics of the micro-holes and these parameters. L2s(54) orthogonal experiments with three processing parameters were designed, as it can get influence relationship between experimental variables through a small amount of experiments. In the OED, the occurrence frequencies of the levels in each column were the same, and the number of occurrences of an ordinal number pair consisting of all of the data in any two columns was the same. This ensured a uniform dispersion of the experimental conditions and eliminated the interference from other factors [18]. As shown in Table 2, a blank column was designated for the error evaluation, and five levels were set for each factor. The laser repetition rate was set to 100 kHz and the blowing pressure was 0.5 MPa. The processing time was 200 s for each micro-hole. -0,5-0.3 -0.1 0.1 0.3 0-5 A'/limi 410 Advances in Production Engineering & Management 14(4) 2019 Femtosecond laser helical drilling of nickel-base single-crystal super-alloy: Effect of machining parameters on geometrical ... Table 2 Factors and levels of orthogonal experiments Factor Level Single pulse energy, B Rotation rate, (r-min-1) Focus downward rate, _(mm-s-1)_ Blank 1 2 3 4 5 60 70 80 90 100 1600 1800 2000 2200 2400 0.005 0.01 0.015 0.02 0.025 3. Results and discussion 3.1 Effect of single pulse energy Fig. 5 shows the variation of the diameters and roundness of the micro-holes with the single-pulse energy. This indicates that, for the exit diameter, there was a leap between 100 |j.J to 110 |j.J in the single-pulse energy from 543.29 |j.m to 951.97 |j.m. When the single-pulse energy is no more than 100 J the shapes of the exits of the micro-holes deviate significantly from a circle, which means that for a better drilling result, a threshold of single-pulse energy is needed when the other related processing parameters are fixed. Fig. 6 shows the SEM photos of the hole exit with single pulse energy of 80 |iJ and 110 J the disparity of roundness is obvious. In fact, a laser with a higher single-pulse energy can produce a high metal vapor pressure that can carry away more material in a molten or liquid state [19]. Furthermore, the laser energy followed a Gaussian distribution, which caused the diameter of the focus on the processed material to increase with increasing pulse energy. During ablation, heat transfer occurred. Although the heat transfer was weak at any instant, its cumulative effect could not be neglected. The scanning path of helical drilling causes the laser to not only process at a certain point, which greatly reduces the heat accumulation and increases the heat release [20]. As the single-pulse energy increased, the entrance and exit diameters increased to maximum values and stayed in the range of about 1050 |j.m and 958 |j.m, and the taper decreased from 0.52 and remained at approximately 0.1. The roundness of the entrance and exit exhibited similar patterns, with diameters of approximately 22 |j.m and 68 |j.m. Single pulse energy C y J) Single pulse energy (y3 Fig. 5 Geometric characteristics change to single pulse energy Fig. 6 Photos of the exit hole with single pulse energy of 80 |J and 110 |J Advances in Production Engineering & Management 14(4) 2019 411 Yin, Wu, Dong, You, Liao 3.2 Effect of rotation rate Fig 7 shows the variation of the diameters, roundness and tapers of the micro-holes with the rotation rate. This indicates that the entrance diameter increased from the top at 1066.68 |j.m, then decreased and remained steady around 1045 |j.m as the rotation rate goes up. The exit diameter remained in the range between 937 |j.m to 952 |j.m. The entrance roundness decreased overall, but the exit roundness increased generally. The rotation rate significantly influenced the exit roundness, because as the rotation rate increased, the number of pulses during a single scanning cycle was reduced. Therefore, the effect of the rotation rate on the processing was essentially due to changes in the number of pulses. The reduction in the number of pulses during one processing cycle resulted in layer-by-layer processing to the bottom, and the corresponding reduction in the amount of erosion caused the exit roundness to deteriorate. The taper remained at approximately 0.1. Compared with other parameters, there was less influence of the rotation rate on the entrance and exit diameters of the micro-holes because the variation of the rotation rate did not affect the number and power of pulse acting on the material. Furthermore, the thermal incubation effect was not evident because the laser spot did not focus on only one area during helical drilling [21]. —■— lin trance roundness —•— lix it roundness 900 1200 1600 1800 2100 Notation speed (r/inir) 000 900 1200 1600 1800 2100 Rotation jpM£ (r/min) Fig. 7 Geometric characteristics change to rotation rate 3.3 Effect of focus downward rate Fig. 8 shows the variation of the parameters of the micro-holes with the downward focus rate. When the downward focus rate increased, the entrance and exit diameters exhibited the opposite trend, and thus, a sharp drop in the taper occurred. The entrance roundness remained between 20 |j.m and 30 |j.m, whereas the exit roundness fluctuated significantly. The downward focus rate was strongly connected to the defocusing amount, which is an important factor influencing the features of the micro-holes [22]. Most of the energy of the laser beam distributed over a certain range of focal depths. The focal plane was where the diameter of the beam waist was at a minimum. Defocusing means the focal plane is not coincident, the deviation distance between which is called the defocusing amount [23]. Fig. 9 shows from left to right illustrations of negative, zero, and positive defocusing. As the focal point moved downward, the negative defocusing amount causes a transition of the phase to occur only by heat transfer on the walls of the micro-holes. A lower ratio of gaseous to molten or liquid states causes the residues to remain bounded at the edges of holes, which enlarged heat-affected zone and decreased the area of the effective heat transfer region. As the focus moved downward, only part of the material on the surface could be melted, whereas the energy of the laser beam concentrated on the focal position moved closer to the back face of the material. Thus, the entrance diameter decreased and the exit diameter increased, causing the taper decreased as the downward focus rate increased. Fig. 10 is the photos of the entrance hole with focus downward rate of 0 and 0.03 mm/s. The difference of the taper is obvious. Influenced by residues, the non-uniform heating contributed to the great fluctuation of the exit roundness. In addition, the appearance of the cone and the change of the taper 412 Advances in Production Engineering & Management 14(4) 2019 Femtosecond laser helical drilling of nickel-base single-crystal super-alloy: Effect of machining parameters on geometrical ... were related to the incident angle of the laser. According to the Fresnel effect [24], when the incident angle of a laser increases, the absorption rate of materials will decrease. During processing, changing the distribution of the laser radiation intensity led to the formation of a cone. When the laser irradiated the micro-holes, the walls could not absorb all the energy, and ablation rate decreased gradually [20]. When the energy density absorbed by the material was lower than the ablation threshold, the ablation process ended. 600 400 60 1 Entrance roundness 'Exil roundness 0.000 0.005 0.010 0.015 0.020 0.025 Speed oï focus down (rtim/s) 0.005 0.010 0.015 0.020 0.02b 6.030 Speed oí focus down Cimn/s) Fig. 8 Geometric characteristics change to focus downward rate Fig. 9 Schematic diagram of the defocusing amount. 'a' is the laser beam; 'b' is the lens and 'c' is the processed material. Fig. 10 Photos of entrance hole with focus downward rate of 0 and 0.03 mm/s Advances in Production Engineering & Management 14(4) 2019 413 Yin, Wu, Dong, You, Liao 3.4 Orthogonal experimental design result and analysis The orthogonal experiment results are shown in Table 3. For training and testing the BPNN, another group of experiments are designed and processed (Table 4). Range and variance analyses [25] are common methods for investigating the influence of each factor on the experiment results. In range analysis, the order of influence of the factors on the results can be seen by comparing the ranges of each factor. There are two parameters in the range analysis, T and R. T is defined as the sum of all of the results in the experiments. denotes the sum of the results of level i in a column. R is defined as the range between the maximum and minimum of Tt for each factor. R = max(Ti) — min(J{) (3) Table 3 Orthogonal experiment result Test A B number Single pulse Rotation rate, energy, |J (r-min-1) C Focus downward rate, [ram's-1] Entrance diameter, |m Exit diameter, |m Entrance roundness, |m Exit roundness, |m 1 60 1600 0.005 1 1032.06 512.46 48.59 56.50 2 60 1800 0.015 3 922.87 661.28 47.30 41.22 3 60 2000 0.025 5 909.87 511.11 36.94 48.41 4 60 2200 0.01 2 1004.80 525.93 33.07 31.10 5 60 2400 0.02 4 911.47 577.78 36.31 52.07 6 70 1600 0.025 2 919.27 604.71 37.91 33.52 7 70 1800 0.01 4 995.20 592.59 28.64 51.60 8 70 2000 0.02 1 937.60 644.44 28.09 41.97 9 70 2200 0.005 3 1036.80 550.61 36.98 61.72 10 70 2400 0.015 5 990.47 723.46 39.24 44.74 11 80 1600 0.02 3 947.20 701.23 32.99 48.25 12 80 1800 0.005 5 1034.00 952.78 36.00 43.44 13 80 2000 0.015 2 988.90 806.98 32.22 50.44 14 80 2200 0.025 4 935.80 617.28 42.48 42.86 15 80 2400 0.01 1 1008.00 850.74 42.39 23.14 16 90 1600 0.015 4 976.00 814.91 48.30 50.25 17 90 1800 0.025 1 921.60 650.94 39.90 25.99 18 90 2000 0.01 3 1017.60 872.73 34.40 29.09 19 90 2200 0.02 5 957.60 708.07 38.78 29.61 20 90 2400 0.005 2 1038.18 960.74 37.26 28.06 21 100 1600 0.01 5 1022.72 915.56 23.60 22.40 22 100 1800 0.02 2 944.40 763.40 31.82 30.98 23 100 2000 0.005 4 1040.72 960.27 32.56 25.74 24 100 2200 0.015 1 989.60 836.60 31.04 42.75 25 100 2400 0.025 3 940.80 718.01 36.08 24.94 Table 4 Supplemental experiments for BPNN modelling Test A B C Entrance Exit Entrance Exit num- Single pulse Rotation rate, Focus downward diameter, diameter, roundness, roundness, ber energy, |J (r-min-1) rate, (mm^s-1] |m |m |m |m B1 60 1600 0.015 915.2 658.59 37.18 52.87 B2 60 1800 0.025 905.12 505.35 41.09 68.56 B3 70 2200 0.02 921.6 654.86 40.71 51.82 B4 80 2400 0.005 1039.98 949.25 38.45 56.19 B5 80 1600 0.015 983.2 777.78 30.05 49.13 B6 90 2000 0.025 918.9 684.44 39.64 27.54 B7 90 1800 0.01 1008.4 884.44 34.48 29.28 B8 100 2200 0.02 942.4 874.93 32.17 38.37 B9 100 2400 0.005 1036.8 969.33 38.88 28.34 414 Advances in Production Engineering & Management 14(4) 2019 Femtosecond laser helical drilling of nickel-base single-crystal super-alloy: Effect of machining parameters on geometrical ... The greater the value of R of the factor is, the greater the influence of the factor is on the result. The range analysis can be used to determine the impact order of the factors intuitively, but the experimental error cannot be estimated. Moreover, the impact magnitudes of the factors cannot be indicated, especially for the OED with more than three levels. Considering the limitation of range analysis, variance analysis is necessary to determine the impact degrees or magnitudes of the factors on the results. The most important parameters in the variance analysis are the sum of squared deviations for the total, for each column, and error. The total sum of squared deviation (SST) and the total freedom (fT) are defined as follows: = ^(x¿-x)2 (4) SST ¿=i fT=n-l (5) The sum of squared deviations for each factor (SSy) and the freedom for each (fj) are defined as follows: SSj m 1 1V ? 1 = ~)Ti2----(6) r ¿—i n ¿=1 fj = m-lQ = l,2.....fc-1) (7) The sum of squared deviations for the experimental error (SSE) and its freedom (fe) are defined as follows: k SSE = SST ~XSSJ (8) 7 = 1 k fe=fT~YJfj (9) 7 = 1 The variances of each factor (MSj) and the one of experimental error (MSE) can be calculated by as follows: SS MSi=J- (10) MSE (11) Je where the F value is used to illustrate the magnitudes of the effects on the factors. This reflects the ratio of the sum of squared deviations of each factor to that of the experimental error. Thus, the F value of each factor is as follows: (12) ] MSE In Eqs. 4 to 12, n = 25, m = 5, k = 4, and r = n/m = 5 for L2s(54) OED. In the F-distribution table, comparing Fj with the value of Fa(fT, fe) can determine whether the factor effect on the results is prominent or not, which is called the F-test. If Fj is larger than Fa, the factor effect on the results is prominent and vice versa. Usually, a = 0.05 is the demarcation of prominence, and a = 0.01 is the demarcation of high prominence [26]. In addition, if the variance of a factor is less than double the variance of the experimental error (MSj < 2 MSE), this factor can be seen as an error and its squared deviation and freedom will be added to those of the error. In the result tables, the factors meeting this criterion are marked with an asterisk '*'. The results of the range and variance analysis are shown in Table 5-12. All of the three factors have significant influences on the hole diameters, while the downward focus rate and the singlepulse energy are the most important factors affecting the entrance and exit diameter results, respectively. For the roundness, the rotation rate was the weakest factor, and the single-pulse energy was the main influencing factor. In the range of experimental data, there was a positive Advances in Production Engineering & Management 14(4) 2019 415 Yin, Wu, Dong, You, Liao correlation between the diameters and the single-pulse energy and a negative correlation between diameters and downward focus rate. The effect of rotation rate on diameters is fluctuant. On the whole, the increase in the single-pulse energy was helpful for reducing roundness, especially for the hole exit The entrance roundness was lower than the exit roundness in general. With a greater pulse energy, more photon energy can be absorbed by the materials. The stronger capability of laser to remove materials makes micro-holes more circular [27]. Fig. 11 represents the exit of the hole processed with single pulse energy at 80 J Rotation rate at 2400 r/min and focus downward rate at 0.01 mm/s. In comparison with the exit hole shown in Fig. 6, the focus moving down is beneficial to the exit hole for its roundness. Defocusing can improve heat transfer through the interior from entrance to exit and material at the back-side can receive the energy from laser more evenly and sufficiently. Table 5 Range analysis for entrance diameters A B C Blank T1 4781.07 4897.25 5181.76 4888.86 T2 4879.34 4818.07 5048.32 4895.55 T3 4913.90 4894.69 4867.84 4865.27 T4 4910.98 4924.60 4698.27 4859.19 T5 4938.24 4888.92 4627.34 4914.66 R 157.17 106.53 554.42 55.37 Order CAB Table 6 Variation analysis for entrance diameters Source of Sum of squares Freedom Variance F-value Critical value Significance variation of F-test (Fa) A 3035.54 4 758.89 22.09 F0.05 = 3.2592 High B 1261.34 4 315.33 9.18 F0.01 = 5.412 High C 43258.14 4 10814.54 314.85 High Blank* 412.18 4 103.05 3.00 Other errors 412.18 12 34.35 Total 47967.21 24 Table 7 Range analysis for exit diameters A B C Blank T1 2788.56 3548.87 3936.86 3495.18 T2 3115.81 3620.99 3757.55 3551.76 T3 3929.01 3795.53 3843.23 3503.86 T4 4007.39 3238.49 3394.92 3562.83 T5 4193.84 3830.73 3102.05 3810.98 R 1405.28 592.24 834.81 315.8 Order ACB Table 8 Variance analysis for exit diameters Source of Sum of squares Freedom Variance F-value Critical value Significance variation of F-test (Fa) A 303899.08 4 75974.78 65.40 F0.05 = 3.2592 High B 44994.62 4 11248.65 9.69 F0.01 = 5.412 High C 97445.99 4 24361.50 20.97 High Blank* 13939.81 4 3484.95 3.00 Other errors 13939.81 12 1161.65 Total 460279.49 24 Table 9 Range analysis for entrance roundness A B C Blank T1 202.21 191.39 191.39 190.01 T2 170.86 183.66 162.10 172.28 T3 186.08 164.21 198.10 187.75 T4 198.64 182.35 167.99 188.29 T5 155.10 191.28 193.31 174.56 R 47.11 27.18 36.00 17.73 Order ACB 416 Advances in Production Engineering & Management 14(4) 2019 Femtosecond laser helical drilling of nickel-base single-crystal super-alloy: Effect of machining parameters on geometrical ... Table 10 Variation analysis for entrance roundness Source of Sum of squares Freedom Variance F-value Critical value Significance variation of F-test (FJ A 309.604 4 77.40 7.97 F0.05 = 3.0069 High B* 98.40 4 24.60 2.53 F0.01 = 4.7726 C 213.18 4 53.30 5.49 High Blank* 56.99 4 14.28 1.49 Other errors 155.39 16 9.71 Total 678.17 24 Table 11 Range analysis for exit roundness A B C Blank T1 229.30 210.92 215.46 190.35 T2 233.55 193.23 157.33 174.10 T3 208.13 195.65 229.40 205.22 T4 163.00 208.04 202.88 222.52 T5 146.81 172.95 175.72 188.60 R 86.74 37.97 48.42 Order ACB Table 12 Variation analysis for exit roundness Source of Sum of squares Freedom Variance F-value Critical value Significance variation of F-test (FJ A 1234.91 4 308.73 10.92 F0.05 = 3.0069 High B* 181.31 4 45.33 1.60 F0.01 = 4.7726 C 689.62 4 172.41 6.10 High Blank* 270.90 4 67.72 2.40 Other errors 452.21 16 28.26 Total 2376.74 24 Fig. 11 Photo of the exit of the hole processed with single pulse energy of 80 |J, rotation rate of 2400 r/min, and focus downward rate of 0.01 mm/s 3.5 Prediction of micro-hole geometry characteristics with back propagation artificial neural network (ANN), modelling with artificial neural network, testing of the model Containing functions of self-learning and self-training, artificial neural networks (ANNs) are suitable for finding relationships between data inputs and outputs, especially for nonlinear data. In general, a neural network can be seen as a mapping relationship between input and output variables [28]. Similar to a biological brain, a neuron is the basic unit of an ANN. One certain neuron in each layer in the network is connected to each neuron in the previous and the next layers. An ANN is a black box that can avoid difficulties in data analysis and modelling. It is suitable for data with uncertainties and a lack of structure. The BPNN is a multilayer feedforward neural network that uses an error back-propagation algorithm [29]. It is widely used in experimental data modelling, even if the mapping relationship seems complex and difficult to formulate [30]. Advances in Production Engineering & Management 14(4) 2019 417 Yin, Wu, Dong, You, Liao Before establishing the BPNN model, experimental data should be normalized to improve the efficiency and sensitivity while training and to prevent overfitting. To ensure that the BPNN model has the ideal extrapolation capability, normalized pre-treatment values should be in the range of 0.2 to 0.8. The data pre-treatment can by calculated as follows: x — min(x) y = (0.8 — 0.2) x-—-^- + 0.2 (13) max(x) — min(x) The reverse data pre-treatment can by calculated as follows: max(x) — min(x) x = (y-0.2)x-08 — 02 +minW (14) In Eqs. 13 and 14, x represents a certain value in the original data of each factor, and y is the corresponding value in the pre-treatment data of each factor, max(x) and min(x) are the maximum and minimum values in original data of each factor, respectively. Of all of the training algorithms used in BPNN, the Levenberg-Marquardt algorithm was an optional one, as it can adjust the network training parameters automatically to ensure the network has real-time updating for its applicable training methods [31]. The common transfer functions used in BPNNs include the hyperbolic tangent sigmoid function (tansig), logarithmic sigmoid function (logsig), and linear function (y = x). Based on previous research [32], using the tansig function as the transfer function in the input-hidden layer and the linear function in hidden-output layer can obtain the best optimization in most BPNN modelling. For a BPNN with there or more layers, the network model can approximate a nonlinear function infinitely as long as there are sufficient neurons [33]. Therefore, the BPNN configuration used in this study contained one input layer, one hidden layer and one output layer. The neural network toolbox in MATLAB was used to establish the BPNN model. It is important to determine the proper number of neurons in the hidden layers. Including too many neurons may lead to the problems like over-fitting and local minima. To ensure the network performance and generalization ability are sufficiently high, the structure of the network should be as compact as possible to achieve high precision. With the experimental formulas YH= o CA-l > k (where k is for the number of data for the input layer, n1 is the number of neurons in hidden layer, and n is the number of neurons in the input layer), n1 = ^n + m + a (where m is the number of neurons in the output layer, a e [1,10]) [34], and n1 = log2 n [35], the number of neurons was selected as 26. As a result, the BPNN with the topological structure of 3-26-4 was chosen. After training the model 49 times, the network converged. Table 13 shows the comparison of the predicted results in this model and the practical experimental data. The total average of the error was 4.06 %. The errors in the diameter were smaller than those of roundness. The factors and scope of the parameters that affected the roundness were more complex than the diameters, which may have caused the roundness values to be more difficult to predict and have a larger deviation than the diameters. The highest error of the result was less than 7 %. Table 13 Comparison of the predicted results in this model and the practical experimental data Parameters Entrance diameter, Exit diameter, Entrance roundness, Exit roundness, |m |m |m |m Test B5 Predicted result 959.89 755.54 33.03 50.96 Actual result 983.2 777.78 30.05 49.13 Error 2.42 % 2.86 % 9.02 % 3.59 % Test B9 Predicted result 1051.33 993.67 37.27 30.21 Actual result 1036.80 969.33 38.88 28.34 Error 1.40 % 2.51% 4.14% 6.60 % 418 Advances in Production Engineering & Management 14(4) 2019 Femtosecond laser helical drilling of nickel-base single-crystal super-alloy: Effect of machining parameters on geometrical ... 4. Conclusion By laser helical drilling experiments, rules for the machining parameters and geometric morphological characteristics of micro-holes were determined in this study. In the range of 80 |iJ to 130 J with the single pulse energy increasing, the diameters increased. The Diameter and roundness of the entrance and exit tended to be stable, while the taper exhibited a slight change. The increase in the rotation rate of the laser can lead to a fluctuation in a range in the values of diameters, but the taper remained almost invariant. The downward focus rate had a remarkable influence on the characteristics of the micro-holes. With the defocusing rate increasing, the diameters changed greatly. The OED shows that when the single-pulse energy was in the range of 60 |iJ to 100 J the focus downward rate had the most significant effect on the entrance diameter. The single-pulse energy effect has the most significant effect on the exit diameter and the roundness. In addition, a BP neural network was used to establish a mapping relationship based on the laser parameters and characteristics of micro-holes with an average error of less than 5 % and the maximum error of less than 7 %. Thus, the results of the laser drilling experiments were predictable, and this process was convenient for optimizing the process parameters, which could be applied to industrial practice. Further research may concentrate on exploring the effect of other laser processing parameters, and different modelling methods for constructing the mapping relationship between the laser parameters and the geometrical morphologies of the micro-holes. Acknowledgement The authors are grateful for the financial support provided by the National Natural Science Foundation of China (grant number 51705440), the Fundamental Research Funds for the Central Universities XMU (grant number 20720180072), the Aeronautical Science Foundation of China (grant number 20170368001), the Shenzhen Fundamental Research Program (grant number JCYJ20170818141303656), and the Natural Science Foundation of Fujian Province, China (grant number 2019J01044). References [1] Huang, H., Yang, L.-M., Liu, J. (2014). Micro-hole drilling and cutting using femtosecond fiber laser, Optical Engineering, Vol. 53, No. 5, Article No. 051513, doi: 10.1117/1.OE.53.5.051513. [2] Rihakova, L., Chmelickova, H. (2017). Laser drilling of alumina ceramics using solid state Nd: YAG laser and QCW fiber laser: Effect of process parameters on the hole geometry, Advances in Production Engineering & Management, Vol. 12, No. 4, 412-420, doi: 10.14743/apem2017.4.268. [3] Shaegh, S.A.M., Pourmand, A., Nabavinia, M., Avci, H., Tamayol, A., Mostafalu, P., Ghavifekr, H.B., Aghdam, E.N., Dokmeci, M.R., Khademhosseini, A., Zhang, Y.S. (2018). Rapid prototyping of whole-thermoplastic microfluidics with built-in microvalves using laser ablation and thermal fusion bonding, Sensors and Actuators B: Chemical, Vol. 255, Part 1, 100-109, doi: 10.1016/j.snb.2017.07.138. [4] Padmanabham, G., Bathe, R. (2018). Laser materials processing for industrial applications, Proceedings of the National Academy of Sciences, India Section A: Physical Sciences, Vol. 88, No. 3, 359-374, doi: 10.1007/s40010-018-0523-5. [5] Ciurana, J., Arias, G., Ozel, T. (2009). Neural network modeling and particle swarm optimization (PSO) of process parameters in pulsed laser micromachining of hardened AISI H13 steel, Materials and Manufacturing Processes, Vol. 24, No. 3, 358-368, doi: 10.1080/10426910802679568. [6] Dubey, A.K., Yadava, V. (2008). Laser beam machining - A review, International Journal of Machine Tools and Manufacture, Vol. 48, No. 6, 609-628, doi: 10.1016/j.iimachtools.2007.10.017. [7] Kamlage, G., Bauer, T., Ostendorf, A., Chichkov, B.N. (2003). Deep drilling of metals by femtosecond laser pulses, Applied Physics A, Vol. 77, No. 2, 307-310, doi: 10.1007/s0Q339-003-2120-x. [8] Zoubir, A., Shah, L., Richardson, K., Richardson, M. (2003). Practical uses of femtosecond laser micro-materials processing, Applied Physics A, Vol. 77, No. 2, 311-315, doi: 10.1007/s00339-003-2121-9. [9] Gruner, A., Schille, J., Loeschner, U. (2016). Experimental study on micro hole drilling using ultrashort pulse laser radiation, Physics Procedia, Vol. 83, 157-166, doi: 10.1016/i.phpro.2016.08.030. [10] Wang, X.C., Zheng, H.Y., Chu, P.L., Tan, J.L., Teh, K.M., Liu, T., Ang, B.C.Y., Tay, G.H. (2010). Femtosecond laser drilling of alumina ceramic substrates, Applied Physics A, Vol. 101, No. 2, 271-278, doi: 10.1007/s00339-010-5816-8. [11] Yang, L., Kong, X., Wang, Y., Ding, Y., Zhang, H., Chi, G. (2016). Laser micro-holes machining technology and its application, Aeronautical Manufacturing Technology, No. 19, 271-278, doi: 10.16080/i.issn1671-833x.2016.19. 032. Advances in Production Engineering & Management 14(4) 2019 419 Yin, Wu, Dong, You, Liao [12] Liu, Y., Zhang, R., Li, W., Wang, J., Yang, X., Cheng, L., Zhang, L. (2018). Effect of machining parameter on femtosecond laser drilling processing on SiC/SiC composites, The International Journal of Advanced Manufacturing Technology, Vol. 96, No. 5-8, 1795-1811, doi: 10.1007/s00170-017-1163-7. [13] Dausinger, F. (2002). Femtosecond technology for precision manufacturing: Fundamental and technical aspects, In: Proceedings of Third International Symposium on Laser Precision Microfabrication, Osaka, Japan, doi: 10.1117/ 12.486506. [14] Abeln, T., Radtke, J., Dausinger, F. (1999). High precision drilling with short-pulsed solid-state lasers, In: Laser Institute of America - Proceedings - LIA, Vol. 88, 195-203, doi: 10.2351/1.5059302. [15] Liao, C., Anderson, W., Antaw, F., Trau, M. (2018). Maskless 3D ablation of precise microhole structures in plastics using femtosecond laser pulses, ACS Applied Materials & Interfaces, Vol. 10, No. 4, 4315-4323, doi: 10.1021/ acsami.7b18029. [16] Kraus, M., Ahmed, M.A., Michalowski, A., Voss, A., Weber, R., Graf, T. (2010). Microdrilling in steel using ultrashort pulsed laser beams with radial and azimuthal polarization, Optics Express, Vol. 18, No. 21, 22305-22313, doi: 10.1364/0E.18.022305. [17] ISO 1101:2017 Geometrical product specifications (GPS) — Geometrical tolerancing — Tolerances of form, orientation, location and run-out, from https://www.iso.org/standard/66777.html. accessed October 27, 2019. [18] Zhu, J., Chew, D.A.S., Lv, S., Wu, W. (2013). Optimization method for building envelope design to minimize carbon emissions of building operational energy consumption using orthogonal experimental design (OED), Habitat International, Vol. 37, 148-154, doi: 10.1016/j.habitatint.2011.12.006. [19] Wang, C., Xue, S., Chen, G., Luan, D., Wang, S., Wang, Y., Wang, S., Liu, J., Wang, Z., Zhang, P. (2018). Influence of laser parameters on micro-hole drilling of Cu50Zr50 amorphous alloys foil, Ferroelectrics, Vol. 523, No. 1, 61-66, doi: 10.1080/00150193.2018.1391557. [20] Fan, N.-N., Xia, Z.-D., Sun, X.-Y., Hu, Y.-W. (2016). Experimental study on stainless steel micro-hole trepanned by femtosecond laser, Laser & Infrared, Vol. 46, No. 10, 1200-1205, doi: 10.3969/j.issn.1001-5078.2016.10.006. [21] Fornaroli, C., Holtkamp, J., Gillner, A. (2013). Laser-beam helical drilling of high quality micro holes, Physics Procedia, Vol. 41, 661-669, doi: 10.1016/j.phpro.2013.03.130. [22] Wang, G.-A., Zhang, Y.-Z., Ni, X.-W., Lu, J. (2007). Effect of deviation distance to focal spot on nanosecond-pulsed-laser drilling rates in air, Chinese Journal of Lasers, Vol. 34, No. 12, 1621-1624. [23] Zou, Z.-Q., Li, J., Hu, L.-Y. (2017). Diameter changing regularity with the laser parameters of nanosecond laser drilling, Optics & Optoelectronic Technology, Vol. 15, No. 5, 58-61. [24] Verbeeck, J., Bertoni, G., Schattschneider, P. (2008). The Fresnel effect of a defocused biprism on the fringes in inelastic holography, Ultramicroscopy, Vol. 108, No. 3, 263-269, doi: 10.1016/j.ultramic.2007.06.007. [25] Wu, X., Leung, D.Y.C. (2011). Optimization of biodiesel production from camelina oil using orthogonal experiment, Applied Energy, Vol. 88, No. 11, 3615-3624, doi: 10.1016/j.apenergy.2011.04.041. [26] Li, X.J., Dong, Y.W., Yin, C.P., Zhao, Q., You, Y.C. (2018). Geometric parameters evolution experiment of hole during femtosecond laser helical drilling, Chinese Journal of Lasers, Vol. 45, No. 5, doi: 10.3788/CIL201845.0502008. [27] Ren, N., Zhang, L., Wang, H., Xia, K., Shi, C. (2017). Orthogonal experiments and variance analysis in Nd:YAG pulsed laser trepanning drilling, Laser & Optoelectronics Progress, Vol. 54, No. 6, doi: 10.3788/LQP54.061408. [28] Dhara, S.K., Kuar, A.S., Mitra, S. (2008). An artificial neural network approach on parametric optimization of laser micro-machining of die-steel, The International Journal of Advanced Manufacturing Technology, Vol. 39, No. 1-2, 39-46, doi: 10.1007/s00170-007-1199-1. [29] Casalino, G. (2018). [INVITED] Computational intelligence for smart laser materials processing, Optics & Laser Technology, Vol. 100, 165-175, doi: 10.1016/j.optlastec.2017.10.011. [30] Majumder, A. (2010). Comparison of ANN with RSM in predicting surface roughness with respect to process parameters in Nd:YAG laser drilling, International Journal of Engineering Science and Technology, Vol. 2, 51755186. [31] Li, M., Wu, H., Wang, Y., Handroos, H., Carbone, G. (2017). Modified Levenberg-Marquardt algorithm for back-propagation neural network training in dynamic model identification of mechanical systems, Journal of Dynamic Systems, Measurement, and Control, Vol. 139, No. 3, Article No. 031012, doi: 10.1115/1.4035010. [32] Guo, Q.-C., He, Z.-F. (2014). Economic forecasting model based on artificial neural network, Computing Technology and Automation, Vol. 33, No. 1, 132-136. [33] Zhang, Y., Gao, X., Katayama, S. (2015). Weld appearance prediction with BP neural network improved by genetic algorithm during disk laser welding, Journal of Manufacturing Systems, Vol. 34, 53-59, doi: 10.1016/j.jmsy.2014. 10.005. [34] Ding, S., Su, C., Yu, J. (2011). An optimizing BP neural network algorithm based on genetic algorithm, Artificial Intelligence Review, Vol. 36, No. 2, 153-162, doi: 10.1007/s10462-011-9208-z. [35] Ai, J.L., Yang, X.Z. (2017). Fault diagnosis of aero-engine based on self-adaptive neural network, Scientia Sinica Technologica, Vol. 48, No. 3, 326-335, doi: 10.1360/N092017-00224. 420 Advances in Production Engineering & Management 14(4) 2019 Advances in Production Engineering & Management Volume 14 | Number 4 | December 2019 | pp 421-434 https://doi.Org/10.14743/apem2019.4.338 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper Advanced risk assessment in reverse supply chain processes: A case study in Republic of Serbia Aleksic, A.a, Runic Ristic, M.b, Komatina, N.a*, Tadic, D.a aUniversity of Kragujevac, Faculty of Engineering, Kragujevac, Serbia bUnion - Nikola Tesla University, Faculty of Management, Sremski Karlovci, Serbia A B S T R A C T A R T I C L E I N F O Management of a reverse supply chain (RSC) often takes place in an uncertain environment, so it is supposed to be analyzed through the proactive approach for avoidance/elimination of risks. Management initiatives based on the assessed risk level and priority of potential failure mode (PFM) should lead to the increase of business effectiveness, the competitive advantage and sustain-ability of the RSC. Therefore, the focus of this research is set to proposing the reliable method that would be user-friendly and suitable for the determination of risk level and priority of PFMs in RSC. Uncertainties related to the severities of Potential Effect(s) of Failure (PEF) and their frequencies', as well as detection of PFMs are described by pre-defined linguistic expressions and modelled by the interval type-2 trapezoidal fuzzy numbers (IT2TrFNs). The assessment of the relative importance of risk factors is set as a fuzzy group decision-making. The weights vector is calculated based on the procedure of fuzzy number comparison. The value of each risk factor at the level of each PFM is assessed through the predefined linguistic expressions modelled by IT2TrFNs. The rank is obtained by modified Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) method. The proposed model is tested on a real-life data from RSC that operates in Serbia. In the domain of practical implications, it may be noticed that the application of the proposed model could decrease the influence of potential causes of failures modes on the overall RSC business activities especially in the terms of strategic management and human resource practices. The novelty of the proposed model may be underlined as it is used for the analysis of different RSC activities and many interconnected issues may be solved by the proposed management measures after conducted analysis. © 2019 CPE, University of Maribor. All rights reserved. Keywords: Reverse supply chain; Risk; Multi-criteria decision analysis; Interval type-2 trapezoidal fuzzy numbers; Fuzzy FMEA framework; Fuzzy TOPSIS *Corresponding author: nkomatina@kg.ac.rs (Komatina, N.) Article history: Received 2 June 2019 Revised 28 November 2019 Accepted 2 December 2019 1. Introduction In the last few decades, the development policies for both, developed and developing countries in the EU is mainly designed in the compliance with the objectives of national sustainable development strategies and the Draft Declaration on Guiding Principles adopted by the EU in 2005. It may be assumed that RSCs play a key role in achieving sustainable development goals. In the scope of activities undertaken as a part of the strategies for achieving sustainable development goals, many crucial issues may arise that need to be solved. Reducing or eliminating those issues that may occur in the recycling processes (RPs) are one of the most important tasks of the managers in the considered economic domain. In the research domains and practice, there are many methods for identification and elimination of potential failures. 421 Aleksic, Runic Ristic, Komatina, Tadic Pervasive changes in the business environment have led to the emergence of a growing number of various types of uncertainties in the RSC's activities. It may be suggested that uncertainties with the greatest influence on the RSC's activities embrace the different aspects of the risk so they may be denoted as causes or Potential Cause(s)/Mechanism(s) of Failure (PCMF). In the literature, the PCMFs are defined at different levels, but they are mostly described in a comprehensive manner at the operational level [1-3]. Consequently, PCMFs may lead to failures. One of the most used methods at the level of industrial processes for the failure analysis is Failure Mode and Effects Analysis (FMEA) [4]. The conventional FMEA is realized through the two steps: (a) each identified failure can be evaluated by three risk factors, which are the severity (S), the occurrences of failure realization (O) and the difficulty of failure detection (D), and (b) identifying the PCMFs that lead to the occurrence of PFMs and determination of management initiatives that can be used to eliminate the influence of the identified PCMFs. Nowadays, the RSC is rather vulnerable and exposed to high levels of risk, so the motivation of this research is to rank the identified PFMs in RSC with respect to S, O, and D simultaneously, as well as their weights. Nevertheless, some researchers believe that the main disadvantages of conventional FMEA are: (a) that the crisp values of risk factors are implemented, and (b) that the relative importance among risk factors are neglected [5]. In order to improve the failure analysis process in many papers, the FMEA is combined with fuzzy sets theory and various fuzzy Multi-Criteria Decision Making Methods (MCDMs) [5, 6]. Utilizing the type-1 fuzzy sets may not be suitable in these cases, so different mathematical tools needed to be improved. The development of mathematical theory, in particular the interval type-2 fuzzy sets (IT2FS), made it possible for imprecisions and uncertainties to be more adequately described in a quantitative sense. The stated facts are aligned with the motivation for this research since the accuracy of the failure analysis may be increased if IT2FS is applied in the scope of modified FMEA. The intention of the obtained data is to provide solid support for the decision-making process in the scope of management, especially in the scope of human resource management [7]. In the scope of research, the modelling of linguistic variables is conducted by an application of the IT2FS [8]. The successfully solved issues within the context of risk management that favors usage of IT2FS can be found in the literature [9]. At the same time, it is worth to mention that the extension of MCDM with IT2FS has been the subject of many researchers' considerations [10]. The main objective of this research is to develop a model that integrates FMEA and fuzzy TOPSIS method with IT2TrFNs for the purpose of proactive analysis and mitigation the PCMFs. According to the obtained rank of PFMs as well as the PCMFs analysis, where PCMFs imply the failure ranked at the first place, the appropriate management initiatives for reduction or elimination of the PCMFs' influence may be defined. These activities may be propagated to the failures' elimination. Lastly, it can be said that the realization of the research goal may be useful in supporting actions regarding national sustainable development strategy. The paper is organized in the following manner. Section 2 presents an overview of the scientific papers from the literature sources. Description of the PFMs that may occur in RP at the level of RSC, handling of possible uncertainties within the S, O, and D, as well as their relative importance, and the proposed algorithm are shown in the Section 3. In the Section 4, the proposed model is illustrated with the real-life data, which originate from one RSC that operates in Central Serbia. Conclusions are presented in the Section 5. The list of abbreviations used in this paper is given in Appendix A. 2. Literature review Improving the RP is one of the most important tasks of each RSC operational management and it can be realized in different ways. Human resource management (HRM) practices aligned with the organizations' strategy can alleviate uncertainties in the RSC's activities, especially in the domain of employees' skills and abilities. HRM activities that are directed towards the enhancement of employees' skills and knowledge results in higher productivity and organizational per- 422 Advances in Production Engineering & Management 14(4) 2019 Advanced risk assessment in reverse supply chain processes: A case study in Republic of Serbia formance [11]. In compliance with the stated, the considered problem is complex and it may be realized through the following steps: (a) identification of PFMs that may occur in RP of one RSC, as well as PCMFs which induce identified PFMs, (b) assessment of the relative importance of S, O, and D, (c) assessment of the values of S, O, and D at the level of each PFM, (d) ranking of the identified PFMs by using the proposed fuzzy TOPSIS, and (e) the determination of suitable measures for reducing or eliminating failures. In the literature, there are many papers where failure analysis is based on FMEA combined with fuzzy sets theory and MCDM [12]. Further, a summary of the literature for both marked sub-problems is provided. There are many approaches for determination of the PCMFs, the PFMs and their PEFs which can be found in the literature. The simplest procedure for identification of the PFMs, the PCMFs and PEFs is a checklist. The effectiveness of this procedure depends on the knowledge, skills and expertise of the decision-makers (DMs). A large number of authors have tried to define the list of PCMFs for a specific group of manufacturing plants, especially when the resources are limited [13, 14]. In addition, the PFMs and their PEFs are identified, with respects to DMs assessment. In some papers, the risk factors values are described by crisp values [15], as it has been proposed in the conventional FMEA [4]. It should be underlined that this assumption is not sufficient in practice since DM cannot express their estimates well enough if they use precise numbers. Some authors suggest that the imprecise numbers could be appropriate for modeling of linguistic expressions that are used to describe the risk factors' values. Development of specific mathematical theories, particularly of fuzzy sets theory [16, 17] has allowed adequate quantitative description for these uncertainties. These uncertainties could be modelled by type-1 fuzzy numbers [18, 19]. Representation of a type-1 fuzzy set is by its membership function whose parameters are of the shape and the location in the universe of discourse. For the determination of the membership function, DMs knowledge or experience is applied. The well-known fact is that the linguistic terms used to describe uncertainty, generally, have a different meaning for different people [8]. Also, determination of the membership function for a larger number of uncertainties is performed under complex an uncertain environment. The concept of a type-2 fuzzy sets demand a large number of complex calculations, consequently, they do not have a wider application in modelling of real uncertainties. The IT2FS may be seen as a special case of a type-2 fuzzy set [8]. The computational processing with the IT2FS is reduced compared to the type-2 fuzzy sets so they are widely employed to solve different decision-making problems [10]. A significant number of uncertainties that exist in various problems could be modelled by the IT2FS [20-22]. In this paper, the risk factors' values are modelled by the interval type trapezoidal fuzzy numbers (IT2TrFNs). Defining the priority of the PFMs by using the procedure proposed in the conventional FMEA [4] has numerous shortcomings. One of the main disadvantages according to the opinion of some researchers is that it neglects the relative importance among risk factors [5]. The PFMs' ranges could be adequately determined by means of FMEA combined with MCDM. In the literature, there are many papers where FMEA is combined with different MCDM with type 1 fuzzy sets [5, 6]. By respecting the assumption that existing uncertainties in risk assessment problem are far more accurately modelled with IT2FS, many authors have expanded the proposed MCDM with IT2FS [10]. Furthermore, special attention is focused on the papers that proposed TOPSIS with IT2FS [20, 23-26]. The comparative analysis of the proposed TOPSIS with IT2TrFN with other similar models is described and presented in Table 1. Taking in account the stated, it can be argued that when the model presented in this paper is compared with the models from the literature, the certain advantages in terms of describing the real problems may be noticed. In the literature, there is a certain number of papers where FMEA is combined with MCDM with IT2FS [23]. In this paper, the ratings of the identified PFMs in FMEA, and their relative importance are modelled by IT2FS. The rank of the identified PFMs with respects to the values and the relative importance is obtained by applying grey relational analysis. In this paper, integrated FMEA and fuzzy TOPSIS with IT2TrFNs have been proposed for the analysis and ranking of PFMs in RPs. Advances in Production Engineering & Management 14(4) 2019 423 Aleksic, Runic Ristic, Komatina, Tadic Table 1 The summarized comparative analysis Chen and Lee, Ghaemi Nasab and Kahraman and Temur et al. Zamri and The proposed 2010 [23] Rostamy- Sari, 2012 [20] 2014 [25] Abdullah, model Malkhalifeh, 2010 2014 [26] [24] The linguistic variables IT2TrFN IT2TrFN IT2TrFN IT2TrFN IT2TrFN IT2TrFN for the relative im- portance of criteria Granularity/relative 7 7 7 7 7 3 importance of criteria Domain of linguistic [0-1] [0-1] [0-1] [0-1] [0-1] [1-5] variables for the rela- tive importance of criteria description Determination of the Fuzzy averag- Fuzzy averaging Ranking of Fuzzy averag- Assessment Fuzzy averag- criteria weights ing method method IT2TrFNs ing method ing method The linguistic variable IT2TrFN IT2TrFN IT2TrFN IT2TrFN IT2TrFN IT2TrFN for the criteria values Granularity/criteria 7 7 7 7 7 7 values Domain of linguistic [0-10] [0-10] [0-10] [0-1] [0-1] [0-1] variable for describing the criteria values PIS/NIS Ranking of Fuzzy averaging The procedure Ranking of Ranking of - IT2TRFNs [27] method; from the con- IT2TRFNs [27] IT2TRFNs Ranking of IT2TRFNs ventional [27] [27] TOPSIS distance Normalized Normalized Euclide- Normalized Normalized - Based on con- Euclidean an distance Euclidean Euclidean ventional distance distance distance TOPSIS and Fuzzy algebra Closeness coefficient Conventional Conventional TOPSIS Conventional Conventional - Based on con- TOPSIS TOPSIS TOPSIS ventional TOPSIS and Fuzzy algebra Rank Conventional Conventional TOPSIS Conventional Conventional Based on the Defuzzification TOPSIS TOPSIS TOPSIS procedure for procedure [28] ranking of and conven- IT2FNs [27] tional TOPSIS The priority of management initiatives could be determined by exact methods [15]. Many authors believe that optimal risk management should be based on the determined risk level [19]. It is assumed that the order of management initiatives realization could be based on the rank of PFMs and analysis of the frequency of PCMFs which lead to the realization of high priority PFMs. In other words, by applying this method, answers to two questions could be obtained: (a) Which impact of PCMFs that needs to be reduced? and (b) How much it needs to be reduced? In this way, the effectiveness of risk management increases, while at the same time costs are reduced, which is further propagated to increase the effectiveness of RSC's business and its competitiveness. 3. Materials and methods Economic development of each state should be should be aligned with its sustainable development. This can be achieved, inter alia, through the continuous realization of the RPs in RSCs. For the activities of the risk assessment in RPs, the decision-making team could be defined in compliance with the needs of the treated enterprise. In the presented research, it is scoped to three members: top manager, manager of the recycling process, and logistics manager. 3.1 Definition of the finite set of PCMFs According to the results in practice, it is perceived that many possible PCMFs may impact one or several PFMs in RSC operating processes. Generally, identified PCMFs may be presented by the set of indices i = {1, ...,i,...,/}, where I present the total number of PCMFs, and the index of each PCMFs is denoted as i,i = 1,...,/. In this paper, PCMFs have been analyzed according to the existing literature [13, 14], and scoped to the operational level in RSC: Demand and supply uncer- 424 Advances in Production Engineering & Management 14(4) 2019 Advanced risk assessment in reverse supply chain processes: A case study in Republic of Serbia tainty (i = 1) - this factor is related to the situations where the unreliable and uncertain resources create supply chain interruption, thus producing risk that can be explained as uncertainty between supply and demand [29]; Failure to select the right suppliers (i = 2) - represents the factor that could lead to downside of the RSC business strategy, its vision for future cooperation, it could increase the RSC lack of the expertise and experience, and it may lead to lower level of competitiveness, quality and reputation [30] as well as other issues in interconnected systems [31]; Lower responsiveness performance (i = 3) - this factor is related to the operative structure of RSC, so it may be assumed that if the RSC operates as a several closed chains it could be less responsive [32]; Inflexibility of supply source (i = 4) - if the RSC is vulnerable to this factor, then it is not able to respond quickly and efficiently to an inconsistent RSC demands and these inflexibilities may arise due to the rate of changes and uncertainties in the business environment [33]; Poor quality or process yield at supply source (i = 5) - this factor is related to the issues that can appear in supply cost, delivery and quality changes, supply base reduction, failure to comply with long-term RSC partners, poor information exchange, poor suppliers technical capabilities, etc. [2]; Coordination complexity/effort (i = 6) - this factor exposes RSC to the most challenging coordination issues, such as: customers' diversity and their different demands, different resources, unanticipated change and level of goal difficulty among RSC members and the customer may lead to an extra coordination burden due to information inaccuracy, different goals of RSC members, or disputes between the partners [34]; Information technology (IT) and information sharing risks (i = 7) - this factor is related to the lack of necessary IT infrastructure and mechanism to capture and propagate information among RSC members in a timely manner [3]; The lack of sustainable knowledge/technology (i = 8) - this factor may arise when there is discrepancy between forecast and actual demand, or lack of sound knowledge and understanding about sustainable technology, operations and methods among partners appear [29]. 3.2 Definition of the finite set of PFMs It may be comprehended that the realization of each PCMFs i,i = 1,...,/, at the level of RSC could lead to the occurrence of one or more PFMs in RPs. These PFMs are identified by DMs with respects to the evidence data, experience and the results of benchmarking analysis. In other words, it is considered that the DMs should define a list of all PFMs that may formally be represented with a set of indices j = {1,...,_/,...,/},_/ = 1,...,]. The total number of identified PFMs is indicated as J, and j,j = 1,...,/ is an index of PFM. For the purpose of research presented in this paper, the PFMs that can be found in operating processes of the treated RSC that is a part of one RSC are: unfair competition impact (j = 1), high level of stock in volatile prices' presence (j = 2), inattention or neglecting working procedures (j = 3), failures of transport vehicles (j = 4), error or discontinuity in information flow (j = 5), failures regarding lack of human skills and knowledge (j = 6), failures induced by natural disasters effects (j = 7), generation of hazardous waste (j = 8), ineffective resource consumption (j = 9), losses induced by inflation and exchange rate differences (j = 10), and failures resulting in hazardous environment work (j = 11). 3.3 Definition of the finite set of risk factors Evaluation and ranking of identified PFMs may be performed in terms of a certain number of risk factors that can be formally presented as a set of risk factors indices = {1, ...,k, ...,K},k = 1,...,K. The total number of risk factors is denoted as K, and k,k = 1,...,K is the index for each risk factor. In the literature, numerous authors suggest that the selection of risk factors for PFMs evaluation in any research domain should be based on FMEA. Respecting this fact, in this paper, identified PFMs of RP are evaluated according to the three risk factors S, O, and D. 3.4 Selection of appropriate linguistic expressions for the assessment of the relative importance of risk factors It is assumed that the determination of the severity, occurrence and detection weights should be stated as the fuzzy group decision-making problem. Each DM could describe the relative im- Advances in Production Engineering & Management 14(4) 2019 425 Aleksic, Runic Ristic, Komatina, Tadic portance of considered risk factors by one of the three defined linguistic expressions which are modelled by the IT2TrFNs: low importance (L) - ((1,1,2,3; 1.1), (1,1,2,2.5; 0.7,0.7)) medium importance (M) - ((1.5,2.5,3.5,4.5; 1,1), (2,2.5,3.5,4; 0.7,0.7)) high importance (H) - ((2,4,5,5; 1.1), (2.5,4,5,5; 0.7,0.7)) The domain values of the IT2TrFNs are defined on the interval [1-5]. The value 1 denotes that risk factor k,k = 1, ...,K has an almost negligible influence in the evaluation of PFMs, and the value 5 means that risk factor k,k = 1, ...,K has a conspicuously extreme big influence in the evaluation of PFMs, respectively. 3.5 Choice of appropriate linguistic expressions for the assessment of the severity and detection of PFMs in RP In general, severities of the manifested PEFs and the detection possibility of the identified PFMs can be determined according to the decision-makers' assessment. Detection of the identified PFMs depends on: (1) the level of security control of RPs and employees in the RSC, and (2) the level of automation of the applied control procedures. There are no defined rules or recommendations on how to determine the number and type of linguistic expressions that DMs should use to describe the values of S and D. The number and type of linguistic expressions depends on the type and the size of the problem as well as on the applied concept of control in the considered RSC. In this paper, DMs expressed their assessments of values S and D at the level of each identified PFM with seven pre-defined linguistic expressions that are modelled by the IT2TrFNs and presented in Table 2. It is worth to mention that defined expressions take into account the conventional FMEA analysis and present state in recycling centers. Table 2 The severity of the PEFs and detection of PFMs during the implementation of the RP Linguistic expres- IT2TrFNs Description of severity Description of detection sions_ Almost no dangers/Certain chance of detection (L1) Low danger/ Very high chance of detection (L2) Low to moderate danger/High chance of detection (L3) Moderate danger/Moderate chance of detection (L4) Moderate to high danger/Remote chance of detection (L5) High danger/Unreliable chance of detection (L6) Extremely high danger/No chance of detection (L7) _out any prior warning ((0,0,0.1,0.25; 1.1), (0,0,0.1,0.2; 0.8,0.8)) Failure has no impact on There are constraints that the recycling process prevent failure ((0,0.15,0.25,0.4; 1.1), (0.05,0.15,0.25,0.35; 0.8,0.8)) There is little or no effect In the recycling process, on the continuity of the the principle of zero fail- recycling process ures is applied or the control is automated ((0.1,0.25,0.35,0.5; 1.1), (0.15,0.25,0.35,0.45; 0.8,0.8)) Failure could result in In the recycling process, minor recycling process the principle of zero fail- problems that can be ures is applier or the overcome with minor control is not automatized modifications to the recycling process ((0.3,0.45,0.55,0.7; 1.1), (0.35,0.45,0.55,0.65; 0.8,0.8)) Failure could result in There is a process for minor customer dissatis- inspections on the sample, faction and/or major but it is not automated recycling process prob- and/or relies on vigilance lems of employers and workers ((0.5,0.65,0.75,0.9; 1.1), (0.55,0.65,0.75,0.85; 0.8,0.8)) Failure could result in a The error can be detected high degree of customer with manual inspection, dissatisfaction and/or but no inspection process major recycling process is in place so that failure problems which require can be detected accidental- reworking ly ((0.6,0.75,0.85,1; 1.1), (0.65,0.75,0.85,0.95; 0.8,0.8)) Failure could result in a The failure can be detected serious recycling process only with a thorough disruption with an inspection, and this is not interruption in service, feasible or cannot be with prior warning readily performed ((0.75,0.9,1,1; 1.1), (0.8,0.9,1,1; 0.8,0.8)) Failure could result in a There is no known ap- breakdown of the total proach for failure detec- recycling process, with- tion 426 Advances in Production Engineering & Management 14(4) 2019 Advanced risk assessment in reverse supply chain processes: A case study in Republic of Serbia The domains of these IT2TrFNs are defined by the interval from 0 to 1. The value 0 denotes that the severity of a PEF is extremely low, and the value 1 denotes that severity of a PEF is extremely high, respectively. Similarly, values are defined in the domain of IT2TrFNs when detection is considered. The value 0, i.e. the value 1 denotes that it is quite certain that the PFM would be detected, or that it is quite certain that there is no possibility of PFM detection, respectively. The above-described procedure for the definition of IT2TrFNs domains, both for severity and for detection, allows a reduction in the complexity and volume of computing since it is not necessary to apply a normalization process nor to take into account the type of these two risk factors according to which PFMs are evaluated. 3.6 Choice of appropriate linguistic expressions for the assessment of the occurrence of PFMs in RP If there is an automated failure detection system, or if there is an up-to-date database of all failures that occur in the considered RSC, then it is possible to calculate the probability of PFM occurrence. In the RSCs that is the subject of research, it is not possible to calculate the probability of occurrence of RPs' PFMs accurately. Hence, in these cases, DMs may assess the possibility of each identified PFM occurrence based on their experience and benchmarking with similar RPs by using pre-defined linguistic expressions. The definition of these linguistic expressions is given according to the assessment of DMs and the relevant literature [18]. The used linguistic expressions are modelled by IT2TrFNs as shown in the Table 3 and they take into account the conventional FMEA analysis and present state in recycling centers. The domains of these IT2TrFNs should be defined into the interval from 0 to 1. The value 0 denotes that the possibility of PFM occurrence is the lowest, and the value 1 denotes that the possibility of PFM occurrence is highest, respectively. In this way, it is not necessary to perform a normalization procedure. Table 3 The occurrence rating scale Linguistic expressions IT2TrFNs Description of severity Low possibility of occurrence (O1) Moderate possibility of occurrence (O2) High possibility of occurrence (O3) Very high possibility of occurrence (O4) Failure is almost inevitable (O5) ((0,0,0.1,0.25; 1.1), (0,0,0.1,0.2; 0.8,0.8)) ((0.1,0.25,0.35,0.5; 1.1), (0.15,0.25,0.35,0.45; 0.8,0.8)) ((0.3,0.45,0.55,0.7; 1.1), (0.35,0.45,0.55,0.65; 0.8,0.8)) ((0.5,0.65,0.75,0.9; 1.1), (0.55,0.65,0.75,0.85; 0.8,0.8)) ((0.75,0.9,1,1; 1.1), (0.8,0.9,1,1; 0.8,0.8)) Failure occurs rarely, or failure occurs about once per year. Failure occurs occasionally, or failure occurs once every 3 months Failure occurs approximately once per month Failure occurs frequently, or failure occurs about once per week. Failure occurs at least once a day, or failure occurs almost every time. 3.7 The proposed fuzzy TOPSIS with the IT2TrFNs The proposed research Algorithm may be realized through a certain number of steps, which are further clarified. Step 1. Fuzzy rating of the relative importance of each considered risk factor is given by each DM. #fce,k = 1.....K, e = 1.....E Step 2. The aggregated relative importance of considered risk factors are calculated by using the fuzzy averaging method: Wi 41' e = l (1) Step 3. The determination of the risk factors weight is based on the ranking values procedure of IT2TrFNs [23]. The likelihood p ,k, k' = 1, ...,K is determined: p (wf = max((1-max(S, 0)),0) (2) where: Advances in Production Engineering & Management 14(4) 2019 427 Aleksic, Runic Ristic, Komatina, Tadic which has the following properties: 0#^)<1 (3) v (#fcu >#fe^) + p (W" >#fcy) = 1 (4) P (f^k — ^fc7) = 0.5 (5) The likelihood p (w^ j ,k,k' = 1,...,K can be considered in a similar way. Determination of the upper fuzzy preference matrix, and the lower fuzzy preference matrix can be presented as: pu = \p and pL = \p (#fcL >#fcL,)l (6) L v L V slKXK The ranking value Rank and the ranking value Rank (w^) are calculated by the proce- dure developed by [27]: K Rank Ront =^.(^-1) • (Ép +2 -1) m The ranking values of IT2TrFN Wk, Rank (wj^ can be calculated by using the following expression: Rank (#fc) = 2 ■ (Rank (#fcu) + Ran/c (9) Step 4. Each severity, the frequency of occurrence, and the possibility of each identified error detection are assessed, xJk,k = 1, ...,K,j = 1,...,/. Step 5. The fuzzy weighted decision matrix [zik]IxK is constructed, where Zjk = wk •xik,k = 1.....K,i = 1.....1. Step 6. The fuzzy Positive Ideal Solution (FPIS), and Fuzzy Negative Ideal Solution (FNIS) zk are defined according to the vertex concept, so that: f+ = ((0,0,0,0; 1.1), (0,0,0,0; 1,1)) and f^ = ((1,1,1,1; 1.1), (1,1,1,1; 1,1)) Step 7. The distance from FPIS, df and FNIS, d~j~ is calculated according to the procedure which is developed in the conventional TOPSIS, so that: K K d+ = ^ l^fc -z]k | and âj \z]k -zk | (10) k=1 k=1 Step 8. The determination of closeness coefficient for each PFM j = 1,...,/: x dr Step 9. The determination of the representative scalar, Ç, of the IT2TrFNs, Ç,,y = 1, ...,/ by using the following procedure [28]: 428 Advances in Production Engineering & Management 14(4) 2019 Advanced risk assessment in reverse supply chain processes: A case study in Republic of Serbia DTraT = - ■ 2 1 (« -O + (a1 a2y+ -a3y + < + (a4 -a^) + (aL -af -a^) + (ßL • a^-a^) . -'a--'""i (12) Step 10. The crisp values, = 1,...,] are sorted in decreasing order. The rank of PFMs in RP is determined according to the obtained rank. At the first place in the rank, there is the PFM that has the greatest impact on the realization of RP. Step 11. The measures for eliminating or decreasing the PCMFs that may lead to PFMs are determined according to the obtained rank and of PFMs and frequency of PCMFs' occurrence. 4. Results and discussion: A case study in Republic of Serbia The testing of the proposed model is realized by using the real-life data from the RSC operating in the region of Central Serbia. In this research, the process of recycling does not take into account the activities out of the recycling center related to processes of purchasing and sale. The input data for the model testing is obtained through the interview with the predefined decisionmaking team of the recycling center. The team is consisted of three members: top manager, manager of the recycling process, and logistics manager. The risk factor values for each identified PFM j = 1,...,/ are assessed and presented in the Table 4 (Step 3 of the proposed Algorithm). The relative importance of risk factors is assessed by DMs (Step 1 of the proposed Algorithm): (k = 1): H, H, H (k = 2): H, M, M (k = 3): M, L, L The aggregated values of risk factors are: = ((2,4,5,5; 1,1), (2.5,4.5,5,5; 0.7,0.7)) Wj = ((2.3,3,4,4.7; 1,1), (2.2,3,4,4.3; 0.7,0.7)) W3 = ((1.2,1.5,2.5,3.5; 1,1), (1.4,1.5,2.5,3; 0.7,0.7)) According to the procedure (Step 4 of the proposed Algorithm), the weighted fuzzy decision matrix is presented in the Table 5. Step 3. The determination of the risk factors' weight is based on the procedure of the ranking values of IT2TrFNs [23]. Determine the likelihood p so that: _ max(4.7 - 5,0) + max(4- 5,0) + max(3-4,0) + max(2.3 - 2,0) + (4.7 - 2) + max(1-1,0) + max(0.7 - 0.7,0) S~ ¡4^ p(vv^ = max - max (27,0)) ,0^ = max(0.6625,0) = 0.6625 In a similar way is calculated the rest of the elements of the upper fuzzy preference matrix, pu, and the lower fuzzy preference matrix pL, so that: pu = 0.5 0.3375 0.1293 0.6625 0.5 0.1200 0.8707 0.8800 0.5 and pu = 0.5 0.2368 0.0360 0.7631 0.5 0.0976 0.9635 0.9024 0.5 Advances in Production Engineering & Management 14(4) 2019 429 Aleksic, Runic Ristic, Komatina, Tadic Table 4 The assessed risk factor values for the recycling process PFM (Potential failure mode) PEF (Potential effect(s) of failure) S PCMF (Potential cause(s)/ _Mechanism(s) of failure) O Current process controls Unfair competition impact (/ = 1) High level of stock in volatile prices' presence G = 2) Inattention or neglecting working procedures (j = 3) Failures of transport vehicles (/ = 4) Error or discontinuity in information flow (/ = 5) Failures regarding lack of human skills and knowledge G = 6) Failures induced by natural disasters effects (j = 7) Generation of hazardous waste (/ = 8) Termination of contracts with L5 customers Financial crisis of the RSC L7 Creation of economic loss due to L5 delay, bad organization and making a scratch Delay in delivery, the impossibility of L3 loading and unloading Disruption of production plans and L2 manufacturing disturbances Injury at work Ineffective resource consumption Economic losses (/ = 9) Losses induced by inflation and exchange rate differences G = 10) Failures resulting in hazardous environment work G = 11)_ Economic losses Injury at work : 1,¿ : ■ 3,i- ■■ 4, i - 2, 4 2,i = 3, 5 = 5,i = 6, = 7, i = 8 2, í ■■ 6, i ■■ l = 6, l = 3, i = 4, 7 L5 i = 6, ¿ = 7, i = ¡ The destruction of products that are L5 i water sensitive Contamination of work and envi- L7 i ronment = 2,i = 3, 5, i 1=5,1= L3 i = 5, ¿ = 6, i = 7 L3 i = 9 L4 i = 5,i = 8 03 Regular check by L2 the top manager O1 Regular check by L5 the financial manager O5 Monitoring by the L1 manager of the recycling process 04 Monitoring by the L2 logistic manager 04 Monitoring by the L2 manager of the recycling process O3 Monitoring by the L3 manager of the recycling process 01 Monitoring by the L3 logistic manager 05 Monitoring by the L1 manager of the recycling process O3 Regular check by L3 the financial manager 02 Regular check by L2 the financial manager 03 Monitoring by the L1 manager of the recycling process Rank (w^ = Rank (W^ = 3-(3 -1 1) ^(0.5 + 0.6625 + 0.8707) +3 - lj = 0.4222 \k=l ' 3 3^(3 — 1) ^(0.5 + 0.7631 + 0.9635) + 3 - 1 j = 0.4544 1 Rank (wj = - • (0.4222 + 0.4544) = 0.4383 « 0.44 The weights of the remaining two risk factors are calculated in a similar way, therefore: Rank (#2) = 0.3631 « 0.36, and Rank (#3) = 0.1986 « 0.2 By using the proposed Algorithm (Step 4 to Step 9) the closeness coefficient for each identified PFM is calculated. The proposed procedure is illustrated by example: ell = ((1,1,1,1; 1,1), (1,1,1,1; 1,1)) — ((0.22,0.29,0.33,0.39; 1,1), (0.24,0.29,0.33,0.37; 0.7,0.7)) + ((1,1,1,1; 1,1), (1,1,1,1; 1,1))-((0.11,0.16,0.20,0.25; 1,1), (0.13,0.16,0.20,0.23; 0.7,0.7)) + ((1,1,1,1; 1,1), (1,1,1,1; 1,1))-((0,0.03,0.05,0.08; 1,1), (0.01,0.03,0.05,0.08; 0.7,0.7)) = ((2.28,2.42,2.52,2.67; 1,1), (2.33,2.42,2.52,2.63; 0.7,0.7)) âî = ((0.22,0.29,0.33,0.39; 1,1), (0.24,0.29,0.33,0.37; 0.7,0.7))-((0,0,0,0; 1,1), (0,0,0,0; 1,1)) + ((0.11,0.16,0.20,0.25; 1,1), (0.13,0.16,0.20,0.23; 0.7,0.7))-((0,0,0,0; 1,1), (0,0,0,0; 1,1)) + ((0,0.03,0.05,0.08; 1,1), (0.01,0.03,0.05,0.08; 0.7,0.7))- ((0,0,0,0; 1,1), (0,0,0,0; 1,1)) = ((0.33,0.48,0.58,0.72; 1,1), (0.38,0.48,0.58,0.68; 0.7,0.7)) The closeness coefficient for PFM (J = 1) is: ((0.33,0.48,0.58,0.72; 1,1), (0.38,0.48,0.58,0.68; 0.7,0.7)) Çi = ((2.61,2.903.1,3.39; 1,1), (2.71,2.90,3.10,3.31; 0.7,0.7)) = ((0.09,0.15,0.20,0.28; 1,1), (0.11,0.15,0.20,0.25; 0.7,0.7)) D 430 Advances in Production Engineering & Management 14(4) 2019 Advanced risk assessment in reverse supply chain processes: A case study in Republic of Serbia The representative scalar, ^ of the IT2TrFNs, is calculated in the following way [28]: 1 ((0.28 - 0.09) + (1-0.15 - 0.09) + (1 • 0.20 - 0.09) =ri-i-+009 + (0.25 - 0.11) + (0.7 • 0.15 - 0.11) + (0.7 • 0.2 - 0.11) ------------- + 0.11 = 0.17 The values of identified PFMs' closeness coefficients are calculated in a similar way. Their values are presented in the Table 5. The rank of the PFMs is obtained by using the proposed Algorithm (Step 10) and given in the Table 5, too. The case study approach examines the effectiveness of the fuzzy logic approach for assessing the product and process-related PFMs within global RSC context. In this way, the management initiatives that should lead to the achievement of the zero failures' principle in the realization processes of RSCs are more precisely defined. The obtained result indicates that the PFM that may endanger the RSC and its operating processes primarily generates the hazardous waste (j = 8) and lower responsiveness performance (j = 3). Taking into the account the input data regarding PCMFs that may originate from PFMs presented in the Table 4, the analysis points two main PCMFs. The analysis of RSC management should tackle firstly ranked PFMs: poor quality or process yield at supply source (j = 5) and lack of sustainable knowledge/technology (j = 8). The poor process yield at supply source (j = 5) may be improved through the several activities. In the first place, the redefinition of supply strategy can be oriented to the enhancement of partnership relations with existing suppliers. Also, the ratio of raw materials or components may be redefined and proposed to existing suppliers. On the other hand, a management team may propose a new model or procedure of supplier selection and introduce new suppliers in compliance with that model. Besides the mentioned measures, the communication with suppliers and partners should be enhanced and continuously improved over time with the development or improvement of existing information systems and communication channels. Internally, company could perform analysis of logistics services' impacts on risk perception by employing suitable computation models [35]. The RSC may cope with the PCMF denoted as a lack of sustainable knowledge/technology (j = 8) through the enhancement of the human resources management aimed at: (a) introduction of contemporary information and communication technologies, (b) implementation of tailored recruitment and selection processes, (c) implementation of long life learning which should consequently lead to improvement of staff effectiveness, (d) enhancement of employees' motivation to learn and improve their skills constantly. These HRM practices should create human resource advantage that will positively influence not only RSC performance but overall business performance. It is worth to mention that other cases ranked at a lower place in the presented model should be also analyzed. This is further propagated to enhance RSC knowledge level, which is the most important resource for competitive advantage achievement in a long-term period. Table 5 The closeness coefficient and rank of PFMs PFMs Rank 1 = 1 ((0.09,0.15,0.20,0.28; 1,1), (0.11,0.15,0.20,0.25; 0.7,0.7)) 0.17 5 1 = 2 ((0.14,0.18,0.22,0.27; 1,1), (0.15,0.18,0.22,0.27; 0.7,0.7)) 0.19 3 1 = 3 ((0.15,0.20,0.25,0.30; 1,1), (0.16,0.20,0.25,0.28; 0.7,0.7)) 0.21 2 } = 4 ((0.06,0.12,0.16,0.24; 1,1), (0.08,0.12,0.16,0.21; 0.7,0.7)) 0.14 6 j = 5 ((0.05,0.11,0.15,0.22; 1,1), (0.07,0.11,0.15,0.20; 0.7,0.7)) 0.13 7-9 I = 6 ((0.10,0.16,0.21,0.29; 1,1), (0.12,0.16,0.21,0.25; 0.7,0.7)) 0.18 4 1 = 7 ((0.07,0.11,0.15,0.22; 1,1), (0.08,0.11,0.15,0.19; 0.7,0.7)) 0.13 7-9 1 = 8 ((0.18,0.23,0.28,0.31; 1,1), (0.21,0.23,0.028,0.30; 0.7,0.7)) 0.24 1 j = 9 ((0.07,0.10,0.15,0.22; 1,1), (0.07,0.10,0.15,0.19; 0.7,0.7)) 0.12 10 j = 10 ((0.02,0.07,0.11,0.18; 1,1), (0.04,0.07,0.11,0.16; 0.7,0.7)) 0.09 11 j = 11 ((0.07,0.12,0.16,0.22; 1,1), (0.08,0.12,0.16,0.20; 0.7,0.7)) 0.13 7-9 Advances in Production Engineering & Management 14(4) 2019 431 Aleksic, Runic Ristic, Komatina, Tadic 5. Conclusion The RSC has a significant influence on the development of every national economy. Effectiveness, competitiveness and long-term sustainability of RSC entities among other things could be increased and continually improved through the application of appropriate risk management activities. The determination of RSC risk management activities should be based on analytic methods, since each solution obtained in an exact way is less encumbered by the subjective decision-makers' perspectives, and therefore it may be considered as more accurate. The main contribution of this paper is the introduction of a model for the purpose of proactive analysis and mitigation the PCMFs of operational risks in an exact manner. The main contribution may be decomposed to the several component contributions: (a) definition of the appropriate linguistic expressions to adequately describe the values of S, O, and D in the domain of RSC (these linguistic terms are defined in compliance with the brainstorming and interviews with enterprise management, quality management documentation and benchmark analysis in the treated economy field), (b) modelling of all existing uncertainties is performed by applying IT2TrFNs, (c) determination of the rank of PFMs is performed by using the proposed method that integrates FMEA, IT2FS, and fuzzy TOPSIS. The proposed method is flexible to the changes in the numbers of PCMFs and PFMs as well as in the shape of membership functions of fuzzy numbers and it can be easily incorporated into the proposed model. The proposed model was tested on real-life data from RSC which operates as an RSC entity in Central Serbia. This paper contributes to both, practice and research. Practical contribution can be addressed to the possibility of decreasing the PCMFs' influence on the RSC business effectiveness and resource consumption. It may be noticed that the model is easy to be utilized, and as such, in the long-term it may be very useful for decision-makers dealing with human resource management initiatives. The main advantage of the presented model presents its convenience to use the risk assessment methodology and to determine RSC management initiatives which should enable avoiding or minimizing the influence of PCMFs. Moreover, the model incorporates human resources management practices that make differentiation by developing the specific pool of human capital that leads to a competitive advantage and sustainable RSC. It can be said that the proposed model describes the considered problem significantly better in comparison to the developed models which can be found in the literature. A wide range of evidence data sometimes cannot be easily found which presents the general limitation of the proposed model. The majority of business entities within RSC perform the data digitalization; it is realistic to assume that many data which exist in the proposed model will exist in the digital form. This will allow easier assessment of the occurrence of failures by applying the advanced statistical methods and neural networks. This would lead, consequently, to the enhancement of the proposed model and the alleviation of some drawbacks of the proposed approach. Future research should cover analysis of RSCs which exist in different industrial domains. This should bring the model validation. Moreover, the development of the software solution could provide significantly more efficient management of the operational risk in SC, which will be the subject of the future research paper. Acknowledgement Research presented in this paper was supported by Ministry of Science and Technological Development of Republic of Serbia, Grant TR-35033, Title: Sustainable development of motor vehicles recycling technology and equipment. References [1] Tang, O., Musa, S.N. (2011). Identifying risk issues and research advancements in supply chain risk management, International Journal of Production Economics, Vol. 133, No. 1, 25-34, doi: 10.1016/i.iipe.2010.06.013. [2] Tummala, V.M.R., Schoenherr, T. (2011). An implementation decision framework for supply chain management: A case study, International Journal of Logistics Systems and Management, Vol. 8, No. 2, 198-213, doi: 10.1504/ IILSM.2011.038603. 432 Advances in Production Engineering & Management 14(4) 2019 Advanced risk assessment in reverse supply chain processes: A case study in Republic of Serbia [3] Dubey, R., Gunasekaran, A., Papadopoulos, T., Childe, S.J., Shibin, K.T., Wamba, S.F. (2017). Sustainable supply chain management: Framework and further research directions, Journal of Cleaner Production, Vol. 142, Part 2, 1119-1130, doi: 10.1016/j.jclepro.2016.03.117. [4] Stamatis, D.H. (2003). Failure mode and effect analysis: FMEA from theory to execution, American Society for Quality, Quality Press, Milwaukee, USA. [5] Adhikary, D.D., Bose, G.K., Bose, D., Mitra, S. (2014). Multi criteria FMECA for coal-fired thermal power plants using COPRAS-G, International Journal of Quality & Reliability Management, Vol. 31, No. 5, 601-614, doi: 10.1108/ IIQRM-04-2013-0068. [6] Liu, H.-C., Liu, L., Liu, N., Mao, L.-X. (2012). Risk evaluation in failure mode and effects analysis with extended VIKOR method under fuzzy environment, Expert Systems with Applications, Vol. 39, No. 17, 12926-12934, doi: 10.1016/j.eswa.2012.05.031. [7] Rousseau, D.M., Barends, E.G.R. (2011). Becoming an evidence-based HR practitioner, Human Resource Management Journal, Vol. 21, No. 3, 221-235, doi: 10.1111/j.1748-8583.2011.00173.x. [8] Mendel, J.M., John, R.I., Liu, F. (2006). Interval type-2 fuzzy logic systems made simple, IEEE Transactions on Fuzzy Systems, Vol. 14, No. 6, 808-821, doi: 10.1109/TFUZZ.2006.879986. [9] Qin, J., Liu, X. (2015). Multi-attribute group decision making using combined ranking value under interval type-2 fuzzy environment, Information Sciences, Vol. 297, 293-315, doi: 10.1016/j.ins.2014.11.022. [10] Celik, E., Gul, M., Aydin, N., Gumus, A.T., Guneri, A.F. (2015). A comprehensive review of multi criteria decision making approaches based on interval type-2 fuzzy sets, Knowledge-Based Systems, Vol. 85, 329-341, doi: 10.1016 /j.knosys.2015.06.004. [11] Batt, R., Colvin, A.J.S. (2011). An employment systems approach to turnover: Human resources practices, quits, dismissals, and performance, Academy of Management Journal, Vol. 54, No. 4, 695-717, doi: 10.5465/amj.2011. 64869448. [12] Liu, H.-C., Liu, L., Liu, N. (2013). Risk evaluation approaches in failure mode and effects analysis: A literature review, Expert Systems with Applications, Vol. 40, No. 2, 828-838, doi: 10.1016/j.eswa.2012.08.010. [13] Song, W., Ming, X., Liu, H.-C. (2017). Identifying critical risk factors of sustainable supply chain management: A rough strength-relation analysis method, Journal of Cleaner Production, Vol. 143, 100-115, doi: 10.1016/j.jclepro. 2016.12.145. [14] Chopra, S., Sodhi, M.S. (2004). Managing risk to avoid supply-chain breakdown, MIT Sloan Management Review, Vol. 46, No. 1, 53-61. [15] Banduka, N., Tadic, D., Macuzic, I., Crnjac, M. (2018). Extended process failure mode and effect analysis (PFMEA) for the automotive industry: The FSQC-PFMEA, Advances in Production Engineering & Management, Vol. 13, No. 2, 206-215, doi: 10.14743/apem2018.2.285. [16] Dubois, D., Prade, H. (1998). An introduction to fuzzy systems, Clinica Chimica Acta, Vol. 270, No. 1, 3-29, doi: 10.1016/S0009-8981(97)00232-5. [17] Zimmermann, H.-J. (2011). Fuzzy set theory - and its applications, Springer Science Business Media, New York, USA. [18] Silva, M.M., de Gusmâo, A.P.H., Poleto, T., e Silva, L.C., Costa, A.P.C.S. (2014). A multidimensional approach to information security risk management using FMEA and fuzzy theory, International Journal of Information Management, Vol. 34, No. 6, 733-740, doi: 10.1016/j.ijinfomgt.2014.07.005. [19] Djapan, M.J., Tadic, D.P., Macuzic, I.D., Dragojovic, P.D. (2015). A new fuzzy model for determining risk level on the workplaces in manufacturing small and medium enterprises, Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability, Vol. 229, No. 5, 456-468, doi: 10.1177/1748006X15581219. [20] Kahraman, C., Sari, I.U. (2012). Multicriteria environmental risk evaluation using type II fuzzy sets, In: Greco, S., Bouchon-Meunier, B., Coletti, G., Fedrizzi, M., Matarazzo, B., Yager, R.R. (eds), Advances in Computational Intelligence. IPMU 2012, Communications in Computer and Information Science, Springer, Berlin, Germany, 449-457, doi: 10.1007/978-3-642-31724-8 47. [21] Liu, H.-C., Li, P., You, J.-X., Chen, Y.-Z. (2015). A novel approach for FMEA: Combination of interval 2-tuple linguistic variables and gray relational analysis, Quality and Reliability Engineering International, Vol. 31, No. 5, 761772, doi: 10.1002/qre.1633. [22] Kiliç, M., Kaya, 1. (2015). Investment project evaluation by a decision making methodology based on type-2 fuzzy sets, Applied Soft Computing, Vol. 27, 399-410, doi: 10.1016/j.asoc.2014.11.028. [23] Chen, S.-M., Lee, L.-W. (2010). Fuzzy multiple attributes group decision-making based on the interval type-2 TOPSIS method, Expert Systems with Applications, Vol. 37, No. 4, 2790-2798, doi: 10.1016/j.eswa.2009.09.012. [24] Ghaemi Nasab, F., Rostami Malkhalifeh, M. (2010). Extension of TOPICS for group DM based on the type two fuzzy positive and negative ideal solutions, International Journal of Industrial Mathematics, Vol. 2, No. 3, 199-213. [25] Temur, G.T., Kaya, T., Kahraman, C. (2014). Facility location selection in reverse logistics using a type-2 fuzzy decision aid method. In: Kahraman, C., Oztayçi, B. (eds), Supply Chain Management Under Fuzziness. Studies in Fuzziness and Soft Computing, Vol. 313, Springer, Berlin, Germany, 591-606, doi: 10.1007/978-3-642-539398 25. [26] Zamri, N., Abdullah, L. (2014). A new qualitative evaluation for an integrated interval type-2 fuzzy TOPSIS and MCGP, In: Herawan, T., Ghazali, R., Deris, M. (eds), Recent Advances on Soft Computing and Data Mining. Advances in Intelligent Systems and Computing, Vol. 287, Springer, Cham, Switzerland, 79-88, doi: 10.1007/978-3-31907692-8 8. [27] Xu, Z.S. (2001). A ranking arithmetic for fuzzy mutual complementary judgment matrices, Journal of Systems Engineering, Vol. 16, No. 4, 311-314. Advances in Production Engineering & Management 14(4) 2019 433 Aleksic, Runic Ristic, Komatina, Tadic [28] Kahraman, C., Oztayçi, B., Sari, I.U., Turanoglu, E. (2014). Fuzzy analytic hierarchy process with interval type-2 fuzzy sets, Knowledge-Based Systems, Vol. 59, 48-57, doi: 10.1016/j.knosys.2014.02.001. [29] Tang, C., Tomlin, B. (2008). The power of flexibility for mitigating supply chain risks, International Journal of Production Economics, Vol. 116, No. 1, 12-27, doi: 10.1016/i.iipe.2008.07.008. [30] Jharkharia, S., Shankar, R. (2007). Selection of logistics service provider: An analytic network process (ANP) approach, Omega, Vol. 35, No. 3, 274-289, doi: 10.1016/i.omega.2005.06.005. [31] Omega, R.S., Noel, V.M., Masbad, J.G., Ocampo, L.A. (2016). Modelling supply risks in interdependent manufacturing systems: A case study, Advances in Production Engineering & Management, Vol. 11, No. 2, 115-125, doi: 10.14743/apem2016.2.214. [32] Simchi-Levi, D., Wei, Y. (2012). Understanding the performance of the long chain and sparse designs in process flexibility, Operations Research, Vol. 60, No. 5, 1125-1141, doi: 10.1287/opre.1120.1081. [33] Sharma, R.K., Kumar, D., Kumar, P. (2005). Systematic failure mode effect analysis (FMEA) using fuzzy linguistic modelling, International Journal of Quality & Reliability Management, Vol. 22, No. 9, 986-1004, doi: 10.1108/ 02656710510625248. [34] Kaur, A., Kanda, A., Deshmukh, S.G. (2008). Supply chain coordination: Perspectives, empirical studies and research directions, International Journal of Production Economics, Vol. 115, No. 2, 316-335, doi: 10.1016/ i.iipe.2008.05.011. [35] Avelar-Sosa, L., García-Alcaraz, J.L., Maldonado-Macías, A.A., Mepa-Muñoz, J.M. (2018). Application of structural equation modelling to analyse the impacts of logistics services on risk perception, agility and customer service level, Advances in Production Engineering & Management, Vol. 13, No. 2, 179-192, doi: 10.14743/apem2018. 2.283. Appendix A List of abbreviations: D Difficulty of failure detection FMEA Failure Mode and Effects Analysis HRM Human resource management IT Information technology IT2FS Interval type-2 fuzzy set IT2TrFN Interval type-2 trapezoidal fuzzy number MCDM Multi-Criteria Decision Making Method O Occurrences of failure realization PCMF Potential Cause(s)/ Mechanism(s) of Failure PEF Potential Effect(s) of Failure PFM Potential failure mode RP Recycling process RSC Reverse supply chain S Severity TOPSIS Technique for Order of Preference by Similarity to Ideal Solution 434 Advances in Production Engineering & Management 14(4) 2019 APEM jowatal Advances in Production Engineering & Management Volume 14 | Number 4 | December 2019 | pp 435-448 https://doi.Org/10.14743/apem2019.4.339 ISSN 1854-625G Journal home: apem-journal.org Original scientific paper A new method for mathematical and simulation modelling interactivity: A case study in flexible job shop scheduling Ojstersek, R.a*, Lalic, D.b, Buchmeister, B.a aUniversity of Maribor, Faculty of Mechanical Engineering, Maribor, Slovenia bUniversity of Novi Sad, Faculty of Technical Sciences, Novi Sad, Serbia A B S T R A C T A R T I C L E I N F O The present study has investigated mathematical and simulation model interactivity for production system scheduling. A mathematical model of a Flexible Job Shop Scheduling Production optimisation problem (FJSSP) was used to evaluate a new evolutionary computation method of multi-objective heuristic Kalman algorithm (MOHKA). Ten Brandimarte and five Kacem benchmarks were applied for evaluation and comparison of MOHKA optimisation results with the Multi-Objective Particle Swarm Optimization algorithm (MOPSO) and Bare-Bones Multi-Objective Particle Swarm Optimization algorithm (BBMOP-SO). Benchmark data sets were divided into three groups, regarding their complexity, from low, middle to high dimensional optimisation problems. The optimisation results of MOHKA show high capability to solve complex multi-objective optimisation problems, especially with real world production systems data. A new robust method is presented of optimisation data interactivity between a mathematical optimisation algorithm and a simulation model. The results show that the presented method can overcome the integrated decision logic of commercial simulation software and transfer the optimisation results into the simulation model. Our interactive method can be used in a variety of production and service companies to ensure an optimised and sustainable cost-time profile. © 2019 CPE, University of Maribor. All rights reserved. Keywords: Flexible job shop scheduling; Mathematical modelling; Simulation modelling; Interactivity; Evolutionary computation; Multi-objective heuristic Kalman algorithm (MOHKA); Multi-objective particle swarm optimization (MOPSO); Bare-bones multi-objective particle swarm optimization algorithm (BBMOPSO) *Corresponding author: robert.ojstersek@um.si (Ojstersek, R.) Article history: Received 15 May 2019 Revised 21 November 2019 Accepted 7 December 2019 1. Introduction In the time of Industry 4.0 [1], which represents highly flexible manufacturing systems, proper scheduling of high-mix low-volume production systems orders is crucial [2]. The high degree of complexity and flexibility of multi-objective optimisation problems increases the problem to achieve optimal scheduling of such production systems significantly [3]. The use of conventional methods [4] and their obtained optimisation results does not achieve the optimally weighted goals of production systems related to appropriate makespan, uniformly high machinery utilisation, financially and timely justified production. For many years, scientists have been using evolutionary computation (EC) methods [5] for the purpose of multi-objective production systems' optimisation [6]. They have been struggling to transfer mathematical models of algorithms into simulation environments to achieve optimised real world production systems indirectly. The complexity of transferring mathematical models into the simulation environments is constrained by the incompatibility of programming environments, integrated decision models within simulation environments, and problems in the reliability of optimisation results' transfers 435 Ojstersek, Lalic, Buchmeister into the simulation environment, in a real world production system. Relevance use and satisfactory optimisation results of EC methods in solving FJSSP are transmitted via the newly proposed algorithms [7, 8] to the real world environment of Industry 4.0 production systems [9, 10]. The complexity of the FJSSP optimisation problem is reflected in the multi-objective nature of the optimisation problem [11]. However, for the transfer and use of new EC methods in a real world environment, the extensive evaluation of algorithms using benchmark and real production datasets is crucial. In this case, the use of simulation modelling is crucial, both in designing the experiment, and in evaluating the optimisation results [12]. The main research problem relates to the compatibility and interactivity between the mathematical model of the EC method optimisation results and the simulation model, which enables the evaluation and application of the EC method in a real world environment [13]. The interactivity between EC methods and the simulation model poses major challenges for researchers, who have not yet studied them thoroughly. In the presented research work, a new method is presented of transferring optimisation results between a mathematical model of MOHKA's EC method [14] and a simulation model in a conventional Simio simulation environment. The presented method provides interactivity between the EC optimisation algorithm and the simulation model, while offering high flexibility of the simulation model and the integration of the EC algorithm optimisation results into the integrated decision logic of the simulation environment [15]. The manuscript is divided into the following sections. Section 2 deals with the mathematical definition of a multi-objective FJSSP optimisation problem. Section 3 presents our own developed method of the multi-objective EC method, called MOHKA. Examples are given of mathematical modelling and experimental testing on five Kacem [16] and ten Brandimarte [17] datasets. Analysis of the MOHKA algorithm optimisation results is performed in comparison with the two existing algorithms MOPSO [18] and BBMOPSO [19]. In this section also the interactivity of the MOHKA mathematical modelling optimisation results and simulation modelling is given. The new proposed block structure of an interactive simulation model method and integrating optimisation results into a simulation model is presented, based on the Kacem and Brandimarte input datasets. Graphically and numerically, the high ability is confirmed of interactivity with optimisation results and the simulation model. Section 3 overviews the new interactive method of using the MOHKA optimisation results in a simulation model, that allows direct use in a real world production system. The high ability to use the proposed method enables further research work on optimally implemented collaborative work places in flexible manufacturing systems. 2. Problem description Production scheduling is defined as decision-making processes that are used on a daily basis in many production and service enterprises [20]. The importance of the decisions taken is, consequently, reflected in the fields of jobs orders, production, transport and distribution of the final products [21]. Production scheduling is the process of optimising, controlling and determination of the limited production system resources (machines, humans, finances etc.). The presented mathematical and simulation modelling method is based on solving an FJSSP multi-objective optimization problem. Given the high degree of difficulty, the following Section presents a notation of abbreviations, a general mathematical description of a multi-objective FJSSP optimization problem, and a mathematical approach how to solve it. 2.1 Notation The notations presented by Graham et al. [22] will be used in the presented paper. i Job (j = 1, ...,n) k Operation (k = 1,..., 0{) n Total number of jobs Oi Number of operations of job ]t Pi,Pij Processing time of job ]t on machine Mj j Machine (J = l,..,m) h Resource (h = 1,..,s) m Total number of machines s Number of limited resources Pik,Pikj Processing time of operation Oik on Mj 436 Advances in Production Engineering & Management 14(4) 2019 A new method for mathematical and simulation modelling interactivity: A case study in flexible job shop scheduling ft Makespan (time required to complete all jobs) /3 Total workload of all machines 2.2 Multi-objective flexible job shop scheduling problem /2 Maximum workload (workload of the most loaded machine) Multi-objective FJSSP is described as: We have n jobs which can be performed on m machines from a set of machines (j = 1,...,m) suitable for carrying out the jobs. The choice of using the machine is made according to the machine occupancy and the suitability of the individual machines to perform the operation. The number of jobs n and number of machines m are given. Each job i has a specific sequence and number of operations 0t. The processing time of the operation Pjk may vary, depending on the machine on which it is performed. For the multi-objective FJSSP some limitations must be made: • One machine can process only one job at a time. • One job can be processed only on one machine ata time. • When the operation starts it cannot be interrupted until the end of the operation, after the completion the next operation can start. • All the jobs and operations have equal priorities at the time zero. • Each machine m is ready at time zero. • Given an operation and the selected machine m, the processing time p^ is fixed. Eqs. 1, 2, and 3 describe three optimisation objectives, as follows: • Makespan (time required to complete all jobs): ft = max {Cj \j = 1,...,n} (1) • Maximum workload (workload of the most loaded machine): /2 = max^^pijkxijk,k = 1,2, ...,m ¿=17=1 Total workload of all machines: n ni m (2) (3) /3 =^^^PijkXijk,k = 1,2,..., m ¿=1j=ik=i where Cj is the completition time of job /¿, and Xijk is a decision variable on which individual machine operation will be processed. In Table 1, we see the FJSSP benchmark example, presented by Kacem et al. [16]. This example has three jobs that must be processed on four machines. The processing time of each operation depends on the machine on which the operation will take place. Fig. 1 shows the optimal solution to the FJSSP optimisation problem presented in Table 1. m5- Mj. My m2- Mi 03,1 02,1 1 2 B 4 5 6 7 8 9 10 t [h] Fig. 1 Gantt chart of optimal solution for a FJSSP optimisation problem Advances in Production Engineering & Management 14(4) 2019 437 Ojstersek, Lalic, Buchmeister Table 1 FJSSP benchmark example Oi,i Mi M2 M3 M4 M5 Ji 0i,i 5 4 2 6 7 02,1 3 5 7 8 4 03,1 5 4 3 2 1 J2 0l,2 9 7 6 7 4 02,2 5 4 7 8 5 03,2 9 1 2 2 4 /3 0l,3 4 7 2 4 8 02,3 5 5 7 5 7 From the optimisation problem defined above, the following limitations can be described: • All machines are available at time t = 0, and any order J can be started at time t = ry. • Only one operation at a time can be performed by a machine. It becomes available to other operations only when the current operation on the machine is completed. • Each operation Oy is the start time ry defined as shown by Eqs. 4 and 5: rU =r/,vl in ri+1J = ru + y^where Yij = min(di,J-,fc) (4) V l 2 represents the number of optimisation parameters, and X represents the feasible set of decision vectors. A set of decision vectors is usually represented by constraint functions. We define the vector-valued objective function as shown in Eq. 7: f:X ^ Rk, f(x) = (/!(x).....fk(x)f. (7) If we want to maximise a function, we can do it by minimising its negative dependence. The element x'eX presents a workable solution or a workable decision. The vector z* =/(x*) eM.k is called the function vector for the feasible solution x*. In multi-objective optimisation, there is usually no viable solution that optimises all of the target functions at the same time. Pareto optimal solutions are solutions that cannot be improved without compromising at least one of the goals of the remaining functions. Using the mathematical notations of Eqs. 8 and 9, we can conclude that the feasible solution x1 e X Pareto is dominated by another solution x2 e X in the case where: /¿(x1) < /¿(x2) for all i e {l, 2, ...,&} and (8) /¿O1) < /¿O2) for at least one j e {l,2, ...,k] (9) is a feasible solution, i*eX and the associated output value /(x*) is Pareto optimal if there is no other solution that dominates it. The Pareto group of optimal solutions is called the Pareto Front. The Pareto Front of multi-objective optimisation problems is limited by two vectors: 438 Advances in Production Engineering & Management 14(4) 2019 A new method for mathematical and simulation modelling interactivity: A case study in flexible job shop scheduling • The nadir vector is defined by Eq. 10: z?ad = sup/i(*) forain = 1,...,fc (10) xex ^ J • The ideal vector is defined by Eq. 11: zideai = inf fi(x) for all i = 1,...,fc (11) The nadir and ideal vector components define the upper and lower bounds for the optimisation functions of Pareto optimal solutions [28]. 3. Results and discussion 3.1 Multi-objective heuristic Kalman algorithm (MOHKA) In order to solve the planning and scheduling problem of FJSSP, researchers use different evolutionary computation methods for solving the multi-objective nature of the problem. The presented research work is based on the use of the Multi-Objective Heuristic Kalman Algorithm (MOHKA), which is based on the mathematical formulation of the Heuristic Kalman Algorithm (HKA), whose operation and use is presented in the cited literature [29]. The positive results of the algorithm's operation and the developed model architecture make it possible to upgrade single-objective optimisation to multi-objective optimisation. MOHKA, the same as HKA, first in the basic interaction, evaluates the solutions using a Gaussian distribution based on the parameters given in the variance and convergence matrix. Based on the definition of the multi-objective Pareto optimal solutions, a limiting function is determined, to ensure that the solutions obtained are within the appropriate limits. The nadir vector zldeal and the ideal vector zldeal, are determined defined by Eqs. 10 and 11. After the evaluation of the suitability of the obtained solutions is completed, the solutions are written to a non-dominant function, which is stored in the data set of non-dominant solutions. The evaluation procedure follows, and the maximum number is selected of the discussed optimisation problem articles . In the final step, an evaluation of the obtained results is performed using a random distribution function. The pseudocode of the algorithm, handling the constraints, updating the non-dominated solution archive, choosing the best samples and pruning the non-dominated solution archive are presented in the literature [14]. Fig. 2 presents a general flow chart ofa MOHKA optimisation algorithm. Fig. 2 Flow chart of MOHKA optimisation algorithm Mathematical modelling of experiments To test the performance of the MOHKA algorithm, we used two of the most well established benchmark data sets for multi-objective optimization of FJSSP. We used five Kacem datasets (Kacem 4x5, Kacem 8x8, Kacem 10x7, Kacem 10x10 and Kacem 15x10) [16] and ten Brandimarte datasets (Mk01 to Mk10) [17]. We divided these benchmark data sets into three groups, according to the complexity of the optimisation problem. Considering the recommenda- Advances in Production Engineering & Management 14(4) 2019 439 Ojstersek, Lalic, Buchmeister tions in the literature [30], we divided the datasets into low, middle and high dimensional optimisation problems. The division of benchmark datasets according to the complexity of the optimisation problems allows a more accurate evaluation of the obtained results, in order to determine the advantages and limitations of the mathematically modelled optimisation algorithm. With the MOHKA algorithm, we optimised three key parameters of a flexible production system: Makespan (MC), total workload of all machines (TW) and maximum workload of an individual machine (MW). The obtained results of the algorithm were compared with the optimal results of the FJSSP problem obtained in the literature [31]. The algorithm was implemented in the MATLAB R2017b software environment, using a PC with an Intel i7 processor with 16 GB of working memory. Mathematical modelling results In order to compare the obtained MOHKA solutions, in addition to the Pareto optimal results, we selected the two currently most advanced algorithms for solving FJSSP optimisation problems: Multi-Objective Particle Swarm Optimization Algorithm (MOPSO) [18] and an improved Bare-Bones Multi-Objective Particle Swarm Optimization (BBMOPSO) algorithm [19]. According to the recommendations [14], the initial parameters of the MOHKA algorithm were set to: N = 300, N = 10, a = 0.3, Na = 100, mr = 0.1, and MaxIter = 3000. The MOPSO and BBMOPSO parameter settings are the same as those in the literature [18, 19]. The optimisation results presented in Tables 3, 4, and 5 represent the optimal solutions of three algorithms, MOHKA, MOPSO, and BBMOPSO, compared to the Pareto optimal solutions. The graphical symbols of the optimisation algorithms on the graphs are presented in Table 2. In order to demonstrate their high ability to solve low, middle and high complex optimisation problems, the following results are presented separately. Numerical and graphical results are given, and comparative comments, advantages and limitations of individual algorithms' comparisons are described. Numerical values MC, TW, MW, presented in Tables 3, 4, 5, are average values of all obtained optimisation results. Table 2 Graphic symbols of the optimisation algorithms Symbol • x + 0 Algorithm MOHKA MOPSO BBMOPSO Pareto optimum • Low dimensional optimisation problems We notice that, with Kacem benchmark datasets, MOHKA solves low-dimensional optimisation problems better. In Table 3, MOHKA achieves Pareto an optimal solution for the Kacem 4x5 dataset, where the MOHKA solutions are located directly in the centre of the Pareto optimal solutions. Compared to MOPSO and BBMOPSO, the Mk01 dataset is also solved better, where the distance between MOHKA optimisation solutions and Pareto optimal solutions are shorter than with the other two algorithms. With the Mk02 dataset, the Pareto optimisation solutions are the closest to the results of the BBMOPSO optimisation algorithm. The comparison between MOHKA and MOPSO shows slightly more optimal solutions given by the MOHKA algorithm. The Kacem 8x8 dataset again demonstrates the equivalence of optimisation results between the MOHKA and MOPSO algorithms. BBMOPSO solved this test data set with an optimally defined solution. BBMOPSO near-optimal optimization results were also presented in the Mk03 dataset, where its solutions are on the Pareto front of optimal solutions. With regard to the MOHKA and MOPSO algorithm solutions, we can claim a similar ability to solve the given data set. In general, we find that the MOHKA algorithm in these low dimensional problems solves the optimization problems as satisfactorily and comparatively as the mentioned comparative algorithm. In two cases, the Kacem 4x5 and Mk01, MOHKA was the most successful, in the others it was comparable to MOPSO and slightly worse than BBMOPSO. 440 Advances in Production Engineering & Management 14(4) 2019 A new method for mathematical and simulation modelling interactivity: A case study in flexible job shop scheduling Table 3 Optimisation results of low dimensional problems Kacem Brandimarte 150 40 140 20 (a) Kacem 4x5 (b) Mk01 (c) Mk02 MC TW MW MC TW MW MC TW MW • 11.5 32 9.5 47.1 161.9 38.4 36.8 157.3 30 + 12 32.3 8.3 52 153.8 43.3 34.5 141 31.5 x 14.7 35.7 9 57 156 48.5 38.3 158 30.8 400 S 200 400 800 200 (d) Kacem 8x8 (e) Mk03 MC TW MW MC TW MW • 25 115.8 19.2 224.3 979.7 218.7 + 16 74 12.5 273 809.2 254.1 x 25.6 101.2 20 263.6 973.1 209.1 • Middle dimensional optimisation problems In the optimisation results of the middle dimensional optimisation problems, presented in Table 4, we find that MOHKA had the largest deviation from the Pareto optimal solutions for the two Kacem datasets in the Pareto graph of solutions Kacem 10x7 and Kacem 10x10. The graphical and numerical results of MOPSO and BBMOPSO proved more optimal optimisation results compared to MOHKA. With the Kacem 10x7 dataset, BBMOPSO achieved the Pareto optimal solution and MOPSO came very close to it. Similarly, the optimisation algorithms solved the Kacem 10x10 dataset, where the difference between MOHKA and MOPSO was smaller, which proves the relevance and comparability of MOHKA with MOPSO. The MOHKA optimisation algorithm solved the Brandimarte middle dimensional optimisation problems better. We can see that, for the two Mk04 and Mk05 datasets, MOHKA solved both data sets most satisfactorily, since the MOHKA algorithm solutions are on the Pareto front of optimal solutions. For the two datasets listed above, BBMOPSO had the most problems, presented on the Mk04 dataset, in which solutions were the furthest from the Pareto optimal solutions. In the case of the Mk05 dataset, MOHKA was comparable to MOPSO. The Mk06 data set was solved relatively well by the MOHKA and BBMOPSO algorithms, with MOPSO farthest from the optimal solution. We concluded that MOHKA performs best in the middle dimensional optimisation algorithms in the Brandimarte datasets, where it dominates the solutions of the other two algorithms. The Kacem datasets were the furthest from the optimal solutions, compared to BBMOPSO and MOPSO. Advances in Production Engineering & Management 14(4) 2019 441 Ojstersek, Lalic, Buchmeister Table 4 Optimisation results of middle dimensional problems Kacem Brandimarte 200 200 300 0 250 TW 670 150 MC (f) Kacem 10x7 (g] Mk04 (h] Mk05 MC TW MW MC TW MW MC TW MW • 16.6 87.4 15.3 85.9 343.4 74.9 189.5 679.3 181.5 + 12 6G.5 11.5 118.1 342.4 1GG.5 2G1.5 676.3 189.9 x 13.5 66 12.5 93.3 35G.1 7G.8 213.5 681.9 189.5 500 200 (i] Kacem 10x10 (j] Mk06 MC TW MW MC TW MW • 16.6 88.4 13 1G4.5 444.9 61.7 + 8.7 42.7 7 123.6 357.1 84.8 x 14.6 63 9.8 123.6 446.8 72.9 • High dimensional optimisation problems MOHKA optimisation algorithm problems in solving middle dimensional Kacem optimisation problems continued with high dimensional optimisation problems, in Table 5, which can be attributed to the construction of Kacem datasets, that are highly susceptible to dropping optimisation algorithms to local minima, to which MOHKA is highly exposed by the KA mathematical structure [14]. Due to the hybridisation and mathematical structure of the algorithm, Kacem 15x10 was best solved by BBMOPSO. We see that all three algorithms were at a relatively short distance from the Pareto optimal solution. Unlike Kacem datasets, MOHKA excelled at Brandimarte datasets, which is more important for solving FJSSP optimisation problems. The structure of the Brandimarte datasets represents real world dataset input of an FJSSP production system. In our case, where we wanted to solve the real world problems of scheduling manufactured systems and establish communication between the optimisation algorithm and the simulation model, this is crucial [32]. The success of solving Brandimarte datasets is paramount in solving FJSSP. We see that MOHKA solved Brandimarte high dimensional optimisation problems Mk08 and Mk09 perfectly for data sets Mk09, and especially Mk08, where we see that MOHKA generated an optimal solution with a high degree of solution delivery reproducibility. Such results demonstrate the robustness and high ability to solve the FJSSP optimisation problem. From presented solutions, we find that MOHKA optimisation results are comparable to the results of MOPSO and BBMOPSO algorithms, especially in Brandimarte datasets, where MOHKA was the most suitable algorithm for solving both medium and high dimensional optimisation problems. MOHKA also demonstrated its competitiveness in low dimensional Kacem datasets. More limitations and deviations from optimal solutions can be detected in medium and high dimensional Kacem datasets. 442 Advances in Production Engineering & Management 14(4) 2019 A new method for mathematical and simulation modelling interactivity: A case study in flexible job shop scheduling Table 5 Optimisation results of high dimensional problems Kacem Brandimarte 300 300 200 TW 600 10(> MC 600 £ 500 400 2550 2500 600 2450 500 700 TW MC (l) Kacem 15x10 (m) Mk07 (n) Mk08 MC TW MW MC TW MW MC TW MW • 24.8 171 21.3 178.2 689.9 169.8 550.1 2506.5 547.5 + 16 91.5 14 189.3 666.2 167.3 597.4 2504 551.9 x 23.4 135.3 19 207 693.7 174.7 599.3 2523.5 545.1 400 200 4000 3000 1000 500 600 400 TW 2000 0 MC TW 1800 200 MC (o) Mk09 (p) MklO MC TW MW MC TW MW • 350.5 2450.9 306.5 358.1 2074.5 226.8 + 534.9 2224.2 383.7 464.2 1888.9 298.8 x 536.1 2448.2 361.2 429 2110.3 258.9 3.2 Simulation modelling In the time of the Industry 4.0, development of simulation tools to optimise production has become increasingly more important. Using advanced simulation methods, we can make manufacturing systems more efficient, financially viable and more competitive in the global market [33]. There are many optimisation problems in production systems that can be solved by the proper use of simulation tools [34]. In our case, we wanted to establish interactivity between the EC optimisation method and the flexible simulation model. The interactive architecture enables communication between the optimisation algorithm and the production system simulation model. Fig. 3 presents the proposed block architecture of a simulation model that enables the interactivity of MOHKA optimisation algorithm results, a simulation model and a real world production system. The interactive architectural model consists of the following phases: • In the first phase, we assign the input data of a real production system or benchmark data sets. The inputs of a real production system must be credible and verifiable. • In the second phase, the implementation of the sequence order optimisation is performed according to the available machines with the MOHKA optimisation algorithm. The algorithm is implemented according to the structure described in Section three. • The third phase transmits the MOHKA optimisation results to the simulation model of a real production system, named the analysed system. Additional optimisation parameters are assigned related to costs, dimensions and setup time. At this stage, we propose the introduction of a new approach for determining the order of operations on an individual machine, which does not depend on the integrated simulation environment decision logic. The approach is presented as follows. Advances in Production Engineering & Management 14(4) 2019 443 Ojstersek, Lalic, Buchmeister • The next phase is the implementation of simulation experiments that allow the calculation of the average multiple iterations' values of the simulation model, and the introduction of simulation scenarios, in order to determine the limitations, changes and proposals for optimisation of the production system. • In the simulation analysis phase, we evaluate the numerical results, and compare the graphical matching of order sequences solutions using MOHKA and the simulation model. Approval of the orders' sequences is confirmed or denied. The importance of data matching ensures the credibility of numerical and simulation results. • The evaluated numerical and simulation results are transferred to a real production system, whereby a change (optimisation) of the production system is enabled, based on an interactive loop. The simulation model was built in the software environment Simio [35], which is a unique software environment for modelling flexible manufacturing or service processes [36]. Simio is based on the use of intelligent objects, and supports both process and object oriented modelling. It is used for discrete and continuous systems. Using the MOHKA algorithm, we solved the FJSSP optimization problem, so we decided to upgrade our existing optimisation results with a suitable simulation model. The following is a new decision logic method that combines the possibilities of testing datasets and optimising real world production systems. The obtained numerical values of the optimisation algorithm, which are given in Table 6, were transferred into the Simio software environment, where the real production system shown in Fig. 4 was modelled. The main advantage of using simulation environments, such as Simio, is the ability to transfer data between the optimisation algorithm, simulation model and real world production system. With the appropriate data transfer method, we can extend the testing and optimisation of production system parameters, in order to extend the numerical model in an optimisation-programming environment that optimises different orders. The simulation model was designed as a flexible modular model built in two parts. The first part of the simulation model is a MOHKA optimisation algorithm that allocates work orders optimally to available machines. In this case, MOHKA optimises three key parameters: Makespan, maximum workload and total workload of all machines. The end result of the MOHKA optimisation algorithm is the optimal allocation of individual work order operations to the available machines. In this case, the order machine execution sequence, their start time, finish time, and machine sequence are obtained. The output results of the MOHKA optimisation algorithm are shown in Table 6. Fig. 3 Block diagram of interactive simulation model 444 Advances in Production Engineering & Management 14(4) 2019 A new method for mathematical and simulation modelling interactivity: A case study in flexible job shop scheduling Fig. 4 Simulation model in Simio Table 6 Optimisation results of MOHKA algorithm transferred to simulation model Order Operation Machine Start time Finish time Machine sequence Ji Oi,i Mi 0 i i 01,2 Mi i 5 2 0i,3 M io 5 s 2 J2 02,1 M io 0 4 i 02,2 M7 4 5 i 02,3 Me 5 s J3 03,1 Ms 0 3 i 03,2 M5 3 5 i 03,3 Ms 5 7 J4 04,i M9 0 i i 04,2 Me i 5 i 04,3 M4 5 6 i J5 05,1 M3 0 4 i 05,2 M3 4 7 2 05,3 M4 7 s 2 The results of the optimisation algorithm shown in Table 6 are transferred automatically to the Simio software environment via the communication interface MATLAB, Excel, Simio [37]. The transmission of optimisation results allows further evaluation of the results using a simulation model of a real world production system. The MOHKA algorithm allocated orders' operations optimally to individual machines, and transferred these data to the Simio simulation environment. When performing simulation experiments in a simulation environment, a limitation occurs in the decision logic of the simulation environment. Simio's simulation environment has an integrated decision logic for order sequencing, which, in a lot of cases, makes it impossible to follow the order sequence given as the solution of the EC optimisation algorithms. In order to eliminate the integrated decision logic of the simulation environment, we introduced a new functional dependency that defines the sequence of execution of individual orders on the assigned machine. In Table 6, the row Machine sequence works by assigning the MOHKA optimisation numerical results of the order sequence and the matching numerical results in Table 6 to the row Machine sequence that defines the sequence of operations according to the MOHKA optimisation solution. Example: According to the MOHKA solution (Fig. 5, top chart), the Oi,i operation is first performed on the Mi machine, followed by the operation Oi,2. The first operation Oi,i to be performed on machine Mi is assigned to this operation consecutive number one, which determines the sequence of execution on the machine. A sequence number two is assigned to operation Oi,2, Advances in Production Engineering & Management 14(4) 2019 445 Ojstersek, Lalic, Buchmeister by the above method. The presented approach defines values for the entire sequence of the Ka-cem 5x10 dataset. The solutions of the approach are shown in Table 6. The smaller the assigned value of the variable Machine sequence is, the sooner it will be executed on the assigned machine. The presented method enables robust, smooth data transfer between the software environments of the optimisation algorithm and the simulation model. With the presented approach, we bypassed the decision logic of the simulation environment, and we could use our own optimisation algorithm to determine the optimum orders' sequence. Fig. 5 shows the Gantt charts' solutions of the optimisation algorithm (top chart) and simulation model (lower chart). The Gantt charts are the same, which proved the high capability of the presented method to integrate the EC method decision logic into the existing simulation environment. The presented approach emphasises the advantage of a modular adaptive design that allows the simulation model to be adapted to a wide range of multi-objective optimisation problems. The simulation model can be adjusted, upgraded or replaced as needed by any part of the two-part structure. Replacing or upgrading individual parts of the presented approach makes the presented method sustainable. M io Mi Ms M7 Ms Ms Mi Ms Mi Mi 0 Mi Ms Mi Me Ms Mi Ms Mi Mi 8 t[h] Crter02 [OfáerOI CMMM OnfcrM 0iw03 M Orí*rCr2 |QW02 I Qfóe-03 ■ OrOw« |oiS«05 OSirOS OrÓcrOS OrwrOl 0rt«01 i -P 1 i i 1 1 0 1234 5678 t[h] Fig. 5 Gantt chart of the optimisation algorithm (top chart) and simulation model results (lower chart) 4. Conclusion The presented research work links our own developed EC method of the MOHKA algorithm with the interactive simulation model. The transfer of optimisation results to a simulation environment can, in many cases, represent limitations in the use of its own decision logic.The presented research results have answered the main research question related to the limitation regarding interactivity of mathematical and simulation modelling. The research problem of a multi-objective FJSSP optimisation problem scheduling was identified initially. Identifying an NP-hard optimisation problem requires the use of advanced EC methods. The proposed MOHKA algorithm shows an example of solving the test Kacem [16] and Brandimarte [17] benchmarks. Fifteen datasets were divided into three difficulty groups. According to the proposal of researchers [25], the division of optimisation problems into three levels of difficulty enables a detailed evaluation of the optimisation results. Optimisation results are evaluated numerically and graphically, and a comparative analysis was performed between MOHKA, MOPSO and BBMOPSO. The advantages and limitations of the individual optimisation results were defined [14], which show 446 Advances in Production Engineering & Management 14(4) 2019 A new method for mathematical and simulation modelling interactivity: A case study in flexible job shop scheduling the characteristics of the individual algorithm with respect to the corresponding mathematical structure. Based on the satisfactory optimisation results, the presented MOHKA method is suitable for optimisation results' interactivity between the mathematical and simulation models. The interactive method shown gives the advantage of transferring optimisation results via a simulation model to a real environment. The importance of transferring the MOHKA algorithm optimisation results to a real world environment demonstrates the high degree of the proposed method's applicability in complex manufacturing systems supported by the Industry 4.0 concept. The graphical results demonstrate the reliability and robustness of the proposed approach, which will be expanded further by modelling and devaluing the flexibility parameter and its dependence on the production cost-time profile. An adequate time justified profile is key in ensuring sustainable production systems. The proposed interactive method represents an advanced, flexible and effective link between mathematical and simulation modelling. The open architectural model allows the extension and application of the method to various optimisation problems for both service and production systems. Further research in the field of Mathematical and Simulation Modelling of flexibility parameters in high-mix low-volume production systems will present the optimum production systems' scheduling importance. Acknowledgement The authors gratefully acknowledge the support of the Slovenian Research Agency (ARRS), Research Core Funding No. P2-0190. References [1] Schwab, K. (2017). The fourth industrial revolution, World Economic Forum, New York, USA. [2] Li, Y., Yao, X., Zhou, J. (2016). Multi-objective optimization of cloud manufacturing service composition with cloud-entropy enhanced genetic algorithm, Strojniški Vestnik - Journal of Mechanical Engineering, Vol. 62, No. 10, 577-590, doi: 10.5545/sv-jme.2016.3545. [3] Chaudhry, I.A., Usman, M. (2017). Integrated process planning and scheduling using genetic algorithms, Tehnički Vjesnik - Technical Gazette, Vol. 24, No. 5, 1401-1409, doi: 10.17559/TV-20151121212910. [4] Huang, J., Suer, G.A. (2015). A dispatching rule-based genetic algorithm for multi-objective job shop scheduling using fuzzy satisfaction levels, Computers & Industrial Engineering, Vol. 86, 29-42, doi: 10.1016/j.cie.2014. 12.001. [5] Mitchell, M. (1998). An introduction to genetic algorithms, MIT press, Cambridge, England. [6] Hao, X., Gen, M., Lin, L., Suer, G.A. (2017). Effective multiobjective EDA for bi-criteria stochastic job-shop scheduling problem, Journal of Intelligent Manufacturing, Vol. 28, No. 3, 833-845, doi: 10.1007/s10845-014-1026-0. [7] Meolic, R., Brezočnik, Z. (2018). Flexible job shop scheduling using zero-suppressed binary decision diagrams, Advances in Production Engineering & Management, Vol. 13, No. 4, 373-388, doi: 10.14743/apem2018.4.297. [8] Wang, H., Wang, W., Sun, H., Cui, Z., Rahnamayan, S., Zeng, S. (2017). A new cuckoo search algorithm with hybrid strategies for flow shop scheduling problems, Soft Computing, Vol. 21, No. 15, 4297-4307, doi: 10.1007/s00500-016-2062-9. [9] Vujica Herzog, N., Buchmeister, B., Beharic, A., Gajsek, B. (2018). Visual and optometric issues with smart glasses in Industry 4.0 working environment, Advances in Production Engineering & Management, Vol. 13, No. 4, 417428, doi: 10.14743/apem2018.4.300. [10] Zhang, H., Liu, S., Moraca, S., Ojstersek, R. (2017). An effective use of hybrid metaheuristics algorithm for job shop scheduling problem, International Journal of Simulation Modelling, Vol. 16, No. 4, 644-657, doi: 10.2507/ IISIMM16(4)7.400. [11] Fu, H.C., Liu, P. (2019). A multi-objective optimization model based on non-dominated sorting genetic algorithm, International Journal of Simulation Modelling, Vol. 18, No. 3, 510-520, doi: 10.2507/IISIMM18(3)C012. [12] Gotlih, J., Brezocnik, M., Balic, J., Karner, T., Razborsek, B., Gotlih, K. (2017). Determination of accuracy contour and optimization of workpiece positioning for robot milling, Advances in Production Engineering & Management, Vol. 12, No. 3, 233-244, doi: 10.14743/apem2017.3.254. [13] Pantic, M., Dordevic, A., Eric, M., Mitrovic, S., Babic, M., Džunic, D., Stefanovic, M. (2018). Application of artificial neural network in biotribological research of dental glass ceramic, Tribology in Industry, Vol. 40, No. 4, 692-701, doi: 10.24874/ti.2018.40.04.15. [14] Ojstersek, R., Zhang, H., Liu, S., Buchmeister, B. (2018). Improved heuristic kalman algorithm for solving multi-objective flexible job shop scheduling problem, Procedia Manufacturing, Vol. 17, 895-902, doi: 10.1016/j.promfg. 2018.10.142. Advances in Production Engineering & Management 14(4) 2019 447 Ojstersek, Lalic, Buchmeister [15] Rupnik, B., Nardin, R., Kramberger, T. (2019). Discrete event simulation of hospital sterilization logistics, Tehnički Vjesnik - Technical Gazette, Vol. 26, No. 5, 1486-1491, doi: 10.17559/TV-20180614102011. [16] Kacem, I., Hammadi, S., Borne, P. (2002). Pareto-optimality approach for flexible job-shop scheduling problems: Hybridization of evolutionary algorithms and fuzzy logic, Mathematics and Computers in Simulation, Vol. 60, No. 3-5, 245-276, doi: 10.1016/S0378-4754(02)00019-8. [17] Mastrolilli, M., Gambardella, L.M. (2000). Effective neighbourhood functions for the flexible job shop problem, Journal of Scheduling, Vol. 3, No. 1, 3-20, doi: 10.1002/(SICI)1099-1425(200001/02)3:1<3::AID-l0S32>3.0.C0:2-Y. [18] Mostaghim, S., Teich, J. (2003). Strategies for finding good local guides in multi-objective particle swarm optimization (MOPSO), In: Proceedings of the 2003 IEEE Swarm Intelligence Symposium, Indianapolis, USA, 26-33, doi: 10.1109/SIS.2003.1202243. [19] Zhang, Y., Gong, D.-W., Ding, Z. (2012). A bare-bones multi-objective particle swarm optimization algorithm for environmental/economic dispatch, Information Sciences, Vol. 192, 213-227, doi: 10.1016/j.ins.2011.06.004. [20] Cajzek, R., Klanšek, U. (2016). Mixed-integer nonlinear programming based optimal time scheduling of construction projects under nonconvex costs, Tehnički Vjesnik - Technical Gazette, Vol. 23, No. 1, 9-18, doi: 10.17559/TV-20140108112928. [21] Becker, C., Scholl, A. (2009). Balancing assembly lines with variable parallel workplaces: Problem definition and effective solution procedure, European Journal of Operational Research, Vol. 199, No. 2, 359-374, doi: 10.1016/ j.ejor.2008.11.051. [22] Graham, R.L., Lawler, E.L., Lenstra, J.K., Kan Rinnooy, A.H.G. (1979). Optimization and approximation in deterministic sequencing and scheduling: A survey, Annals of Discrete Mathematics, Vol. 5, 287-326, doi: 10.1016/S0167-5060(08)70356-X. [23] Lin, L., Gen, M. (2018). Hybrid evolutionary optimisation with learning for production scheduling: State-of-the-art survey on algorithms and applications, International Journal of Production Research, Vol. 56, No. 1-2, 193223, doi: 10.1080/00207543.2018.1437288. [24] Miettinen, K.M. (2012). Nonlinear multiobjective optimization, Springer, New York, USA. [25] Deb, K., Agrawal, S., Pratap, A., Meyarivan, T. (2000). A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II, In: Schoenauer, M., Deb, K., Rudolph, G., Yao, X., Lutton, E., Merelo, J.J., Schwefel, H.-P. (eds), Parallel Problem Solving from Nature PPSN VI, Lecture Notes in Computer Science, Vol. 1917, Springer, Berlin, Germany, 849-858, doi: 10.1007/3-540-45356-3 83. [26] Deb, K., Jain, H. (2014). An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: Solving problems with box constraints, IEEE Transactions on Evolutionary Computation, Vol. 18, No. 4, 577-601, doi: 10.1109/TEVC.2013.2281535. [27] Hwang, C.-L., Masud, A.S.M. (1979). Multiple objective decision making - Methods and applications: A state-of-the-art survey, Springer-Verlag, Berlin, Heidelberg, New York, doi: 10.1007/978-3-642-45511-7. [28] Hassanzadeh, H.R., Rouhani, M. (2010). A multi-objective gravitational search algorithm, In: 2nd International Conference on Computational Intelligence, Liverpool, United Kingdom, 7-12, doi: 10.1109/CICSyN.2010.32. [29] Ojstersek, R., Zhang, H., Palcic, I., Buchmeister, B. (2017). Use of heuristic Kalman algorithm for JSSP, In: Andraš. A. (ed.), XVII International scientific conference on industrial systems, Novi Sad, Faculty of Technical Sciences, Department for Industrial Engineering and Management, Novi Sad, Serbia, 72-77. [30] Demir, Y., Kur§at I§leyen, S. (2013). Evaluation of mathematical models for flexible job-shop scheduling problems, Applied Mathematical Modelling, Vol. 37, No. 3, 977-988, doi: 10.1016/j.apm.2012.03.020. [31] Chiang, T.-C., Lin, H.-J. (2013). A simple and effective evolutionary algorithm for multiobjective flexible job shop scheduling, International Journal of Production Economics, Vol. 141, No. 1, 87-98, doi: 10.1016/j.ijpe.2012. 03.034. [32] Prester, J., Buchmeister, B., Palčič, I. (2018). Effects of advanced manufacturing technologies on manufacturing company performance, Strojniški Vestnik - Journal of Mechanical Engineering, Vol. 64, No. 12, 763-771, doi: 10.5545/sv-jme.2018.5476. [33] Lu, Y. (2017). Industry 4.0: A survey on technologies, applications and open research issues, Journal of Industrial Information Integration, Vol. 6, 1-10, doi: 10.1016/j.jii.2017.04.005. [34] Ojstersek, R., Buchmeister, B. (2017). Use of simulation software environments for the purpose of production optimization, In: Proceedings of the 28th DAAAM International Symposium on Intelligent Manufacturing and Automation, Zadar, Croatia, 750-758, doi: 10.2507/28th.daaam.proceedings.106. [35] Joines, J.A., Roberts, S.D. (2015). Simulation modeling with SIMIO: A workbook, 4th Edition, CreateSpace Independent Publishing Platform, Sewickley, USA. [36] Yang, W., Takakuwa, S. (2017). Simulation-based dynamic shop floor scheduling for a flexible manufacturing system in the industry 4.0 environment, In: Proceedings of the 2017 Winter Simulation Conference, Las Vegas, USA, 3908-3916, doi: 10.1109/WSC.2017.8248101. [37] Dehghanimohammadabadi, M., Keyser, T.K. (2017). Intelligent simulation: Integration of SIMIO and MATLAB to deploy decision support systems to simulation environment, Simulation Modelling Practice and Theory, Vol. 71, 45-60, doi: 10.1016/j.simpat.2016.08.007. 448 Advances in Production Engineering & Management 14(4) 2019 Advances in Production Engineering & Management Volume 14 | Number 4 | December 2019 | pp 449-460 https://doi.Org/10.14743/apem2019.4.340 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper Surface features of chromium alloyed carbon steel specimens after salt-spray tests in NaCl solution Varga, G.a, Torok, T.b, Felho, C.a*, Orosz-Szirmai, G.b, Rez, I.c institute of Manufacturing Science, University of Miskolc, Miskolc, Hungary institute of Metallurgy, University of Miskolc, Miskolc, Hungary cStarters E-Components Generators Automotive Hungary Ltd, Miskolc, Hungary A B S T R A C T A R T I C L E I N F O The most common corrosion testing procedures use flat test coupons fabricated from a given steel material. However, workpieces that undergo machining and finishing, especially those with complex geometry, may be more subject to surface degradation on their curved surfaces. In the long run, this may adversely affect the smooth operation of the tested component. This study investigates surface features of machined and finished chromium alloyed steel specimens with rather complex geometries. Changes in several important surface features (topography, roughness, cylindricity, chemical degradation rate and corrosion products) were all measured periodically (after 24, 96 and 192 hours) on several sets of manufactured carbon steel planetary axle specimens exposed to a chemically aggressive medium (aqueous NaCl spray) at 35 °C. In the neutral salt spray (NSS) testing cabinet the cylindrical parts of the chromium alloyed carbon steel shafts showed quite severe and uneven chemical degradation, with the formation of several iron oxide-hydroxide products (rust) observed together with some chromium compounds. After careful removal of the relatively loose corrosion products, the exposed bare shaft surface's geometrical changes showed steady (close to linear) increase in the roughness values throughout the duration of the corrosion tests. This chemical attack caused significant changes in the surface topography as well. It was found that the average values of roughness parameters after the 24hour test were about two and a half times higher than the original values, while they increased by four-fold in the 96-hour test and by approximately eightfold in the 192-hour test. Furthermore, it was found that the values of the 3D roughness parameters (Sa, Sz, Sq) are on average twice that of their 2D counterparts (Ra, Rz, Rq) on the corroded surfaces. Circularity and roundness error data showed a similar increase with salt-spray test time duration. © 2019 CPE, University of Maribor. All rights reserved. Keywords: Surface features; Surface topography; Roundness error; Cylindricity deviation; Corrosion; Surface roughness; Carbon steel; Chromium alloyed steel; Salt-spray test; NaCl solution *Corresponding author: csaba.felho@uni-miskolc.hu (Felho, C.) Article history: Received 15 February 2019 Revised 23 September 2019 Accepted 25 September 2019 1. Introduction Many rotating machine parts like steel shafts (e.g. planetary axles) built into different motor vehicles operating with narrow gaps between the shafts and their counterparts. As such machine elements should run for years without any major surface degradation (like wear and/or corrosion) causing the risk of operational malfunction, they should be properly designed, built, assembled and operated, paying attention to the initial surface finishing as well as to the ever-changing actual surface state of the given steel rods (like rotating shafts) with elapsing time of operation. In this respect, the chemical impact of the environment can also be important, as car- 449 Varga, Török, Felho, Orosz-Szirmai, Réz bon steels are prone to corrosion and surface corrosion products may also adversely influence the operation of such steel machine elements. In the literature, the relationship between corrosion and surface roughness is bidirectional. One part of the studies deals with the effect of the initial surface roughness of the part exposed to corrosion on the behaviour of corrosion and investigates the corrosion process [1-9]. The effect of initial surface roughness on corrosion has been investigated on workpieces of various material grades, such as AISI 304 stainless steel [2], X70 high strength steel [3], IN718 superal-loy [4], magnesium [5], magnesium alloy [6], aluminium alloy 6061 [7] and alloy 690TT, with large wt% content of nickel and chromium [8]. The other group of studies analyses how corrosion exacerbates the surface roughness of the test component, such as in [9]. The relationship between the surface roughness and corrosion resistance of ferrite stainless steel 21Cr was identified by Lee et al. in [10]. Kostadin et al. [11] have investigated the role of chilled air-cooling in corrosion resistance during turning of mar-tensitic stainless steel X20Cr13. Toloei et al. [12] showed that the smoother the surface is, the more resistive it is against corrosion, as the lower roughness acts as a better barrier to penetration of the aggressive electrolyte into the metal substrate. The surface roughness of steel pipes exposed to the corrosive effects of the marine environment was investigated by Gathimba et al. [13]. They also concluded that a relationship exists between surface roughness and corrosion degradation, and only surface roughness height parameters are correlated strongly with corrosion degradation. This research project was designed to follow the changing surface condition of a given steel rod specimen (part of an electric motor) after keeping it in a salt-spray cabinet for different periods of time. Afterwards both the surface topography and the rate of corrosion were determined. The surface topography features investigated in the research were the surface roughness, circularity and cylindricity errors. The surface corrosion products were also analysed in order to be able to estimate their probable harmful effect during operation of such electrical motor elements in chemically highly aggressive/corrosive environments. 2. Materials and methods The steel rod specimens were received from an automotive part manufacturer producing also electric motors, whose unused shafts were removed and sent to the testing laboratories. The chemical elementary composition of the specimens was analysed by a GD OES spectrometer type Profiler 2 using certified reference steel samples/etalons, and the measured composition is given in Table 1. As can be seen from Table 1, this type of steel also contains chromium in a relatively high percentage (> 1 %) in addition to the other common components of such a widely used C45 type round bar steel grade. A simplified drawing of the tested part of the given planetary axle is shown in Figs. 1 and 2. The surface finishing of the manufactured rods was made by grinding to roughness values of Rz4 and Rz6.3 (Fig. 2). Table 1 Elementary composition of the steel rod specimens C Mn Si Ni Cr Cu Mo Ti V P S Al wt% wt% wt% wt% wt% wt% wt% wt% wt% wt% wt% wt% 0.411 0.764 0.175 0.021 1.025 0.005 0.010 0.009 0.017 0.005 0.012 0.021 LO S 047 ' ' | s . 014.6 . - - I40 o 5 ---- --------------------- 13 50 10 fi . 2 fi 5.7 104.5 ,10.5 450 Fig. 1 Simplified drawing of the investigated planetary axle Advances in Production Engineering & Management 14(4) 2019 Surface features of chromium alloyed carbon steel specimens after salt-spray tests in NaCl solution 2.1 Specimen preparation For determining the initial surface topography of the steel rod specimens, their original surface conditions were maintained as received (they were kept in a desiccator) and only a very light physical cleaning (de-dusting) preparation was applied before commencing the surface roughness measurements. Afterwards, the specimens were transferred to the corrosion testing laboratory, where all of the preparation steps were administered according to the ISO 9227 standard for neutral salt spray testing procedures. For the selective removal of the corrosion products, the ISO 9227 recommendation was applied, i.e. the oxides-hydroxides (rust) of the corroded specimen surface were dissolved in an inhibited (with hexamethylenetetramine) HCl solution, playing special attention to avoid over-pickling, i.e. partial dissolution of the bare steel after chemical solubilisation of the corrosion products. Before measuring the topography data, the specimens were also thoroughly cleaned and then kept in the desiccator to avoid any re-oxidation and/or surface post-contamination. 2.2 Surface characterisation: Surface roughness Surfaces of the samples were analysed by an Altisurf© 520 three-dimensional surface topography measuring machine. This equipment has three different measuring heads: an inductive gauge, an LH-G32 laser sensor and a CL2 confocal chromatic sensor (CCS) with an MG140 magnifier. The most accurate measuring head, the CSS sensor was used in the current investigations, as it can measure surface roughness of the order of 5nm [14]. Three measurement positions were defined on each section investigated (A-C, Fig. 2), and the measurements were performed in three angular positions at 120° from each other. The three angular positions were marked on the face of the rod. Thus, 18 profiles and 18 surfaces were recorded for each surface. Surface 1 has a roughness specification of Rz4, while Surface 2 has a more permissive tolerance of Rz6.3. The investigated surface roughness parameters were mainly the prescribed maximum height of the roughness profile Rz and its 3D counterpart Sz, and some common roughness parameters as Ra, Rq and Sa, Sq in 3D. The 3D roughness parameters were included in order to be able to characterize the surface more precisely. The The surface roughness measurement and evaluation parameters were selected according to ISO 4288:1996 and ISO 4287:1997 for profiles and ISO 25178-2:2012 for 3D measurements. Thus, the evaluation length was 4 mm for profiles and the cut-off length was Ac = 0.8 mm. Rz4 V ■-íf-É E E Rz6.3/ ■ » is Surface 1 Surface 2 Fig. 2 Surface roughness measurement positions on the rod parts of the tested planetary axle 2.3 Surface characterisation: Circularity and cylindricity errors Circular and cylindricity error measurements were made on the Talyrond 365 Circular and Shape Error Measuring Equipment. During the measurements we determined the range of cylindricity and circularity error characteristics to be tested. For cylindricity deviations, these are: CYLt (peak to valley cylindricity deviation), CYLp (peak to reference cylindricity deviation), and CYLv (reference to valley cylindricity deviation) [15]. In the case of circularity errors, RONt (roundness total), RONp (roundness peak), and RONv (roundness valley) values were deter- Advances in Production Engineering & Management 14(4) 2019 451 Varga, Török, Felho, Orosz-Szirmai, Réz mined. On the Surface 1 cylindrical surface, both cylindricity and circularity measurements were carried out along circles perpendicular to the axis of the cylinder corresponding to the A, B and C points selected. On the Surface 2 cylindrical surface the cylindricity error was measured at 5 places, of which only three (A, B, and C) places are shown here (Fig. 2). To evaluate the measurement data, Taylor Hobson [iltra software for the Taylor Hobson's Talyrond 365 gauging device was used. 2.4 Corrosion testing in neutral salt spray In a neutral salt spray (NSS) testing chamber the specimens are exposed to a highly corrosive environment where tiny and fresh solution droplets (salty fog) continuously arrive to and wet the solid surface, maintaining a more or less coherent aqueous NaCl solution film containing dissolved oxygen as well. In such an artificial environment, which is somewhat similar to natural windy situations with precipitating salty water droplets close to the sea, carbon steel specimens start corroding according to the well-known electrochemical reaction mechanism that is thoroughly described, for example, in a recent review article of Alcántara et al. [16]. Namely, in the aqueous salt solution containing dissolved oxygen, there will be migrating tiny local anode and cathode sites with iron dissolution (oxidation to iron ions) and oxygen reduction, respectively. Then the iron(II) and iron(III) cations can react with the anions (first of all OH- anions) present nearby with the formation of soluble and non-soluble oxide-hydroxides and other corrosion products depending on the actual local electrochemical reactions and physical transport circumstances (diffusion, solution flows, etc.). In the NSS chamber operated at the Corrosion Testing Centre of Starters E-Components Generators Automotive Hungary Ltd., the specimens were placed parallel to each other (Fig. 3), and the operational parameters listed in Table 2 were used. Fig. 3 Specimen (steel planetary axles) positions in the NSS chamber equipped with air flow operating aqueous salt spray nozzles Table 2 Operational parameters of the NSS testing cycles and subsequent surface studies Test cycle times NSS solution contact Surface examination Topography measurements 24 h 35 °C, continuous After rust removal Roughness and cylindricity 96 h 35 °C, continuous After rust removal Roughness and cylindricity 192 h 35 °C, continuous After rust removal Roughness and cylindricity 2.5 Corrosion products testing phase analysis As is often observed during NSS testing of such grade of steels, the specimens started rusting quite fast and after 96 hours exposure time much of the loosely and unevenly formed corrosion products could be relatively easily removed from the cylindrical surface of the corroded steel axles. The rust scrapings were homogenized, and the phase composition was determined by the 452 Advances in Production Engineering & Management 14(4) 2019 Surface features of chromium alloyed carbon steel specimens after salt-spray tests in NaCl solution X-ray powder diffractometer Bruker D8 Advance operated with the following parameters: Co Ka radiation, 40 kV cathode ray tube voltage and 40 mA current. The measured diffraction peaks were then identified, and the corresponding mineralogical phases are also shown on the XRD spectrum (Fig. 11). 3. Results and discussion 3.1 Surface topography changes in roughness Representative segments are indicated on the cylindrical surface of the 96-hour sample where the change in surface roughness can be traced on both on the initial and the corrosion tested surfaces. Surface 1 of Sample No. 3 at measurement position B and angular position 1 can be seen on Fig. 4. It can be observed that the respective roughness parameters show an approximately four-fold increase after the testing. The investigated roughness values for different salt spray test durations are summarized in Table 3, and the values are visualized in Figs. 5 and 6. The roughness values increase remarkably as the test duration time becomes longer (Fig. 6). Several roughness parameters were evaluated during the surface roughness investigations of the corroded and then handled (i.e. cleaned and rust-freed) surfaces, namely the arithmetical means roughness Ra, the maximal roughness Rz, the root mean square roughness Rq and their 3D counterparts. It can be concluded from the measurement results that all 3D roughness values are always higher. This is because the corroded surface itself is very inhomogeneous, as can be seen in Fig. 4. Here, Fig 4a shows the initial surface (before the corrosion test), which has a typical ground surface topography: the surface has a regular pattern, which comes from the grinding wheel structure and from the applied processing parameters. However, after the corrosion tests and rust removal procedure, the character of the surface significantly changes: the regular pattern has transformed into the irregular surface shown in Fig. 4b because of the previously mentioned surface inhomogeneity. However, the Rz values and their 3D counterpart Sz are of the greatest interest during this research, as these values were defined directly on the drawing as surface quality requirements. Therefore, only the Rz and Sz values are included in Figs. 5 and 6, but the nature of the change in roughness values with the salt spray test cycle time is very similar for the other roughness characteristics as well. Based on the data it can be stated that the values of the roughness parameters are on average two and a half times higher (2.5) after the 24-hour test compared to their original values; they increased approximately four-fold (4.3) after 96 h and more than eight-fold (8.7) after 192 hours of salt spray testing. However, it is observed that the values of the 3D roughness parameters (Sa, Sz, Sq) are on average more than twice as high as their 2D equivalent (Ra, Rz, Rq) on the corroded surfaces. One possible reason for this can be that the corrosion appears crater-like at certain points on the surface and therefore has a high surface inhomogeneity, which is better detected by 3D roughness characteristics. Table 3 Roughness values of the two surfaces for different test durations Surface 1 Surface 2 Initial 24 h 96 h 192 h Initial 24 h 96 h 192 h Ra 0.587 1.019 1.418 3.658 0.631 1.424 2.167 2.930 Rz 5.845 8.893 11.387 23.780 6.053 10.980 16.223 22.280 Rq 0.755 1.333 1.807 4.714 0.808 1.777 2.767 3.938 Sa 0.648 1.574 3.690 9.998 0.666 2.403 4.317 7.348 Sz 6.739 18.467 34.633 74.620 6.773 26.500 42.867 71.840 Sq 0.831 2.172 4.843 12.886 0.851 3.243 5.600 9.378 Advances in Production Engineering & Management 14(4) 2019 453 Varga, Török, Felho, Orosz-Szirmai, Réz |jm 4 2 0 -2 -4 Ra 0.602 |im Rz 5.97 |im Rq 0.777 |im 4 mm [Jin 10 -10 -20 Ra Rz Rq 2.29 |im 14.4 |im 2.78 |im 4 mm X = 2 mm Y = 2 mm Z = 6.78 Mm X = 2 mm Y = 2 mm Z = 3S.9 pm (a) before (b) after 96 h salt spray test and rust removal Fig. 4 2D and 3D roughness images and values for Sample No. 3, Surface 1, Measurement position B, Angular position 1 80 70 60 I 50 X 40 £ 30 20 10 .......I.I.I Before After Before After Before After Before After Before After Before After Surface 1 Surface 2 Surface 1 Surface 2 Surface 1 Surface 2 Sample 1 Sample 3 Sample 5 Rz Sz Fig. 5 Changes in Ra, Rq and Sa, Sq roughness parameters for the different samples and surfaces before and after NSS test 80 70 60 - 50 40 g 30 20 10 0 I 24 96 Saltspray testcycletime, h Fig. 6 Changes in Rz and Sz roughness parameters with salt spray test cycle time 454 Advances in Production Engineering & Management 14(4) 2019 Surface features of chromium alloyed carbon steel specimens after salt-spray tests in NaCl solution The quantification of the deviation change is facilitated by the introduction of a deterioration factor: where: P Parameter (e.g. Rz, Sz, CYLp, CYLv, CYLt, RONp, RONv, RONt) t Time (e.g. 24 h, 96 h, 192 h) i Refers to the initial state Pi means the value of the parameter at the initial (before the salt spray test) state. Table 4 Values of deterioration of Rz and Sz roughness parameter values with duration of salt spray test DRz,t, % Dsz,t, % ti = 24 h 79.324 286.974 t2 = 96 h 213.947 525.978 t3 = 192 h 263.874 949.065 3.2 Changing of circularity errors and cylindricity deviations Fig. 7 displays the change in circularity error and the cylindricity deviations for the initial and the tested surface in the cylindrical surface area of the 96-hour sample for the representative segments shown in Subsection 3.1. The data from Table 5 shows that the longer the duration of the salt test was, the greater the value of the circularity error and the cylindricity deviation after the rust removal. This is observable on Fig. 8 as well, where the results after each test intervals are compared. Fig. 8 shows changes of RONt and CYLt parameters for the 96-hour salt spray test time for Surface 1 (A) and Surface 2 (C) for the three specimens. Similar results can be found on the other surfaces of the examined samples. It is apparent from Fig. 8 that the measured CYLt peak to valley cylindricity deviation and the measured RONt roundness total error show very similar values, both for Surface 1 and for Surface 2. In addition, for the 24-hour and 192-hour salt spray test the deviations for Surface 1 were the greater, so further investigations are detailed only for Surface 1. Figs. 9 and 10 illustrate the relationship between the three cylindricity deviation parameters for Surface 1, as well as the three circularity error parameters and the salt spray test cycle time. It is clear that with the increase in salt spray test time, the nature of increase in both the cylindricity deviation and the circularity error is almost identical. Applying Equation (1), the deterioration factor can be calculated for CYLp, CYLt, RONp, RONt as well (Table 6). Table 5 Cylindricity deviation and circularity error values of the two surfaces for different test durations Sample 1 Sample 3 Sample 5 Surface 1 Surface 2 Surface 1 Surface 2 Surface 1 Surface 2 Before After Before After Before After Before After Before After Before After (Initial) 24 h (Initial) 24 h (Initial) 96 h (Initial) 96 h (Initial) 192 h (Initial) 192 h RONt 8.40 18.37 3.91 15.18 3.37 44.63 3.63 48.31 5.93 70.46 4.70 69.94 RONp 6.97 5.44 2.43 4.36 1.92 15.13 1.81 16.91 4.32 24.46 2.85 27.02 RONv 1.44 12.97 1.48 10.82 1.46 29.51 1.82 31.60 1.61 38.72 1.85 42.91 CYLt 4.65 33.11 6.28 20.57 4.13 55.64 5.23 64.02 8.60 90.33 31.47 76.57 CYLp 3.13 12.39 4.02 4.80 2.41 17.12 2.64 18.85 6.44 33.54 28.76 25.17 CYLv 1.52 20.72 2.26 15.77 1.72 38.56 2.60 42.17 2.16 56.79 2.71 51.39 Advances in Production Engineering & Management 14(4) 2019 455 Varga, Török, Felho, Orosz-Szirmai, Réz Before (initial state) RONt 3.06 |im RONp 1.69 |im RONv 1.37 |im After 96 h salt spray test and rust removal RONt 49.91 |im RONp 14.14 |im RONv 35.77 |im CYLt 4.13 |im CYLp 2.41 |im CYLv 1.72 |im CYLt 55.64 |im CYLp 17.12 |im CYLv 38.56 |im Runout (|jm) . 3.151— ■ Z Hägfit (mm) 147.003 :--1 136.998 C.11^ Eccentricity Fig. 7 Circularity and cylindricity deviation images and values for Sample No. 3, Surface 1 (Measurement position B for the circularity figures) ■ J L. 1 . ll _ - i - j 100 90 80 E 70 a 60 £ 50 £ 40 uj 30 20 10 0 Before After Before After Before After Before After Before After Before After Surface 1 Surface 2 Surface 1 Surface 2 Surface 1 Surface 2 Sample 1 Sampled Sample 5 ■ RONt ■ CYLt Fig. 8 Change in RONt roundness total error and CYLt peak to valley cylindricity deviation parameters for the 96-hour salt spray cycle time for various samples 456 Advances in Production Engineering & Management 14(4) 2019 Surface features of chromium alloyed carbon steel specimens after salt-spray tests in NaCl solution Table 6 Values of deterioration of cylindricity deviation (CYLp and CYLt) and circularity error (RONp and RONt) values with duration of salt spray test DcYLp,t, % DcYLt,t, % DRONp,t, % DronU, % ti = 24 h 210.53 471.85 23.64 211.36 t2 = 96 h 329.07 860.97 178.13 656.44 t3 = 192 h 740.60 1460.10 455.91 1094.24 100 Initial 24 h 96 h 192 h Salt spray test cycle time CYLp CYLv CYLt Fig. 9 Relation between CYLp, CYLv and CYLt cylindrical errors and salt spray cycle time, Surface 1 80 Initial 24 h 96 h 192 h Salt spray test cycle time RONp RONv RONt Fig. 10 Relation between RONp, RONv and RONt circular error parameters and salt spray cycle time, Surface 1 We found that the value of the deterioration factor for CYLt becomes approximately 3.1 times higher when the salt spray test time was 192 hours related to the value which resulted at 24 hours. This ratio is approximately 5.2 when examining the RONt circularity error parameter. 3.3 Surface corrosion products (Fe-Cr-O-OH phases) Smooth and reliable operation of any machine part requires special attention also to the risk of chemical surface degradation, and hence deformation due to corrosion. The changing rate of corrosion with exposure time and the many different forms of corrosion products - even in the "simple" case of iron base alloys - increases the importance of testing actual examples closely related to industrial applications of steels. The tested type of common carbon steel was additionally alloyed with chromium, an element whose chemical affinity to oxygen is much higher than that of iron. For the simplest cases, i.e. after calculating [17] the relevant standard formation Gibbs energies (AG8) of the two metal trioxides at 25 °C, it was found that AGe(Cr2Os) « -1053 kJ/mol and AG8(Fe2O3) « -741 kJ/mol, so it is not surprising that oxide components containing chromium were also found amongst the corrosion products. While studying the influence of chromium on the flow-accelerated corrosion behaviour of four different low alloy steels (C « 0.08 %, and Cr « 0, 0.5, 2, and 5 %, respectively) in 3.5 % NaCl solution, Jiang et al. [18] observed no significant difference in weight loss between the carbon steel and Cr-containing steels in stagnant solution; however, with flowing (0...2.4 m/s) solutions the carbon steel exhibited the highest weight loss. Among the corrosion products they could identify the following major phases: FeOOH/Fe(OH)3 , y-FeOOH, a-FeOOH and FeCr2O4 by Raman spectroscopy, while XPS analysis revealed the presence of other chromium(III) oxides and hydroxides (Cr(OH)3, CrOOH) in Advances in Production Engineering & Management 14(4) 2019 457 Varga, Török, Felho, Orosz-Szirmai, Réz addition to the iron chromite (FeCr2O4). They also mention some kind of chromium enrichment of the corrosion products, which is in line with our EDS results observed at the more extensively rusted spots of the tested Cr-containing carbon steel axle. At these spots a stronger adhesion of the Cr-containing corrosion products can also be assumed. In the salt spray environment (i.e. in the NSS chamber) quite many insoluble corrosion products containing iron were formed together with some encapsulated water soluble components (like NaCl), which could be detected and identified by the applied x-ray diffraction technique (Fig. 11). 1 1 Wt % At % ■ l B I Wt % At % ■ 2 1 Wt % At % O 1.44 4.84 O 25.03 53.54 O 14.75 36.88 Cl 0.34 0.52 Cl 0.28 0.27 Cl 4.43 5.00 Cr 1.36 1.40 Cr 1.69 1.45 Cr 3.50 2.70 Mn 1.23 1.20 Mn 0.47 0.29 Mn 1.09 0.79 Fe 95.63 92.03 Fe 72.17 44.22 Fe 76.24 54.63 ♦ Maghemite Fe^HC^ • Hematite FejO^ A Magnetite Fe304 « Iron chloride FeClj ■ Halite NaCl O Goethite a-FeO(OH) Bernalite Fe(OH)3 A Lepidocrocite y-FeO(OH) © Molyzite FeCI3 ft Iron oxide chloride 0) Chromite FeCr203 2-Theta ■ Scale Fig. 11 SEM images and EDS microprobe analysis taken at points 1A, 1B, and 2 (above) and the chemical compounds (phases) observed by X-ray diffraction at the corroded specimen surface kept in NSS testing chamber for 96 hours 458 Advances in Production Engineering & Management 14(4) 2019 Surface features of chromium alloyed carbon steel specimens after salt-spray tests in NaCl solution Distribution of the corrosion products formed on the cylindrical surface is rather uneven and the presence of chloride was detected by means of EDS microprobe analysis. In the X-ray diffraction peaks iron (III) chloride and -oxychloride could be clearly identified together with some NaCl(s) inclusions (halite crystals). This reflects well the widely stated contribution and strong influence of the chloride ions in the corrosion mechanism of many metals in an aqueous environment in contact with air [19]. At least part of the highly water-soluble metal chlorides - often formed as chloride intermediates during a rather complex surface oxidation process of an iron base alloy -will migrate and can also be washed away from the corroding interface (anodic local areas) by moving streams of the liquid environment This effect can also reduce somewhat the overall growth rate of the often only loosely adhered non-soluble corrosion products of the otherwise highly porous rust layer(s). After careful removal of the adhered rust from the corroded specimens, the exposed remaining topography of the bare metal also revealed the uneven formation of rust with varied chemical composition. As to the roughness of the surface, it can be stated that with the increase in the duration of the corrosion tests, the roughness of the corroded surface increases. However, the roughness of the corroded surface can be measured only after the proper corrosion removal procedure, which has a significant effect on the surface topography as well. Therefore, the careful selection of the proper corrosion removal procedure is important. 4. Conclusions Based on our extensive and detailed experimental studies the following major conclusions could be drawn: • Neutral salt spray (NSS) corrosion testing of a precision machined cylindrical chromium alloyed steel shaft revealed uneven surface degradation. • Hence, in addition to identification of the surface corrosion products (listed in Fig. 11), experiments revealed that it is very important to determine also a detailed surface topography for the type of machine components due to their rather complex nature and the uneven degradation mechanism of the tested environment. • It was also observed that the initial roughness of the surface affects its corrosion behaviour during NSS testing. On the basis of the measured data it can be stated that the values of the roughness parameters were two and a half times higher (2.5) after the 24-hour test; those were increased four-fold (4.3) after the 96-hour test and more than eight-fold after 192 hours (8.7). It was observed that the values of the 3D roughness parameters (Sa, Sz, Sq) are on average twice as high as their 2D equivalents (Ra, Rz, Rq) on the corroded surfaces. This is because corrosion appears as small craters at certain points on the surface and therefore the surface has high inhomogeneity, which is better detected by 3D roughness characteristics. • Moreover, the NSS medium affects not only the roughness parameters but the form errors of the steel shaft as well. It was determined that the value of the deterioration factor for CYLt becomes approximately 3.1 times higher when the salt spray test time is 192 hours related to the value which resulted at 24 hours. This ratio is approximately 5.2 when examining the RONt circularity error parameter. • The investigations presented in the paper may help automotive companies to become more aware of the nature of corrosion on surfaces and to produce more corrosion-resistant components. Acknowledgement This research was supported by the European Union and the Hungarian State, co-financed by the European Regional Development Fund in the framework of the GINOP-2.3.4-15-2016-00004 project, aimed to promote the cooperation between the higher education and the industry. Advances in Production Engineering & Management 14(4) 2019 459 Varga, Török, Felho, Orosz-Szirmai, Réz References [1] Li, Y., Cheng, Y.F. (2016). Effect of surface finishing on early-stage corrosion of a carbon steel studied by electrochemical and atomic force microscope characterizations, Applied Surface Science, Vol. 366, 95-103, doi: 10.1016/ j.apsusc.2016.01.081. [2] Bajt Leban, M., Mikyska, C., Kosec, T., Markoli, B., Kovac, J. (2014). The effect of surface roughness on the corrosion properties of type AISI 304 stainless steel in diluted NaCl and urban rain solution, Journal of Materials Engineering and Performance, Vol. 23, No. 5, 1695-1702, doi: 10.1007/s11665-014-0940-9. [3] Xu, M., Zhang, Q., Yang, X., Wang, Z., Liu, J., Li, Z. (2016). Impact of surface roughness and humidity on X70 steel corrosion in supercritical CO2 mixture with SO2, H2O, and O2, The Journal of Supercritical Fluids, Vol. 107, 286297, doi: 10.1016/j.supflu.2015.09.017. [4] Pradhan, D., Mahobia, G.S., Chattopadhyay, K., Singh, V. (2018). Effect of surface roughness on corrosion behavior of the superalloy IN718 in simulated marine environment, Journal of Alloys and Compounds, Vol. 740, 250-263, doi: 10.1016/i.iallcom.2018.01.042. [5] Yayoglu, Y.E. (2016). Corrosion characteristics of magnesium under varying surface roughness conditions, Graduate theses and dissertations, University of South Florida, USA. [6] Walter, R., Kannan, M.B. (2011). Influence of surface roughness on the corrosion behaviour of magnesium alloy, Materials & Design, Vol. 32, No. 4, 2350-2354, doi: 10.1016/i.matdes.2010.12.016. [7] Almansour, A., Azizi, M., Jesri, A.M., Entakly, S. (2015). Effect of surface roughness on corrosion behavior of aluminum alloy 6061 in salt solution (3.5%NaCl), International Journal of Academic Scientific Research, Vol. 3, No. 4, 37-45. [8] Seo, M.J., Shim, H.-S., Kim, K.M., Hong S.-I., Hur, D.H. (2014). Influence of surface roughness on the corrosion behavior of Alloy 690TT in PWR primary water, Nuclear Engineering and Design, Vol. 280, 62-68, doi: 10.1016/ i.nucengdes.2014.08.023. [9] Hagen, C.M.H., Hognestad, A., Knudsen, O.0., Sorby, K. (2019). The effect of surface roughness on corrosion resistance of machined and epoxy coated steel, Progress in Organic Coatings, Vol. 130, 17-23, doi: 10.1016/ i.porgcoat.2019.01.030. [10] Lee, S.M., Lee, W.G., Kim, Y.H., Jang, H. (2012). Surface roughness and the corrosion resistance of 21Cr ferritic stainless steel, Corrosion Science, Vol. 63, 404-409, doi: 10.1016/i.corsci.2012.06.031. [11] Kostadin, T., Cukor, G., Jakovlievic, S. (2017). Analysis of corrosion resistance when turning martensitic stainless steel X20Cr13 under chilled air-cooling, Advances in Production Engineering & Management, Vol. 12, No. 2, 105114, doi: 10.14743/apem2017.2.243. [12] Toloei, A., Stoilov, V., Northwood, D. (2013). The relationship between surface roughness and corrosion, In: Proceedings of the ASME 2013 International Mechanical Engineering Congress and Exposition, San Diego, California, USA, doi: 10.1115/IMECE2013-65498. [13] Gathimba, N., Kitane, Y., Yoshida, T., Itoh, Y. (2019). Surface roughness characteristics of corroded steel pipe piles exposed to marine environment, Construction and Building Materials, Vol. 203, 267-281, doi: 10.1016/ i.conbuildmat.2019.01.092. [14] Rishikesan, V., Samuel, G.L. (2014). Evaluation of surface profile parameters of a machined surface using confo-cal displacement sensor, Procedia Materials Science, Vol. 5, 1385-1391, doi: 10.1016/i.mspro.2014.07.456. [15] Whitehouse, D.J. (2011). Handbook of Surface and Nanometrology, 2nd edition, CRC Press, Boca Raton, USA, doi: 10.1201/b10415. [16] Alcántara, J., de la Fuente, D., Chico, B., Simancas, J., Diaz, I., Morcillo, M. (2017). Marine atmospheric corrosion of carbon steel: A Review, Materials, Vol. 10, No. 4, 406, doi: 10.3390/ma10040406. [17] Roine, A. (2002). Outokumpu HSC Chemistry, Version 5.1, Chemical reaction and equilibrium software with extensive thermochemical database, Outokumpu Research Oy, Finland. [18] Jiang, S., Chai, F., Su, H., Yang, C. (2017). Influence of chromium on the flow-accelerated corrosion behaviour of low alloy steels in 3.5% NaCl solution, Corrosion Science, Vol. 123, 217-227, doi: 10.1016/i.corsci.2017.04.024. [19] Verma, C., Ebenso, E.E., Quraishi, M.A. (2017). Corrosion inhibitors for ferrous and non-ferrous metals and alloys in ionic sodium chloride solutions: A review, Journal of Molecular Liquids, Vol. 248, 927-942, doi: 10.1016/ i.molliq.2017.10.094. 460 Advances in Production Engineering & Management 14(4) 2019 Advances in Production Engineering & Management Volume 14 | Number 4 | December 2019 | pp 461-471 https://doi.Org/10.14743/apem2019.4.341 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper A new approach for product quality prediction of complex equipment by grey system theory: A case study of cutting tools for CNC machine tool Pang, J.H.a, Zhao, H.a, Qin, F.F.b*, Xue, X.B.a, Yuan, K.Y.a aCollege of Mechanical and Electrical Engineering, Wenzhou University, Wenzhou, P.R. China bSchool of Information Science and Engineering, Wenzhou Business College, Wenzhou, P.R. China A B S T R A C T A R T I C L E I N F O To compete in total global market, product quality has attracted the attention of manufacturers as an important mean of product differentiation. As effective product quality prediction method is the key technology for quality control system, a new prediction model and calculation method inspired by the grey system theory is proposed in this paper. Our practical evaluation shows that the quality of complex equipment was improved. Firstly, a new method of grey forecasting model for complex equipment was proposed, and the principle and method of grey predictive model with several variables were introduced. Secondly, this article discussed grey system theory model and showed how to use it in the forecasting process. Then, the quality prediction model and method using grey theory were set up with quality characteristics of cutting tools for Computer Numerical Control (CNC) machine tool. Finally, analysis of the test system showed that the applied predicting model and method were feasible and effective. This new method is also applicable to predict product quality of other complex electromechanical products which are composed a number of systems and subsystems. © 2019 CPE, University of Maribor. All rights reserved. Keywords: Quality control; Computer numerical control (CNC); Machine tool; Quality prediction; Grey system theory *Corresponding author: E-mail: helloqinfeifei8@sohu.com (Qin, F.F.) Article history: Received 30 September 2019 Revised 9 December 2019 Accepted 13 December 2019 1. Introduction Improving design and manufacturing quality are the reasons for the development of quality standard evaluation process of complex equipment production in the world. Complex equipment quality prediction is an important means for understanding the developmental tendency of electronic products. Product quality prediction using performance characteristics represents a difficult controllable variable in complex engineering environments. Due to quality characteristics, quality control, production technology and the technical conditions in engineering system practice, it is possible for errors to occur during product development and manufacturing processes [1]. The quality prediction phase is hence an important part of the manufacturing process, where emphasis is placed on managing resources and controlling operations to optimize costs and quality [2]. Therefore, a new quality predicting model and calculation method for complex equipment based on quality characteristics of engineering system should be identified for improving product quality under all foreseeable situations. In this paper, a new model and calculation method has been perfected to predict production quality for complex equipment before actual implementation in manufacturing process [3]. Using analysis and prediction methods to do quality management and prediction work has become 461 Pang, Zhao, Qin, Xue, Yuan general knowledge in many fields [4]. Boosting quality prediction in manufacturing enterprises and advancing their quality management level have become the basic direction of complex equipment manufacturing industry [5]. This also helps product designers in their design decision-making by providing sufficient feedback information on capability and level of manufacturing quality [6]. Based on numerous survey results, this paper analyzed the current status of quality prediction in the manufacturing industry and proposed new quality prediction model and calculation method using grey system theory [7]. The application of grey system theory in complex equipment quality prediction is very important and useful. In this paper, we developed a fundamental model providing supports on product quality prediction of complex equipment associated with a particular sampling policy for a manufacturing system [8]. Real-time and accurate short-term quality forecasting had become a critical matter in intelligent engineering systems. In this paper, the current complex equipment quality prediction methods were systematically summarized, and the main contents of grey system were elaborated. In addition, the grey system theory model was used for complex equipment quality prediction. Finally, an ensemble data assimilation procedure was developed and evaluated for the quality prediction model. Based on quality prediction models and methods have been discussed in many literatures and a lot of effective techniques and algorithms in the last 20 years. A methodology for sorting raw materials into homogenous groups with constant and optimized processing is presented to the most important factors causing unstable end-product quality by the fuzzy-c-means algorithm [9]. This paper has presented a simple data-based method in which measurements of other process variables are related to product quality, which were predicted recursively based on the measurements of reactor cooling rate [10]. A reliable model of a process both for the steady-state and unsteady-state regimes was established to predict the evolution of product composition with reliability of prediction and model reduction the neural model [11]. A novel offline modeling for product quality prediction of mineral processing was presented by using the least-squares support vector machine [12]. This paper has presented augmented reality and virtual reality as a successful tool in quality and defects management in construction industry [13]. A modeling approach to address the end-of-batch product quality prediction problem for batch processes was developed by using the three-way data set and the data set through the variable direction [14]. 20-lumped kinetics model and a linear partial least squares model was applied to the on-line aromatics yield prediction for a commercial continuous catalyst regeneration plat forming process [15]. An online product quality prediction method was developed based on offline clustering and online recognition of the operating modes a nonlinear dimension reduction method and multi-mode quality prediction method [16]. In order to predict software quality using a plain learner or a boosting algorithm which incorporates sampling, this paper presented the proposed techniques to several groups of datasets from two real-world software systems [17]. A novel online final product quality prediction scheme was proposed in this paper for the improvement of quality prediction in multi-phase batch processes, which explored the different effects of process variables in different phases on final product quality [18]. A developed model was used to predict final quality of products in a wide range of mineral content and temperature treatment data using principal component analysis and the influence of certain minerals [19]. This paper deals with the prediction of wafers quality in the semiconductor industry based on the pattern recognition principle by using a historical data of health indicators [20]. A survey of regularized linear regression methods using feature reduction and variable selection methods was presented to predict the water quality using the production equipment data and regression parameter optimization [21]. The potential of model predictive control on an offshore production unit starting was presented to control the gas-lift and ensure quality specifications of products of primary processing of petroleum through computer simulation [22]. A new methodology that utilizes a stochastic model to capture this complex relationship for better prediction that utilizes a linear model to represent the impact of tool wear was approved to capture the impact of quality degradation [23]. A novel Multi-Phase Support Vector Regression based soft sensor model for online quality prediction of glutamate concentration was presented by using the proposed soft sensor model for online product quality prediction [24]. An improved maximum like- 462 Advances in Production Engineering & Management 14(4) 2019 A new approach for product quality prediction of complex equipment by grey system theory: A case study of cutting tools ... lihood estimation method was presented to reliability modelling of CNC machine tools with enlarging censored failure data [25]. On the basis of the reviewed documents, this paper has proposed the product quality prediction model and method of complex equipment via grey system theory, and tested and revised the method by the empirical method. The rest of this paper is organized as follows. Section 2 discusses the proposed quality prediction model of the research, and calculation method via grey system theory. Section 3 discusses and analyzes the results of experiments conducted in tooling system of CNC machine tool. Finally, some useful conclusions and future work are summarized in Section 4. 2. Materials and methods 2.1 Proposed quality prediction model Aimed at process prediction in dynamic systems, a quality prediction model that meets the criteria for forecasting mechanism of time-varying quality characteristics is proposed in this article. Grey system theory is suitable for dynamic systems and is popularly applied in the prediction realm [26]. In this paper, a complex equipment quality prediction model was set up using the grey system theory which provided quality assurance of engineering systems. Here, an application example on the manufacturing process of complex equipment was used to validate the effectiveness of the prediction process and method. The key to a successful prediction of manufacturing processes by mathematical method is the selection of an appropriate model and quality characteristics of engineering system [27]. A viable solution of prediction system will be essential based on technology of information sharing and transparency in different departments and business units [28]. Based on properties of complex equipment manufacturing process, the prediction model in quality control based on gray system theory is illustrated in Fig. 1. A nonparametric prediction algorithm for the forecasting of recurrent random events in complex engineering systems is expressed in the proposed quality prediction model. In the former stage, the main task of product quality prediction is feature selection and data preprocessing. In addition, analysis and calculation are made through the grey system theory method, which formed the basis for system and quantitative analysis. Grey system theory is an effective method applicable to quality prediction when samples are small [29]. Then, on the basis of analyzing engineering system effluent quality prediction, we put forward a novel predictive method by grey system model GM (1, 1). The proposed forecast is a short-term quality prediction operation and is important for fault prevention of engineering machineries and prolonging their service life. The proposed quality prediction model is a fine prediction mode and method, which features automatic calculation and better reasonability and practicability. On the basis of above study, the proposed quality prediction model based on grey system theory is shown in Fig. 2. Input Performance char acteristics Process parameter variations Processing inconsistencies Pait stiffness and stability Parts accuracy Coupling between parts y Prediction process and method (Grey system theory) Output Reliability Stability Precision of retention Fig. 1 Quality prediction process and method Advances in Production Engineering & Management 14(4) 2019 463 Pang, Zhao, Qin, Xue, Yuan Fig. 2 Proposed quality prediction model via grey system theory As seen in the Fig. 2, the complex equipment quality prediction is an intelligent system based on related research in this study. In some ways, the prediction system greatly improved the quality forecasting and delivered a better decision backing for quality management of the manufacturing industry [30]. 2.2 Calculation method via grey system theory Grey system theory was firstly developed in 1989 and has been used as a qualitative and quantitative analysis method in a system when available information is uncertain and incomplete [31]. This system has developed rapidly in recent years [32] and was successfully used in the field of quality assessment, risk forecasting, information decision, etc. Generally speaking, grey system theory includes grey generating, grey relational analysis, grey forecasting, grey decision making and grey control [33]. Grey predicting is an important part of grey system theory. Grey forecasting is an important part of the grey system theory when information is imperfect [34]. The grey system model GM (1, 1) has been successfully applied to many fields of engineering system; moreover, it has been used in prediction realm [35]. Grey prediction model GM (1, 1) is applied to fit with the data sequence of performance characteristics in engineering system. Based on the grey system theory and method, a prediction model was devised to predict the quality level in the future was built. The relationship between the vibration and product quality prediction of complex equipment is analyzed with grey system theory and a new appraisal method for product quality prediction combining the qualitative and quantitative analysis by using the grey theory is composed here. The process analysis show that the grey theory is of optimal prediction effect for product quality prediction of complex equipment. The grey model GM (1, 1) of grey system theory for product quality is put forward, and it provides scientific basis for the quantitative prediction result. The influence degree of the different variables for product quality prediction of complex equipment was analyzed by using grey model GM (1, 1). According to the prediction result of manufacturing complex equipment, we analyzed and predicted the product quality in the molding process of grey system theory. Due to currently shortage of grey system theory quality predicting, an adaptive quality forecasting method of grey model GM (1, 1) is presented here. The differential equation of GM (1, 1) is defined as: ^+aX™=u (1) dt where X^ is defined as an accumulating sequence, and a, u are estimate parameters, which can be estimated by using least square method. 464 Advances in Production Engineering & Management 14(4) 2019 A new approach for product quality prediction of complex equipment by grey system theory: A case study of cutting tools ... (1) Construction of original sequence based on once accumulation method. The original sequence is defined in this paper as the following: x(o) = (j(o)(l), j(o)(2), xm(3),...,X<0>(n)) (2) where i = 1,2, —, n. Therefore, the original sequence can be generated by using once accumulation in the following equation: *(1)(o = (3) (2) Evaluation of parameters a and u Based on the grey prediction model GM (1, 1), the numerical values of parameters a,u can be calculated by using the generalized least squares method. The matrix and some constants are defined as: -[^«(1) + X<1>(2)]/2 1 -[^«(2) + X<1>(3)]/2 1 A = -[X^in- 1) + *<1>(n)]/2 1. (4) where the uncertain matrix A is norm bound and time-varying. The original sequence is defined as: yji = [z(0)(2),Z(0)(3),^,Z(°)(n)]T (5) Thus, an equation to identify parameters from experimental data can be expressed as: a= [a] = (ATA)"1ATyn (6) (3) Construction of the GM (1, 1) mode The GM (1, 1) model was constructed by using interactions and dependencies between these elements by the following forecast formula: £«(¿ + 1) = - (7) where (1) = (1), (i) = tfW (i) - tfW (i - 1), i = 2,3, -, n. (4) Verification of the accuracy of the prediction model Using the grey theory, our model accuracy can hence be compared, enhance verified with actual product quality feedback. By solving these equations, statistical expressions of expectation and mean square covariance about eigenvalues is given as following: 50 = jS2J(n - 1) (8) where S0 represents mean square covariance of the prediction model system to be minimized in the next step: 502 = i;=1[*(o)(o-*(o)]2 (9) *(0) = l"=1*(0)(0/n, i = 2,3,-,n (10) (5) Computation of mean-square deviation This simplifies the situation in the series = — where the mean-square de- viation was calculated. In general, the formula of the mean-square deviation is shown as followed: 51 = V512/(n-1) (11) 512 = Vn p)(0-,-C°)]2 (12) '¿=1 £-(0)= (13) Advances in Production Engineering & Management 14(4) 2019 465 Pang, Zhao, Qin, Xue, Yuan Variance analysis revealed that the proposed method in this paper is more effective than traditional iterative methods. To test the exactness and robustness of the improved algorithm, a variance ratio is used here as a measure of efficiency as followed: c = S1/S0 (14) The small probability of error can be computed in the same way: p = |£(o)(0 -£<0)| < 0.6745 x S0 (15) (6) Determining the accuracy grade of prediction model In Table 1, we categorized the prediction accuracy of different model by its small error probability and variance ratio. Accordingly, the prediction model GM (1, 1) is in the Accuracy grade I category with its reliable forecast result. Table 1 Grade specification of prediction accuracy Accuracy grade Precision grade Values of small error probability Values of variance ratio I Very good > 0.95 < 0.35 II Good > 0.80 < 0.50 m General >0.70 <0.65 IV Unqualified < 0.70 > 0.65 (7) Forecasting using obtained prediction model Computation results were compared with experimental findings, verifying that this model can accurately predict product quality. If the resultant test is found to be acceptable, existing data can be transformed into numerical prediction for future data, which is calculated as: tf<0>(n + 1) =£W(n+1)-£W(n) (16) + 2) = £(1)(n + 2) - £(1)(n + 1) (17) where + 1) and + 2) are established by comparing predicted values with actual data. The predicted values should be arranged according to the given test set order. The forecast data of grey prediction GM (1, 1) were used as the input variable with real life results as output. 3. A case study The application of grey system are introduced in this case study, the quality prediction problem of cutting tools for CNC machine tool of five axes simultaneous moving are also described. The grey system modeling is applied on evaluation of the structural and mechanical failure of CNC machine tool with the consideration of quality characteristics. To ensure better forecast accuracy, many satisfactory forecasting methods are assessed as the base of product quality prediction model of CNC machine tool. The cutting tools system is an important bridge that links the machine tool and the workpiece. The quality of cutting tools strongly influences the machining accuracy of the workpiece. This paper expounds the value of multiple quality characteristics to cutting tools in China manufacturing enterprises. The operation system of cutting tools for CNC machine tool is shown in Fig. 3. The quality characteristics of cutting tools for CNC machine tool are obtained through static analysis and laboratory tests. The extraction of quality characteristics of cutting tools has great importance during the normal operation of CNC machine tool. As a result of this processing, the quality prediction is more useful for human and machine perception in further manufacturing process tasks. In this paper, we review the present research and key techniques for process quality, including the preprocessing methods, feature extraction and classifier design methods. Based on the proposed prediction GM (1, 1) model of grey system, the mechanical properties and quality characteristics of cutting tools for CNC machine tool in this article are listed in Table 2. 466 Advances in Production Engineering & Management 14(4) 2019 A new approach for product quality prediction of complex equipment by grey system theory: A case study of cutting tools ... Fig. 3 Operation system of cutting tools for CNC machine tool Table 2 Mechanical properties and quality characteristics of cutting tools Number MTBF Hardness Bending strength Fracture toughness Elasticity modulus hour HRA MPa MPa- m1/2 GPa 1 921.56 93.6 650 7.52 508 2 928.12 93.5 660 7.58 512 3 932.87 94.2 655 8.21 522 4 938.43 94.3 720 8.53 530 5 942.35 94.9 745 8.43 541 6 948.24 94.7 750 8.53 528 Moreover, the focus of this case study is on the quality prediction method of cutting tools of CNC machine tool based on the five quality characteristics, including mean time between failures (MTBF), hardness, bending strength, fracture toughness and elasticity modulus. The quality criterion of MTBF is used widely for qualitative assessment. One effective way of describing strength of cutting tools is in terms of hardness. Additionally, the bending strength of cutting tools was also explained. The fracture toughness of material was adequate to ensure operation quality of the cutting tools. The elasticity modulus of cutting tools can be obtained by using experimental testing with analysis and inference methods. Table 2 compares quality characteristics of cutting tools based on key aspect of hardness, bending strength, fracture toughness and elastic modulus that affect the performance of machine tools. From the original sequence derived previously based on grey system GM (1, 1), we can get: ^(0) = (x(°)(1), *(°>(2),-, *W(6)) = (921.56,928.12,932.87,938.43,942.35,948.24) Based on Eq. 3, the original sequence can be computed by using once accumulation as summarized below: = (ZW(1),^(1)(2),--,^(1)(6)) = (921.56,1849.68,2782.55,3720.98,4663.33,5611.57) Then, the uncertain matrix A is given by Eq. 4 in the following: -1385.62 1 -2316.12 1 -3251.77 1 -4192.16 1 -5137.45 1- Advances in Production Engineering & Management 14(4) 2019 467 Pang, Zhao, Qin, Xue, Yuan Based on Eq. 5, the original sequence drawn from matrix A and the data from Table 2 are represented as shown below: "928.12" 932.87 yn= 938.43 942.35 -948.24- Then, the numerical values can be mixed and matched to create the appropriate combination for grey system theory. That is the key to computation in practice. The calculation process from this simplified theory of this new method is hence reported as followed. Ta - A'A r61825863.1394 L -16283.1050 (atA)-i = [ 16283.1050] 5.0000 J 0.00000011 0.00037016] 0.00037016 1.40545883Í On the other hand, a structurally simple equation to identify parameters from experimental data can be obtained by using Eq. 6: -0.00530] L920.73885J As mentioned above, the relationship between a and u is the following: - = -58.71186551 a So, the grey system GM (1, 1) is obtained based on Eq. 7: £(1)(i + 1) = 174615.08891=1__ 935.2617 S0 = JS$/(n-1) = 9.728 The relation between the number of sampling points and remainder error is presented in this paper. Then, the choice is between original value and a processed value. The results of the computations by using grey system are tabulated in Table 3. Table 3 Predicted value residual error of prediction model Number Original value Predicted value Residual error Relative error No. X(°\i) (%%) 1 921.56 921.5600 0 0 2 928.12 928.0817 0.038338 0.0041% 3 932.87 933.0144 -0.14443 -0.0155% 4 938.43 937.9734 0.456586 0.0487% 5 942.35 942.9588 -0.60876 -0.0646% 6 948.24 947.9706 0.269405 0.0284% From the comparison sheet, the main factor influencing accuracy of product quality forecasting is the uncertainty of original value for CNC machine tool. GM (1, 1) model method can meet the demand of forecasting while with less original data at high accuracy. To counter the fluctuating case of original data, the forecasting precision can be achieved by using the model of GM (1, 1) of grey system theory. The forecast result of high precision was got in forecast of quality prediction of CNC machine tool via GM (1,1) model by means of selecting revision of residual error, which was shown in Fig. 4. 468 Advances in Production Engineering & Management 14(4) 2019 A new approach for product quality prediction of complex equipment by grey system theory: A case study of cutting tools ... 940 838 936 934 932 930 923 926 924 922 920 1 • 1 • 1 1 1 — ■ -- Original value J* - • Predicted value X ( ft / - / / - 10 Fig. 4 Comparison chart of original value and predicted value In addition, the formula for the mean-square deviation is given by: £<°) =i2-=1£(0)(i) = 0.0019, = V^i/(6 - 1) = 0.3671 In the next step, the variance ratio can be given by Eqs. 14 and 15 in the following formula: c = S1/S0 = 0.3671/9.7281 = 0.0377, 0.6475 xS0 = 0.6475 X 9.7281 = 6.2989 Then, the small probability of error can be computed through the following: p = |£W(i) - ei°)| < 6.2989 |£(0)(1) - £<°)| = 0.0019 < 6.2989, |£^°)(2) - £<0>| = 0.0365 < 6.2989 |£W(3) - £<°)| = 0.1463 < 6.2989, |e(°)(4) - £(0)| = 0.4547 < 6.2989 |£(0)(5) - £(0)| = 0.6106 < 6.2989, |ei°)(6) - £(0)| = 0.2675 < 6.2989 Based on the variance ratio and the data from Table 1, it can be concluded that these forecasting models are acceptable, c = 0.0377 < 0.35, p = 1. Hence the model precision was found to be very good and can be used for forecast. Thus, the results of the comparison of original data and prediction data are displayed in Fig. 5. The serial test predictive value of four factors to MTBF was 953.01 hour. And the experiments show that the complex equipment quality prediction system based on more sensors and grey system theory model is efficacious. Although this paper focuses on CNC machine tool, the proposed methods apply to other complex equipment as well. Piediction data Original data 905 910 915 920 925 930 935 940 945 950 955 MTBF (Hour) Fig. 5 Contrast relation between original data and prediction data Advances in Production Engineering & Management 14(4) 2019 469 Pang, Zhao, Qin, Xue, Yuan 4. Conclusion The quality prediction technique and the numerical prediction method were studied in this article. A new algorithm based on grey system theory in the forecasting course is discussed and evaluated. We simulated the prediction process as an intelligence system for complex equipment as well as took quality characteristics prediction factors as input data and the prediction quantities as output data. This paper also compiled a microcomputer program which deals with the process prediction by using grey system theory. Self-adaptive prediction was implemented to promote accuracy of predicted results by adjusting input parameters of the engineering system. In order to improve the prediction precision, grey system theory model and method were used in these prediction processes. The proposed quality prediction model was used to represent past and current quality characteristics of engineering system and to predict future results of the manufacturing process. It is hoped to elevate the prediction technology and quality of complex equipment to better serve the manufacturing enterprises. This method can also be used in the field of other complex equipment quality prediction and quality target prediction. In conclusion, this paper presents calculations based on the data of CNC machine tool using grey forecast model. The results of the case study show that the quality level of complex equipment is automated deployment and utility computing. Verification results also showed that the quality prediction model could give satisfactory quality prediction accuracy of complex equipment. Acknowledgement This work was supported by the National Natural Science Foundation of China (No. 71671130), and the Provincial Natural Science Foundation, Zhejiang, China (No. LY20G010014). References [1] Farahani, A., Tohidi, H., Shoja, A. (2019). An integrated optimization of quality control chart parameters and preventive maintenance using Markov chain, Advances in Production Engineering & Management, Vol. 14, No. 1, 5-14, doi: 10.14743/apem2019.1.307. [2] Bhandari, S.H., Deshpande, S.M. (2011). Analysis of engineered surfaces for product quality monitoring, International Journal of Computers and Applications, Vol. 33, No. 4, 284-292, doi: 10.2316/lournal.202.2011.4. 202-2670. [3] Çiflikli, C., Kahya-Ozyirmidokuz, E. (2012). Enhancing product quality of a process, Industrial Management & Data Systems, Vol. 112, No. 8, 1181-1200, doi: 10.1108/02635571211264618. [4] Brahmachary, T.K., Ahmed, S., Mia, M.S. (2018). Health, safety and quality management practices in construction sector: A case study, Journal of System and Management Sciences, Vol. 8, No. 2, 47-64. [5] Davrajh, S., Bright, G. (2013). Advanced quality management system for product families in mass customization and reconfigurable manufacturing, Assembly Automation, Vol. 33, No. 2, 127-138, doi: 10.1108/0144515131130 6636. [6] Brombacher, A., Hopma, E., Ittoo, A., Lu, Y., Luyk, I., Maruster, L., Ribeiro, J., Weijters, T., Wortmann, H. (2012). Improving product quality and reliability with customer experience data, Quality and Reliability Engineering International, Vol. 28, No. 8, 873-886, doi: 10.1002/qre.1277. [7] Amanna, A., Price, M.J., Thamvichai, R. (2011). Grey systems theory applications to wireless communications, Analog Integrated Circuits and Signal Processing, Vol. 69, No. 2-3, 259-269, doi: 10.1007/s10470-011-9719-1. [8] Sosa, M., Mihm, J., Browning, T. (2011). Degree distribution and quality in complex engineered systems, Journal of Mechanical Design, Vol. 133, No. 10, Article ID 101008, doi: 10.1115/1.4004973. [9] Berget, I., Nœs, T. (2002). Optimal sorting of raw materials, based on the predicted end-product quality, Quality Engineering, Vol. 14, No. 3, 459-478, doi: 10.1081/QEN-120001883. [10] Pan, Y., Lee, J.H. (2003). Recursive data-based prediction and control of product quality for a PMMA batch process, Chemical Engineering Science, Vol. 58, No. 14, 3215-3221, doi: 10.1016/S0009-2509(03)00190-8. [11] Chetouani, Y. (2009). Model-order reduction based on artificial neural networks for accurate prediction of the product quality in a distillation column, International Journal of Automation and Control, Vol. 3, No. 4, 332-351, doi: 10.1504/IIAAC.2009.026780. [12] Ding, J., Chai, T., Wang, H. (2011). Offline modeling for product quality prediction of mineral processing using modeling error PDF shaping and entropy minimization, IEEE Transactions on Neural Networks, Vol. 22, No. 3, 408-419, doi: 10.1109/tNn.2010.2102362. [13] Ahmed, S., Hossain, M.M., Hoque, M.I. (2017). A brief discussion on augmented reality and virtual reality in construction industry, Journal of System and Management Sciences, Vol. 7, No. 3, 1-33. 470 Advances in Production Engineering & Management 14(4) 2019 A new approach for product quality prediction of complex equipment by grey system theory: A case study of cutting tools ... [14] Ge, Z., Song, Z., Gao, F. (2012). Statistical prediction of product quality in batch processes with limited batch-cycle data, Industrial & Engineering Chemistry Research, Vol. 51, No. 35, 11409-11416, doi: 10.1021/ie202554r. [15] Hou, W., Zhao, L., Zhang, L., Wang, G. (2013). On-line production quality prediction for a commercial naphtha catalytic reforming process, Information Technology Journal, Vol. 12, No. 18, 4553-4560, doi: 10.3923/itj.2013. 4553.4560. [16] Zhao, F., Lu, N., Yang, Y. (2013). Product quality prediction method for injection molding process based on operating mode recognition, Huagong Xuebao/CIESC Journal, Vol. 64, No. 7, 2526-2534, doi: 10.3969/j.issn. 0438-1157.2013.07.030. [17] Gao, K., Khoshgoftaar, T.M., Wald, R. (2014). The use of under- and oversampling within ensemble feature selection and classification for software quality prediction, International Journal of Reliability, Quality and Safety Engineering, Vol. 21, No. 1, doi: 10.1142/S0218539314500041. [18] Tang, X., Li, Y., Guo, J., Xie, Z. (2014). Final quality prediction for multi-phase batch process based on phase cumulative product quality model, Transactions of the Institute of Measurement and Control, Vol. 36, No. 5, 696708, doi: 10.1177/0142331213501688. [19] Arsenovic, M., Pezo, L., Stankovic, S., Radojevic, Z. (2015). Factor space differentiation of brick clays according to mineral content: Prediction of final brick product quality, Applied Clay Science, Vol. 115, 108-114, doi: 10.1016/ j.clay.2015.07.030. [20] Melhem, M., Ananou, B., Djeziri, M., Ouladsine, M. Pinaton, J. (2015). Prediction of the Wafer quality with respect to the production equipments data, IFAC-PapersOnLine, Vol. 48, No. 21, 78-84, doi: 10.1016/j.ifacol.2015.09.508. [21] Melhem, M., Ananou, B., Ouladsine, M., Pinaton, J. (2016). Regression methods for predicting the product's quality in the semiconductor manufacturing process, IFAC-PapersOnLine, Vol. 49, No. 12, 83-88, doi: 10.1016/ j.ifacol.2016.07.554. [22] Ribeiro, C.H.P., Miyoshi, S.C., Secchi, A.R., Bhaya, A. (2016). Model predictive control with quality requirements on petroleum production platforms, Journal of Petroleum Science and Engineering, Vol. 137, 10-21, doi: 10.1016/ j.petrol.2015.11.004. [23] Hao, L., Bian, L., Gebraeel, N., Shi, J. (2017). Residual life prediction of multistage manufacturing processes with interaction between tool wear and product quality degradation, IEEE Transactions on Automation Science and Engineering, Vol. 14, No. 2, 1211-1224, doi: 10.1109/TASE.2015.2513208. [24] Zheng, R., Pan, F. (2017). Multi-phase support vector regression soft sensor for online product quality prediction in glutamate fermentation process, American Journal of Biochemistry and Biotechnology, Vol. 13, No. 2, 90-98, doi: 10.3844/ajbbsp.2017.90.98. [25] Yang, Z., Zhu, D., Chen, C., Tian, H., Guo, J., Li, S. (2018). Reliability modelling of CNC machine tools based on the improved maximum likelihood estimation method, Mathematical Problems in Engineering, Vol. 2018, Article ID 4260508, doi: 10.1155/2018/4260508. [26] Bai, Y., Wang, P., Xie, J. (2014). Optimization of urban water supply schemes based on grey system theory, International Journal of Control and Automation, Vol. 7, No. 9, 239-246, doi: 10.14257/ijca.2014.7.9.20. [27] Longstaff, A.P., Fletcher, S., Parkinson, S., Myers, A. (2013). The role of measurement and modelling of machine tools in improving product quality, International Journal of Metrology and Quality Engineering, Vol. 4, No. 3, 177184, doi: 10.1051/ijmqe/2013054. [28] Nazifa, T.H., Ramachandran, K.K. (2018). Exploring the role of information sharing in supply chain management: A case study, Journal of System and Management Sciences, Vol. 8, No. 4, 13-37. [29] Chuang, T.-F., Chang, Y.-H. (2014). Comparison of physical characteristics between Rana latouchtii and Rana adenopleura using grey system theory and artificial neural network, Ecological Engineering, Vol. 68, 223-232, doi: 10.1016/j.ecoleng.2014.03.038. [30] Bettayeb, B., Bassetto, S.J., Sahnoun, M. (2014). Quality control planning to prevent excessive scrap production, Journal of Manufacturing Systems, Vol. 33, No. 3, 400-411, doi: 10.1016/j.jmsy.2014.01.001. [31] Lee, Y.-C., Wu, C.-H., Tsai, S.-B. (2014). Grey system theory and fuzzy time series forecasting for the growth of green electronic materials, International Journal of Production Research, Vol. 52, No. 10, 2931-2945, doi: 10.1080 /00207543.2013.857057. [32] Sang, G., Shi, K., Liu, Z., Gao, L. (2014). Missing data imputation based on grey system theory, International Journal of Hybrid Information Technology, Vol. 7, No. 2, 347-356, doi: 10.14257/ijhit.2014.7.2.30. [33] Zhou, Z., Sang, N. (2012). Image edge detection algorithm based on grey system theory, International Journal of Digital Content Technology and its Applications, Vol. 6, No. 23, 782-789. [34] Shao, Z. (2012). The evaluation on green textile production system based on grey theory, International Journal of Advancements in Computing Technology, Vol. 4, No. 2, 65-71. [35] Tserng, H.P., Ngo, T.L., Chen, P.C., Tran, L.Q. (2015). A grey system theory-based default prediction model for construction firms, Computer-Aided Civil and Infrastructure Engineering, Vol. 30, No. 2, 120-134, doi: 10.1111/ mice.12074. Advances in Production Engineering & Management 14(4) 2019 471 Advances in Production Engineering & Management Volume 14 | Number 4 | December 2019 | pp 472-482 https://doi.Org/10.14743/apem2019.4.342 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper Effect of the manufacturer quality inspection policy on the supply chain decision-making and profits Hu, H.a*, Wu, Q.a, Zhang, Z.a, Han, S.a aSchool of Economics and Management, Yanshan University, Qinghuangdao, P.R. China A B S T R A C T A R T I C L E I N F O Due to competitive pressure and information asymmetry, manufacturers will produce quality inspection avoidance behaviour to gain short-term economic benefits, but this behaviour affects the ultimate quality and safety of the product. This paper studies the two-echelon supply chain consisting of a manufacturer and a retailer, and analyses whether the manufacturer's quality inspection avoidance behaviour model is considered or not. This paper discusses the impact of quality inspection level, quality loss cost, product repair cost, product return rate on the profit and optimal decision-making behaviour of both actors of the supply chain. It is found that when the manufacturer's quality inspection avoidance level is high, the increase of retailer' quality inspection effort level, manufacturer's internal failure cost, consumer product return rate and retailer' external quality loss cost will lead to the decrease of manufacturer's quality effort level instead of increasing. Finally, the numerical study is given to verify the above conclusion, and analysed the influence of different parameters on the optimal decision and supply chain actors profits. © 2019 CPE, University of Maribor. All rights reserved. Keywords: Supply chain; Decision-making; Quality inspection policy; Quality inspection avoidance; Incentive mechanism; Product return; Profit *Corresponding author: huhaiju@ysu.edu.cn (Hu, H.) Article history: Received 10 July 2019 Revised 6 December 2019 Accepted 12 December 2019 1. Introduction Product quality is an important factor that affects consumers' purchase intention. Nowadays, more and more enterprises are outsourcing product design and manufacturing activities to manufacturers, and the price competition among enterprises is increasingly fierce. The huge competitive pressure requires manufacturers to further reduce production cost, which aggravates the incentive for manufacturers to conduct production with adulterated quality. The most commonly used measure of quality management and control in enterprises is quality inspection. But quality inspection has two major drawbacks: (a) they are costly, and (b) they are prone to lead to quality inspection evasion motivation. In recent years, there have been many safety accidents caused by product quality problems, such as the 2008 Sanlu milk powder illegal addition of melamine incident, the 2016 Samsung Note7 explosion, and Toyota recall. In this case, it is particularly important for companies to compensate for returned products in order to win back consumers [1]. A large number of facts show that some manufacturers conceal quality information and successfully avoid the inspection of buyers, industries and even countries. In order to control the quality of the products in the supply chain, it is necessary to coordinate the supply chain so that every participating enterprise in the supply chain is willing to strive to provide qualified products for customers. Jeong et al. [2] believes that the improvement of quality is conducive to the development of long-term cooperative relations among supply chain members. In contrast to information asymmetry, information symmetry can benefit upstream suppliers 472 Effect of the manufacturer quality inspection policy on the supply chain decision-making and profits [3]. In order to achieve supply chain coordination, various incentive strategies and supply chain contracts are generally considered to reduce the occurrence of such quality and safety incidents. The remainder of the study is proceeds as follows. Section 2 discusses the relevant literature. Section 3 develops our model. We analysed whether the manufacturer's quality inspection avoidance behaviour model is considered or not. In Section 4, incentive mechanism is put forward to influence the quality decision of manufacturer. Section 5 verifies our basic model with results and numerical illustrations. Section 6 present conclusions and limitations. 2. Literature review Three streams of literature are relevant with our research: (a) research on moral hazard caused by quality control strategies in the supply chain, (b) research on how to develop quality control contract and design quality incentive mechanism in supply chain, (c) research on how to implement quality risk control in the design of supply chain quality contract. Gao at al. [4] established the decision-making control model of moral hazard when they studied the cost sharing of internal and external quality losses. Fan at al. [5] considered the issue of corporate social responsibility, and believed that social responsibility is conducive to improving product quality and product demand, while making supply chain enterprises more profitable. Liu at al. [6] believes that CSR cost-sharing can maximize the profits of the whole supply chain and improve the competitiveness of supply chain members. Starbird [7], in the process of studying the design of supply chain quality contract, proposed how to design supply chain punishment, reward and inspection strategies, and put forward the quality incentive mechanism accordingly. Based on the incentive function of punishment and reward, Zhang at al. [8] proposed a cooperative mechanism combining order quantity, penalty cost and bonus. Chao at al. [9] considered the problem of information asymmetry under the cost-sharing contract of recall, and proposed how to conduct quality improvement incentive and contract design in product recall. Lee at al. [10] studied that buyback and revenue-sharing contracts could not achieve supply chain coordination under the uncertain quality of manufacturers and imperfect detection mechanism of retailers, and put forward quality compensation contracts. Hu at al. [11] think that the initiative recall of unqualified products and compensation can positively affect customers' purchase intention in the perspective of service recovery. Huang at al. [12] believes that offering money-back guarantee will help improve quality, demand, consumer surplus and retailers' profits, while reducing the return rate of defective products. McWilliams [13] pointed out that both high-quality products and low-quality products have the power to provide refund guarantee, and the refund guarantee can improve the profit level of low-quality enterprises. Zhang [14] believes that consumer characteristics (naive or sophisticated) exert subtle influence on the seller's strategy choice between partial returns and full returns. Giannoccaro and Pontrandolfo [15] proposed that revenue sharing contracts had a coordinating effect on supply chain members, but Xiao at al. [16] found that revenue sharing contracts were more beneficial to the upstream of the supply chain. Yan at al. [17] take the quality level of products provided by suppliers as reference standard. It is found that suppliers and retailers will voluntarily cooperate in accordance with contracts to achieve balanced results when retailers pay same attention to their own earnings and fairness. Zhu at al. [18] studied the buyer's determination of the inspection level based on the investment level of the manufacturer's production process, and the results showed that when the manufacturer's investment level increased, the buyer's quality inspection level decreased significantly. Babich and Tang [19] dealt with the adulteration of product quality through the delayed payment mechanism and the detection mechanism, and found that the detection mechanism could not completely prevent the adulteration of suppliers' products, while the delayed payment method could. Cao at al. [20] analysed the influence of product quality level on supply chain actors from the perspective of pricing and guarantee time, and found that manufacturers should set high wholesale price and long guarantee time for producing high-quality products. When consumers' sensitivity to warranty periods falls below the threshold, manufacturers of low-quality products can make more profit than their competitors. Cao at al. [21] compared the influence of deferred payment mechanism, inspection mechanism and traceability Advances in Production Engineering & Management 14(4) 2019 473 Hu, Wu, Zhang, Han mechanism on supplier adulteration. It was found that increasing the buyer's product liability cost could effectively reduce supplier's willingness to adulterate. By reducing inspection cost and traceability cost, the buyer could better restrain supplier adulteration behaviour. On the basis of the existing research, this paper considers the manufacturer's quality inspection evasive behaviour, discusses the quality decision and profit change of the upstream and downstream members of the supply chain in this case, and further expands the research on the relationship between manufacturer's quality management behaviour and supply chain performance. 3. Model formulation and notation A manufacturer M and a retailer R form the traditional two-echelon supply chain system. The retailer purchases a certain number of products from the manufacturer, and then sells them to consumers in the market for economic profit. In order to ensure the delivery quality of final products, retailer conducts quality inspection and incentive on manufacturer to avoid the cost of quality loss caused by potential quality risks. In the case of limited retailer' quality inspection level and information asymmetry, manufacturer will invest a certain cost to avoid retailer' quality inspection, so as to obtain greater profits in short term. The notation used in this paper is given in Table 1. Supposing the manufacturer's quality effort level is x, the manufacturer's evasive inspection effort level h, and the retailer's quality test effort level y, where x, h and y e (0,1) correspond to the product quality qualification rate, the success rate of evasive inspection, and the success rate of quality inspection, respectively. Hypothesis the inspection of retailer can correctly identify the products without defects, but there is a certain probability of accepting the defective products. If the quality inspection is passed, the retailer will purchase a certain number of products from the manufacturer at the wholesale price w and sell them to the customer at the market price p. If not, the product will be repaired or reprocessed by the manufacturer and then delivered to the retailer. The reprocessed product will pass the quality inspection with probability 1 [4], and the resulting internal failure cost is m. The manufacturer's unit production cost is c. In addition, the retailer's quality inspection level is limited. When customers find low-quality products and return them to the retailer, the retailer will face external quality loss cost s and product return rate t. Let ¿ax2, -¿by2, and ^eh2 respectively be the manufacturer's quality effort cost function, the retailer's quality inspection cost function, and the manufacturer's quality inspection avoidance cost function. The three cost functions are strictly convex. The higher the level of product quality, detection and evading detection, the higher the cost it will pay. Table 1 Model parameters and decision variables Notation Description X Manufacturer's quality effort level y Retailer's quality inspection effort level h The evading level of manufacturer's quality inspection m Internal failure cost per unit product t The possibility of return when the customer receives a nonconforming product P The retail price offered by the retailer to the customer w The wholesale price charged by a manufacturer to a retailer c Manufacturer's unit production cost a Manufacturer of quality coefficient b Retailer coefficient of quality testing effort e Manufacturer detection avoidance effort coefficient s Cost of quality penalty faced by retailer 474 Advances in Production Engineering & Management 14(4) 2019 Effect of the manufacturer quality inspection policy on the supply chain decision-making and profits Research basic hypothesis are: (1) Both the manufacturer and the retailer are rational decision-makers, risk-neutral, and each pursues the maximization of expected profits. (2) When the products produced by the manufacturer are qualified, there is no "rejection" phenomenon in the retailer' quality inspection. (3) s > p > w > c. When the product is returned, the retailer will not only return p to the consumers, but also need some extra expenses, such as logistics cost and economic compensation. (4) a> e. If the evading effort coefficient of the manufacturer is less than the quality effort coefficient of the manufacturer, the manufacturer will have the motive of evading the quality inspection. (5) m m. When the cost coefficient of manufacturer's quality effort is larger than the cost of manufacturer's quality loss, the manufacturer will have the motive of avoiding quality inspection. 3.1 Without considering manufacturer's quality inspection avoidance behaviour (model 1) In terms of the above model notations and assumption, the profits of the manufacturer and the retailer, nM and nR are defined as follows: nM(x) = w — c — (1 — x)ym — ^ax2 (1) n*(y) = p-w-(1-y)(1-x)ts-±by2 (2) (1 — x)ym is that the product has not passed the retailer's quality inspection and returned to the manufacturer for repair, so that the qualified rate is 1, resulting in internal failure cost per unit product; (1 — y)(1—x)ts is external penalty imposed by a retailer who fails to detect a defective product and then flows to the market and is returned by the consumer. We could have the optimal quality strategy of the manufacturer as x = R(y) = and the retailer as y = R(x) = By solving ^g^ = 0, = 0, we have manufacturer and retailer optimal strategy (3) mst+ab * ast r A -> y =--(4) J mst+ab Proposition 1: Without considering the manufacturer's quality inspection avoidance behaviour: The manufacturer's optimal quality effort level x increases with the retailer's inspection effort level y. With the increase of internal failure cost, external penalty and return rate, the manufacturer's optimal quality level will be improved. 3.2 Considering the manufacturer's quality inspection avoidance behaviour (model 2) The model assumes that the manufacturer has the motivation to evade quality inspection, and may cheat, forge or bribe the retailer in the process of quality inspection, expecting to invest with low quality cost and obtain excess economic profit. In this model, the manufacturer determines its quality effort level and quality inspection avoidance level, and the retailer determines the quality inspection level of the manufacturer. The behaviour choices of both actors are not understood by the other actor, and they make independent decisions with the goal of maximizing their own economic profits. Then, the expected profits of manufacturer and retailer are n'„(x,h) and n^(y). n^(x,h) = w — c — (1— x)y(1 — h)m — -eh2 - -ax2 (5) n^(y) = p-w-(1-y + yh)(1-x)ts--by2 (6) Advances in Production Engineering & Management 14(4) 2019 475 Hu, Wu, Zhang, Han d2nf (V) Because of —^— = —b<0, n^(x) has the optimal solution. From Eq. 5 to find the Hessian matrix about x and h for the manufacturer's profit function nj^, we can get H[n^(x, ft)] = r —a —ymi . . \—ym —e J. When ae — y¿m¿>0, there are the optimal solutions of Eq. 5. By solving = o, dnMÍx'h~) _q, dn^O) = 0, we have manufacturer and retailer option oh ay mal strategy. x = (7) h = ^^ (8) y= (9) The combination of the optimal strategy for (x, h, y) , can achieve their maximum benefit. 4. Retailer's quality incentives to manufacturer This chapter mainly considers the problem of quality incentive in the case of quality evading. Based on the consideration of manufacturer's quality inspection evading behaviour, it analyses several common measures in actual quality management, and whether there is the possibility of failure when the manufacturer has the motive of quality inspection evading. 4.1 incentive effect of quality inspection the level y Generally speaking, quality inspection can effectively improve the quality of the manufacturer's products. Denote = , £^ = q. We have ox oh ay Denote ^^ =0. We have oy = eym-f"f (10) ^ ae-y2m2 v 7 — me(y2m2-2amy+ae) (ae-y2m2)2 ( Theorem 1: When y e (0,-—-——), manufacturer's quality effort level increases with the re- (Q-Qg \ -—-,1 I, manufacturer's quality effort level decreases with the retailer's inspection level. Quality inspection will always improve the quality of the manufacturer's products in model 1, but in model 2, if the quality inspection exceeds a threshold, the manufacturer is more willing to use evasive methods. Denote ^f^ =0, ^^ =0. We have oxoy ahoy dx _ em(l-2ft) (12) dy ae-y2m2 dh _ am-2m2y(l-h) rt■3^ dy ae-y2m2 Proposition 2: If h >0.5, then <0, ^ >0. dy dy Because of a> e, manufacturer' quality evasion motivation will always be greater than its quality effort level. When the manufacturer has high quality avoidance motivation, the retailer will force the manufacturer to reduce the level of quality efforts and invest more in quality avoidance. Manufacturer can ease the pressure on quality inspection and repair by avoiding methods, but retailer will face more returns of unqualified products, seriously reducing its expected profits. 476 Advances in Production Engineering & Management 14(4) 2019 Effect of the manufacturer quality inspection policy on the supply chain decision-making and profits 4.2 Incentive effect of internal failure quality cost m The internal failure quality cost is the repair of nonconforming products, less than the external penalty cost flowing to the market, and effectively improves the retailer's expected revenue. Denote £I^=o, ^=0. We have dxdm dhdm dydm dm am K J am v ' dh dx e—=-—ym+(1-x) — m + (1-x)y (14) dm dmy dm y K dm dm v J dm dm Denote —amstx — 5bm2y2 — ehmst + amst + emst + abe = U, U > 0 , we have: dh _ aby(l-2x) dm U dx _ bey{l-2h) dm U dy asty(x-l)+esty{h-l)+4bmy3 Kdm U (15) Proposition 3: If x>1/2, then |^<0;If x<1/2, f^>0. If £>1/2, then |^<0. If dm dm dm h <1/2, then >0. dm If the quality level of the product is high, the increase in internal failure cost will weaken the manufacturer's motivation to avoid detection, and manufacturer is more willing to increase the investment in product quality. Otherwise, the greater the loss of internal quality, the stronger the motivation to avoid quality inspection. When manufacturer has strong quality avoidance motivation, internal failure cost will force manufacturer to reduce the level of quality efforts and invest more in quality inspection avoidance. Owing to x < h and a> e, the increase of quality avoidance is larger than that of quality effort at the beginning. If h < 1 /2, the increase of internal failure quality cost will reduce the manufacturer's profits, but the repair cost is greater than the investment of quality avoidance cost. 4.3 Incentive effect of unqualified product return rate t The parameter t is the percentage of unqualified products accepted by the retailer but returned by the consumers. As the defective products return rate increase, expected benefit of the seller is being lost in the form of compensation [22]. Denote ^f^ =0, ^f^ =0, ^ =0. We have: dxdt dhdt dydt dy aes(l-x-h) ~T = —----(16) dt U 1 J dx _ ems(l-x)(ft-l)(2ft-l) dt ~ U ( ) Proposition 4: If x + h >1, then ^ <0. If x + h <1, then ^ > 0. If h < 1/2, then ^ >0. If dt dt dt h>V2, £1, then ^ <0. If x +ft <1, then ^ > 0. If h< 1/2, then ^ >0. If * as as as h> 1/2, then — <0. ' ds When the manufacturer has strong quality avoidance motivation, the increase of external quality penalty cost will force the manufacturer to reduce the level of quality efforts and invest more in quality inspection avoidance. Quality inspection avoidance will replace the manufacturer's quality efforts. If ft <1 — x, retailer's quality inspection effort level increases with external quality penalty cost. 5. Results and discussion of numerical simulation In this section, numerical simulation will be used to compare the two situations between not considering manufacturer's quality evading behaviour and considering manufacturer's quality evading behaviour, and to analyse the influence of quality inspection level, internal failure cost, external quality penalty cost and product return rate on the quality decision-making of supply chain actors. x1,y1, M1 and R1 are respectively the manufacturer's product quality effort level, the retailer's quality inspection level, the manufacturer's expected profit and the retailer's expected profit without considering the manufacturer's quality inspection avoidance iour. x2, y2, M2 and R2 are respectively the manufacturer's product quality effort level, the retailer's quality inspection level, the manufacturer's expected benefits and the retailer's expected benefits considering manufacturer's quality inspection avoidance behaviour. We use the following parameters to illustrate this model, p = 35, c = 8, w = 18, a = 9, e = 8, b = 8, and its results are presented in Figs. 1 to 8. 5.1 Optimal decision and expected benefits with respect to y Here, let m = 7.5 , t = 0.5 , s = 36. The calculation result is shown in the following Figs. 1 and 2. Without considering manufacturer quality inspection avoidance behaviour, the quality effort level of the manufacturer increases with y. The manufacturer's expected benefits decrease with y. The retailer's expected benefits increase at first and then decreases with y. When considering manufacturer's quality avoidance behaviour, the level of manufacturer's quality inspection avoidance increases with y. When the evading level of manufacturer's quality inspection is greater than 0.5, the manufacturer's quality effort decreases with y, so manufacturer's quality inspection evading input will take the place of manufacturer's quality effort input. The expected benefits of the manufacturer decrease with y. The retailer's expected revenue is convex. In equilibrium, no matter how y changes, the level of manufacturer's quality effort in model 1 is greater than that in model 2. Under model 2, the manufacturer's profit is improved, but the retailer's profit is far less than model 1. 478 Advances in Production Engineering & Management 14(4) 2019 Effect of the manufacturer quality inspection policy on the supply chain decision-making and profits |-.vi.....h| I-Ml.....El m: R3 Fig. 1 Optimal decision with respect to y Fig. 2 Expected benefits with respect to y 5.2 Optimal decision and expected benefits with respect to m Here, let t = 0.5, s = 36, me [5,10], and the calculation result is shown in the Figs. 3 and 4. Without considering manufacturer' quality inspection avoidance behaviour, the manufacturer quality effort level increases with m. The retailer quality inspection level decreases with m. At the same time, the manufacturer's expected benefits decrease with m, and the retailer expected benefits increase with m, but the impact is not significant. When considering manufacturer's quality inspection behaviour, the retailer's quality inspection level decreases with m. Owing to x < 0.5, the evasion level of quality inspection of manufacturer increases with m. Manufacturer' quality effort level is convex function. When h > 0.5, the manufacturer's quality effort level decreases with m. The manufacturer and retailer's expected revenue decreases with m. In equilibrium, no matter how y changes, the level of manufacturer's quality effort in model 1 is greater than that in model 2. Although the level of retailer' quality inspection decreased with m, model 1 was more than model 2 at the beginning, and later less than model 2. Internal failure quality cost is negative with manufacturer's profit, and positive with retailer's profit. Under model 2, the manufacturer's profit is slightly improved, but the retailer's profit is far less than model 1. ■Ml.....El--M2 - Fig. 3 Optimal decision with respect to m Fig. 4 Expected benefits with respect to m Advances in Production Engineering & Management 14(4) 2019 479 Hu, Wu, Zhang, Han 5.3 Optimal decision and expected benefits with respect to t Here, let m = 7.5, s = 36, t e [0.1,0.8], and the calculation result is shown in the Figs. 5 and 6. Without considering manufacturer quality inspection avoidance behaviour, manufacturer quality effort level increases with t. Retailer's quality inspection level increases with t. Manufacturer and retailer's expected profits decrease with t, but when t > 0.3, the downward trend is gentle. When considering manufacturer quality inspection avoidance behaviour, the manufacturers' quality inspection avoidance level increases with t. Manufacturer' quality effort level is convex function. When h> 0.5, the manufacturer's quality effort level decreases with t. Owing to h + x < 1, retailer quality inspection efforts increase with t. At the same time, the manufacturer and the retailer expected the profit to decrease with t. In equilibrium, no matter how t changes, the level of manufacturer's quality effort in model 1 is greater than that in model 2. Although the level of retailer' quality inspection increased with t, model 1 was more than model 2 at the beginning, and later less than model 2. Unqualified product return rate is negative with manufacturer and retailer's profit. If the return rate is large enough, the retailer will face losses. D.2 0.3 0.4 O.i 0.6 DJ I [—— xl..... vl h x]--v] — -h-x]| Fig. 5 Optimal decision with respect to t -Ml..... El--M2 Fig. 6 Expected benefits with respect to t 5.4 Optimal decision and expected benefits with respect to s Here, let m = 7.5, t = 0.5, se [30,40], and the calculation result is shown in the Figs. 7 and 8 show that without considering manufacturer quality inspection avoidance behaviour, manufacturer quality effort level increases with s, retailer quality inspection level increases with s. Manufacturer and retailer expected profits decrease with s, but the downward trend is gentle. When considering manufacturer's quality inspection avoidance behaviour, the evasive level of manufacturer's quality inspection increases with s. Manufacturer' quality effort level is convex function. When h > 0.5, the manufacturer's quality effort level decreases with s. Since h + x < 1, retailer' quality inspection efforts increase with s. At the same time, the manufacturer and retailer's expected profit decrease with t. In equilibrium, no matter how s changes, the level of manufacturer's quality effort in model 1 is greater than that in model 2. In the two models, the optimal quality inspection level of the retailer is very similar. External quality penalty cost is negative with manufacturer and retailer's profit. In model 1, increasing external quality loss cost will improve the manufacturer's quality level, and the retailer can offset some negative impact through the improvement of quality assurance ability. 480 Advances in Production Engineering & Management 14(4) 2019 Effect of the manufacturer quality inspection policy on the supply chain decision-making and profits K 1..... vl h k]--V] — -li-x: I Ml.....R1 M2 R2 Fig. 7 Optimal decision with respect to s Fig. 8 Expected benefits with respect to s 6. Conclusion This paper studies the two-echelon supply chain system composed of manufacturer and retailer. In this system, the information asymmetry between manufacturer and retailer cannot eliminate the opportunistic behaviour of manufacturer due to the pursuit of short-term economic benefits. Comparative analysis of the level of quality inspection, internal quality loss, external quality loss, defective product return rate to actors of the supply chain optimal decision and expected benefits, it is concluded that the manufacturer' evading of quality inspection will weaken their motivation of quality effort and have a significant negative impact on retailer' revenue, which is not beneficial to coordination and cooperation of the supply chain from the long run. If the evading level of manufacturer's quality inspection is more than 0.5, manufacturer's product quality effort level will decrease with external quality penalty, product return rate, internal failure cost. The increase of retailer' investment in quality inspection will lead to the reduction of manufacturer' efforts in product quality and the further enhancement of inspection avoidance incentives. If the product quality level provided by the manufacturer is high, the level of quality inspection avoidance of the manufacturer will decrease with the increase of internal quality loss. If h + x < 1, retailer' investment in quality inspection will increase with product return rate and external quality penalty. Information sharing is a leading capability [23]. The retailer can invest in this capability to design a contract to achieve sustainable performance, rather than in quality inspection that would incur more costs. The development of supply chain actors should be based on the harmonious and win-win business circle. For the convenience of calculation and research, this paper only considers the traditional two-echelon supply chain model. In actual management of the enterprise, a complete supply chain is usually made up of multiple enterprise node, and affected by external environment. The quality decision-making of supply chain is relatively complex, and it is necessary to further study. In addition, although the model construction has drawn valuable conclusions, it is not very clear about the path and mechanism of manufacturer's evasive behaviour on quality and performance. Acknowledgement This work is supported by National Natural Science Foundation of China (No. 71704151), Social Science Foundation of Hebei province (HB19GL006), Research funding of Hebei Key Research Institute of Humanities and Social Sciences at Universities (JJJD1801). Advances in Production Engineering & Management 14(4) 2019 481 Hu, Wu, Zhang, Han References [1] Yoo, S.H. (2014). Product quality and return policy in a supply chain under risk aversion of a supplier, International Journal of Production Economics, Vol. 154, 146-155, doi: 10.1016/i.iipe.2014.04.012. [2] Jeong, E.-B., Park, G.-W., Yoo, S.H. (2019). Incentive mechanism for sustainable improvement in a supply chain, Sustainability, Vol. 11, No. 13, doi: 10.3390/su11133508. [3] Liu, Y., Li, J., Quan, B.-T., Yang, J.-B. (2019). Decision analysis and coordination of two-stage supply chain considering cost information asymmetry of corporate social responsibility, Journal of Cleaner Production, Vol. 228, 10731087, doi: 10.1016/i.iclepro.2019.04.247. [4] Gao, C., Cheng, T.C.E., Shen, H., Xu, L. (2016). Incentives for quality improvement efforts coordination in supply chains with partial cost allocation contract, International Journal of Production Research, Vol. 54, No. 20, 62166231, doi: 10.1080/00207543.2016.1191691. [5] Fan, J., Iang, X., Ni, D. (2019). A study on corporate social responsibility and product quality in supply chains under different channel power structures, Chinese Journal of Management, Vol. 16, No. 5, 754-764. [6] Liu, Y., Quan, B.-T., Li, J., Forrest, J.Y.-L. (2018). A supply chain coordination mechanism with cost sharing of corporate social responsibility, Sustainability, Vol. 10, No. 4, doi: 10.3390/su10041227. [7] Starbird, S.A. (2001). Penalties, rewards, and inspection: Provisions for quality in supply chain contracts, Journal of the Operational Research Society, Vol. 52, No. 1, 109-115, doi: 10.1057/palgrave.iors.2601048. [8] Zhang, C.-H., Ren, J.-Y., Yu, H.-B. (2006). Supply chain collaboration mechanism based on penalty and bonus under asymmetric information, Chinese Journal of Management Science, Vol. 14, No. 3, 32-37. [9] Chao, G.H., Iravani, S.M.R., Savaskan, R.C. (2009). Quality improvement incentives and product recall cost sharing contracts, Management Science, Vol. 55, No. 7, 1122-1138, doi: 10.1287/mnsc.1090.1008. [10] Lee, C.H., Rhee, B.-D., Cheng, T.C.E. (2013). Quality uncertainty and quality-compensation contract for supply chain coordination, European Journal of Operational Research, Vol. 228, No. 3, 582-591, doi: 10.1016/i.eior.2013. 02.027. [11] Hu, H., Diebarni, R., Zhao, X., Xiao, L., Flynn, B. (2017). Effect of different food recall strategies on consumers' reaction to different recall norms: A comparative study, Industrial Management & Data Systems, Vol. 117, No. 9, 2045-2063, doi: 10.1108/IMDS-10-2016-0464. [12] Huang, F., Song, H.-M., Yang, H., Yang, J.-M., Wang, Q.-L., Wu, J.-W. (2019). The impact of money back guarantees on quality and service of supply chain products, Industrial Engineering and Management, Vol. 24, No. 3. [13] McWilliams, B. (2012). Money-back guarantees: Helping the low-quality retailer, Management Science, Vol. 58, No. 8, 1521-1524, doi: 10.1287/mnsc.1110.1497. [14] Zhang, J. (2014). Quality improvement or perception enhancement? The role of consumer behavior and returns policy, Journal of System and Management Sciences, Vol. 4, No. 3, 56-80. [15] Giannoccaro, I., Pontrandolfo, P. (2004). Supply chain coordination by revenue sharing contracts, International Journal of Production Economics, Vol. 89, No. 2, 131-139, doi: 10.1016/s0925-5273(03)00047-1. [16] Xiao, D., Pan, K.-W. (2012). Quality control and coordination mechanism in supply chain based on revenue sharing contract, Chinese Journal of Management Science, Vol. 20, No. 4, 67-73. [17] Yan, F., Liang, G.-Q., Liu, X., Nin, L. (2018). Supply chain contract coordination considering supplier's quality investment under fairness preference, Operations Research and Management Science, Vol. 27, No. 3, 50-58. [18] Zhu, L.-L., Yu, T., Xia, T.S. (2013). Product quality control contract model in a two-echelon supply chain, Chinese Journal of Management Science, Vol. 21, No. 1, 71-79. [19] Babich, V., Tang, C.S. (2012). Managing opportunistic supplier product adulteration: Deferred payments, inspection, and combined mechanisms, Manufacturing & Service Operations Management, Vol. 14, No. 2, 301-314, doi: 10.1287/msom.1110.0366. [20] Cao, K., He, P. (2018). Price and warranty competition in a supply chain with a common retailer, INFOR: Information Systems and Operational Research, Vol. 56, No. 2, 225-246, doi: 10.1080/03155986.2017.1363590. [21] Cao, Y., Hu, H.-L., Wan, G.-Y. (2017). Game analysis and mechanism choices of supplier product adulteration behavior, Operations Research and Management Science, Vol. 26, No. 7, 54-63. [22] Kumar, V., Ekwall, D., Wang, L. (2016). Supply chain strategies for quality inspection under a customer return policy: A game theoretical approach, Entropy, Vol. 18, No. 12, doi: 10.3390/e18120440. [23] Nazifa, T.H., Ramachandran, K.K. (2019). Information sharing in supply chain management: A case study between the cooperative partners in manufacturing industry, Journal of System and Management Sciences, Vol. 9, No. 1, 1947. 482 Advances in Production Engineering & Management 14(4) 2019 APEM jowatal Advances in Production Engineering & Management Volume 14 | Number 4 | December 2019 | pp 483-493 https://doi.Org/10.14743/apem2019.4.343 ISSN 1854-625G Journal home: apem-journal.org Original scientific paper Hybrid fuzzy multi-attribute decision making model for evaluation of advanced digital technologies in manufacturing: Industry 4.0 perspective Medic, N.a*, Anišic, Z.a, Lalic, B.a, Marjanovic, U.a, Brezocnik, M.b aUniversity of Novi Sad, Faculty of Technical Sciences, Novi Sad, Serbia bUniversity of Maribor, Faculty of Mechanical Engineering, Maribor, Slovenia A B S T R A C T A R T I C L E I N F O Manufacturing is currently at a turning point from mass production to customized production. The implementation of the Industry 4.0 concept, leading to automation and digitalization of manufacturing processes, is therefore considered vital for companies that aim to follow emerging trends in production. Research in this field is primarily focused on companies from developed countries, while companies from transition countries have difficulties to adapt to new business environment. The aim of this paper is to evaluate the use of advanced digital technologies in manufacturing companies from transition countries (i.e. Serbia) in the context of Industry 4.0. To address this problem, an evaluation method based on Fuzzy Analytic Hierarchy Process (FAHP) and Preference Ranking Organization Method for Enrichment Evaluations (PROMETHEE) is proposed. FAHP was used to determine criteria weights as an input for PROMETHEE method which was then used to evaluate advanced digital technologies. For this purpose, the dataset from the European Manufacturing Survey gathered in 2018 from Serbian manufacturing companies is used. The results of this empirical research revealed that production planning and scheduling, digital exchange of data with suppliers/customers, and production control systems play vital role for manufacturers in the context of industry 4.0. These results could serve to manufacturers for their strategic orientation and decision making. © 2019 CPE, University of Maribor. All rights reserved. Keywords: Industry 4.0; Manufacturing; Digitalization; Advanced technologies; Multi-attribute decision making (MADM); Fuzzy analytic hierarchy process (FAHP); PROMETHEE method *Corresponding author: medic.nenad@uns.ac.rs (Medic, N.) Article history: Received 10 May 2018 Revised 12 December 2019 Accepted 15 December 2019 1. Introduction Ever since the beginning of industrialization, technological improvements have led to paradigm shifts which are called industrial revolutions [1]. The fourth industrial revolution (i.e. Industry 4.0) is triggered by the introduction of emerging technologies (e.g., Internet of things, wireless sensor networks, big data, cloud computing, embedded system, and mobile Internet) into the manufacturing environment [2]. The process of introducing Industry 4.0 in manufacturing companies should include the following types of integration [3]: • Horizontal integration through value networks to facilitate inter-corporation collaboration, • Vertical integration of hierarchical subsystems inside a factory to create a flexible and reconfigurable manufacturing system, • End-to-end engineering integration across the entire value chain to support product customization. 483 Medic, Anisic, Lalic, Marjanovic, Brezocnik Manufacturers that follow these trends should be able to produce customized and small-lot products efficiently and profitably. In order to achieve these standards, advanced digital technologies have become the focus of the research related to Industry 4.0 as they are considered as one of the main enablers of Industry 4.0 [4]. Having this in mind, the "smart factory" is recognized as one of the key features of Industry 4.0 [5]. The smart factory includes following advanced digital technologies: • Mobile/wireless devices for programming and operation of equipment and machinery [6], • Digital solutions in production (e.g. tablets, smartphones) [6], • Software for production planning and scheduling (e.g. ERP) [7], • Digital exchange of product/process data with suppliers/customers (e.g. supply chain management) [8], • Near real-time production control system (e.g. systems of centralized operating and machine data acquisition) [9], • Systems for automation and management of internal logistics (e.g. RFID) [1], • Product-lifecycle-management-systems [10], • Virtual reality or simulation [11]. Research related to Industry 4.0 is primarily conducted in manufacturing companies from developed countries, since this concept is developed in leading manufacturing economies of the world [12]. The aim of this research is to evaluate the use of advanced digital technologies in manufacturing companies in the context of Industry 4.0 in transition countries (i.e. Serbia). This evaluation includes a comparison of the aforementioned advanced digital technologies based on a set of criteria. For this purpose, Multi-Criteria Decision Making (MCDM) methods should be used. MCDM problems can be classified into two main categories: Multi-Attribute Decision Making (MADM) and Multi-Objective Decision Making (MODM). MADM is more appropriate for discrete problems associated with evaluation or ranging of predetermined and limited number of alternatives using a set of criteria. MODM methods are suitable for continuous problems of design or planning, with the aim of achieving aspired goals within given constraints [13]. Since the main concern of this research is to evaluate the use of advanced digital technologies in manufacturing companies, MADM methods will be used, as they are designed to deal with this kind of problems. MADM methods have emerged as a common tool in research related to manufacturing that involves evaluation procedures. Recently, hybrid MADM methods that combine different MADM methods have become increasingly present in literature. From the range of individual tools, only Analytic Hierarchy Process (AHP) is used more than hybrid MADM methods [14]. Furthermore, hybrid Fuzzy MADM (FMADM) methods are becoming more and more utilized in research. In most cases Fuzzy AHP (FAHP) was combined with other methods (i.e. TOPSIS, VIKOR, and PRO-METHEE) [15]. TOPSIS and VIKOR are compromise ranking methods proposed for determining the most preferred alternative based on the closeness to the ideal solution. PROMETHEE method is an outranking method which is based on the pairwise comparison in order to determine the dominance among alternatives [13]. For the evaluation of the use of advanced digital technologies in manufacturing companies it is more important to determine the dominance among alternatives by comparing them to each other, rather than focusing on finding out which of the alternatives is the closest to the ideal solution. Therefore, the PROMETHEE method seems to be more suitable for this research. Similar approach was proposed for selection of organizational innovations in manufacturing companies [16]. Furthermore, the literature review revealed that FAHP [17] and PROMETHEE [18] are primarily used in the research related to manufacturing sector. In the PROMETHEE method, it is assumed that the decision maker is able to appropriately weight the criteria, as there are no specific guidelines for this procedure. Therefore, it is usually combined with AHP, since it is recommended that PROMETHEE should be strengthened with the ideas of AHP in the phase of determining criteria weights [19]. Furthermore, fuzzy logic was introduced in the procedure of determining criteria weights with AHP to reduce vagueness and uncertainty of the decision-makers' judgement [20]. 484 Advances in Production Engineering & Management 14(4) 2019 Hybrid fuzzy multi-attribute decision making model for evaluation of advanced digital technologies in manufacturing ... In this paper, a hybrid FMADM method combining FAHP and PROMETHEE was employed to evaluate the use of advanced digital technologies in manufacturing companies in the context of Industry 4.0. More specifically, the main contribution of this paper is using a hybrid FMADM method combining FAHP and PROMETHEE to evaluate advanced digital technologies in manufacturing companies from transitional countries (i.e. Serbia) that contribute the most to the production principles of Industry 4.0. In this way, the research related to advanced digital technologies in the context of Industry 4.0 will be extended to transitional economies, since current research in this field is typically conducted in manufacturing companies from developed countries. The remainder of the paper is structured as follows. Section 2 describes the materials, methods and data that were used in this research, while Section 3 presents the research results and discussion. Finally, Section 4 contains the conclusion, including the identified limitations of the study and suggestions for further research. 2. Materials, methods, and data This work proposes a hybrid FMADM model for evaluating the use of advanced digital technologies in manufacturing companies. More specifically, advanced digital technologies are evaluated in terms of their contribution to the production principles of Industry 4.0. For this purpose, FAHP and PROMETHEE were used. FAHP was applied to determine criteria weights, while PROMETHEE was used for the evaluation of advanced digital technologies. The procedure of the proposed model is presented in Fig. 1. The AHP method was developed by Saaty [21]. It is based on pairwise comparison using a nine-point scale. The use of crisp numbers for pairwise comparison in traditional AHP is considered insufficient and imprecise due to the vagueness and uncertainty of the decision-makers' judgment [22]. In addition, the opinion of the decision makers is usually expressed in linguistic Fig. 1 General model for the evaluation of advanced digital technologies Advances in Production Engineering & Management 14(4) 2019 485 Medic, Anisic, Lalic, Marjanovic, Brezocnik form. As a result, fuzzy logic was introduced into pairwise comparison process of AHP to reduce this deficiency, as it is designed to deal with the problems concerning subjective uncertainty. Fuzzy set theory is based on the idea that the elements have a degree of membership in a fuzzy set [23]. Fuzzy membership functions (i.e. fuzzy numbers) that featured most often in fuzzy logic are the following: monotonic, triangular, and trapezoidal [24]. Triangular fuzzy numbers (TFNs) are the most utilized in FMADM studies, due to their suitability to the nature of experts' linguistic evaluations [25]. A TFN denoted as a = (I, m, u) where I < m < u, has the triangular-type membership function as in Eq. 1: 0, x< 1or x>u>\ x — I P-ai^) = \m-l' l