ISSN 1854-6250 APEM journal Advances in Production Engineering & Management Volume 13 | Number 3 | September 2018 Pfl Published by CPE apem-journal.org University of Maribor 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 deski@apem-journal.org Janez Gotlih desk2@apem-journal.org Website Technical Editor Lucija Brezocnik lucija.brezocnik@um.si 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 Lanndon A.Ocampo, University of the Philippines Cebu, 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 © 2018 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 13 | Number 3 | September 2018 | pp 233-368 Contents Scope and topics 236 Sustainable manufacturing - An overview and a conceptual framework for continuous 237 transformation and competitiveness Hussain, S.; Jahanzaib, M. A hybrid grey cuckoo search algorithm for job-shop scheduling problems under 254 fuzzy conditions Yang, F.; Ye, C.M.; Shi, M.H. Design, finite element analysis (FEA), and fabrication of custom titanium alloy cranial 267 implant using electron beam melting additive manufacturing Ameen, W.; Al-Ahmari, A.; Mohammed, M.K.; Abdulhameed, O.; Umer, U.; Moiduddin, K. Dynamic integration of process planning and scheduling using a discrete particle swarm 279 optimization algorithm Yu, M.R.; Yang, B.; Chen, Y. An integral algorithm for instantaneous uncut chip thickness measuring in 297 the milling process Li, Y.; Yang, Z.J.; Chen, C.; Song, Y.X.; Zhang, J.J.; Du, D.W. Change-point estimation for repairable systems combining bootstrap control charts and 307 clustering analysis: Performance analysis and a case study Yang, Z.J.; Du, X.J.; Chen, F.; Chen, C.H.; Tian, H.L.; He, J.L. Multi-objective optimization for delivering perishable products with mixed time windows 321 Wang, X.P.; Wang, M.; Ruan, J.H.; Li, Y. Game theoretic analysis of supply chain based on mean-variance approach under 333 cap-and-trade policy He, L.; Zhang, X.; Wang, Q.P.; Hu, C.L. Decision-making strategies in supply chain management with a waste-averse and 345 stockout-averse manufacturer Jian, M.; Wang, Y.L. Genetic programming method for modelling of cup height in deep drawing process 358 Gusel, L.; Boskovic, V.; Domitner, J.; Ficko, M.; Brezocnik, M. Calendar of events 366 Notes for contributors 367 Journal homepage: apem-journal.org ISSN 1854-6250 (print) ISSN 1855-6531 (on-line) ©2018 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-Integrate d M anufacturing Cutting and Forming Processes Decision Support Systems 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 Advances in Production Engineering & Management Volume 13 | Number 3 | September 2018 | pp 237-253 https://doi.Org/10.14743/apem2018.3.287 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper Sustainable manufacturing - An overview and a conceptual framework for continuous transformation and competitiveness Hussain, S.a*, Jahanzaib, M.a industrial Engineering Department, University of Engineering and Technology, Taxila, Pakistan A B S T R A C T A R T I C L E I N F O Sustainable manufacturing is advancing amidst the changing context and emerging paradigms. Business success and longterm sustainability today depends upon transforming continuously, and maintaining competitiveness through building context specific capabilities. To address these needs, this research presented a comprehensive overview of sustainable manufacturing, in first part, contributed by numerous aspects in the field. The second part proposes a conceptual framework comprising three interconnected elements: 'Ideal', 'Strategy', 'Architecture'. The ideal is best depicted as an exploration and choice context. Synthesis of stakeholders' desires and systematic discovery of opportunities, in context of larger containing systems, manifestes into desired attributes and characteristics of products, technology and enterprise system that are to be approached continuously. The strategy element is a match and transformation context. Strategic planning, focussed on continuous identification and building of capabilities, evolves into a broader concept of the business which enhances the firm's capacity to adapt to changing contexts and meet its proposed ends. The architecture is a function and execution context at the operational level. It combines capabilities, organization and operational structure, and value creation processes to perfom desired function in context of an agreed upon business concept. A systemic and potentially viable approach, embedded with specific capabilities to integrate sustainability into the core of a manufacturing business, is thus proposed which sets this research work apart. Contextual interrogation gives it a new way of constructing the big picture of the issue, i.e. sustainable manufacturing. This simplistic scheme is supplemented and guided by multi-aspect research in sustainable manufacturing, circular economy, capabilities, strategy and transformation, and systems thinking. The proposed framework is expected to fulfill the key needs of enterprises in sustainable manufacturing business, i.e. to transform and maintain competitiveness in a fast paced environment. © 2018 CPE, University of Maribor. All rights reserved. Keywords: Sustainable manufacturing (SM); Sustainability; Circular economy (CE); Strategy; Architecture; Capabilities; Systems thinking (ST) *Corresponding author: Sajj2u@yahoo.com (Hussain, S.) Article history: Received 17 November 2017 Revised 8 August 2018 Accepted 24 August 2018 1. Introduction Manufacturing is an indispensable feature of the modern society owing to its vital strength in terms of societal service, and pivotal role in contributing towards national and global economy. However, at the same time, these human activity systems collectively consume a considerable amount of resources, and generate wastes and emmissions thus affecting the environment and society [1]. It has been largely realized that natural resources and regenerative capacity of the ecosystem is not infinite. Current production and consumption patterns are being seen as largely responsible to assure a foreseeable future [2]. Consequently, manufacturing is evolving, within the framework of triple bottom line 3BL (environment, economy and society), for pursuing the big picture of sustainable development [3, 4]. ENVIRONMENT DOMAIN ECONOMIC DOMAIN SOCIAL DOMAIN Environmental issues (climate change, global warming, pollution, ozone layer depletion) Eco-system concerns (waste and emissions, landfill, extinction of species, bio-diversity, restoration and conservation) Policies, regulations and directives Environmental performance Impact assessment (frameworks, e.g. LCA) Design philosophies (DfE, bio-mimicry, cradle to cradle etc.) Resource conservation (minerals, energy, water) Material and energy tracking (frameworks, e.g. materials flow analysis, energy flow analysis Strategies (waste minimization, resource efficiency, eco-efficiency) Economic paradigm (e.g. Circular Economy) Socioeconomic trends Production and consumption patterns Policies, business strategy, business model Economic growth and development Economic advantage, competitive advantage Financial performance (value added, stakeholders' return, profit, potential financial benefits etc.) Economic performance (products/operations/facilities) Product lifecycle costs (frameworks, e.g. LCC) Economic impact (products, manufacturing chain) Societal needs and values, issues and trends, lifestyle and culture Social interaction and collaboration, feedback and inputs, perspective Policies, regulations and directives Organization and social behavior Social performance, responsibility, reputation Social value, social benefits (local, national, global) Product benefits (customer, consumer, society) Social equity, standard of living, quality of life Health and safety , working conditions, employment opportunities, education and training, community well-being, Social impact (frameworks, e.g. SLCA) Fig. 1 Sustainability aspects relevant to manufacturing business Sustainable Manufacturing (SM from here on) involves simultaneous consideration of multiple factors and their interaction across three dimensions. Interconnectivity and co-dependence of firms with others in the entire value network even inceases the breadth and overall complexity of the issue (SM), attributed to increased plurality and high degree of interdependence between elements [5-7]. Dealing with such a complex system necessitates a holistic understanding of the issue. Systems Thinking (ST from here on) is being recognized as one of the key competencies to address sustainability problems systemically across multiple domains and at different scales [8]. ST promotes holism which provides theoretical awareness of the system and avoids unintended consequences of sub-optimization [9, 10]. Sustainability, viewed as a system property, implies analysis and creation of sustainable value from the perspective of larger containing systems [11]. Yet, balancing the needs and benefits for 3BL requires robust methods incorporating an integratd analysis and evaluation of interaction among system elements. Nevertheless, exploring through vast literature, various aspects relevant to the issue (SM) are presented in three domains in Fig. 1. Manufacturing enterprises, among others, are faced with an increasing complexity, uncertainty, change and diversity. There are a number of factors whose interplay is supposed to give emergence to this complex and diverse milieu. Naming few are; constantly shifting external environment, growing environmental and social concerns, inter-connectedness of industrial problems, faster learning and development needs, and transformation in societal needs etc. Additionally, there are concerns on rapidly depleting natural resources, and increasing price and volatility of raw materials and commodities [12]. Moreover, it has been realized that efficient manufacturing processes, technological innovation, optimization abilities, and other success factors in context of industrial production (e.g. economy of scale) are insufficient to ameliorate the overall situation [13]. Further, it has been learnt that gradual improvements at the level of product or process are inadequate to embrace long term business success and sustainability. Sustainability frameworks and approaches, in relation to manufacturing, are generally aimed at development and implementation of sustainability strategy, sustainable business models, design and innovation capabilities, performance improvement and assessment, and decision making tools etc. Many of the proposed frameworks emphasize systemic management of sustainability concepts and requirements in context of broader systems, and a deep involvement of stakeholders. The proposed frameworks and approaches demonstrate some utility in different industrial aspects. But they lack, generally, in managing the issue (SM) holistically towards continuous transformation and competitiveness in manufacturing enterprises. Moreover, despite widespread literature and enormous research on the issue, yet there is need to explore approaches which are potentially viable amidst the changing contexts to ensure long term sustainability. Nevertheless, researchers are learning about the capabilities of various frameworks that can bring about a change at the level of entire enterprise as a system. Numerous propositions can further increase the speed of learning and implementation. SM businesses have to seek ways to contribute to sustainable development objectives effectively, now and in the future, while advancing amidst the changing contexts. This needs continuous support and guidance that could systematically enable the enterprise's capability to meet its desired sustainability ends. Successful transformation requires more relevant and focused designs, and a mind-set change in relation to existential purpose of the firm, its stakeholders, society and environment A framework, however, may be established that serves as a guideline for the enterprise to create and re-create itself, to fit the proposed ends. The conceptual framewok proposed in this research is an abstract concept structured by three interacting elements: 'Ideal', 'Strategy', and 'Architecture'. The 'Ideal' or ideal intent of the enterprise is chosen by stakeholders through exploration of opportunities within the contextual space. This proposition may be sussessfully realized provided the enterprise system is technologically and operationally capable enough and the one that can learn and adapt effectively. The 'Strategy' or strategic intent is formulated in context of ideal intent. It represents a mix of manufacturing and business level capabilities whose development and strategic deployment is potentially vital to gain competitive advantage and capture the proposed ends at large. Strategic planning and decision making is aimed towards evolving an innovative business concept that is operationally viable. The 'Architecture' at down stream converts this concept into reality. It formally describes the overall business architecture framework for SM including its functional aspects, structural issues, and value creation processes and activities. This simplistic scheme may herald long term business success and sustainability provided the concept is operationalized in a systematic and meaningful manner. The framework design and function is complemented by emerging concepts and approaches in ST, SM research and Circular Economy (CE from here on), capabilities, and various frameworks for strategy and architectrure. Contextual interogation gives it a new way of constructing the big picture of the issue, i.e SM. Rest of the paper unfolds this scheme in quite a natural order. After the conceptual framework has been introduced in this section, the next section provides a comprehensive overview of the issue (SM) in various aspects and from multiple perspectives, to construct its contextual position. Section 3 describes proposed framework and its individual elements in detail. The relationship of proposed conceptual approach with social environment and level of economic development in a country along with its operationalization and main obstacles pertinent to a specific environment are discussed in Section 4. Finally, major findings from our discussion are drawn in Section 5, which concludes this research. 2. Sustainable manufacturing: An overview SM has been received with even greater recognition in recent years as a comprehensive strategy for improving the economic, social and environmental performance of a firm. New management concepts (product lifecycle management and assessment, eco-effectiveness etc) , innovative technologies, and design and engineering approaches have been introduced to incorporate sus-tainability in manufacturing [14]. Research frameworks encompass numerous aspects at varying scale (enterprise, supply chain, manufacturing cell and operations etc.). Adopting a system design approach, Tonelli et al. [15] proposed a framework for sustainable industrial system involving supply chain design, sustainability performance measurement and management, and organisational change . Using DMAIC (Define, Measure, Analyze, Improve, Control) approach, Koho et al. [16] proposed a concept aimed at supporting Finnish manufacturing companies in realizing sustainable production while emphasizing measurement as basis of improvement. Subic et al. [17] developed a framework focused on development and assessment of relevant capabilities across the entire supply chain. Zhang et al. [7] yearned for a framework, in line with LCA and SLCA, to integrate sustainability assessment with decision making methodology that assists in process planning at the production cell level. For sustainability in machining operations, Masood et al. [18] proposed a milling technique with the aim to enhance cutting tool life. Hassine et al. [19] presented a case study on sustainable turning operation utilizing a specifically developed multi objective optimization technique. An emerging research interest is observed on creating value through sustainable industrial production systems (Sustainable Supply Chains SSCs and Closed Loop Supply Chains CLSCs). Govindan et al. [20] highlighted that reverse logistics (RL) and CLSC is now a revenue opportunity instead of a cost-minimization approach. Schenkel et al. [21] segregated value creation into four types of values (economic, environmental, information and customer) and presented an insight on maximizing the value by coordinating the forward and reverse supply chains. Ageron et al. [22] recognized that sustainability issues when integrated in SCM, as an intended strategy, contribute to business in terms of competitive advantage. Gopal Krishnan et al. [4] argued the necessity of devising sustainable operations along the value chain and approaching the 3BL in an integrated manner. Galal et al. [23] studied sustainability, mainly emission reduction, across an agri-food supply chain by developing a discrete-event simulation model. From the system's perspective, SM needs to be analyzed in context of larger containing systems [11]. Understanding of external interactions becomes significant becacause of a broad range of factors across the product life cycle and their connectivity with broader systems (supply chain, market, and ecosystem) [24]. Fiksel [25] developed a triple value 3V model by dividing the physical world into three types of interacting systems (industrial, societal, and environmental) ,wherein the intricate dynamic linkages and flows of value among them may provide an architecture for system modelling. On the concepts of complexity theory, Moldavska [26] proposed a systemic model whereby a system contributes to not only its own sustainability, defined in terms of attractors, but also to that of larger containing system. From an operational perspective, Zhang et al. [7] proposed definitions and assessment framework for SM based on ST. Moldavska et al. [27] proposed a framework for sustainability assessment system for manufacturing organizations , guided by ST. A sustainable product is just as effective and at least the same quality as the product it is replacing while produced using fewer resources and better materials. In addition to product's attributes, SM focuses on how the product is made, with the aim to minimize its impact. In a lifecy-cle view, it includes input materials, product design, production and distribution processes, and end of life (EoL) management. Within Lifecycle Thinking framework, manufacturers look beyond the supply chain to maximize the product performance throughout the life cycle. Lifecycle Engineering (LCE) is aimed at designing products that meet requirements for quality, cost, manufac-turability, and consumer appeal, while at the same time minimising of environmental impacts. Product life cycle, if analysed systematically, may drive advancements in product and process technology, and other improvement strategies. Sharing and managing of product data, information and knowledge across this entire lifecycle forms the essence of Product Lifecycle Management (PLM) which seeks to integrate stakeholders, processes, resources and information in a product-centric knowledge environment, to make informed product decisions [28]. SM Principles include; using less energy and materials, using renewable energy, using fewer hazardous materials and toxic chemicals, closing the loop of resource flow, and designing products that can be easily and economically reused or re-manufactured and consume less energy during use [15]. Similar principles have been proposed by other researchers towards achieving sustainable ends. Alayon et al. [29] referred to 9 principles of sustainable production relevant to products, energy and material; economic performance; worker and community development; and natural environment. Esmaeilian et al. [30] reffered to seven main steps proposed by Ray Anderson, to become what he called a 'Prototypical Company of the 21st Century'. For long-term sustainability, economic growth needs to be decoupled from environmental constraints through technological changes and innovation in industrial systems [31]. CE and C2C are relevant frameworks here. Based on principles of Green Engineering, they approach sustainability from an eco-effective design perspective. Principles together with relatively new economic paradigm of CE, and manufacturing paradigms (e.g. C2C, industry 4.0), may create sustainable competitive advantage alongside a positive impact on environment and society. Value creation in CE, as popularized by EMF, is based principally on elimination of waste and use of renewable energy during the entire lifecycle. Resource yields are optimised by maximizing the number of technical cycles (repair, reuse, or remanufacturing) and the time spent in each of the multiple use cycles of product. Increased inflow of virgin materials into the economy may be substituted by cascading or diversifying the reuse of products across value chain. Moreover, the use of uncontaminated material streams (pure inputs) can increase collection and redistribution efficiency, and maintain material quality thus improving product longevity and material productivity. Furthermore, the principle of dematerialization(minimizing comparative materials use) is at play during processes of design, making and distribution of products to ensure re-useability with lesser changes, faster returns, and higher potential savings [12]. Contrary to the linear production pattern (take, make, waste), this regenerative mode is manifested in higher resource utility and optimization while minimising supply risks and enhancing of natural capital [32]. C2C design protocol, based on the strategy of designing ecologically intelligent products and processes, offers a good premise for CE. It redefines the problem to address its source, and shifts the end of pipe liabilities to product design. A design of industrial systems on this concept necessitates flow of materials in safe, regenerative, closed-loop cycles. To achieve an optimal effectiveness close to the natural system, inputs and outputs should be safe and beneficial [33, 34]. Industry 4.0 is an evolutionary industrial paradigm centered on achieving quality, efficiency and flexibility at the same time. Businesses need to converge their capabilities, to match with this era of fast paced and broad scale disruptions [35]. The manufacturing systems in 'Industry 4.0' utilize information and communication technology (ICT from here on) infrastructure in an internet of things to deliver smart products. Cross-linked value creation modules embedded with ICT and Cyber-Physical Systems (CPS) enable efficient allocation of resources , decentralized decision-making, and efficient work design and organization [36]. SM practices are aimed at reduction in each of energy, water, waste and emissions, improved product quality, enhanced corporate image and market competitiveness, better access to new markets, reduced costs, and increased profits. Practices range from very simple process improvements to innovative product design [3], and large investments in new technologies for cleaner production. Creation of wildlife habitats, promotion of community engagement and other external practices e.g. product stewardship can help integrate stakeholders' views into the business operations while contributing to a more sustainable world [37]. Firms primarily focus on conserving resource use and reducing physical waste which often widens to lifecycle sustain-ability through adoption of suitable policies and tools whereby an initial focus on internal operations grows to involve entire supply chain. Integrating SM concepts and strategies into the core of business requires long-term decision making [29] and radically rethinking the existing business models. In a system's view, activities across all business functions need to be focussed and coordinated instead of isolated improvements in either of product design, manufacturing technology or product end-of-life management. SM strategies can be broadly categorized under the two themes or dimensions [14]; Eco-efficiency is aimed at increasing economic development through lowering of resource intensity throughout the life cycle [38] while Eco-effectivity is directly tied to the reduction of environmental or social burdens [14]. Strategies take on a different perspective at each level of enterprise. At the manufacturing process level, the focus is on improving efficiency and quality, optimization of process parameters, and reduction in processing cost. At the supply chain level, consideration may be given to the use of renewable energy across the chain, recirculation of products and materials into product lifecycles, and effective management of raw materials and services for an improvement in social and environmental impacts. At the enterprise level, the focus is on policy analysis, sustainable business development, and achievement of business objectives like productivity, competitiveness and profitability. Firms may struggle at the onset of their efforts towards SM. Yet many firms embrace an early success. Their efforts and the manner in which certain intiatives were taken could be mimiced and modelled for greater understanding. Moreover, challenges encountered by others in the course of successful transformation may be helpful from the perspective of re-using the relevant experiences and knowledge. Manufacturing businesses need support and guidance if they are to contribute to sustainable ends. Challenges and barriers faced need to be clarified. Future challenges include competition in terms of both inputs (resources) and outputs (products and services) , diversity of customer base, rapidly changing consumer lifestyles and preferences, increasing business risks and uncertainties, adherence to regulations, avoiding use of technologies that have adverse impacts on environment, and focussing on total economic, social, and environmental cost of a product instead of low cost production [39]. The financial aspect of sustainability initiatives can be an enormous obstacle to the overall progress. Capital investment decisions need to be based on cost benefit justification because technologies required in reducing plant emissions to zero might be both complex and costly. Other commonly encountered barriers include time and resource constraints, general business concerns, and corporate culture. Aligning internal company culture with sustainability objectives in day-to-day business operations is not an easy. Manufacturing businesses need fundamental changes in their missions, and significant improvements in ways of doing business [16]. Existing strategies, procedures and standards might need a radical transformation to adapt to this new role. Transformation might need a redesign of the organization which may pose difficulties in several aspects e.g. defining enterprise structure, boundaries, decision making processes, and resolving conflicting goals etc. [40]. Manufacturing businesses today are coped with challenges of increasing product variety and shortening product life cycles [5, 41]. While large and diversified product-mix can capture additional business opportunities, added economic benefits and enhanced consumer value, it also poses increased complexity and unpredictability to the manufacturing system. In case of SM, the concerns are even larger in terms of managing variety across the entire product life cycle. Manufacturing system's flexibility and adaptability is a critical requirement to accommodate changes brought about due to management of multiple generations of a variety of products. From the economic perspective, effective resource utilization, and integration of product design with manufacturing processes will have to be ensured to attain sustainable competitive advantage in these environments [5]. Sustainability assessment, performance measurement and reporting is a challenging task in manufacturing organizations [7, 27] due to the complexity and dynamics of manufacturing operations and diverse theoretical perspectives on sustainability requirements, often normative in nature [24]. WBCSD, GRI, OECD, ISO, UNEP etc. have developed standards and guidelines, yet there are difficulties in several aspects. GRI's general as well as sector-specific assessment and reporting guidelines need customization for their befitting to the needs of a firm and its stakeholders. Life Cycle Thinking (LCT), a more popular framework, integrates three elements-Lifecycle Assessment LCA, Lifecycle Costing LCC, and Social Lifecycle Assessment SLCA- which together provide a methodological framework for Life Cycle Sustainability Assessment (LCSA) [42, 43]. Although an LCA study has certain limitations and challenges, it is a robust and most widely used framework for evaluation of potential impacts of a product during entire life cycle [7]. However, SLCA and LCC methods are relatively less developed and lack scientific consensus - too few tools and fully developed databases are available - hence rarely conducted in process industries [43]. Despite huge amount of research in sustaiability assessment, there is an absence of comprehensive, universally accepted and widely applicable metrics. Organizations find limitations in terms of interpretation and implementation capabilities relevant to a variety of tools, indicators, and metrics on the list. Organizations are even recognizing the importance of developing sustainability tools and metrics specific to their unique business environment [27] . The real value of assessment lies in attaining of sustainability objectives and providing valuable input to decision-makers rather than external reporting only [27, 44]. To avoid sub-optimization, there is a need to develop frameworks and tools, from a unified and holistic perspective, capable to assess and monitor the entire organization as a system instead of focussing on individual areas and functions. Moreover, the value added through such an endeavour should be in balance with the resources required [27]. 3. Proposed conceptual framework Continuous renewal of enterprise intent and translating it into an implementation plan is a contested problem [45]. Understanding the nature and purpose of enterprise, as a starting point, may guide in establishing requirements from a functional perspective. An enterprise as a system consists of structure and processes which collectively perform economic function, characterized by products and services delivered in context of an environment [46, 47]. The principal objective of an enterprise is the 'development' that serves to contribute to the development of itself and the larger containing systems. Greater this capacitation, better the firms can attain their objectives [48]. The framework proposed in this research is based on enterprise's desire to contribute to the sustainability of society, environment and economy, and stay compettitive in relation to changing contextual position of SM. As presented in Fig. 2, the framework as a concept system contains three functional elements or subsystems interrelated with each other: 'Ideal', 'Strategy', 'Architecture', and is it self part of a larger whole, i.e the enterprise. The ideal is best depicted as an exploration and choice context. Systematic discovery of opportunities within all important business dimensions and a coevolved proposition of products and enterprise system are what un-derlye a holistic and clarified vision embedded with stakeholders' desires. The enhanced capacity to sense opportunities manifests into an improved understanging of what satisfies stakeholders, and which attributes (of product, system) can lead competitive advantage? The ideal intent thus evolved serves as a reference to be approached continuously throughout the transformation cycle. The strategy element is best represented as a match and transformation context. Strategic planning is aimed at identification of capabilities corresponding with opportunities; formulation of policies and guidelines to enble an effective and efficient use of enterprise's means; setting of business objectives and goals; and commitment of resources and assets. The desirable consequence of this planning process is a broader concept of the business. Continuous identification and development of capabilities enhances capacity to adapt to changing contexts. The effectiveness is reflected in the speed and degree to which specific capabilities can be created and matched with opportunities. The architecture, at the operational level, is a function and execution context. It combines capabilities, organization structure and business processes to perfom desired functions in context of an agreed upon business concept. The context of the issue (SM) and hence the ideal intent keep changing. In response, the strategic planning reformulates the business concept and renews relevant capabilities. The architecture reconfigures capabilities and assets to adapt to changing business needs. This continuous process harnesses the capabilities and collective learning of the enterprise as a whole under dynamic business requirements. By operationalizing the individual elements and managing interaction among them, a manufacturing enterprise can be transformed continuously and systematically in relation to its desired sustainability ends. In subsequent part individual elements are explained in terms of their relationship with each other and with larger systems. Various aspects of SM are discussed alongside as deemed relevant to the nature and scope of each element. Fig. 2 Structural representation of framework 3.1 Ideal Sustainable products should not only be safe for environment but also create opportunities for life-systems, in addition to satisfying the needs of business and society. Determining what the environment (life system) actually wants, to attain its own ends, is beyond scientific deliberation, and seems unattainable. Nevertheless, businesses can approach this limit by consistently turning social needs and desires to environmental advantage specifically, rather merely creating opportunities that could positively impact [49]. This could serve as an ideal to be approached endlessly. Apple has recently pronounced to move its entire supply chain to closed loop and renewable energy in the wake of undeniable climate change and depleting earth resources. Though it seems unattainable currently as the report says, "It sounds crazy, but we're working on it", Apple has embarked on its journey by launching projects related to capabilities enhancement and technology development [50]. An SM enterprise delivers in context of a hierarchy of larger systems: industry (economic system); marketplace and society (social system); and eco-system or environment. Stakeholders' interests can significantly drive manufacturing towards sustainability and competitiveness [51]. This whole system's view helps understand the way a manufacturing enterprise interacts with its surroundings, throughout the entire life cycle of products. More pragmatically, however, the ideal can be chosen in relation to opportunities within present and future context of SM, continuously changing business environment, and SM principles, practices and industry callenges. Because external environment keeps changing, it necessitates proposing of products and technology in a larger context, i.e. emerging social and environmental needs, and trends and changes in business environment (e.g industry drivers, technological innovation, product lines and competitive priorities etc.). This concept is presented in Fig. 3. Through this systemic inquiry, stakeholders' desires can be translated into full range of potential business opportunities, and desired properties and boundaries of enterprise system. To validate the content of a proposed ideal, interactive planning (IP) suggests two means: 1) inter-subjectivity and consensus among stakeholders with regard to a proposal, and 2) experimental feasibility [52]. Creating a vision is a key requirement for successful organizational change, consistent with Ackoff's notion of an idealized design. Re-design and transformation are desired consequences of a thoughtfully synthesized vision. External vision (insights, beliefs, and assumptions) is based on interpretation of all the factors in the context that have an impact on firm's operations and performance. Context can be viewed as a field of related parts [53], constructed for an issue, say SM. Descriptions of contextual pattern within this socio-ecological field form the present context. In regards to future, scenarios are chosen that describe alternative, comparably plausible lines of changing overall field conditions. Although there is no way of designating the one future that is most likely to emerge [54], scenarios can facilitate cognition of multiple possible futures in concrete detail and guide the context specific design of new business concepts or innovating around existing ones. Business environment is a sort of design space, a context in which to conceive or adapt the business. In fast-paced competitive environments, assumptions about how market forces, industry forces and emerging trends unfold give us the design space to develop potential business options [55]. The boradened scope of activities in SM, and large number of actors along the value chain necessitate all concerned to have a voice in decisions pertaining to products, technology, market, and enterprise system at large. Guided by System's methods relevant to stakeholders' engagement, and Ackoff's interactive planning (IP) and his view of the organization as social systemic in nature, the process of co-creating the ideal intent can be made more valuable. Firm's stakeholders and their participation are centric to this process. We need to understand stakeholders' needs and interests, and pregressively combine them. By establishing a generative context for the engagement of diverse perspectives of multiple stakeholders, a synergistic proposition can be evolved that is meaningful in a context of use [56]. An important insight is that no participative process can include every possible perspective, i.e. comprehension is impossible. In an interconnected world, knowing about any situation has limits or boundaries [57]. In relation to system boundary, we need to determine what is in our control, who (people and issues) is to include, how to justify exclusions, and how to address marginalization. 3.2 Strategy An enterprise's quest towards its ultimate ends is corresponding with the capabilities possessed by its ends-seeking system. Strategic intent of the firm is sought mainly to convey the purpose of the enterprise and how it is expected to perform in relation to those ends. Strategic objectives set the direction for designing business solutions and shaping operational structure [58]. Policies govern the selection of means and instruments by which the objectives are to be pursued, and principles are formulations of values to be preserved in such selections [48]. Principles provide guidance in making strategically consistent choices from the perspective of sustainability. Capabilities, derived from strategic principles, are a core concept to communicate and realize strategic intent. Building of capabilities has emerged as a new approach under strategic planning [59] to maintain competitive advantage and attain long-term business sustainability in a turbulent environment. A firm needs specific capabilities to enhance its capacity to create sustainable value and adapt to changing circumstances faster than competitors. The strategic planning process, using ideal intent as as a set of inputs, yields strategic intent manifested in an agreed upon business concept. This process is presented in Fig. 4. Capabilities identified in this process are built alongside to enable meaningful execution of this business concept. Capabilities are superior business processes, represented as a set of knowledge, skills, and competences [60]. The capability of organization as a whole is a high-level routine (or collection of routines) which combines tasks from different business functions into a particular output or solution of significant importance to the firm [61]. Capabilities not only provide a platform for strategy but also impact growth and profitability which in turn provide means for further strengthing of them. For the manufacturing sector this view of strategy has important bearings on organization's skill base, and dynamic decision making especially on infrastructural aspects [59]. The principles, policies, Vision Opportunities Proposed (products & system) Strategic planning Transformation Sustainability— Capabilities Competitiveness Business concept Objectives Fig. 4 Evolving strategic intent Advances in Production Engineering & Management 13(3) 2018 and systems of actions, underlying the infrastructural choices (Organization, manufacturing planning and control, quality, new product introduction and human resources), govern and reinforce the capabilities and resources that affect a firm's operations. Strategic importance of building capabilities is being realized beyond mere competitive reasons due to a shift in strategy focus from decreasing costs to increasing added value [62]. To embrace sustainability as a new business opportunity in manufacturing, a firm's key management function is to bridge the gap between strategic intent and required capabilities. Capabilities need to be developed in all of the important business dimensions (product, market, technology etc.) to capture entire range of opportunities. The future of manufacturing industry is characterized by capabilities in skills development, CE, technology and innovation, supply chain integration and digitization. Pre-production (R&D, product design) capabilities include those to de-velope innovative product concepts and its supporting technology. Within production segment, overall capabilities can be added to by undertaking fabrication of increasingly complex and innovative products, flexibile operations and optimized process planning, and improving resource efficiency through new fabrication methods. In post-production segment, capabilities are required in terms of efficiency and reliability to perform marketing, distribution, and EoL activities [35]. Both the individuals' learning capacity (competence) and capabilities of the organization are central to innovating the business in more sustainable ways and accommodate new needs of sustainability. Table 1 illustrates these capabilities along with salient requirements to attain them. Table 1 Capabilities for sustainable manufacturing Capability Salient requirements Integrating sustainability into the core of business Identification and evaluation of opportunities Adopt a business perspective, and transform the traditional business mindset. Innovation in business to turn environmental and social concerns to competitive advantage without incurring increased costs. Articulate and maintain a holistic view of sustainable value aimed at proposing a balanced value and distribution of economic costs and benefits among all actors. Judicious commitment of financial resources to investment opportunities. Create differentiated and innovative business models and enhance capacity to create, adjust and replace business models._ Gathering and filtering of information from both inside and out side the enterprise, and from core as well as the periphery of relevant business environment. Comprehending the likely evolution of industry, technology, and marketplace. Establish search and discovery methods and procedures, and allocate resources. Enhance interpretive skills of individuals and learning capacity of the organization. Aassessment of relevant opportunities through decision rules and informed conjectures._ Education and training focussed mainly on internalization of sustainability concepts, creative and innovation abilities, functional capabilities and life-cycle-thinking. Building of cross-domain skills (technology, engineering, electronics, robotics, computational and computer sciences) to perform manufacturing functions in industries of future. Cost-effective, innovative and exploratory teaching methods e.g. in-person training and coaching, experiential environments and digital methods. Adopting a stuctured approach for capability alignment and assessment, i.e. tools, methods, procedures, quantitative targets, credible metrics etc._ Building of individual skills and competence Culture adaptation, transition processes, and strategic direction to appreciate difficulties. Technological capabilities for an improved reverse logistic set-up aimed at high value recovery. Strategic use of ICT to establish an efficient collaboration and knowledge sharing. Capitalizing on information value to accelerate innovation in design, and identify more opportunities for additional value creation across the chain. Product life cycle planning with a clear focus on product reuse/re-manufacturing. Understanding materials (formulation, quality, innovation etc.), and a broader integration of material science into the business case and across value chain. Operationaliza-tion of CE Effective management of design process within a lifecycle thinking framework. Capture and reuse the knowledge acquired throughout the lifecycle to modify existing designs. Understanding the criticality of decisions at early stages of design in view of lifecycle costs. Use fewer materials, renewable materials, non-hazardous materials, and substitute materials. Evaluate materials against chosen criteria e.g. service life, recoverability, separateability etc. Design products that consume lesser resources during production, packing, distribution, recovery, reuse, re-manufacture and recycling. Evaluate competing product designs for price and ecological impact. Design products that are re-usable, re-manufacturable, recyclable, and bio-degradable._ Product design & development Table 1 (Continuation) Evaluate reusability, and remanufacturability of design alternatives via pre-determined criteria e.g. technological and economic feasibility, and environmental benefits. LCA, cost-benefit analysis and financial modelling to reveal potential designs for incorporation. Incorporate design for sustainability (DfX) and other strategies e.g. design for upgrade, design for disassembly, design for modularity, and design for reliability etc._ Production and delivery Innovation and flexibility in production and process technology. Choose and develop production processes that are efficient, resource conserving, processing of minimum material, and less polluting. Understand technology, energy and material as three interrelated components, vital to improving resource efficiency. Effective deployment and use of innovative technologies (e.g. ICT) for process control, energy saving, production efficiency, and safety etc. Efficient transportation and logistics to reduce environmental impact and cost of distribution._ End of life (EoL) management Establishment of reverse logistics activities, and their alignment with after market and remarketting businesses. User friendly and cost efficient collection and segregation depending on product quality, and effective segmentation to support dynamic recovery and economics of value capture. Capitalize on knowledge acquired during lifecyle to facilitate EoL processes. Devise processes that enable product re-usability with minimum possible resources to gain best environmental and economic advantages. Leverage external resources to increase the throughput and performance of EoL processes e.g. outsourcing the core acquisition, integrating the supply chain with aftermarket divisions, and accomplishing economies of scale in re-manufacturing by sharing resources with others._ Standards Establish norms and rules that lead to a significant reduction of activities and transaction costs. Strategic use of standards to add to the firm's capabilities and manufacturing base, re-engineering of the technology and processes, and redesign of products. Adoption of international standards towards a faster knowledge acquisition and technological transition. Establishing uniform standards along the value chain for standardization of product-related information, complete life cycle assessment of products, sustainability performance evaluation of all upstream and downstream partners, and optimization of end product sustainability. Efficient collaboration and data exchange, distribution of right information in a fast and reliable way, and holistic management of complex tasks through the use of communication standards._ Allignment and collaboration Extending capabilities in terms of exploiting additional resources through broad-based external search and close interaction and collaboration with other enterprises, entities and institutions in the business environment and society. Strategic involvement of customers in product design and development with the aim to include their perspective in value proposition, discover new opportunities and business models, develop more successful and radical value creation solutions, reduce product costs, anticipate the potential for applying new technologies, and retain customers longer. Extending capabilities in terms of exploiting additional resources or shifting certain activities (e.g. procurement, product design etc.). Establishing strategic partnering with supplier to capture benefits of economies of scale by outsourcing to contract manufacturers and large suppliers, assimiliating the capability of suppliers in terms of those capabilities that are closely linked and mutually reinforcing, capitalizing on rapid innovations by component suppliers ahead of the competition, improving supplier operations, and optimizing supplier-manufacturer interface to reduce material costs and gain a better insight on environmental impact reduction opportunities along the supply chain. Leveraging benefits of co-specialization and complementary innovations towards differentiated product offerings and unique cost savings, improved performance, sustainable competitive advantage, and enhanced entrepreneurial capacities._ 3.3 Architecture The architecture serves the purpose of realizing an agreed upon business concept, formulated as part of firm's strategic intent It links business strategy with business roles and processes, and structures the responsibility over business activities along the value chain [46]. It determines effectiveness in attaining firm's objectives keeping in view their strategic impact. From a functional perspective, it manages interaction among structure, processes and activities, and resources and capabilities to deliver sustainable value to all stakeholders. This concept is presented in Fig. 5. Within its functional framework, it utilizes a common and preferably customized knowledge base that allows manipulation of information about different business aspects from a variety of perspectives. • Business concept Objectives • Economic • Environmental • Social Fig. 5 Architecture framework for sustainable value creation In system's view, architecture design involves a holistic understanding of how function, possible structures and processes should interact and communicate with each other to produce desired outputs, within a larger context. The planning, learning, and control function is thus an integral element here, responsible to evaluate the effectiveness of organization as a whole [10]. In functional design process, a comprehensive review of the strategic intent is conducted first to understand business priorities, relationship among objectives and transformation requirements. BA role is then determined in relation to business needs including required capabilities and core activities. The detailed business analysis next establishes a clear view of the business in terms of its organization and value chain structure. The final solution dictates initiatives (e.g. investment analysis, new product rollout, capability outsourcing etc.) and measurement criteria relevant to set objectives. This process, however, is ongoing and continues to redeem the once finalzed structure, processes, activities and even the business concept. Analogous to Checkland's [63] abstract concept of 'system as an adaptive whole', an architecture responds to the changing needs and requirements, interactively manages the value creation factors, and establishes a fit between problem space and proposed solution. The desirable emergent properties characterizing such a system could be e.g. a unique and differentiated business value or improved sustaina-bility performance etc. Environmental and social concerns today, changing market needs, technological innovations, and capabilities based competition are driving manufacturing businesses towards sustainability. Structural flexibility on the other hand must complement the prevalent dimension of the competition; based on either of new product offerings, market dominance or technological leadership, towards an efficient and conflict free internal environment. A multidimensional mode of structure - capable of producing different outcomes in the same or different environments - can interactively manage the interdependent business dimensions thus eliminating the need for a change of emphasis from one business orientation to another [10, 48]. [45] Strategic intent is interdependent with firm's operational structure which establishes relations between people, process, technology, managerial aspects, and information. The operational structure must ensure that wealth creation, cost reduction, and environmental and social perfromance are improved simultaneously to support the sustainability strategy. Processes are primarily concerned with how the actual output(s) is to be produced through an integrated chain of activities performed by different groups in an organization, with a strong interface among them. Processes are technologically driven and subject to continuous change and improvement. A firm must comprehend the relevance of technologies to emerging competitive challenges, understand the flow and interface between system elements, and acquire operational knowledge of processes. In addition, these technological processes should maintain compatibility with already in place organizational processes that are concerned with creating integration, alignment and synergy among the parts of an organization [10]. A firm's core competence (unique resources and strengths) drives its core processes. In case of SM, the core processes in the chain include design and development of new products, production and delivery of products, and EOL management. Salient requirements in terms of capability of these processes have been illustrated in Table 1. Each process must be designed to achieve its priorities. Operational effectiveness of these processes may be assessed with key metrics corresponding to those key activities in the value chain that influence the critical dimensions of sustainability performance. In addition to key activities, the firm needs to perform other essential activities e.g. to train its workforce, to develop new capabilities, and to link new activities to the existing activity system [64]. The nature of products and services, and the activities performed to deliver them collectively determine the structure of value chain [46]. The choices on content, structure and governance of activities shape the architecture in terms of where these are to be performed and what kind of resources are required to perform them? The entire collection of activities can be managed conveniently by aggregating them at different levels (major functions at top-level and specific activities at lowest level) [64], or segregating them according to their core, supporting or peripheral nature. 4. Discussion The proposed conceptual approach rests on some key tenets, elaborated in section 3, which collectively could enable long-term sustainability (ability to transform and maintain competitiveness) for manufacturing businesses; among them greater capacitation of businesses to contribute to sustainable development, understanding stakeholders' needs and intersts, a coevolved vision and business proposition, enhanced capacity to discover opportunities systematically, learning capacity, competence and adptability of individuals and organization, building of context specific capabilities in all core segments of the business and their strategic use towards creating added value, innovative solutions and faster responce to changing market needs and technological challenges, and holistic understanding of problems and solutions in a larger context to improve economic, environmental and social perfromance simultaneously . The existing social environment and level of economic development of a country could be relevant/ significant in view of implementation of proposed concept across its entire manufacturing sector. Specific strategies and policies that consider both economic growth and social development are vital to nurture an environment for sustainable business development. Leading a competitive advantage demands educated and competent manpower. Quality of education has a direct impact on skill levels of individuals. Equal opportunities for education, training, employment and health increase productivity and growth. Exploration and capturing of new business opportunities is a knowledge intensive and enterpreneurial endeavour. National policies for matching skills with business needs could add to the ability of businesses to transform and create new value. Strong financial means in a country and better access to those, along with other supporting conditions, may encourage entrepreneurship and pursuits of competitiveness and innovativeness. National spending on research and development could pave the way for strong collaborative networks in a society which foster co-creation. The proposed conceptual framework serves as a guideline in enabling continuous transformation and competitiveness of a manufacturing enterprise. However, it is not an implementation methodology or a procedure in itself that explicilty defines any preference as to how to get there. It will be instantiated using a chosen or developed methodology; conceptual frameworks, generally, do not include or imply any process model for implementation. Nevertheless, useful insights have been drawn from relevant research in reference to enterprise. First, formalization or standardization of proposed framework's elements, in context of specific manufacturing environment, is critical for its consistency and common understanding across management levels and common people. An explicit management procedure may be established to enable systematic implementation and execution through relevant tools and techniques. Potential challenges and any unforeseen conditions or obstacles, posed by the dynamics of the external or internal environment, may be identified. Building a robust culture of sustainability is vital to successful implementation of proposed concept. Sustainability education and skills in the domain of change management and organizational behaviour are critical to deployment process. Top level management must lead the implementation process to elevate its status. Employees who understand day-to-day business routines may be encouraged to get them involved in generating innovative ideas. Implementation objectives should be defined. Moreover, it needs to be assured that specified initiatives and activities contribute to transformation goals throughout the course of im- plementation. Communicating to the employees and other stakeholders that why the company has chosen this continuous transformation path, i.e. "ideal—strategy—architecture" is essential. Pertinent here is to frame and relate stakeholders' apprehensions in an understandable form. Company wide communiqué, at the bare minimum, may include business needs, implementation issues and requirements, and expected benefits of such a transformative change or redesign. Despite above recommendations, an average manufacturing environment may yet face multiple challenges on the way to become a sustainable competitive firm in accordance with transformation path proposed in this research. Capturing opportunities for sustainable value creation and building and evaluation of capabilities (Table 1) are no less than a paradigm shift for such environments. The proposed concept is centred on sustainability being a core component of the business (strategy) along with its continual renewal to gain long-term success while average manufacturers most often target short-term economic gains. Nevertheless, a deliberate support from government and society, in various aspects, could play a significant part to bring them at par with proposed requirements through a systematic incremental approach. Identification and evaluation of opportunities (idealization stage) could be a knowledge intensive and costly affair posing resource constraints. This process may be facilitated through a network of manufacturers, researchers and domain experts thereby helping traditional manufacturers seek alternate products and new markets to create competitive advantage. Experts in cultural assessment and change management can assist in developing and promoting a company culture needed for integrating sustainability into the everyday business. Government can sponsor participation mechanisms during idealization and strategic planning stage that could enable businesses to create new value innovatively with existing resources instead of huge investments in new technology. Evaluation of resource efficiency and lifecycle impacts needs specialized skills wherein provisions can be made for free of cost availability of such services. All such initiatives, at significantly low cost to businesses, could motivate them towards implementation of proposed concept. 5. Conclusion Despite enormous research on SM, rarely the studies have focussed on continuous business transformation according to stakeholders' desires while maintaining competitiveness in a changing context. This research tried to fill this gap by first recognizing the need for a holistic and systemic representation of the issue along with a potentially viable approach embedded with capability to integrate sustainability concepts in manufacturing. A comprehensive overview of the issue, contributed by numerous aspects, provided an overall design context. A holistic and integrated process of planning and realization was aspired in relation to sustainability goals and a systematic transformation. Ultimately, framework was designed to provide a simplistic and progressive approach that could produce consistently the sustainable outcomes, and meet proposed ends. The three elements that jointly constitute this approach, i.e. 'Ideal', 'Strategy', and 'Architecture', have been presented from their functional perspective and contextual relation with each other and larger context A systemic platform is thus established by constructing and operationalizing the individual elements. The dependency and relationship of elements with each other drives the conclusion that when executed in isolation, each element could still create value but limited to its context of use and initiating choices.The holistic view can maximize the cumulative sustainable value by managing the interaction among its elements, and maintaining the traceability from proposed ideal to strategic intent, and finally to architecture which has implications for functional effectiveness of the enterprise. Transformation can be more effectively achieved and repeated by accomplishing what is necessary and in concert with stakehold-ers'desires within each of the three elements. The enhanced capacity of enterprise to satisfy these requirements is the key to attain proposed sustainability ends. Continuous identification and re-newal of capabilities build this capacity to exploit opportunities for sustainable value creation. Embracing this conceptual framework and its systematic operationalization will drive continuous transformation of a SM business in a fast paced environment. It also seems fair to conclude that proposed framework can play a critical role in maintaining competitiveness since it suggests several capabilities whereby challenges can be addressed through meeting salient requirements underlying each capability. However, the framework is not yet embedded with methodology, procedures or tools for implementation, and there is room for further comprehension in terms of implementation strategy. Further research is needed in actual settings to confirm whether it can achieve the claimed ojectives in diverse manufacturing environments. The experience so gained will demonstrate its conceptualized role of integrating sustainability requirements and addressing challenges, pertinent to manufacturing businesses. References [1] Duflou, J.R., Sutherland, J.W., Dornfeld, D., Herrmann, C., Jeswiet, J., Kara, S., Hauschild, M., Kellens, K. (2012). Towards energy and resource efficient manufacturing: A processes and systems approach, C1RP Annals, Vol. 61, No. 2, 587-609, doi: 10.1016/j.cirp.2012.05.002. [2] Bi, Z. (2011). Revisiting system paradigms from the viewpoint of manufacturing sustainability, Sustainability, Vol. 3, No. 9, 1323-1340, doi: 10.3390/su3091323. [3] Garetti, M., Taisch, M. (2012). Sustainable manufacturing: Trends and research challenges, Production Planning & Control, Vol. 23, No. 2-3, 83-104, doi: 10.1080/09537287.2011.591619. [4] Gopalakrishnan, K., Yusuf, Y.Y., Musa, A., Abubakar, T., Ambursa, H.M. (2012). Sustainable supply chain management: A case study of British aerospace (BAE) systems, International Journal of Production Economics, Vol. 140, No. 1, 193-203, doi: 10.1016/j.ijpe.2012.01.003. [5] ElMaraghy, H., Schuh, G., ElMaraghy, W., Piller, F., Schonsleben, P., Tseng, M., Bernard, A. (2013). Product variety management, C1RP Annals, Vol. 62, No. 2, 629-652, doi: 10.1016/j.cirp.2013.05.007. [6] Jawahir, I.S., Bradley, R. (2016). Technological elements of circular economy and the principles of 6R-based closed-loop material flow in sustainable manufacturing, Procedia C1RP, Vol. 40, 103-108, doi: 10.1016/j.procir. 2016.01.067. [7] Zhang, H., Calvo-Amodio, J., Haapala, K.R. (2015). Establishing foundational concepts for sustainable manufacturing systems assessment through systems thinking, International Journal of Strategic Engineering Asset Management, Vol. 2, No. 3, 249-269, doi: 10.1504/IISEAM.2015.072124. [8] Lonngren, J., Svanstrom, M. (2016), Systems thinking for dealing with wicked sustainability problems: Beyond functionalist approaches, In: Leal Filho, W., Nesbit, S. (eds.), New developments in engineering education for sustainable development, Springer, Cham, Switzerland, 151-160, doi: 10.1007/978-3-319-32933-8 14. [9] Jackson, M.C. (2006). Creative holism: A critical systems approach to complex problem situations, Systems Research and Behavioral Science, Vol. 23, No. 5, 647-657, doi: 10.1002/sres.799. [10] Gharajedaghi, J. (2011), Systems thinking: Managing chaos and complexity: A platform for designing business architecture, Third edition, Elsevier, Burlington, USA. [11] Gaziulusoy, A.I. (2015). A critical review of approaches available for design and innovation teams through the perspective of sustainability science and system innovation theories, Journal of Cleaner Production, Vol. 107, 366-377, doi: 10.1016/j.jclepro.2015.01.012. [12] World Economic Forum (prepared in collaboration with Ellen MacArthur Foundation and McKinsey & Company). Towards the circular economy: Accelerating the scale-up across global supply chains, from http://www3.weforum.org/docs/WEF ENV TowardsCircularEconomy Report 2014.pdf.accessed August 23, 2016. [13] Teece, D.J. (2007). Explicating dynamic capabilities: The nature and microfoundations of (sustainable) enterprise performance, Strategic Management Journal, Vol. 28, No. 13, 1319-1350, doi: 10.1002/smj.64. [14] Oertwig, N., Galeitzke, M., Schmieg, H.-G., Kohl, H., Jochem, R., Orth, R., Knothe, T. (2017). Integration of sustainability into the corporate strategy, In: Stark, R., Seliger, G., Bonvoisin, J. (eds.), Sustainable manufacturing: Challenges, solutions and implementation perspectives, Springer, Cham, Switzerland, doi: 10.1007/978-3-31948514-0. [15] Tonelli, F., Evans, S., Taticchi, P. (2013). Industrial sustainability: Challenges, perspectives, actions, International Journal of Business Innovation and Research, Vol. 7, No. 2, 143-163, doi: 10.1504/IjBiR.2013.052576. [16] Koho, M., Tapaninaho, M., Heilala, J., Torvinen, S. (2015). Towards a concept for realizing sustainability in the manufacturing industry, Journal of Industrial and Production Engineering, Vol. 32, No. 1, 12-22, doi: 10.1080/ 21681015.2014.1000402. [17] Subic, A., Shabani, B., Hedayati, M., Crossin, E. (2012). Capability framework for sustainable manufacturing of sports apparel and footwear, Sustainability, Vol. 4, No. 9, 2127-2145, doi: 10.3390/su4092127. [18] Masood, I., Jahanzaib, M., Haider, A. (2016). Tool wear and cost evaluation of face milling grade 5 titanium alloy for sustainable machining, Advances in Production Engineering & Management, Vol. 11, No. 3, doi: 10.14743/ apem2016.3.224. [19] Hassine, H., Barkallah, M., Bellacicco, A., Louati, J., Riviere, A., Haddar, M. (2015). Multi objective optimization for sustainable manufacturing, application in turning, International Journal of Simulation Modelling, Vol. 14, No. 1, 98-109, doi: 10.2507/IJSIMM14( 1)9.292. [20] Govindan, K., Soleimani, H., Kannan, D. (2015). Reverse logistics and closed-loop supply chain: A comprehensive review to explore the future, European Journal of Operational Research, Vol. 240, No. 3, 603-626, doi: 10.1016/ j.ejor.2014.07.012. [21] Schenkel, M., Caniels, M.C.J., Krikke, H., van der Laan, E. (2015). Understanding value creation in closed loop supply chains - Past findings and future directions, Journal of Manufacturing Systems, Vol. 37, Part 3, 729-745, doi: 10.1016/i.imsv.2015.04.009. [22] Ageron, B., Gunasekaran, A., Spalanzani, A. (2012). Sustainable supply management: An empirical study, International Journal of Production Economics, Vol. 140, No. 1, 168-182, doi: 10.1016/i.iipe.2011.04.007. [23] Galal, N.M., El-Kilany, K.S. (2016). Sustainable agri-food supply chain with uncertain demand and lead time, International Journal of Simulation Modelling, Vol. 15, No. 3, 485-496, doi: 10.2507/IISIMM15(3)8.350. [24] Fiksel, J. (2003). Designing resilient, sustainable systems, Environmental Science & Technology, Vol. 37, No. 23, 5330-5339, doi: 10.1021/es0344819. [25] Fiksel, J. (2012). A systems view of sustainability: The triple value model, Environmental Development, Vol. 2, 138-141, doi: 10.1016/i.envdev.2012.03.015. [26] Moldavska, A. (2016). Model-based sustainability assessment - An enabler for transition to sustainable manufacturing, Procedia CIRP, Vol. 48, 413-418, doi: 10.1016/i.procir.2016.04.059. [27] Moldavska, A., Welo, T. (2016). Development of manufacturing sustainability assessment using systems thinking, Sustainability, Vol. 8, No. 1, doi: 10.3390/su8010005. [28] Terzi, S., Bouras, A., Dutta, D., Garetti, M., Kiritsis, D. (2010). Product lifecycle management - From its history to its new role, International Journal of Product Lifecycle Management, Vol. 4, No. 4, 360-389, doi: 10.1504/IIPLM. 2010.036489. [29] Alayón, C., Safsten, K., Johansson, G. (2016). Conceptual sustainable production principles in practice: Do they reflect what companies do?, Journal of Cleaner Production, Vol. 141, 693-701, doi: 10.1016/i.iclepro.2016.09.079. [30] Esmaeilian, B., Behdad, S., Wang, B. (2016). The evolution and future of manufacturing: A review, Journal of Manufacturing Systems, Vol. 39, 79-100, doi: 10.1016/i.imsy.2016.03.001. [31] Haas, W., Krausmann, F., Wiedenhofer, D., Heinz, M. (2015). How circular is the global economy?: An assessment of material flows, waste production, and recycling in the European union and the world in 2005, Journal of Industrial Ecology, Vol. 19, No. 5, 765-777, doi: 10.1111/iiec.12244. [32] Ellen MacArthur Foundation. Towards a circular economy: Business rationale for an accelerated transition, from https://www.ellenmacarthurfoundation.org/assets/downloads/TCE Ellen-MacArthur-Foundation 9-Dec-2015.pdf, accessed October 7, 2016. [33] McDonough, W., Braungart, M., Anastas, P.T., Zimmerman, J.B. (2003). Peer reviewed: Applying the principles of green engineering to cradle-to-cradle design, Environmental Science & Technology, Vol. 37, No. 23, 434-441, doi: 10.1021/es0326322. [34] McDonough, W., Braungart, M. (2002). Design for the triple top line: New tools for sustainable commerce, Corporate Environmental Strategy, Vol. 9, No. 3, 251-258, doi: 10.1016/S1066-7938(02)00069-6. [35] World Economic Forum. Manufacturing our future: Cases on the future of manufacturing, from http://www3. weforum.org/docs/GAC16 The Future of Manufacturing report.pdf. accessed May 2, 2017. [36] Stock, T., Seliger, G. (2016). Opportunities of sustainable manufacturing in industry 4.0, Procedia CIRP, Vol. 40, 536-541, doi: 10.1016/i.procir.2016.01.129. [37] Miroshnychenko, I., Barontini, R., Testa, F. (2017). Green practices and financial performance: A global outlook, Journal of Cleaner Production, Vol. 147, 340-351, doi: 10.1016/i.iclepro.2017.01.058. [38] Abdul Rashid, S.H., Evans, S., Longhurst, P. (2008). A comparison of four sustainable manufacturing strategies, International Journal of Sustainable Engineering, Vol. 1, No. 3, 214-229, doi: 10.1080/19397030802513836. [39] Chatha, K., Butt, I. (2015). Themes of study in manufacturing strategy literature, International Journal of Operations & Production Management, Vol. 35, No. 4, 604-698, doi: 10.1108/IJ0PM-07-2013-0328. [40] Kappelman, L.A., Zachman, J.A. (2013). The enterprise and its architecture: Ontology & challenges, Journal of Computer Information Systems, Vol. 53, No. 4, 87-95, doi: 10.1080/08874417.2013.11645654. [41] van Iwaarden, J., van der Wiele, T. (2012). The effects of increasing product variety and shortening product life cycles on the use of quality management systems, International Journal of Quality & Reliability Management, Vol. 29, No. 5, 470-500, doi: 10.1108/02656711211230481. [42] Halog, A., Manik, Y. (2011). Advancing integrated systems modelling framework for life cycle sustainability assessment, Sustainability, Vol. 3, No. 2, 469-499, doi: 10.3390/su3020469. [43] Friedrich-Schiller-University Jena. Roadmap for sustainability assessment in european process industries, Technical report, from https://www.researchgate.net/publication/301823430, accessed 13 May 2016. [44] Zhang, H., Haapala, K.R. (2015). Integrating sustainable manufacturing assessment into decision making for a production work cell, Journal of Cleaner Production, Vol. 105, 52-63, doi: 10.1016/i.iclepro.2014.01.038. [45] Hendrickx, H.H.M. (2015). Business architect: A critical role in enterprise transformation, Journal of Enterprise Transformation, Vol. 5, No. 1, 1-29, doi: 10.1080/19488289.2014.893933. [46] University of Cambridge. State-of-practice in business modelling and value-networks, emphasising potential future models that could deliver sustainable value, from http://www.sustainvalue.eu/publications/D2 1 Final Rev1 0 web.pdf, accessed November 10, 2016. [47] Grassl, W. (2012). Business models of social enterprise: A design approach to hybridity, ACRN Journal of Entre-preneurship Perspectives, Vol. 1, No. 1, 37-60. [48] Ackoff, R.L. (1999), Re-creating the corporation: A design of organizations for the 21st century, Oxford University Press, USA. [49] Mathews, F. (2011). Towards a deeper philosophy of biomimicry, Organization & Environment, Vol. 24, No. 4, 364-387, doi: 10.1177/1086026611425689. [50] Apple Inc. Environmental responsibility report, from https://images.apple.com/environment/pdf/Apple Environ mental Responsibility Report 2017.pdf accessed May 15, 2017. [51] Ocampo, L.A., Clark, E.E., Tanudtanud, K.V.G., Ocampo, C.O.V., Impas Sr, C.G., Vergara, V.G., Pastoril, J., Tordillo, J.A.S. (2015). An integrated sustainable manufacturing strategy framework using fuzzy analytic network process, Advances in Production Engineering & Management, Vol. 10, No. 3, 125-139, doi: 10.14743/apem2015.3.197. [52] Haftor, D.M. (2011). An evaluation of R.L. Ackoff s interactive planning: A case-based approach, Systemic Practice and Action Research, Vol. 24, No. 4, 355-377, doi: 10.1007/s11213-010-9188-y. [53] Rhyne, R.F. (1994). Contextual discipline: Its essentiality within social-systems analysis, Technological Forecasting and Social Change, Vol. 47, No. 3, 277-292, doi: 10.1016/0040-1625(94)90069-8. [54] Amer, M., Daim, T.U., Jetter, A. (2013). A review of scenario planning, Futures, Vol. 46, 23-40, doi: 10.1016/j. futures.2012.10.003. [55] Osterwalder, A., Pigneur, Y. (2010), Business model generation: A handbook for visionaries, game changers, and challengers, John Wiley & Sons, Hoboken, New Jersey, USA. [56] Midgley, G. How can co-creation produce better outcomes for people and communities?, from https://www. researchgate.net/publication/303752603, accessed July 24, 2016, doi: 10.13140/RG.2.1.2373.3360. [57] Midgley, G. (2000), Systemic intervention: Philosophy,methodology and practice, Springer Science, New York, USA, doi: 10.1007/978-1-4615-4201-8. [58] Ackoff, R.L. (1990). Strategy, Systems Practice, Vol. 3, No. 6, 521-524, doi: 10.1007/BF01059636. [59] Hayes, R.H., Pisano, G.P. (1996). Manufacturing strategy: At the intersection of two paradigm shifts, Production and Operations Management, Vol. 5, No. 1, 25-41, doi: 10.1111/j.1937-5956.1996.tb00383.x. [60] Bilge, P., Seliger, G., Badurdeen, F., Jawahir, I.S. (2016). A novel framework for achieving sustainable value creation through industrial engineering principles, Procedia C1RP, Vol. 40, 516-523, doi: 10.1016/j.procir. 2016.01.126. [61] Mills, J., Platts, K., Bourne, M. (2003). Competence and resource architectures, International Journal of Operations & Production Management, Vol. 23, No. 9, 977-994, doi: 10.1108/01443570310491738. [62] World Economic Forum. The future of manufacturing: Driving capabilities, enabling investments, from http://www3.weforum.org/docs/Media/GAC14/Future of Manufacturing Driving Capabilities.pdf,accessed May 2, 2017. [63] Checkland, P. (2012). Four conditions for serious systems thinking and action, Systems Research and Behavioral Science, Vol. 29, No. 5, 465-469, doi: 10.1002/sres.2158. [64] Zott, C., Amit, R. (2010). Business model design: An activity system perspective, Long Range Planning, Vol. 43, No. 2-3, 216-226, doi: 10.1016/j.lrp.2009.07.004. APEM jowatal Advances in Production Engineering & Management Volume 13 | Number 3 | September 2018 | pp 254-266 https://doi.Org/10.14743/apem2018.3.288 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper A hybrid grey cuckoo search algorithm for job-shop scheduling problems under fuzzy conditions Yang, F. , Ye, C.M.a, Shi, M.H.a aBusiness School, University of Shanghai for Science and Technology, Shanghai, P.R. China bCollege of Management, Henan University of Traditional Chinese Medicine, Zhengzhou, Henan, P.R. China A B S T R A C T A R T I C L E I N F O This paper aims to acquire the precise makespan or delivery period in jobshop scheduling (JSP) under fuzzy conditions. To this end, the author designed a grey scheduling model and a hybrid grey cuckoo search (HGCS) algorithm in the following steps. Firstly, three- and four-parameter interval grey numbers were introduced to depict the fuzzy makespan and delivery period, respectively; then, the possibility measure and necessity measure were defined, and the tardiness credibility index was proposed to estimate the probability of job tardiness. After that, a grey mixed integer programming model was developed to minimize the mean tardiness credibility, and the HGCS was proposed to solve the model. Finally, simulations were conducted on the classical example of 6(3) x 6. The results show that the proposed algorithm outperformed the basic cuckoo search. The research findings shed new light on the JSP under fuzzy conditions. © 2018 CPE, University of Maribor. All rights reserved. Keywords: Job-shop scheduling problem (JSP); Grey scheduling; Fuzzy condition; Cuckoo search (CS); Credibility; Possibility measure; Necessity measure *Corresponding author: yangfeng1126@126.com (Yang, F.) Article history: Received 22 May 2018 Revised 29 July 2018 Accepted 22 August 2018 1. Introduction Job-shop scheduling (JSP) is a common NP-hard scheduling problem in modern industry. The problem is about the processing of a number of jobs on several machines. Each of these jobs needs to go through multiple processes and can run on several machines. The constraints of the JSP are as follows: the jobs must be treated in a feasible sequence of processes and follow the same order of priority; the same job should be processed only once on the same machine; the same machine can only process one job at a time. Under the above constraints, the JSP aims to optimize performance indices like the minimum makespan or the minimum tardiness through rational arrangement of the process sequence and start time. Many algorithms have been developed to solve classical JSPs. For instance, Zhang et al. [1] proposed a hybrid meta-heuristics algorithm based on shuffled frog leaping algorithm (SFLA), intelligent water drops algorithm (IWD) and path relinking (PR) algorithm. Modrak et al. [2] converted an m-machine problem into a 2-machine problem, and created a Johnson algorithm to solve it. Chaudhry et al. [3] presented a genetic algorithm with process planning and scheduling. Considering sequence flexibility, this algorithm solves classical JSPs using the spreadsheet of domain independence. Huang et al. [4] suggested solving classical JSPs by sequence flexibility and minimized the makespan with an improved genetic algorithm. In classical JSPs, the production is scheduled under the ideal condition, that is, all factors are determined. In real-world production scheduling, however, both makespan and delivery period are highly indeterminate owing to uncertainties like machine breakdown, unfavourable environment or wrong decisions. To reflect the exact production situation, it is necessary to treat the makespan and delivery period as fuzzy variables. Recent years has seen much attention being paid to the JSP under fuzzy conditions. The most common solution is fuzzy scheduling, that is, describing the fuzzy makespan and delivery period with triangular and trapezoidal fuzzy numbers. For instance, Ishii et al. [5] proposed the concept of fuzzy delivery period after exploring various jobs, and then investigated single- and double-machine open-shop scheduling under fuzzy delivery period. Noori-Darvish et al. [6] created a two-objective fuzzy programming model for open-shop scheduling, whose sequence relies on the job preparation time, fuzzy makespan and fuzzy delivery period. Focusing on two-machine scheduling, Gharehgozli et al. [7] put forward a fuzzy mixed integer programming model for scheduling sequence dependent on the job preparation time and fuzzy delivery period. Simeunovic et al. [8] built up a workforce scheduling model based on artificial neural network. For triangular and trapezoidal fuzzy numbers, both sides are monotonous linear functions. Therefore, linear, uniform variations can be expected for the probability function of the fuzzy number value on both sides of the most likely value (i.e. the window of the most likely makespan and the most likely delivery period). Nevertheless, the probability functions of makespan and delivery period are not necessarily linear in actual production. The possible linear and nonlinear functions of the two indices are shown in Fig. 1. Considering the fuzziness of makespan and delivery period, a scheduling constraint should be expressed as an interval rather than a fixed value. Moreover, the ideal makespan and delivery period in the classical JSPs are often the most likely values in fuzzy situations; the probability functions for the deviation from the most likely values are not necessarily linear ones. These features are consistent with those of interval grey numbers. In grey system theory, an interval grey number has different probabilities in taking different values. This unique feature has been widely discussed in the academia. For example, reference [9] introduces the interval grey number into multi-criterion optimization, with the aim to solve various fuzzy decision-making problems in the real world. Considering the impact of fuzzy grey number on prediction results, reference [10] converts the interval grey number sequence into a real number sequence and establishes a prediction model of interval grey number. Based on the definition of the basic interval number, references [11-18] put forward the three-parameter interval number and a multi-attribute decision-making method. Among them, reference [11] discusses the grey target decision-making with known maximum probability of attribute values: the index weight was determined through subjective or objective weighting, a comprehensive optimization model was established for index weight, and the attribute values were identified by three-parameter interval numbers. In light of the above, this paper establishes a grey JSP model, which treats fuzzy makespan and delivery period as three- and four-parameter interval grey numbers, respectively, and defines the tardiness credibility index of jobs based on the possibility measure and necessity measure. Then, the fuzzy tardiness was determined as the optimization target and a grey mixed integer programming model was established. Finally, the proposed model was solved by an intelligent algorithm. Mik(t) Mi(rf) Mik(t) Mi(rf) t d t Z¡ d (a) (b) (c) (d) Fig. 1 (a) and (b): Linear probability functions, (c) and (d): Nonlinear probability functions 2. Grey job shop scheduling problem 2.1 Problem description In actual scheduling, the job processing is more or less affected by some fuzzy factors. Thus, the scheduling process should be measured by probabilities. Here, the fuzzy makespan and delivery period are expressed as three- and four-parameter interval grey numbers, and the scheduling target is to optimize the mean tardiness credibility of the job. To realize the goal, it is necessary to review the definition of grey number. Definition 1 [20] A number whose value falls in an interval rather than at a fixed point is known as the grey number and can be denoted as ®. A three-parameter interval grey number can be expressed as a(®) e[a",a*,a+] (0 < a~ a+ Definition 2 [19] A four-parameter interval grey number can be expressed as e [a~,a*",a*+,a+] (0 < a~ | -[(a*+)2 - (a+)2]"1, a*+ a+ According to the above description, the grey JSP can be described as: • Let P = { Pi,P2,—,Pn) be a set of n jobs, with pt being the i-th job (i = 1,2, ...,n). • Let M = { m1,m2,...,mm} b e a set of m machines, with m.j being the y'-th machine ( = 1,2, .„,m). • Let OP = { op1,op2,...,opn} be a set of process sequences, with opi = { opil,opi2,., opim] being the process sequence of job and opik being the machine number of the i-th job in the /i-th process. If opik = 0, the k-th process of the i-th job is not implemented (k = 1,2, ...,m). • Let T = {ail(0),ai2(0),., aim(0)} be the set of makespans of job pt, with aik() e\T[k,T*k,T*k~\, whereTj*^, T[kand Tik are respectively the most likely, the shortest and the longest makespans of the /c-thprocess. The function describes the possibility of completing /i-thprocess of the i-th job at time t. • Let = (d[r,d*~,d*+,d+) be a four-parameter interval grey number indicating the delivery period of the i-th job. The interval (d*",d*+) is the most likely delivery period of the i-th job. The other constraints are the same with those in the classical JSP. Because the delivery period can accurately reflect the market changes, the ultimate target of our scheduling is to find a feasible scheduling plan that minimizes the mean tardiness of the job. 2.2 Grey number operator Suppose all jobs to the processed on the /c-th machine are allocated into one sequence, in which the i-th job is preceded by the job sequence 1,2,...,/-1. Thus, the makespan of the i-th job equals the total makespan of all jobs in the preceding sequence. To obtain the total makespan, the addition of two three-interval grey numbers should be defined as follows [21]. Definition 3 Let aik{= [a",a*,a+] and aJk(= Then, we have: aik(®) + a]k(®) = [«" +P~,a* +/3*,a+ + £+] According to Definition 3, the grey makespan a^®) of the i-th job can be obtained as: «¿(0) = lq=iaqk(®), q = 1,2,...,n (3) Where aqk(is the grey makespan of the q-th job on the k-th machine. Before obtaining the job tardiness, it is also necessary to define the reduction operator of four-parameter interval grey numbers. Theorem 1 For the membership function (2) of four-parameter interval number, there exist x1 e[ca ,a* ], x2 e[a*+,a^"] and x3 = (x2 — x2)2 that satisfy: =Vdj(.X 2) (ai )2 -0/)2 -x. (a-)* -(a/)l -(a*-)2 + «) Proof: If mDí(xi) =Mdj(x 2) , then (ai ) [(a* > - (a¿ )' r X. «F -(<)' 1-J-l 2 Solving this equation, we have: x2 = ([(a/+)2 _(a/) (ai )2~xi 1 i 2 Hence, ■ (af)2-(a|)2-x; 1 1 2 3 (aí )2~x¡ H(a} 1 1 T+W)2 1 J - 1 1 1 Hence, there exist x1, x2 and x3 = (x2 — x2)2 that satisfy: 1 MD¡(xi) =Mdj(X 2) = (ai )2 "O/)2 -x. Q.E.D. Definition 4 Let ¿¿(®) = [«( ,at* ,a[*+,al;f] and bj(®) = [aj be two four-parame- ter interval grey numbers. Through reduction operation, another four-parameter interval grey number can be obtained as ¿¿(®) — 6/(®) = [«r — aj~,a?~ — o*~,a*+ — a*+— af]. Whereas makespan and delivery period are respectively three- and four-parameter interval grey numbers, the three-parameter interval grey number should be treated as a special type of four-parameter interval grey number, that is, aifc(g) e[aifc a*fc,a^]o aifc(g) e[aik a*ik,a*ik, a^]. By definition 4, the grey tardiness of the i-th job in the k-th machine can be obtained as Til®) = ai(g) -¿¡(g) = -a«+,arfc- -«r+,«ifc ~«r,aîk ~ai > Definition 5 Let ¿¿(g) = [a[r,a[*_,al*+,a^] and 6j(g) = [aj~,a*~,a*+,af] be two four-parameter interval grey numbers. If af >aj~, a*~ >af~, a*+ >o*+and a* then ¿¿(g) >bj(g). Similarly, if af 0) and Nec(x > 0) respectively stand for the possibility and necessity of x > 0. Definition 7 Let ^(g) = (ai~,a*~,a*+,a+) be the grey tardiness of the i-th job. Then, the tardiness credibility of the i-th job can be defined as the weighted sum of possibility measures and necessity measures and denoted as Con^x > 0). Therefore, Coni(x> 0) = SPos(x> 0) + (1 — 8)Nec(x > 0)(0 < S < 1). If 0 < S < 1, the following can be derived from Definition 7: If tardiness credibility of the i-th job is Coni(x> 0) = 1, then Nec(x> 0) = 1, i.e. the job production must be delayed; If Coni(x> 0) = 0, then Pos(x> 0) = 0, i.e. the production of the i-th job cannot be delayed. The coefficient S is determined by various factors, such as the decision-maker and the extent of the tardiness. Its value is negatively correlated with the conservativeness of the decision. The following properties can be obtained according to the definition of credibility: Property 1 If the grey tardiness of the i-th job is Tt(g) = (a[~,a*~,a*+,a[+)(ai" 0) of the job can be expressed as: Con¿(x> 0) = i 0, < <0 Sat + *+ , < >0, a*+ <0 aí ~ai S, «r >0, a*" <0 a*~-5a~[ *-, ai " >0, a¡~ <0 1, af >0 (4) Let aik(®),Di(®) and ^(g) be the grey makespan, the grey delivery period and the grey tardiness of the i-th job in the JSP, respectively. Suppose the feasible scheduling plans of all jobs are allocated to the set FS and the objective function is denoted as /(x). For a known scheduling plan x e FS, it is possible to establish the following grey mixed integer programming model: rninf(x) = I?=1 U --isL + a*28 + + oJmljk/n (5) xeS V "i "i "i "i J s.t. Zfc=iZ"=irnikJ = 1, i = 1,2, ...,n (6) £f=1 mikj <1, k = 1,2,...,m;j = 1,2,...,n (7) a-ikiS.®) = 2r=iS;|Liajfe(®)mijfe i,ji = 1,2, ..-,n; k = 1,2, ...,m (8) Ttm = aikj.m -Did®) i,jt = 1,2.....n (9) mikj = 0,1 i,j = 1,2,...,n;k = 1,2,...,m (10) ^ = 0,1 q 6 {1,2,3,4} (11) The goal of this model is to obtain the optimal scheduling plan that minimizes the value of the objective function /(x), i.e. the mean tardiness credibility of all jobs. Eq. 6 specifies that each job must occupy only one position in the process sequence on each machine; Eq. 7 stipulates that each position in the process sequence on each machine can only be occupied by one job; Eq. 8 defines aikj.(m), the grey makespan of the i-th job at the y'rth position of the process sequence on the k-th machine; Eq. 9 defines the grey tardiness T^®) of the i-th job; Eq. 10 defines the indicator variable mikj (mikj = 1 if the i-th job occupies the y'-th position in the process sequence on the k-th machine, and mikj = 0 if otherwise); Eq. 11 ensures that = 1 if a)q <0 and o>q = 0 if otherwise. Theorem 2 If the grey tardiness of the i-th job is T¡(®) = ,a*~ ,a*+ ®). Proof. Since the tardiness credibility coefficient S 6 [0,1], it is clear that 0 < 8af/(af — a*+)< S< (a*~ — 8al)/(a*~ —af) < 1 for the piecewise function of Eq. 4. Thus, we have 0 < 5a[+/(a+ -<+)< S< «" -Sap/Cap -a(")<1. It can be seen that T¡(®) — aj(®) = (ai~ — o*+,a*+ — a*",«* — af) (a¿~ — ^ a(", a*~ ~ct*+ ~aJ~ — and a* — ctf ConT.m_a.m. Q.E.D. Theorem 3 In a scheduling sequence S, if several jobs exist in the process sequence of the k-th machine, the grey makespan and delivery period of the i-th job are aik(®) and D¿(®), respectively, and the grey makespan and delivery period of the j-th are a/fc(®) and Dj(®), respectively. Note that aifc(®) = (a("fe,a*fe-,a*+,a+), a,-fc(®) = (a/fe,a*;-,a^k+,aJ+fe) , i,j = 1,2, ...,r and k = 1,2, ...,m; r is the number of jobs in the process sequence of the k-th machine; i and j are respectively the positions of the i-th and j-th jobs in the process sequence of the k-th machine. If A(®) 7^(0). It can be seen from Theorem 2 that the tardiness credibility Con^^ >ConT.0y Then, the grey tardiness 7^/(0) of the i-th job in scheduling sequence S' and that 7)(0) of the y'-th job in scheduling sequence S can be expressed as: 77(0) = I^iiawfc(0) + a,-fc(0) + + aifc(0) -£»¿(0) (16) 3} (0) = iL'ii awk (0) + aifc (0) + avfc (0) + o/fc (0) -Dj (0) (17) Since D,-(0) >£¿(0), we have T/(0) >7}(0) according to Eqs. 16 and 17. It can be seen from Theorem 2 that the tardiness credibility ConTt0 >ConT.(0y Therefore, Eqs. 12 and 13 show that the sum of tardiness credibility of the process sequence of the k-th machine is fk(S') >//¿(5). This means the scheduling sequence after exchanging the i-th job and the y'-th job is no better than that before the exchange. Thus, the y'-th job must be processed after the i-th job. Q.E.D. 3. Design of hybrid grey cuckoo search 3.1 Cuckoo search Artificial intelligence algorithm is usually used to solve job shop scheduling models. Artificial intelligence algorithms include swarm intelligence optimization algorithm, decision support algorithm [22] [23], data fusion algorithm [24] and machine learning algorithm [25], etc. Xu et al. [26] proposed a bat algorithm to solve the problem of two-flexible job shop scheduling. Cuckoo search (CS) [27] is a heuristic search algorithm developed by Yang and Deb in 2009. It was inspired by the obligate brood parasitism of some cuckoo species by laying their eggs in the nests of other host birds. The cuckoo algorithm has better performance than particle swarm optimization [28, 29]. Since it was established, the CS has been extensively studied by scholars around the world. In 2010, Yang and Deb applied the CS in the multi-objective solution problem [30]. In 2011, Valian et al. improved the CS based on feedforward neural network [31]. In the same year, Walton et al. proposed an improved Brown free gradient CS [32]. The basic idea of the CS comes from the breeding behaviour of cuckoo and the Levy flight mode of birds. The Levy flight is a random walk in which the step-lengths have a probability distribution that is heavy-tailed. When defined as a walk in a space of dimension greater than one, the steps made are in isotropic random directions. Considering its ability to avoid the local optimum trap and achieve good global optimization, the CS is adopted to solve the grey mixed integer programming model in this paper. The CS is suitable for small scale experiments rather than optimization problems with heavy computing load. Relying solely on random walk, the search for optimal solution consumes lots of computing power, and cannot guarantee the fast convergence of the algorithm. To improve the computing performance, it is necessary to integrate the features of the problem into the CS in a particular situation. Below are the steps to create the hybrid grey cuckoo search (HGCS). 3.2 Code design The CS is not directly applicable to the coding problem of discrete processes. Thus, the coding rule here is based on continuous processes. Suppose the position of each nest represents a feasible scheduling solution. For n jobs and m machines, the code of the feasible solution can be described by the locations of nxm nests. The code stands for a process sequence, in which each job must appear exactly m times. Taking a4-job 3-machine problem as an example, if the position code is 132143123442, the corresponding job processing sequence should be (J-L1, /31, /2,i, /1,2, /4,1, /3,2, /1,3,/2,2, /3,3, /4,2, /4,3, /2,3), with Ji,j being the y'-th process of the i. In other words, the first job and the third job should be processed successively on the first machine; then, the sec- ond job should be processed on the first machine, the first job should be processed on the second machine, and so on; Finally, the second job should be processed on the third machine. 3.3 Generation of initial population In the CS, the nest location must be initialized at the start However, the initial population is of poor quality because it is generated randomly without using any knowledge related to the problem to be solved. To improve the quality, the priority rules for minimum makespan and earliest delivery period in the JSP were introduced to develop a grey heuristic algorithm for initializing the locations of some nests. The specific steps are as follows. For the four-parameter interval grey number w(®) = (a~,a*",a*+,a+), let sign(w(®)^ = S(a~ + a*~) + (1 — 5)(a*+ + a+) be the grey mark of w(®), and S be the tardiness credibility coefficient. Then, sign(aik(®)) = S(a^k + a*fc) + (1 — 5)(a*fc + a*k) holds for the grey makespan aik(®) = (.aik,aik,atk), and sign^D^®)) = S(dl + + (1 — 5)(d*+ + holds for the grey delivery period £¿(0) = (d[;_,d*~,d*+,d[:f). Thus, the minimum makespan and earliest delivery period can be expressed as Min{sign{aik{®))} and Min{sign{Di{®))}, respectively. Then, the said grey heuristic can be established through the following steps: Step 1: Set up a set of grey markers for jobs, machines and delivery period: Pk = {sign(alk(®)),sign(a2k(®)), ...,sign(ank(®))}(k 6 [1,m]), Mt = {sign(ai1(®)),sign(ai2(®)),...,sign(aim(®))}(l 6 [1,n]) and Dq = {sign(Di(®)),sign(D2(®)), ...,sign(Dn(®))}. M( =Ml-sign(alk,(®j) Step 2: Let sign^Dq(®)j = Min(Dq), sign(aik>(®)) = Min(M{) and sign(al>k(®)) = Min(P{), that is, assign the i'-th job to the k'-th machine. Then, we have Pk = Pk — sign^al>k(®)^, Dq =Dq — sign^D^ (®)), with i = 1,2,..., n and k = 1,2,..., m. Step 3: Repeat Step 2 until the jobs are assigned to each machine. To diversify nest locations in the initial population, a third of the entire population was generated through the above steps, while the remaining parts were still generated randomly. The basic parameters of the CS were configured as: the dimension of the search space is m, the number of nests is n, the probability that the host bird discovers the egg laid by the cuckoo is Pa, and the maximum number of iterations is MaxT. After that, the initial nest locations were converted into process sequence through process coding, the value of objective function was computed corresponding to each nest location, and the optimal nest location was determined. Here, the goal is to calculate minf(x),x 6 5, i.e. the minimum mean tardiness credibility of all jobs in the fuzzy JSP. 3.4 Generation of candidate population First, the optimal nest location of the previous generation was preserved, and the fitness function fitness = Q — f(x), x 6 S of each nest location was calculated, with Q being a suitable positive number and f(x) being the objective function. Then, the candidate population was generated by exploring paths and locations through Levy flight. The direction of Levy flight is arbitrary, and the step length obeys a heavy-tailed probability distribution. To sum up, the candidate population can be generated by: xf+1 =xf + a ©Levy(A), i = 1,2, ...n (18) where xf is the location of the i-th nest in the t-th generation; a is the step length control coefficient, which is usually 1; © is the point-to-point multiplication; Levy{X) is the path of the random search of Levy flight. The direction of Levy flight obeys uniform distribution. 3.5 Optimal selection The locations of all nests in the current population were improved by Theorem 3, that is, the delivery periods of the nest locations were compared. If the grey delivery period £¿(0) < Dj(®) and aik(0) — D¿(0) 8 6 4 2 °0.1 1 10 100 1000 3000 Particle Size (|jm) Fig. 7 Powder particle size distribution measured by laser diffraction technique Particle size distribution Fig. 8 Scanning electron microscope image of Ti6Al4V ELI powder particles Astigmatism lens Focus lens Fig. 9 Schematic of EBM process Fig. 10 EBM setup from ARCAM Bonded powder Support structures Fig. 11 Powder recovery system The entire built process takes place under high vacuum of 0.1 Pa to 0.001 Pa. After the completion of built, fabricated parts are cooled down to room temperature under helium environment to prevent oxidation. The build envelope with fabricated custom implant is then taken to the powder recovery system (PRS) to blast the sintered powder using compressed air as shown in Fig. 11. The support structures are then removed and the thickness of the implant is measured using the screw gauge. The thickness was found to be approximately 0.52 mm as shown in Fig. 12. The weight of the fabricated implant is measured using digital weighing machine. The weight was found to be close to the weight of the bone portion that was removed from the skull, assuming the bone density to be p =210 kg m-3 [39]. The details of the fabricated implant and the equivalent bone portion are listed in Table 5. To evaluate the accuracy of the fabricated implant, it is first scanned using ViuScan hand scanner as show in Fig. 13. Then Geomagic Qualify software is used for the 3D comparison. The Geomagic quality is one of the commonly used techniques to graphically represent the surface deviation between two models. The EBM fabricated cranial implant was taken as test model and the original CT scan was taken as reference model. Using this software, the scanned EBM fabricated model was superimposed onto the reference CAD model with the Best Fit Alignment function using 1500 points. The differences were analyzed as positive and negative deviations. A positive deviation occurs if the test (scanned) model is larger than the reference (Original CT scan) model and a negative deviation if it is smaller. The mean deviations and root mean square of the deviations (RMSD) between test and reference model were calculated [41]. In order to demonstrate the magnitude, location and direction of the discrepancies between them, a color coded model was developed as shown in Fig. 14. Fig. 12 (a) The fabricated implant with support structure, (b) Implant after the removal of support Table 5 Weight details of the fabricated and the equivalent bone portion Material Thickness (mm) Volume (mm3) Weight (g) Bone portion Varying thickness: 3.18-5.76 mm 149114.285 31.314 Ti6Al4V ELI implant 0.52 9758.140 41.574 (a) (b) Fig. 13 (a) Implant scanning, (b) Scanned implant model V ¥ w 60529 j 5 0945 ! 4 1362 ' i 3 1778 22194 12610 _ m* [ -1.261o[ j •2.21941 -3 1 7781 -<4.13621 -5 09451 -6 0529 ! / Fig. 14 The 3D comparison result Fig. 15 Fabricated custom implant: (a) Top view, (b) Side view, (c) Front view The results showed that the average positive and negative deviations are 0.321 mm and -0.3515 mm respectively with standard deviation 0.3246 mm and the RMSD deviation is 0.5380 mm. The fit of the fabricated implant was evaluated by fixing it on the 3D printed plastic skull model as shown in Fig. 15. The visual inspection from all the sides shows that fabricated implant successfully achieved a good fit and aesthetics. 5. Conclusion In this study, a thin solid light weight custom cranial implant is designed, analyzed and fabricated using EBM. The craniofacial implant was designed based upon the patient CT scan images using MIMICS® and 3-MATIC®. Finite element analysis was carried out in order to assess the quality of the developed implant design under different loading conditions. The Ti6Al4V ELI implant of 0.5 mm thickness was then fabricated using EBM technology. The fitness of the fabricated implant was tested on the plastic prototype of the skull. The results showed the successful fabrication of 0.5 mm thick custom cranial implant for skull defect reconstruction via EBM technology while maintaining the structural strength requirements for reconstruction of skull defect . The main advantage of this approach for the skull defect reconstruction is the customizability, flexibility and less lead time for implant fabrication. However, this technology is still expensive and require significant amount of validation and testing before the implant is finally used for reconstruction of the defect. Acknowledgements The authors are grateful to the Deanship of Scientific Research, king Saud University for funding through Vice Deanship of Scientific Research Chairs. References [1] Gibson, I., Rosen, D.W., Stucker, B. (2010). Additive manufacturing technologies: Rapid prototyping to direct digital manufacturing, Springer, New York, USA, doi: 10.1007/978-1-4419-1120-9. [2] Morgan, L., Lindhe, U., Harrysson, O. (2003). Rapid manufacturing with electron beam melting (EBM) - A manufacturing revolution?, In: Proceedings of Annual International Solid Freeform Fabrication Symposium - An Additive Manufacturing Conference, Austin, Texas, USA, 433-438. [3] Berce, P., Chezan, H., Balc, N. (2005). The application of rapid prototyping technologies for manufacturing the custom implants, In: Proceedings of 8th ESAFORM Conference on Material Forming, Cluj-Napoca, Romania, 27-29. [4] Christensen, A., Kircher, R., Lippincott, A. (2008). Qualification of electron beam melted (EBM) Ti6Al4V-ELI for orthopaedic applications, In: Proceedings of MPMD, Materials and Processes for Medical Devices Conference, Palm Desert, USA, 48-53. [5] Li, X., Wang, C., Zhang, W., Li, Y. (2009). Fabrication and characterization of porous Ti6Al4V parts for biomedical applications using electron beam melting process, Materials Letters, Vol. 63, No. 3-4, 403-405, doi: 10.1016/j.matlet. 2008.10.065. [6] Abdel-Haleem, A.K., Nouby, R., Taghian, M. (2011). The use of the rib grafts in head and neck reconstruction, Egyptian Journal of Ear, Nose, Throat and Allied Sciences, Vol. 12, No. 2, 89-98, doi: 10.1016/j.ejenta.2011.08.004. [7] He, Y., Ye, M., Wang, C. (2006). A method in the design and fabrication of exact-fit customized implant based on sectional medical images and rapid prototyping technology, The International Journal of Advanced Manufacturing Technology, Vol. 28, No. 5-6, 504-508, doi: 10.1007/s00170-004-2406-y. [8] Durham, S.R., McComb, J.G., Levy, M.L. (2003). Correction of large (> 25 cm2) cranial defects with "reinforced" hy-droxyapatite cement: Technique and complications, Neurosurgery, Vol. 52, No. 4, 842-845, doi: 10.1227/01.NEU. 0000054220.01290.8E. [9] Ducic, Y. (2002). Titanium mesh and hydroxyapatite cement cranioplasty: A report of 20 cases, Journal of oral and maxillofacial surgery, Vol. 60, No. 3, 272-276, doi: 10.1053/joms.2002.30575. [10] Engstrand, T., Kihlstrom, L., Lundgren, K., Trobos, M., Engqvist, H., Thomsen, P. (2015). Bioceramic implant induces bone healing of cranial defects, Plastic and Reconstructive Surgery - Global Open, Vol. 3, No. 8, doi: 10.1097/G0X.000 0000000000467. [11] Charnley, G., Judet, T., Pirou, P., Garreau de Loubresse, C. (2000). Titanium femoral component fixation and experience with a cemented titanium prosthesis, In: Learmonth, I.D. (ed.), Interfaces in total hip arthroplasty, Springer, London, UK, 3-10, doi: 10.1007/978-1-4471-0477-3 1. [12] Moiduddin, K., Darwish, S., Al-Ahmari, A., ElWatidy, S., Mohammad, A., Ameen, W. (2017). Structural and mechanical characterization of custom design cranial implant created using additive manufacturing, Electronic Journal of Biotechnology, Vol. 29, 22-31, doi: 10.1016/j.ejbt.2017.06.005. [13] Niinomi, M. (2008). Mechanical biocompatibilities of titanium alloys for biomedical applications, Journal of the Mechanical Behavior of Biomedical Materials, Vol. 1, No. 1, 30-42, doi: 10.1016/j.jmbbm.2007.07.001. [14] Schipper, J., Ridder, G.J., Spetzger, U., Teszler, C.B., Fradis, M., Maier, W. (2004). Individual prefabricated titanium implants and titanium mesh in skull base reconstructive surgery. A report of cases, European Archives of Oto-Rhino-Laryngology and Head & Neck, Vol. 261, No. 5, 282-290, doi: 10.1007/s00405-003-0686-8. [15] Mazzoli, A., Germani, M., Raffaeli, R. (2009). Direct fabrication through electron beam melting technology of custom cranial implants designed in a PHANToM-based haptic environment, Materials & Design, Vol. 30, No. 8, 3186-3192, doi: 10.1016/j.matdes.2008.11.013. [16] Goyal, S., Goyal, M. K. (2014). Restoration of large cranial defect for cranioplasty with alloplastic cranial implant material: A case report, The Journal of Indian Prosthodontic Society, Vol. 14, No. 2, 191-194, doi: 10.1007/s13191-012-0185-y. [17] Lieger, O., Richards, R., Liu, M., Lloyd, T. (2010). Computer-assisted design and manufacture of implants in the late reconstruction of extensive orbital fractures, Archives of Facial Plastic Surgery, Vol. 12, No. 3, 186-191, doi: 10.1001/archfacial.2010.26. [18] Joffe, J.M., Nicoll, S.R., Richards, R., Linney, A.D., Harris, M. (1999). Validation of computer-assisted manufacture of titanium plates for cranioplasty, International Journal of Oral and Maxillofacial Surgery, Vol. 28, No. 4, 309-313, doi: 10.1016/S0901-5027(99)80165-9. [19] da Costa, D.D., Pedrini, H., Bazan, O. (2009). Direct milling of polymethylmethacrylate for cranioplasty applications, The International Journal of Advanced Manufacturing Technology, Vol. 45, No. 3-4, 318-325, doi: 10.1007/s00170-009-1978-y. [20] da Costa, D.D., Lajarin, S.F. (2012). Comparison of cranioplasty implants produced by machining and by casting in a gypsum mold, The International Journal of Advanced Manufacturing Technology, Vol. 58, No. 1-4, 1-8, doi: 10.1007/s00170-011-3388-1. [21] Saldarriaga, J.F.I., Velez, S.C., Posada, M.D.A.C., Henao, I.E.B.B., Valencia, M.E.C.A.T. (2011). Design and manufacturing of a custom skull implant, American Journal of Engineering and Applied Sciences, Vol. 4, No. 1, 169-174, doi: 10.3844/ajeassp.2011.169.174. [22] Dujovne, A.R., Bobyn, J.D., Krygier, J.J., Miller, J.E., Brooks, C.E. (1993). Mechanical compatibility of noncemented hip prostheses with the human femur, The Journal of Arthroplasty, Vol. 8, No. 1, 7-22, doi: 10.1016/S0883-5403 (06)80102-6. [23] Parthasarathy, J. (2014). 3D modeling, custom implants and its future perspectives in craniofacial surgery, Annals of Maxillofacial Surgery, Vol. 4, No. 1, 9-18, doi: 10.4103/2231-0746.133065. [24] Cho, H.R., Roh, T.S., Shim, K.W., Kim, Y.O., Lew, D.H., Yun, I.S. (2015). Skull reconstruction with custom made three-dimensional titanium implant, Archives of Craniofacial Surgery, Vol. 16, No. 1, 11-16, doi: 10.7181/acfs.2015. 16.1.11. [25] Murr, L.E., Gaytan, S.M., Martinez, E., Medina, F., Wicker, R.B. (2012). Next generation orthopaedic implants by additive manufacturing using electron beam melting, International Journal of Biomaterials, Vol. 2012, doi: 10.1155/ 2012/245727. [26] Al-Ahmari, A., Nasr, E.A., Moiduddin, K., Alkindi, M., Kamrani, A. (2015). Patient specific mandibular implant for maxillofacial surgery using additive manufacturing, In: Proceedings of 2015 International Conference on Industrial Engineering and Operations Management (IEOM), Dubai, UAE, 1-7, doi: 10.1109/IE0M.2015.7093788. [27] Mazzoli, A., Germani, M., Raffaeli, R. (2009). Direct fabrication through electron beam melting technology of custom cranial implants designed in a PHANToM-based haptic environment, Materials & Design, Vol. 30, No. 8, 3186-3192, doi: 10.1016/j.matdes.2008.11.013. [28] Li, X., Wang, C., Zhang, W., Li, Y. (2009). Fabrication and characterization of porous Ti6Al4V parts for biomedical applications using electron beam melting process, Materials Letters, Vol. 63, No. 3-4, 403-405, doi: 10.1016/j.matlet. 2008.10.065. [29] Saldarriaga, J.F.I., Velez, S.C., Posada, M.D.A.C., Henao, I.E.B.B., Valencia, M.E.C.A.T. (2011). Design and manufacturing of a custom skull implant, American Journal of Engineering and Applied Sciences, Vol. 4, No. 1, 169-174, doi: 10.3844/ajeassp.2011.169.174. [30] Syam, W.P., Al-Shehri, H.A., Al-Ahmari, A.M., Al-Wazzan, K.A., Mannan, M.A. (2012). Preliminary fabrication of thin-wall structure of Ti6Al4V for dental restoration by electron beam melting, Rapid Prototyping Journal, Vol. 18, No. 3, 230-240, doi: 10.1108/13552541211218180. [31] Al-Ahmari, A., Nasr, E.A., Moiduddin, K., Anwar, S., Al Kindi, M., Kamrani, A. (2015). A comparative study on the customized design of mandibular reconstruction plates using finite element method, Advances in Mechanical Engineering, Vol. 7, No. 7, doi: 10.1177/1687814015593890. [32] Tsouknidas, A., Maropoulos, S., Savvakis, S., Michailidis, N. (2011). FEM assisted evaluation of PMMA and Ti6Al4V as materials for cranioplasty resulting mechanical behaviour and the neurocranial protection, Bio-Medical Materials and Engineering, Vol. 21, No. 3, 139-147, doi: 10.3233/BME-2011-0663. [33] Allen, T., Goodwill, S., Haake, S. (2008). Experimental validation of a tennis ball finite-element model for different temperatures, In: Estivalet, M., Brisson, P. (eds.), The Engineering of Sport 7, Springer Nature, Switzerland, 125-133, doi: 10.1007/978-2-287-09411-8 15. [34] Chamrad, J., Marcian, P., Borak, L., Wolff, J. (2016). Finite element analysis of cranial implant, In: Proceedings of 32nd conference with international participation computational mechanics, Zelezna Ruda, Czech Republic, doi: 10.13140/ RG.2.2.19870.13124. [35] Parle, D., Ambulgekar, A., Gaikwad, K. (2015). Patient specific design and validation of cranial implants using FEA, In: Altair Technology Conference, Bengaluru, India, 1-9. [36] Ridwan-Pramana, A., Marcian, P., Borak, L., Narra, N., Forouzanfar, T., Wolff, J. (2017). Finite element analysis of 6 large PMMA skull reconstructions: A multi-criteria evaluation approach, Plos One, Vol. 12, No. 6, doi: 10.1371/ journal.pone.0179325. [37] Koike, M., Martinez, K., Guo, L., Chahine, G., Kovacevic, R., Okabe, T. (2011). Evaluation of titanium alloy fabricated using electron beam melting system for dental applications, Journal of Materials Processing Technology, Vol. 211, No. 8, 1400-1408, doi: 10.1016/j.jmatprotec.2011.03.013. [38] Karlsson, J. (2015J. Optimization of electron beam melting for production of small components in biocompatible titanium grades, Doctoral thesis, Acta universitatis upsaliensis. [39] Margulies, S.S., Thibault, K.L. (2000). Infant skull and suture properties: Measurements and implications for mechanisms of pediatric brain injury, Journal of Biomechanical Engineering, Vol. 122, No. 4, 364-371, doi: 10.1115/1. 1287160. [40] Ameen, W., Al-Ahmari, A., Mohammed, M.K., Mian, S.H. (2018). Manufacturability of overhanging holes using electron beam melting, Metals, Vol. 8, No. 6, 1-24, doi: 10.3390/met8060397. [41] Al-Ekrish, A.A., Alfadda, S.A., Ameen, W., Hormann, R., Puelacher, W., Widmann, G. (2018). Accuracy of computer-aided design models of the jaws produced using ultra-low MDCT doses and ASIR and MBIR, International Journal of Computer Assisted Radiology and Surgery, 1-8, doi: 10.1007/s11548-018-1809-4. [42] Colic, K., Sedmak, A., Legweel, K., Milosevic, M., Mitrovic, N., Miskovic, Z., Hloch, S. (2017). Experimental and numerical research of mechanical behaviour of titanium alloy hip implant, Tehnicki vjesnik - Technical Gazette, Vol. 24, No. 3, 709-713, doi: 10.17559/TV-20160219132016. [43] Ramakrishnaiah, R., Al Kheraif, A.A., Mohammad, A., Divakar, D.D., Kotha, S.B., Celur, S.L., Hashem, M.I., Vallittu, P.K., Rehman, I.U. (2017). Preliminary fabrication and characterization of electron beam melted Ti-6Al-4V customized dental implant, Saudi Journal of Biological Sciences, Vol. 24, No. 4, 787-796, doi: 10.1016/j.sjbs.2016.05.001. Advances in Production Engineering & Management Volume 13 | Number 3 | September 2018 | pp 279-296 https://doi.Org/10.14743/apem2018.3.290 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper Dynamic integration of process planning and scheduling using a discrete particle swarm optimization algorithm Yu, M.R.a*, Yang, B.a, Chen, Y.b aSchool of Mechatronic Engineering, North University of China, Taiyuan, P.R. China bNorth Automatic Control Technology Institute, Taiyuan, P.R. China A B S T R A C T A R T I C L E I N F O Because of the inherent relationship between process planning and scheduling, integration of process planning and scheduling (IPPS) provides a new path for further improvements of these two activities. Therefore, a novel two-phase IPPS approach is put forward in this paper. In the new method, the preplanning phase generates a process network for each job with consideration of the static shop floor status. After that, the final planning phase simultaneously creates the process plan of each job and the scheduling plan according to the current shop floor status. Based on the modified definition of IPPS and the proposed mathematical model, the IPPS problem and the dynamic IPPS problem can be solved together. Furthermore, a discrete particle swarm optimization (DPSO) algorithm is proposed to solve the IPPS optimization problem. In the DPSO algorithm, the particles update their positions by crossing with their own historical best positions (pbests) and the global best position of the population (gbest). In order to avoid local convergence, an external archive is introduced to keep more than one elite, and the gbest of each particle is randomly selected from the external archive. Furthermore, mutation operation is introduced to enhance the local search ability of DPSO algorithm. Finally, some comparative results are given to verify the efficiency and effectiveness of the proposed IPPS method and the DPSO algorithm as well as the dynamic IPPS method. © 2018 CPE, University of Maribor. All rights reserved. Keywords: Process planning; Scheduling; Dynamic integration; Mathematical model; Optimization; Discrete particle swarm optimization (DPSO) *Corresponding author: yumingrang@163.com (Yu, M.R.) Article history: Received 4 December 2017 Revised 21 August 2018 Accepted 24 August 2018 1. Introduction In the common manufacturing systems, which transform the raw material or semi-finished product into final product using kinds of machines, some preparation activities such as materials, tools, process plans, scheduling plans and so on. In these activities, process planning and scheduling are two crucial functions which are usually carried out sequentially. In other words, the process plans, which are the outcomes of process planning, are transferred into the scheduling system, which assigns operations to specified machines at appropriate moments according to the precedence relations in the process plans, shop floor status, scheduling criteria, etc. Obviously, the effectiveness of the scheduling results should be strongly dependent on the process plans. However, in the past years, these two activities are often executed sequentially, in which the scheduling plans are generated separately according to the fix process plans obtained by the process planning. However, when considering the disturbances in manufacturing process such as arrival of urgent jobs, due date change and machine breakdown, the traditional method, in which process planning and scheduling are separately treated, seems to be inadequate for following reasons [1]: • Traditionally, the process plans are generated by process planners under some ideal assumptions, for example, the resources in the shop floor are always available, which is unrealistic in a real manufacturing environment; • The conventional process planning methods provide deterministic process plans to scheduling system, which ignores the possibility of improvement for scheduling with alternative process plans; • Because of the time delay between process planning phase and scheduling phase, even if the dynamic shop floor status is considered in advance, it may change greatly when the schedule plan is executed, thus the generated optimal process plans may become suboptimal or even invalid; • In traditional way, the process planning generates optimal plans with the consideration of process plan criterion, and the scheduling works in a similar manner, conflicts may appear in such a separate way. To overcome these shortages, Chryssolouris et al. [2] first proposed the concept of integration of process planning and scheduling (IPPS). When process planning and scheduling are integrated, the manufacturing system can respond promptly to disruptions. As a result, the resource utilization and the productivity will be improved. In the following sections, some literature reviews on IPPS are given in Section 2, and Section 3 describes the materials and methods. Then the discrete particle swarm optimization algorithm is presented in Section 3, which is followed by the application of DPSO in optimizing the proposed IPPS problem in Section 4, which is followed by some experiment results in Section 5 and conclusions in Section 6. 2. Literature review Ever since the concept of IPPS was proposed by Chryssolouris et al. [2], a lot of research papers on IPPS could be found. Based on the descriptions in some existing publications on IPPS [1, 3-6], the IPPS methods can be divided into three categories: non-linear process planning (NLPP), closed loop process planning (CLPP) and distributed process planning (DPP). These IPPS approaches and related contributions are described below. NLPP: The NLPP method attempts to generate as many process plans for each part as possible under ideal shop floor status with consideration of the process flexibility, sequence flexibility [7] and operation flexibility. Here, the process flexibility represents the possibility of machining a specified feature using alternative processes or series of processes, hereafter, use macro-level plan to stand for process or a serial of processes for machining a feature. And the sequence flexibility stands for possibility of different sequences of the selected macro-level plans. Finally, the operation flexibility means that a specified process may be arranged to different machines. After that, different priority levels are assigned to these process plans according to the optimization objectives of process planning, e.g. total machining cost/time, etc. And then, these process plans will be tested in the scheduling system according to their priority levels. CLPP: In the CLPP method, the process plans are always feasible with respect to the current scheduling environment because they are generated just before the jobs are released to the shop floor, in the other words, the status of each resource in the shop floor is known to the process planner. Some of the CLPP approaches first generate off-line plans according to the static shop floor status, and then make some necessary on-line refinement on the basis of the availability of resources on the shop floor at the subsequent scheduling phase. DPP: The DPP method performs both process planning and scheduling simultaneously in a hierarchic manner. In most of the DPP methods, they can be divided into two phases, which are preplanning and final planning. Some other DPP approaches may have three phases: preplan- ning, pair planning and final planning. Similarly, at the preplanning phase, the process planning function recognizes the features and the feature relationships of each job, and determines the macro-level plan for each feature. Meanwhile, the scheduling function estimates the required machine capabilities. Then, at the final planning phase (which can be further divided into pair planning and final planning in some DPP methods), both the final process plan of each part and the scheduling plan are generated simultaneously, in this way, the integration occurs at the moment that the machining processes and the available resources are matched. Both the NLPP and CLPP methods are one-way information flows as shown in Fig. 1, in which process planning is still executed before scheduling. In this way, the objectives of process planning and scheduling may be conflicted. Therefore, Both NLPP and CLPP are incapable of finding the global optimal solution. With this in mind, the DPP method, which performs the technical and capacity-related planning tasks simultaneously, may be the only one IPPS method that integrates the functions of process planning and scheduling. However, in most of the DPP methods [8, 9], the macro-level plan of each feature is determined in the preplanning stage; and the selected macro-level plans for all features are used for matching the operation capabilities of the available resources in the final planning stage. It means that the flexibility of the final planning phase is reduced because of ignoring the process flexibility. As mentioned before, optimization algorithm is another research focus on IPPS because of its NP-hard characteristic [11, 12]. Kinds of optimization algorithms, such as genetic algorithm (GA) [13], particle swarm optimization (PSO) algorithm [8], etc. have been used for optimizing the IPPS problems. However, most of these algorithms are designed for one of the three existing IPPS methods. Therefore, the performances of those optimization algorithms may be restricted by the corresponding IPPS methods. As is known to all, the basic purpose of IPPS is handling the inevitable disturbances in the manufacturing processes by treating process planning and scheduling as a whole. However, even in the newest publications on IPPS, dynamic IPPS is still less mentioned in the literatures for its complexity. For example, Zhang and Wong [14] proposed an object-coding genetic algorithm for IPPS, and Wang etal. [15] addressed a systematic approach for optimizing the IPPS problems in sustainable machining. The models in these two papers are still designed for static IPPS, in which the results of IPPS are final process plans and schedules. It is worth mentioning that even the optimal solution is obtained using integrated method, this solution may also be deteriorated in the execution stage because of the inevitable disturbances. Therefore, dynamic IPPS is necessary for further improving the performances of the manufacturing system. To sum up, an effective IPPS method and the corresponding optimization algorithm are still necessary to be developed to get full integration of process planning and scheduling. Therefore, this paper proposes a novel IPPS method to overcome the shortages of the three existing IPPS methods. Based on the modified definition of IPPS, the mathematical model of the proposed IPPS method is presented, which is also suitable for dynamic IPPS method. Finally, a discrete particle swarm optimization is developed for optimizing the proposed IPPS problem. Fig. 1 Information flows of: (a] NLPP and (b) CLPP [10] 3. Materials and methods 3.1 The proposed IPPS method As shown in Fig. 2, the proposed IPPS method is a two-phase approach, which also contains preplanning and final planning. Be different from the existing DPP methods, at the preplanning phase, instead of determining the macro-level plan for each feature, the process network, which is widely used in representing the flexibility of process plan [16, 17], of each part to be machined is generated according to the static resources. A process network types of nodes, which are starting node, intermediate node and ending node, more details for the process network can be found in the reference of Ho et al. [17]. Then the process networks of all the parts to be machined are stored and will be used in the final planning phase. With consideration of optimization criteria of both process planning and scheduling, an optimal plan, which includes the process plans of all parts to be machined and the scheduling plan, is generated according to the current shop floor status. Of course, an effective optimization algorithm is necessary because of the complexity of the problem. After that, the optimal plan is released to the shop floor immediately. Suppose the makespan of the optimal plan is T and current time is to. Because disturbances such as machine breakdown, rush order and order cancellation may occur during the time interval to to to + T, a dynamic IPPS method, in which both the process plans of the jobs and the schedule are changeable, is presented to handle the disruptions. Production plan Feature-based design data Static resources PP knowledge 4=3 ¡Dynamic resources Optimization algorithms 0=S Preplanning Process networks Final planning Process Scheduling plans plan Dynamic shop floor status Update the process networks Manufacturing Dynamic IPPS execution monitoring Fig. 2 The structure of the proposed IPPS method (PP stands for process planning) 3.2 Mathematical model of the proposed IPPS method In order to solve the IPPS and the dynamic IPPS in a uniform manner, the definition of IPPS [18] is modified as follows. Given a set of n jobs, at the time t, the number of the unfinished features of each part is Nf, the number of the available machine at time t is Mt. considering the flexibility of process, sequence and operation, determine the macro-level plan for each unfinished feature and the corresponding machine (machine set) as well as the sequences of the operations on each machine. Meanwhile, the precedence constraints in each part are satisfied and some optimization objectives can be achieved. Based on the proposed IPPS method and the modified definition of IPPS, the mathematical model of the proposed IPPS method is founded. In this paper, the optimization objectives of process planning and scheduling are minimum machining time and minimum makespan, separately. Firstly, some reasonable assumptions and notations are given below: • Jobs are independent. Job preemption is not allowed and each machine can handle only one job at any moment. • The different operations of the same job cannot be machined simultaneously. • The setup time and the transport time are negligible or included in the processing time. Mt the total number of the available machines at the time t; N the total number of jobs to be machined at the time t; Nif the number of the unfinished features of the z'-th job at the time t; Nijo the number of alternative macro-level plans of the j-th feature of z'-th job; Ozjk the k-th alternative macro-level plan of the j-th feature of z'-th job; Mzjki the l-th alternative machine (machine set) of the macro-level plan Ok tijki the machining time of macro-level plan Ozjk on machine l; cZjkl the completion time of macro-level plan Oijl on machine k; (1 the k-th alternative macro-level the k-th alternative macro-level plan is selected for y'-th feature of i-th job othervise the Z-th machine is selected for the macro-level plan Otjki othervise f1 the Z-th machine is selected for the io Objectives: • Minimizing the total machining time: N Nif Ni]o Mt tíjkl XYijkl x%ijk (1) ¿ = 1j = ik = l 1 = 1 • Minimizing the makespan: F2 = Min{Makespan = Max{cijkl xXijk xy¿jfeí} (2) Subject to: • The j-th feature of the i-th job can only choose one macro-level plan from its alternatives: • The macro-level plan Oijk can only choose one machine (machine set) from its alternatives: Other constraints can be found in Li et al. [19]. It is worth mentioning that the numbers of the available machines and the jobs to be processed as well as the unfinished features of each job are varied as time goes by. Therefore, this mathematical model is also suitable for the dynamic IPPS method. And these kinds of information of dynamic shop floor status are provided by the manufacturing execution system. 3.3 Particle swarm optimization algorithm The proposed IPPS method involves in four tasks, which are selection of macro-level plan for each unfinished feature, selection of machine (machine set) for each selected macro-level plan, sequencing the selected macro-level plans and determining the starting time of each process on its corresponding machine. It is a typical NP-hard problem and much more complex to solve than that of existing IPPS methods, which are partially involved in the four referred tasks. Particle swarm optimization (PSO) algorithm has been widely used in many optimization problems since it was proposed in 1995 [20]. However, the continuous characteristic of PSO algorithm restricts its applications in combinatory optimizations problems. Some researchers [21, 22] used modified encoding methods to transform the scheduling problems into continuous versions, and then the PSO can be used to solve these scheduling problems. However, they are incapable of solving the optimization problem of the proposed IPPS method since its complexity. Therefore, this paper proposes a discrete particle swarm optimization (DPSO) algorithm for the proposed IPPS method. Standard PSO algorithm In most of the researches, the standard PSO refers to a modified PSO algorithm [23], in which an inertia weight is introduced to improve the optimizing ability. In the PSO algorithm, a swarm, which contains a number of particles with certain positions and velocities, is initialized randomly. Then each particle dynamically adjusts its velocity and position according to its own and the population's experiences. This can be explained as follows: (3) (4) = 'uVfa + c1r1(pbestid -Xfd) + c2r2(gbestd -Xfd) Vdmin ^dmax ^dmin + V¿td+1 «Xj^ «X, X, V, dmax f^id Vid+1 x. dmax (5) (6) yt+1 N Admax is the velocity of the d-th dimension of the i-th particle at the time t, and it means the distance to be traveled in the d-th dimension from its current position. w is the inertial weight used for regulating the trade-off between the global exploration and local exploration abilities of the swarm. Ci and C2 are the acceleration constants, riand r2 are two random functions within the range of [0,1]. pbestid represents the d-th dimension of the z-th particle's best historical position, andgbestd represents the d-th dimension of the swarm's global best position. The particles update their positions and velocities according to Eq. 5 and Eq. 6, and the pbest and gbest are also updated according to the fitness function. Thus all the particles move to the global best solution to finish the search process. Discrete PSO algorithm It can be found that the position and the velocity in standard PSO are both continuous variables from Eq. 5 and Eq. 6, which are meaningless for the IPPS optimization problems. This paper proposes a discrete version of PSO, in which the particles update their positions according to the following equation: xr = x< t+i gbest (7) Xf+1 and Xf are the positions of the z-th particle in current and the next iteration, Pfbest is the historical best position of the z-th particle, and Pg^est is the best position of the population in the tth iteration. ® stands for the crossover operator. Eq. 7 shows that the particles adjust their positions according to their own and the population's experiences, which maintains the advantages of the standard PSO algorithm. The discrete particle swarm optimization algorithm can be explained as follow according to Fig. 3: • The particleZf and its pbest P[bestare considered as two parents Pi and P2; Oi and O2 are two offspring of them. Select the better one of these two offspring for the next crossover, e.g. O2 is selected; • O2 andP^best, which is the gbest of Xf, are treated as two new parents P[and P2; and 0'21 are two offspring of them. Select the better one of these two offspring to be used as the new position of Xf, that is Xt t+i 1 2 3 2 3 1 3 1 PLcj 2 1 3 3 1 1 3 2 OI 1 2 3 3 3 1 2 1 1 o2(P{) 2 1 3 2 1 1 3 3 Pgbesl (•P2) 1 3 1 2 3 3 1 2 r br;+1 (o,) 1 1 3 2 1 2 3 3 02 2 3 1 2 3 3 1 1 Fig. 3 Illustration of discrete particle swarm optimization algorithm Generally, the gbest of the swarm is unique. Under this condition, each particle will cross with that gbest, which may lead to local convergence. For this reason, this paper introduces an external archive to keep the elites, which are the first Nea (Nea the size of the external archive) parti- cles with lowest fitness values, generated in each iteration. Then each particle randomly selects an individual from the external archive as its gbest in each iteration, and the members of the external archive are also updated iteratively. 4. The application of discrete particle swarm optimization (DPSO) in optimizing the proposed IPPS problem When applying the DPSO algorithm to optimize the IPPS problem, some necessary explanations and definitions should be given in advance. The information of a process network can also be listed in a table. Table 1 lists the information contained in four process networks of four jobs, in which '999' means that machine cannot finish the process, and Oi2n-Oi2i2 means the macro-level plan O121 contains two processes. It is worth mentioning that different macro-level plans of the same feature may have different machining time on the same machine. For example, the machining time of Oiii and O112 on machine Mi is 3 and 5 units, separately. The machine priority level for a macro-level plan is defined as follows. Machine priority level. For a macro-level plan, the machine priority level (1>2>3>4>5, etc.) is determined by the machining time of the machine or machine set. The longer the machining time of a machine is, the higher its machine priority level will be. And if the machining time of two machines (machine sets) for the same macro-level plan are equal, the bigger the number of the machine is, the higher its machine priority level will be. Table 1 Information contained in the four process networks of four jobs Jobs Features Alternative Alternative machines and machining time macro-level plans Ml M2 M3 M4 Ms Oiii 999 3 4 999 4 O112 5 999 3 4 999 O113 3 4 999 999 4 O1211-O1212 2/999 3/4 999/3 999/4 999/999 O1221-O1222 3/2 999/3 4/999 4/3 3/999 Oi3i 4 5 999 3 999 O132 999 999 3 4 5 Oi33 3 4 3 5 3 O211 6 999 5 999 999 O212 999 6 999 6 5 O2211-O2212 3/4 999/2 4/3 999/999 3/999 O231 5 3 999 4 999 O232 4 999 4 5 5 O233 4 4 5 999 3 O241 5 3 4 999 4 O242 4 999 5 999 6 O3111-O3112-O3113 999/3/999 4/5/3 4/4/6 999/4/5 5/999/4 O321 4 999 5 3 4 O322 5 5 4 6 4 O33i 4 6 4 5 999 O332 3 5 3 4 999 O333 999 4 999 3 5 O3411-O3412 999/4 4/5 5/3 3/999 999/999 O342 10 10 9 999 999 O411 4 999 5 999 4 O421 3 4 999 3 5 O422 5 3 4 999 999 O423 3 4 999 4 3 O4311-O4312-O4313 5/3/999 4/5/4 5/999/6 999/5/999 4/6/5 O441 5 6 5 999 5 O442 6 999 5 6 4 O4431-O4432 999/4 2/3 4/999 999/999 3/4 J4 Fii (before F12) F12 Fi3 F21 F22 F23 (before F24) F24 F31 F32 F33 (before F34) F34 F41 F42 (before F43) F43 F44 Ji J3 Table 2 Machine priority levels of the alternative machines Jobs Features Alternative Machine priority level(1>2>3>4>5) macro-level plans 1 2 3 4 5 Ji F11 O111 M2(3) M3(4) Ms(4) Mi(999) M4(999) (before O112 Ma(3) M4(4) Mi(5) M2(999) Ms(999) F12) O113 Mi(3) M2(4) Ms(4) M3(999) M4(999) F12 O1211-O1212 Mi+Ma(5) Mi+M2(6) Mi+M4(6) M2+M3(7) M2(7) O1221-O1222 Mi(5) Ms+Mi(5) M5+M2(6) M4(7) M4+M2(7) F13 O131 M4(3) Mi(4) M2(5) M3(999) Ms(999) O132 Ma(3) M4(4) Ms(5) Mi(999) M2(999) O133 Mi(3) M3(3) Ms(3) M2(4) M4(5) J2 F21 O211 M3(5) Mi(6) M2(999) M4(999) Ms(999) O212 Ms(5) M2(6) M4(6) Mi(999) M3(999) F22 O2211-O2212 Mi+M2(5) MS+M2(5) Mi+M3(6) Mi(7) M3(7) F23 O231 M2(3) M4(4) Mi(5) M3(999) Ms(999) (before O232 Mi(4) M3(4) M4(5) Ms(5) M2(999) F24) O233 Ms(3) Mi(4) M2(4) M3(5) M4(999) F24 O241 M2(3) M3(4) Ms(4) Mi(5) M4(999) O242 Mi(4) M3(5) Ms(6) M2(999) M4(999) J3 F31 O3111-O3112-O3113 Mi+M2(i0) M2+M4(ii) M2(i2) Mi+Ms(i2) M3(i4) F32 O321 M4(3) Mi(4) Ms(4) M3(5) M2(999) O322 M3(4) Ms(4) Mi(5) M2(5) M4(6) F33 O331 Mi(4) M3(4) M4(5) M2(6) Ms(999) (before O332 Mi(3) M3(3) M4(4) M2(5) Ms(999) F34) O333 M4(3) M2(3) Ms(5) Mi(999) M3(999) F34 O3411-O3412 M4+M3(6) M4+Mi(7) M2+M3(7) M3(8) M2(9) O342 M3(9) Mi(i0) M2(i0) M4(999) Ms(999) J4 F41 O411 Mi(4) Ms(4) M3(4) M2(999) M4(999) F42 O421 Mi(3) M4(3) M2(4) Ms(5) M3(999) (before O422 M2(3) M3(4) Mi(5) M4(999) Ms(999) F43) O423 Mi(3) Ms(3) M2(4) M4(4) M3(999) F43 O4311-O4312-O4313 Mi+M2(i2) M2(i3) Mi+Ms(i3) Mi+M3(i4) Ms(i5) F44 O441 Ms(4) M3(5) Mi(6) M4(6) M2(999) O442 Mi(5) M3(5) Ms(5) M2(6) M4(999) O4431-O4432 M2(5) M2+Mi(6) M2+MS(6) Ms(7) Ms+Mi(7) Based on the definition of machine priority level, the machine priority levels of the machines shown in Table 1 can be obtained as shown in Table 2. As referred before, the process network is generated based on static shop floor status in the preplanning phase, thus in the final planning phase, the status of the shop floor as well as the machine priority levels should be updated. 4.1 Encoding and decoding Encoding refers to problem mapping, which translates a problem to a particle. For the jobs listed in Table 2, a possible particle which contains three sections is shown in Fig. 4. The feature string, whose length equals to the total number of the features of all jobs, represents the manufacturing sequence of the features. For example, '1, 2, 3' represent the features of job 1, and '4, 5, 6, 7' represent the features of job 2, etc. The MLP string stands for the selected macro-level plans for corresponding features, e.g. the second 2 in the MLP string corresponds to the number 9 in the feature string. It means that the second feature of the third job (F32) has selected its second alternative macro-level plan (O322). The machine string represents the machine priority levels of the selected machines, e.g. the macro-level plan referred before (O322) has selected the first prior machine (M3) from its alternatives because the value of the corresponding position in the machine string is 1. Feature string 3 1 2 4 9 10 12 15 8 11 6 7 5 13 14 MLP string 1 3 2 1 2 3 1 2 1 2 1 2 1 2 1 Machine string 2 1 2 2 1 1 2 2 1 2 2 1 2 3 2 Fig. 4 A particle for the four jobs listed in Table 2 (MLP stands for macro-level plan) Decoding means solution generation, which translates a particle to a solution of the optimization problem. In the proposed IPPS method, the solution contains the process plan of each job and the scheduling plan. Based on the encoding method referred before, the decoding process is described as follows (take the particle shown in Fig. 4 for example): Step 1: Translate the feature string to the features machining sequence of each job, for example, the four jobs' feature machining sequence are '3, 1, 2', '4, 6, 7, 5', '9, 10, 8, 11' and '12, 15, 13, 14', separately; meanwhile, the machining sequence of all the features is also determined, that is F13-F11-F12-F21-F32-F33-F41-F44-F31-F34-F23-F24-F22- F42-F43; Step 2: Translate the MLP string to the selected macro-level plans for corresponding features; Step 3: Translate the machine string to the selected machines for corresponding macro-level plans; Step 4: Combine the results of steps 1, 2 and 3 to generate the process plan for each job. Step 5: According to the results of steps 1, 2 and 3, the features machining sequence, selected macro-level plans and corresponding machines are known, and then the decoding method proposed by Zhang et al. [24] is introduced to decode the particle to an active scheduling plan. 4.2 Fitness function In this paper, the optimizing objectives of process planning and scheduling are minimum total machining time and minimum makespan, respectively. Although these two objectives have the same dimension, their orders of magnitude may vary enormously. Thus the normalization of the data is still necessary. Let Di(Ti, Mi) (1 < i < N) stand for the assessment of the solution corresponding to the i-th particle, in which N is the number of the particles in the swarm. Ti and Mi represent the total machining time and the makespan, respectively. The procedure of normalization can be explained as follows: T.' = Ti-mun{Tj,lO24l(M2) O311l(M2)-O3112(Ml)-O3113(M2) -O32l(M4)-O333(M4)-O341l(M4) -O3412(M3) O44l(M5)-O41l(M5)-O423(M5) -O431l(M2)-O4312(Ml)-O4313(M2) 71 O111(M2)-O1211(M1)-O1212(M3)- O131(M4) O21l(M3)-O23l(M4)-O221l(Ms)- O2212(M2>O24l(M5) O3111(M2)-O3112(M1)-O3113(M2)-O321(M4)-O331(M3)-O3411(M4)- O3412(M3) O421(M4)-O441(Ms)-O411(Ms)- O4311(M2)-O4312(M1)-O4313(M2) 75 O113(M1)-O132(M3)-O1211(M1) -O1212(M3) O211(M3)-O231(M4)-O2211(M5) -O2212(M2)-O241(M5) O3111(M2)-O3112(M1)-O3113(M2) -O321(M4)-O333(M4)-O3411(M4) -O3412(M3) O442(M5)-O411(M5)-O423(M5) -O4311(M2)-O4312(M1)-O4313(M2) 73 J3.2 Jl.2 Jl.3 J 1.4 J4.5 El I I J3.4 I J3.5 I J3.6 J4.6 J2.4 27 m, m2 m3 m4 m5 J3.2 ju ■ 23 J3.1 1 1.1 J3.3 j4.4 j2.4 j4.6 h j3.5 [ht j4.1 j2.2 j3.4 j jl.4 j j3.6 j [4.2 U.3 J2.3 J2.5 1 1 ~jü~1 pü7| Q JíT] [^3] QI J4.6 j s s J3.6 1 j4.2 i. J2.3 J2.5 5 10 15 20 25 (a) Makespan = 27 30 0 5 10 15 20 (b) Makespan = 23 5 10 15 20 (c) Makespan = 22 Fig. 9 Gantt charts of the schedules (a) NLPP and CLPP, (b) DPP and (c) the proposed IPPS method 5.3 Case 3 As referred before, based on the modified definition and the mathematical model of IPPS, dynamic scheduling can be extended to dynamic IPPS for effectively handling the disruptions occurred in the execution process. In this section, two typical dynamic events, which are machine breakdown and arrival of urgent jobs, are designed to verify the effectiveness of the dynamic IPPS method. Machine breakdown The optimal scheduling plan obtained by the proposed IPPS method is released to the shop floor. Suppose that at the time t = 5 the machine M4 is broken-down with a relatively long repair time. Under this condition, traditionally, reactive scheduling methods often transfer the affected oper- ations to their alternative machines for instead, which is the so called route changing rescheduling [35]. In the route changing rescheduling method, neither the alternative macro-level plan of each feature nor the operation sequence of each job is unchangeable. However, in the dynamic IPPS method, the process plans are also changeable, which may improve the flexibility of rescheduling. The comparative results are shown in Table 6 (process plans) and Fig. 10 (scheduling plans), respectively. From the results we can find that both the total machining time and makespan are reduced by the dynamic IPPS method. Table 6 Process plans of the route changing rescheduling strategy and the dynamic IPPS method Job Route changing rescheduling Dynamic IPPS method Ji J2 J3 J4 T Ol13(Ml)-Ol32(M3)-Ol21l(Ml)-Ol212(M3) O21l(M3)-O23l(M2)-O221l(M5) -O2212(M3)-O24l(M2) O311l(M2)-O3112(Ml)-O3113(M2)-O32l(Ml)-O333(M2)-O341l(M2)-O3412(M3) O442(M5)-O41l(M5)-O423(M5)-O431l(M2) -O4312(Ml)-O4313(M5) 76 Ol13(Ml)-Ol32(M3)-Ol21l(Ml)-Ol222(M3) O21l(M3)-O233(M5)-O24l(M2) -O221l(M5)-O2212(M2) O311l(M2)-O3112(Ml)-O3113(M2)-O322(M3) -O341l(M2)-O3412(M3) -O332(M3) O442(M5)-O41l(M5)-O423(Ml)-O431l(M5) -O4312(Ml)-O4313(M5) 74 m, m2 m3 m4 m5 M4 breakdown I H7| [J7] ^ h M4 breakdown J2.1 J2.4 I a ■J4.2 J4.3 m, m2 m3 m4 m5 1 n«i J4.3 J J1.3 J J4.5 J J3.1~| 1 J.3.3 J2.3 J3.5 J2.5 ! J2.1 Jl.2 J3.4 J 1.4 J3.6 J3.7 i 1 ! J4.1 j J4.2 J2.2 J4.4 J2.4 J4.6 : 1 5 10 15 20 (a) Makespan = 26 5 10 15 20 (b] Makespan = 24 Fig. 10 Gantt chart of route changing rescheduling method (a) and the dynamic IPPS method (b) Arrival of urgent jobs The optimal scheduling plan obtained by the proposed IPPS method is released to the shop floor. Suppose two urgent jobs are released to the shop floor at the moment t = 5. The information of these two jobs is listed in Table 7. Generally, the urgent jobs should be completed as soon as possible. The comparative results are shown in Fig. 11, in which the dynamic IPPS method also outperforms the route changing rescheduling method. From the above comparative results, it is obvious that the dynamic IPPS method is much more suitable than other rescheduling methods because it can increase the flexibility of rescheduling. Of course, the robustness and the stability of the dynamic IPPS method should be further studied. Arrival of urgent jobs 28 Arrival of urgent jobs 25 m, m2 m3 m4 •J3.2 J5.1 J2.3 J1.3 m, J2.2 J3.3 J4.4 J4.5 J4.6 Jó.6 Jfi-l J3.4 Jó.3 J3.5 J3.6 Jl.4 ■J4.2 J5.2 J2.5 0 25 m2 m., m4 ms 30 0 1 Jl.l Í3.2 J5.I J 1.3 |Jó.3 J4.5 J2.5 J5.5 J3.1 j J2.2 J3.3 J4.4 Jö.4 Jó.5 Jó.6 J2.1 J6.1 J5.2 J5.4 J3.7 Jl.2 Jó.2 J3.4 J3.5 J3.6 Jl.4 J4.I 'J4.2 J4.3 1 J2.3 J5.3 1 J4.6 1 20 5 10 15 20 (a] Makespan = 28 Fig. 11 Gantt charts of the route changing rescheduling method (a) and the dynamic IPPS method (b) 5 10 15 (b) Makespan = 25 25 Table 7 Information of the two urgent jobs Job Feature Alternative macro-level plan Alternative machines and machining time Mi M2 M3 M4 Ms J5 F51 O511 4 6 999 5 4 (before F52) Osi2 5 999 3 4 999 O513 999 4 5 999 6 F52 Os21 5 6 999 5 4 O522 999 5 4 5 999 F53 Os31 3 5 4 999 4 O532 4 999 3 4 999 O533 999 4 999 5 3 F54 O5411-O5412 4/3 999/4 4/999 5/4 999/999 Os42 999 6 7 999 7 J6 F61 O611 4 999 5 4 999 O612 5 4 999 6 5 F62 O621 999 5 4 6 5 O622 4 3 999 5 4 F63 O631 5 4 5 4 999 (before F64) O632 4 999 3 999 5 O633 4 4 5 6 999 F64 O6411-O6412-O6413 3/4/999 999/4/3 4/999/5 5/4/4 999/999/4 m, m2 m3 m4 m5 Arrival of urgent jobs I 28 Arrival of urgent jobs 25 J3.1 j3, Js.l J2.3 Jl.3 h .4 h.s J5.5 j J2.2 1 J3.3 J44 J4.5 J4.6 Í6.6 J2.1 h.2 h A J3.7 J. Je.i J3.4 J,3 ,3, J3.6 Jl.4 '4.1 I 'h.2 J5 h. J2.5 0 m, m2 m3 m4 m5 l'o 0 1 Jl.l {3.2 J5.1 Jl.3 h.3 J4.5 J2.5 J5.5 J3.1 J2.2 J3.3 J4.4 h A Jó.5 - Jö.6 J2.1 Jó.l ^sT] J2.4 J5.4 J3.7 Jl.2 h.2 J3.4 J3.5 J3.6 Jl.4 J4.I J4.2 J4.3 1 J2.3 J5.3 1 J4.6 1 20 10 15 20 25 30 0 5 10 15 (a] Makespan = 28 (b] Makespan = 25 Fig. 11 Gantt charts of the route changing rescheduling method (a) and the dynamic IPPS method (b) 25 6. Conclusion This paper proposes a novel IPPS method which has combined the advantages of the three existing IPPS methods. Based on the modified definition of IPPS, the mathematical model of the proposed IPPS method is presented. This model is also suitable for the dynamic IPPS method, which can greatly improve the flexibility of rescheduling. Furthermore, a discrete version of particle swarm optimization algorithm is put forward to solve the optimization problem of the proposed IPPS method. In the DPSO algorithm, the particles update their positions by crossing with their own historical best positions (pbests) and the global best position of the population (gbest). Therefore, the continuous PSO algorithm can be used for the combinatory optimization problems such as optimization of the IPPS problems. In order to avoid local convergence, an external archive is introduced to keep more than one elite, and the gbest of each particle is randomly selected from the members of the external archive. Finally, the mutation operation is introduced to enhance the local search ability of the DPSO algorithm. As shown in the experiments and results of case 1, the proposed DPSO algorithm stands out from several typical intelligent optimization algorithms for its efficiency and effectiveness. Besides, the PDSO is designed for the proposed IPPS problem, which involves in four tasks as referred before, and the process planning, scheduling or the three existing IPPS methods are partially of the proposed IPPS problem, therefore, the DPSO algorithm could also be used to solve those problems. In this work, we select the total machining time to judge the process plan and the makespan for scheduling plan, to simplify, we have combined these two objectives into one weighted function. However, other objectives such as total machining cost, total tardiness and utilization of machines may also be considered to complete the proposed IPPS method. It is no doubt that the Pareto-based optimization approaches will make this work more perfect, which will be studied in our future work. Furthermore, the robustness and the stability of the dynamic IPPS method should be further studied as well. Acknowledgements This work was supported by Shanxi Foundation Research Projects for Application (201701D221146) and Natural Science Foundation of North University of China (Grant No. 2017003). References [1] Phanden, R.K., Jain, A., Verma, R. (2011). Integration of process planning and scheduling: A state-of-the-art review, International Journal of Computer Integrated Manufacturing, Vol. 24, No. 6, 517-534, doi: 10.1080/09511 92x.2011.562543. [2] Chryssolouris, G., Chan, S., Cobb, W. (1984). Decision making on the factory floor: An integrated approach to process planning and scheduling, Robotics and Computer-Integrated Manufacturing, Vol. 1, No. 3-4, 315-319, doi: 10.1016/0736-5845(84)90020-6. [3] Yu, M., Zhang, Y., Chen, K., Zhang, D. (2015). Integration of process planning and scheduling using a hybrid GA/PSO algorithm, The International Journal of Advanced Manufacturing Technology, Vol. 78, No. 1-4, 583-592, doi: 10.1007/s00170-014-6669-7. [4] Tan, W., Khoshnevis, B. (2000). Integration of process planning and scheduling - A review, Journal of Intelligent Manufacturing, Vol. 11, No. 1, 51-63, doi: 10.1023/a:1008952024606. [5] Li, X., Gao, L., Zhang, C., Shao, X. (2010). A review on integrated process planning and scheduling, International Journal of Manufacturing Research, Vol. 5, No. 2, 161-180, doi: 10.1504/I|Mr.2010.031630. [6] 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/IJSIMM 16(4)7.400. [7] Huang, X.W., Zhao, X.Y., Ma, X.L. (2014). An improved genetic algorithm for job-shop scheduling problem with process sequence flexibility, International Journal of Simulation Modelling, Vol. 13, No. 4, 510-522, doi: 10.2507/ IISIMM13(4)C020. [8] Guo, Y.W., Li, W.D., Mileham, A.R., Owen, G.W. (2009). Optimisation of integrated process planning and scheduling using a particle swarm optimisation approach, International Journal of Production Research, Vol. 47, No. 14, 3775-3796, doi: 10.1080/00207540701827905. [9] Leung, C.W., Wong, T.N., Mak, K.L., Fung, R.Y.K. (2010). Integrated process planning and scheduling by an agent-based ant colony optimization, Computers & Industrial Engineering, Vol. 59, No. 1, 166-180, doi: 10.1016/j.cie. 2009.09.003. [10] Zhang, H.-C., Merchant, M.E. (1993). IPPM - A prototype to integrate process planning and job shop scheduling functions, CIRP Annals, Vol. 42, No. 1, 513-518, doi: 10.1016/S0007-8506(07)62498-6. [11] Dai, M., Tang, D., Xu, Y., Li, W. (2014). Energy-aware integrated process planning and scheduling for job shops, Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture, Vol. 229, No. 1, 13-26, doi: 10.1177/0954405414553069. [12] Galzina, V., Lujic, R., Saric, T. (2012). Adaptive fuzzy particle swarm optimization for flow-shop scheduling problem, Tehnicki vjesnik - Tehnical Gazette, Vol. 19, No. 1, 151-157. [13] Chaudhry, I.A., Usman, M. (2017). Integrated process planning and scheduling using genetic algorithms, Tehnicki vjesnik - Technical Gazette, Vol. 24, No. 5, 1401-1409, doi: 10.17559/TV-20151121212910. [14] Zhang, L., Wong, T.N. (2015). An object-coding genetic algorithm for integrated process planning and scheduling, European Journal of Operational Research, Vol. 244, No. 2, 434-444, doi: 10.1016/j.ejor.2015.01.032. [15] Wang, S., Lu, X., Li, X.X., Li, W.D. (2015). A systematic approach of process planning and scheduling optimization for sustainable machining, Journal of Cleaner Production, Vol. 87, No. 914-929, doi: 10.1016/j.jclepro.2014. 10.008. [16] Kim, Y.K., Park, K., Ko, J. (2003). A symbiotic evolutionary algorithm for the integration of process planning and job shop scheduling, Computers & Operations Research, Vol. 30, No. 8, 1151-1171, doi: 10.1016/S0305-0548(02) 00063-1. [17] Ho, Y.-C., Moodie, C.L. (1996). Solving cell formation problems in a manufacturing environment with flexible processing and routeing capabilities, International Journal of Production Research, Vol. 34, No. 10, 2901-2923, doi: 10.1080/00207549608905065. [18] Guo, Y.W., Li, W.D., Mileham, A.R., Owen, G.W. (2009). Applications of particle swarm optimisation in integrated process planning and scheduling, Robotics and Computer-Integrated Manufacturing, Vol. 25, No. 2, 280-288, doi: 10.1016/j.rcim.2007.12.002. [19] Li, X., Gao, L., Shao, X., Zhang, C., Wang, C. (2010). Mathematical modeling and evolutionary algorithm-based approach for integrated process planning and scheduling, Computers & Operations Research, Vol. 37, No. 4, 656667, doi: 10.1016/j.cor.2009.06.008. [20] Kennedy, J., Eberhart, R. (1995). Particle swarm optimization, In: Proceedings of ICNN'95 - International Conference on Neural Networks, Perth, Australia, 1942-1948, doi: 10.1109/ICNN.1995.488968. [21] Lin, T.-L., Horng, S.-J., Kao, T.-W., Chen, Y.-H., Run, R.-S., Chen, R.-J., Lai, J.-L., Kuo, I.-H. (2010). An efficient jobshop scheduling algorithm based on particle swarm optimization, Expert Systems with Applications, Vol. 37, No. 3, 2629-2636, doi: 10.1016/j.eswa.2009.08.015. [22] Sha, D.Y., Hsu, C.-Y. (2006). A hybrid particle swarm optimization for job shop scheduling problem, Computers & Industrial Engineering, Vol. 51, No. 4, 791-808, doi: 10.1016/j.cie.2006.09.002. [23] Shi, Y., Eberhart, R. (1998). A modified particle swarm optimizer, In: 1998 IEEE International Conference on Evolutionary Computation Proceedings. IEEE World Congress on Computational Intelligence, Anchorage, USA, 6973, doi: 10.1109/ICEC.1998.699146. [24] Zhang, G., Gao, L., Shi, Y. (2011). An effective genetic algorithm for the flexible job-shop scheduling problem, Expert Systems with Applications, Vol. 38, No. 4, 3563-3573, doi: 10.1016/j.eswa.2010.08.145. [25] Li, W.D., Ong, S.K., Nee, A.Y.C. (2004). Optimization of process plans using a constraint-based tabu search approach, International Journal of Production Research, Vol. 42, No. 10, 1955-1985, doi: 10.1080/0020754031 0001652897. [26] 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. [27] Kacem, I., Hammadi, S., Borne, P. (2002). Approach by localization and multiobjective evolutionary optimization for flexible job-shop scheduling problems, IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews), Vol. 32, No. 1, 1-13, doi: 10.1109/tsmcc.2002.1009117. [28] Brandimarte, P. (1993). Routing and scheduling in a flexible job shop by tabu search, Annals of Operations Research, Vol. 41, No. 3, 157-183, doi: 10.1007/bf02023073. [29] Li, J.-Q., Pan, Q.-K., Gao, K.-Z. (2011). Pareto-based discrete artificial bee colony algorithm for multi-objective flexible job shop scheduling problems, The International Journal of Advanced Manufacturing Technology, Vol. 55, No. 9-12, 1159-1169, doi: 10.1007/s00170-010-3140-2. [30] Wang, X., Gao, L., Zhang, C., Shao, X. (2010). A multi-objective genetic algorithm based on immune and entropy principle for flexible job-shop scheduling problem, The International Journal of Advanced Manufacturing Technology, Vol. 51, No. 5-8, 757-767, doi: 10.1007/s00170-010-2642-2. [31] Xia, W., Wu, Z. (2005). An effective hybrid optimization approach for multi-objective flexible job-shop scheduling problems, Computers & Industrial Engineering, Vol. 48, No. 2, 409-425, doi: 10.1016/j.cie.2005.01.018. [32] Moslehi, G., Mahnam, M. (2011). A Pareto approach to multi-objective flexible job-shop scheduling problem using particle swarm optimization and local search, International Journal of Production Economics, Vol. 129, No. 1, 14-22, doi: 10.1016/j.ijpe.2010.08.004. [33] Bagheri, A., Zandieh, M., Mahdavi, I., Yazdani, M. (2010). An artificial immune algorithm for the flexible job-shop scheduling problem, Future Generation Computer Systems, Vol. 26, No. 4, 533-541, doi: 10.1016/j.future.2009. 10.004. [34] Xing, L.-N., Chen, Y.-W., Yang, K.-W. (2009). An efficient search method for multi-objective flexible job shop scheduling problems, Journal of Intelligent Manufacturing, Vol. 20, No. 3, 283-293, doi: 10.1007/s10845-008-0216-z. [35] He, W., Sun, D.-H. (2013). Scheduling flexible job shop problem subject to machine breakdown with route changing and right-shift strategies, The International Journal of Advanced Manufacturing Technology, Vol. 66, No. 1-4, 501-514, doi: 10.1007/s00170-012-4344-4. APEM jowatal Advances in Production Engineering & Management Volume 13 | Number 3 | September 2018 | pp 297-306 https://doi.Org/10.14743/apem2018.3.291 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper An integral algorithm for instantaneous uncut chip thickness measuring in the milling process Li, Y.a, Yang, Z.J.a, Chen, C.b' *, Song, Y.X.a, Zhang, J.J.a, Du, D.W.c aCollege of Mechanical Science and Engineering, Jilin University, Changchun, Jilin, P.R. China bCollege of Engineering, China Agricultural University, Beijing, P.R. China institution of Oceangraphic Instrumentation, Shandong Academy of Science, Qingdao, P.R. China A B S T R A C T A R T I C L E I N F O Instantaneous uncut chip thickness (IUCT) calculation is an essential work for dynamic cutting force prediction accurately in milling process. This study presents an integral algorithm in polar coordinate system for measuring the thickness of transient uncut chip. The milling trajectory, cycloidal motion, is adopted in the formulation. Both milling continuity and cutter run-out are also considered in this model. The developed model offers a methodology for calculating the IUCT precisely. Furthermore, a series of simulations are carried out under different processing parameters. The results suggest that increasing both the feed per tooth and number of teeth can surge the width of IUCT slightly, but decrease with smaller cutter radius. The milling force simulations are validated by the experiment results measured in the reference and compared with classical approximate method, showing the proposed IUCT model providing good applications in instantaneous milling force predictions. © 2018 CPE, University of Maribor. All rights reserved. Keywords: Milling; Instantaneous uncut chip thickness; Dynamic cutting forces; Integral algorithm *Corresponding author: chenchao2018@cau.edu.cn (Chen, C.) Article history: Received 7 December 2017 Revised 29 July 2018 Accepted 24 August 2018 1. Introduction Milling is an important processing method in the field of machining. Milling force that originates in the cutter-workpiece interface is an important factor that affects vibrations, milling efficiency, surface quality and the chatter stability of CNC tools [1-3]. Various analytical prediction models for cutting force have been established to improve the accuracy and reliability of milling force [4-6]. These models differ considerably, but the cutting force and IUCT are closely related [7, 8]. Besides, the force prediction accuracy is also important for monitoring of work conditions, such as tool life, surface roughness, even the machining stability of the milling process [9, 10]. Martel-lotti was first to use an approximated formula to calculate chip thickness. The simplified tooth path is considered as circular and lack of tooth eccentricity or tooth run-out [11]. The IUCT expression of the tooth path is Tt(a) = frN •sin a, where frN is feed per tooth and a is the instantaneous angular position of the tooth, which is widely used and researched [12, 13]. S. Spiewak then proposed an improved model of calculating IUCT in milling process. He facilitated stepwise and orderly increases in model sophistication until a desirable level of performance was achieved. The application results indicates that the speed, accuracy and reliability of monitoring and control potentially have been all improved, while the calculation process is rather complex and time consuming [14]. Another classical model is the one proposed by Li et al. in 2001. Instead of using a numerical method, they calculated the thickness of instantaneous uncut chip with Taylor's series by analyzing true tooth trajectories [15]. In Kumanchik's model, an analytic expression for uncut chip thickness in milling process was formulated while considering a series of impact factors, such as the cycloidal motion of teeth, run-out, and uneven teeth spacing [16]. Some people then expressed interest in the model through making improvements or modifications in the calculations. The main contribution of Ge Song and his co-authors is their investigation on the accurate positions of geometric points in the profile of IUCT with different level of cut width. This study also adopted the iterative algorithm to improve the accuracy of the thickness of uncut chip at the chip cross-section [17]. Chip thickness in circular interpolation was examined in 2008. Sai et al. found an opposing trend between chip thickness and the radius of circular trajectory. The regions or movements of the cutter should be also considered [18]. N. Grossi improved the trochoidal motion formula by introducing run-out values by two variables- the distance (d) between geometric centers of two adjacent cutters and the instantaneous cutter angle a [19]. However, the present methods for measuring IUCT ignored the true path shape of the cutter tooth. Most used models entail a circular tool-path approximation evenly. As mentioned earlier, cutting force is primarily a function of IUCT in an analytical cutting force approach. These simplified models could not achieve satisfied accuracy definitely. Besides, Milling is a continuous process and the approximate circumferential model considers the IUCT as discrete line that connects the tool center to the current tooth's cutting edge which fails to describe the transient chip formation and also limits the application scope of formulas. This study introduces an integral algorithm in polar coordinate system for measuring the thickness of transient uncut chip. The milling cutter trajectory, cycloidal motion, is adopted in the novel formula, in which the IUCT is associated with feed rate per teeth, number of teeth and cutter radius. Both the effect of milling continuity and the cutter run-out are also considered in this method. Furthermore, different processing parameters are used to obtain a series of IUCT. The simulation milling forces are validated by the experiments and the classical method, showing the advantages and accuracy of using the new proposed method. The nomenclature is given in the Appendix A. 2. Milling tooth path and equation Researchers investigated the path of cutter tip during milling. Most of them assumed the tooth path of cutting points as circular [14]. However, the true path shape of the cutter should be considered as a trochoid, which is a more reasonable approach [1, 15]. Fig. 1(a) describes movement in the direction of a straight line of the milling cutter. The radius of the cutter is assumed as r, where in the teeth are evenly distributed along the circumference. When the milling cutter is in clockwise rotation, the total teeth are labeled as 1,2,3,...,i,...,N. The speed of motor spindle is given as n (r/min). The feed rate of the workpiece is labeled as ty (mm/min). Thus, milling feed rate can be calculated using the Eqs. 1 and 2 [20]: vf=frN^N^n (1) (2) n = TfD (a) (b) Fig. 1 Model for milling cutter rotation: (a) milling along the liner direction; (b) milling tooth trajectory in Cartesian coordinate system Variable frN is feed distance per tooth for one cycling, v is the rate of cutting, D is cutter diameter, which is twice the radius of the milling tool. The reasonable parameters can be given according to machining experience. Cutter movement consists of its own rotation and rectilinear motion along the cutting direction. The feed distance of the tool center varies with the rotation angle of the cutter. Thus, feed distance for per angle /rad(mm/rad) can be deduced from Eqs. 1 and 2: frad=(3) As shown in Fig. 1(b), the initial center of the tool is labeled as C0, which coincides with the origin of the Cartesian coordinate system. The direction of the Z-axis is opposite to the feed path of the workpiece and the 7-axis is normal to the Z-axis. The Z-axis is perpendicular to the XOY coordinate plane and has the same direction as the cut depth. The equations represent the trajectories of arbitrary cutting point Q, namely i tooth. The tip of i + 1 tooth can be expressed as in [15]: Xi = frzVi + rsin(<^ -at) (4) Vi =rcos(q)i - ait) si Js1 where angle a is measured with its polar axis along the direction of X-axis positive. The above equation is transformed into connect the instantaneous angle of cutter^ and the area of uncut chip thickness, where (15) a = n n 2~Vi, 0 i o -0.05 -0.1 frt0.2r5N4 frt0.2r8N4 " - '' \ \ "V S 0 100 120 140 160 180 Angle Fig. 8 Effect of different radius of milling cutter on uncut chip thickness (frt = 0.2mm/tooth, r = 5 mm and 8 mm, N = 4) 6. Conclusion A method based on polar integral algorithm for calculating the thickness of transient uncut chip thickness is proposed for the cutting force predictions in milling process. The analytical equation is also suitable for engineering application. The main findings of the research: • True milling trajectory considered as trochoid motion was established in the polar coordinates. Both milling continuity and milling trajectory are considered when developing a new model for IUCT calculations. The new developed model offers a methodology for calculating the IUCT precisely. • The milling forces obtained by experiment in reference [21] and simulations with the developed IUCT model are compared in the paper. The integral model proves to be more reasonable in terms of the location where the peak appears. Results indicate that the outcomes followed by the polar integral algorithm are close to the reality. • Different process parameters are simulated and compared according to the results calculated by proposed model. The comparisons suggest that increasing both the feed per tooth and number of teeth can surge the width of IUCT slightly, but decrease with smaller cutter radius. Acknowledgement The authors would like to express their acknowledgment to National Science and Technology Major Project (2016ZX04004007), China Postdoctoral Science Foundation (2017M621203) and Science and Technology Project in Education Department of Jilin Province (JJKH20180133KJ). References [1] Omar, O.E.E.K., El-Wardany, T., Ng, E., Elbestawi, M.A. (2007). An improved cutting force and surface topography prediction model in end milling, International Journal of Machine Tools and Manufacture, Vol. 47, No. 7-8, 12631275, doi: 10.1016/j.ijmachtools.2006.08.021. [2] Matsubara, A., Ibaraki, S. (2009). Monitoring and control of cutting forces in machining processes: A review, Monitoring and Control of Cutting Forces in Machining Processes, Vol. 3, No. 4, 445-456, doi: 10.20965/ijat.2009. p0445. [3] Mwinuka, T.E., Mgwatu, M.I. (2015). Tool selection for rough and finish CNC milling operations based on tool-path generation and machining optimisation, Advances in Production Engineering & Management, Vol. 10, No. 1, 18-26, doi: 10.14743/apem2015.1.189. [4] Milfelner, M., Cus, F., Balic, J. (2005). An overview of data acquisition system for cutting force measuring and optimization in milling, Journal of Materials Processing Technology, Vol. 164-165, 1281-1288, doi: 10.1016/ j.jmatprotec.2005.02.146. [5] Liu, X.-W., Cheng, K., Webb, D., Luo, X.-C. (2002). Improved dynamic cutting force model in peripheral milling. Part I: Theoretical model and simulation, The International Journal of Advanced Manufacturing Technology, Vol. 20, No. 9, 631-638, doi: 10.1007/s001700200200. [6] Ikua, B.W., Tanaka, H., Obata, F., Sakamoto, S. (2001). Prediction of cutting forces and machining error in ball end milling of curved surfaces -I theoretical analysis, Precision Engineering, Vol. 25, No. 4, 266-273, doi: 10.1016/ S0141-6359(01)00077-0. [7] Qu, S., Zhao, J., Wang, T., Tian, F. (2015). Improved method to predict cutting force in end milling considering cutting process dynamics, The International Journal of Advanced Manufacturing Technology, Vol. 78, No. 9-12, 1501-1510, doi: 10.1007/s00170-014-6731-5. [8] Han, Z., Jin, H., Fu, H. (2015). Cutting force prediction models of metal machining processes: A review, In: Proceedings of 2015 International Conference on Estimation, Detection and Information Fusion (ICEDIF), Harbin, China, 323-328, doi: 10.1109/ICEDIF.2015.7280216. [9] Gradisek, J., Kalveram, M., Weinert, K. (2004). Mechanistic identification of specific force coefficients for a general end mill, International Journal of Machine Tools and Manufacture, Vol. 44, No. 4, 401-414, doi: 10.1016/ j.iimachtools.2003.10.001. [10] Saric, T., Simunovic, G., Simunovic, K. (2013). Use of neural networks in prediction and simulation of steel surface roughness, International Journal of Simulation Modelling, Vol. 12, No. 4, 225-236, doi: 10.2507/IISIMM12 (4)2.241. [11] Martellotti, M.E. (1945). An analysis of the milling process, Part II: Down milling, Transactions of ASME, Vol. 67, No. 4, 233-251. [12] Gonzalo, O., Beristain, J., Jauregi, H., Sanz, C. (2010). A method for the identification of the specific force coefficients for mechanistic milling simulation, International Journal of Machine Tools and Manufacture, Vol. 50, No. 9, 765-774, doi: 10.1016/i.iimachtools.2010.05.009. [13] Budak, E., Altinta§, Y., Armarego, E.J.A. (1996). Prediction of milling force coefficients from orthogonal cutting data, Journal of Manufacturing Science and Engineering, Vol. 118, No. 2, 216-224, doi: 10.1115/1.2831014. [14] Spiewak, S. (1995). An improved model of the chip thickness in milling, CIRPAnnals, Vol. 44, No. 1, 39-42, doi: 10.1016/s0007-8506(07)62271-9. [15] Li, H.Z., Liu, K., Li, X.P. (2001). A new method for determining the undeformed chip thickness in milling,Journal of Materials Processing Technology, Vol. 113, No. 1-3, 378-384, doi: 10.1016/S0924-0136(01)00586-6. [16] Kumanchik, L.M., Schmitz, T.L. (2007). Improved analytical chip thickness model for milling, Precision Engineering, Vol. 31, No. 3, 317-324, doi: 10.1016/i.precisioneng.2006.12.001. [17] Song, G., Li, J., Sun, J. (2013). Approach for modeling accurate undeformed chip thickness in milling operation, The International Journal of Advanced Manufacturing Technology, Vol. 68, No. 5-8, 1429-1439, doi: 10.1007/ s00170-013-4932-y. [18] Sai, L., Bouzid, W., Zghal, A. (2008). Chip thickness analysis for different tool motions: For adaptive feed rate, Journal of Materials Processing Technology, Vol. 204, No. 1-3, 213-220, doi: 10.1016/i.imatprotec.2007.11.094. [19] Grossi, N., Sallese, L., Scippa, A., Campatelli, G. (2015). Speed-varying cutting force coefficient identification in milling, Precision Engineering, Vol. 42, 321-334, doi: 10.1016/i.precisioneng.2015.04.006. [20] Fu, G. (2000). Basic knowledge of metal cutting, In: Shi, Q.F, Wang, L. (eds.), Metal cutting manual, (Third edition), Shanghai Science and Technology Press, Shanghai, 1-31, (in Chinese). [21] Altintas, Y. (2012). Manufacturing automation: Metal cutting mechanics, machine tool vibrations, and CNC design, (Second edition), Cambridge University Press, New York, USA, doi: 10.1017/CB09780511843723. [22] Moufki, A., Dudzinski, D., Le Coz, G. (2015). Prediction of cutting forces from an analytical model of oblique cutting, application to peripheral milling of Ti-6Al-4V alloy, The International Journal of Advanced Manufacturing Technology, Vol. 81, No. 1-4, 615-626, doi: 10.1007/s00170-015-7018-1. [23] Armarego, E.J.A., Deshpande, N.P. (1993). Force prediction models and CAD/CAM software for helical tooth milling processes. II. Peripheral milling operations, International Journal of Production Research, Vol. 31, No. 10, 2319-2336, doi: 10.1080/00207549308956860. Appendix A Nomenclature Tt(a) IUCT calculated by traditional model

0,1 > 0,0 > 0 (1) where A is the intensity parameter, ft is the shape parameter, and t is the operating time. Following the NHPP, the number of failures during the time interval (t1,t2] equals to W(t1,t2) = E(N(t2)-N(t1)) = f a)(pL)dpL (2) Jti The same can be inferred that W(0,t) = E(N(t)) represents the expected number of failures through time interval [0,t], thereby the cumulative failure intensity function is given by W(t) = i = At? (3) Jo Then the corresponding reliability function is derived as R(t) = (4) As well as the cumulative density function F(t) = 1- R(t) = 1 - e~XtP (5) And the probability density function = = w(t)K(t) = XptP~ie~ktP (6) at However, a single NHPP model can only illustrate the situation that failure intensity and failure time are strictly monotonic, which is unable to describe the non-monotonic trends of different life stages. Therefore, a sectional model of multiple NHPPs is required to fit the bathtub-shaped curve in order to describe the rules of failures in different failure periods, which could help to obtain the change point. 2.1 Two sectional NHPP model Experts have developed various models to describe complex and diverse data distributions [22], most of them are constructed assuming that systems are non-repairable or 'repair as new'. It is not appropriate for the concern in our case. Based on this, we propose a sectional model involving two NHPPs, representing the early failure period and the random failure period, respectively. 0<^o (7) (w2(t) = A2fat^~\ t0 1, the failure intensity increases as t progresses, it means that the failure rate will rise with system ageing. When fi = 1, the failure intensity is a constant. Accordingly, the cumulative failure intensity, which is also known as the cumulative failure number, has a change pattern seen in Fig. 1. Through further mathematical derivation, the conclusion agrees with the trend of failure intensity, it indicates that the growth rate of cumulative failure number will first reduce with time in early failure period, when it comes to random failure period, the growth rate will maintain a steady state, after that it will continue to rise during wear-out failure period. In view of our proposed sectional model, it is aimed at describing the changing trend of the early failure period and the random failure period. Through setting the constraints that failure intensity and cumulative failure intensity are both continuous at t0, the accurate change point t0 could be obtained. which is expressed as (M1(t0) = ù)2(t0) IWi(to) = W2(to) Substitute Eq. 7 and Eq. 8 in Eq. 9 (9) Ai ftic/1"1^^2"1 A^1 = -htop2 (10) Since the cumulative failure intensity function is much in evidence to be continuous at t0 with the sectional two NHPP modeling, the change point t0 can be obtained as Mir7 <") At this point, the bathtub-shaped failure intensity and its change point can be achieved with our proposed sectional model, which can be adopted to assist the failure process monitoring. Burn-in Useful life Wear out Time to failure Time to failure Fig. 1 Bathtub-shaped curves of the failure intensity and corresponding cumulative failure number 3. Change-point estimation combined SPC with sequential clustering For the purpose of estimating the exact change point, a specific bootstrap control chart to monitor the performance of repairable systems is established. The objective is to track possible transition points by detecting all out-of-control signals during the failure process. These candidate points help to classify data patterns into the two phases of early failure period and random failure period, respectively. Then, the clustering techniques are adopted to eliminate interference data within out-of-control signals in order to identify the optimum transition point. Finally, the accurate change point can be estimated by fitting the two sections of our proposed failure intensity model with the segmented observation data. The proposed approach integrates the advantages of SPC method and clustering analysis. Considering the characteristics of SPC, two possible states are first predetermined as in-control and out-of-control, it will simplify the clustering model with a known number of clusters. Meanwhile, the data is monitored in time series, the order of points are preserved for the iteration to search for the optimal clusters. At last, cluster settings will only be examined when an out-of-control signal is detected. The number of iterations is dramatically reduced to save computing time. Therefore, the SPC method is beneficial for improving efficiency for change-point estimation. Additionally, the intervention of sequential clustering approach is for benefits of lower interference and better accuracy, because the SPC chart in this study is required to be more sensitive to sustained shifts and gradual changing trends rather than sudden changes or random noises. 3.1 The bootstrap control chart Traditional control charts [23], like standard Shewhart charts, can only be strictly applied to the normal distribution by monitoring the shifts of mean and variance. As to those modified CUSUM charts and EWMA charts [24], they are always set up for some specific distributions. In our case, the failure process is modelled with a sectional NHPP, it's difficult to find the applicable control charts and corresponding methods to establish its control limits. Therefore, a bootstrap control chart with no limits of any specific distribution is developed. The particular SPC chart is established based on Monte Carlo simulation [25], which can be implemented as the following steps shown in Fig. 2. Step 1: Determine the stable values of parameter As and They can be defined by experience or MLEs with observations collected in Phase I, which represents a stable state in SPC theory. Step 2: Generate bootstrap random variables til, ti2, —,tin from the NHPP model with parameters As and in size n. n equals to the data volume of Phase II, which is the subsequent monitoring process based on the estimated control limits. The method from [26] is improved to generate random variables of failure times til,ti2, ■■■,tin following the NHPP sequentially. Step 3: Calculate the MLEs from til, ti2, ■..,tin to obtain and ^ of the NHPP models. Step 4: Obtain the p-quantile Wpi with Wpi = (—(1/2;) In p)1/ where Wpi represents 100 p-th percentile of interest. Step 5: Repeat the steps from Step.2 to Step.4 for B(B > 1000) times to obtain B groups of p-quantiles. Order them from the smallest to the largest as Wpl < Wp2 < ••• < WpB. Then the centre line (CL) is yielded as a mathematical expectation of acquired data, the lower control limit (LCL) is the i-th p -quantile Wpi, and the upper control limit (UCL) is the (B — j)-th value as where i = aB/2, representing that there're i estimators be- yond the control limits. The false alarm risk a indicates the probability when the system is diagnosed to be out of control while it's actually in control. Fig. 2 Procedure to establish the bootstrap control chart After control limits are obtained, the bootstrap control chart can be established through online monitoring. Time between failures (TBF) observations are plotted after each failure, and process shifts can be detected under some suitable runs rules. Take an instance as shown in Fig. 3, we assume there are 200 observations in total, the first 100 points in Phase I are supposed to be in control. Then the control limits, including UCL, LCL and CL, are calculated following the steps from Step 1 to Step 5. When it's converted to Phase II, the other 100 points are plotted on the established control chart, which should have been detected to be out of control as we expected. However, the fact in Fig. 3 shows that several observations (red solid points) in Phase I are beyond the control limits, meanwhile, a large number of observations (black solid points) in Phase II are within the control limits. It infers that the proposed control chart has a problem of random interference, the accuracy of detection will be reduced and it would be hard to determine which signal is the optimal transition point. In view of this problem, the sequential clustering approach is introduced to combine with the proposed SPC method. 3.2 The sequential clustering approach The purpose of introducing clustering approach is to remove the disturbance in SPC and extract the optimal point f, it can detect whether systems convert to another failure period and divide observations in two phases. In the proposed sequential clustering approach, two possible clusters are defined: for the i-th out-of-control signal Ot, all observations before Ot are classified to the in-control cluster, and all observations after Ot are automatically classified to the out-of-control cluster. In the interest of obtaining the optimal transition point t, first of all, the clustering model with applicable validity indices is developed. Then, the objective function of each cluster setting according to time series is examined. At last, the optimal transition point t is obtained by optimising the objective function. Two validity indices and the objective function from [16] are illustrated and improved as follows. Clusters within variation A cluster within variation Vwithin is expressed as the distance between observations and their cluster centres, given by Cin of in-control cluster centre and Cout of out-of-control cluster centre, as shown in Fig. 4. In consideration of Phase I and Phas e I I based on the specific bootstrap control chart, Cin in Phase I equals to the in-control mean. Cin in Phase II is substituted as the expectation of the change-point model, which is also defined as the CL of established control chart. Therefore, Vwithin in Phase I is r n Vwithin =J^(Xi~Cin)2 + ^ (12) i=l i=T+l where Xt is the average value of i-th subgroup of observations, and c-v£i c - y -ÍÍ in / _ ' °out / ¿—i T n — ¿=1 ¿=T+1 (13) What's more, VWithin in Phase II is K L It, within —CL)2 + ^ -CoutY i=l i=T+l (14) Clusters between variation A cluster between variation is expressed as the distance between cluster centres and the total cluster centre of all observations CT, as shown in Fig. 4. Similar with the situation of clusters within variation Vwithin, CT is substituted as the CL of bootstrap control chart in Phase II, so as to obtain the expression of Vbetween in Phase I as ^between ~T(Çin ^V)2 + (n T)(C0Ut CT)2 where CT = Xt/n and Vbetween in Phase II is ^between = T(.Cin~CL)2 + (n — T )(C0Ut — CL)2 (15) (16) UCL LCL 0 20 40 60 80 100 120 140 160 180 200 Sample number Fig. 3 An example of established bootstrap control chart UCL 100 120 140 160 180 200 Failure number 1st to (t-1)th sample as in-control cluster, t th to n th sample as out-of-control cluster Calculate the clusters validity indices Calculate Clusters Total Variation as the objective function The sequential clustering approach No Fig. 4 Illustration of clusters within and between variation Identify the optimal clustering setting by minimizing the objective function Fit the proposed sectional NHPP model with in- and out-of- control observations Obtain the accurate change point \ \through maximum likelihood estimation mle ! Fig. 5 Change-point estimation procedure The objective function Clusters total variation Vtotal is defined as the objective function. It integrates the within and between variations for the purpose of valuating different cluster settings more completely and getting the optimal clusters more efficiently. The objective function has the expression as Vtotal =V\vithin + ^between [17] For the classification in Phase I, it is derived from Eq. 12 and Eq. 15 r n Vtotal = _Ci")2 + Z & ~Cout)2 -r(Cin-CT)2 + (n- T)(C0Ut ~CT)2 [18] i=l i=T+l What's more, for the classification in Phas e I I, Eq. 14 and Eq. 16 is combined as r n Vtotal = -CL)2 + ^ & -Couty -r(Cin-CL)2 + (n- r)(C0Ut -CL)2 [19] i=l i=T+l To find out the optimal transition point t and assign observations to in- and out-of-control clusters, the proposed objective function should be minimised. f = argmin Vtotai(r) [20] 3.3 Change-point estimation procedure The overall structure of proposed change-point estimation procedure is summarised as Fig. 5. First of all, the bootstrap SPC chart is constructed,possible combinations of in-and out-of-control clusters are classified by each out-of-control signal. Then, the sequential clustering approach is adopted, proposed validity indices are examined sequentially for each clustering setting and the objective function is optimised to obtain the best assignment of observations. Finally, the best assignment of observations is used to fit the proposed two sectional NHPP model, the more precise change-point estimator could be calculated through maximum likelihood estimation. As for the part of maximum likelihood estimation, let f(t1, t2, ..., t^,... ,tn) denotes the joint probability density of failure times 0 < t1 < ••• < t? < ••• < tn, where n is the random number of failures. Among them,tx, t2,..., belong to the optimal in-control cluster, and t?, ti+1, ...,tn belong to the optimal out-of-control cluster. The joint probability density can be constructed with a failure intensity u>(t) in Eq. 7 and a cumulative failure intensity W(t) in Eq. 8. iaift)*rr t^e-^, o). tQk is the start time of vehicle k from the distribution centre, k = 1,2,..., |K"|. • Mixed time window penalty costs To improve the quality of distribution services, customers require the distributor arrives within specified time windows; if not, they need to pay the corresponding wait or delay costs. The quality of fresh goods is sensitive to time. Any vehicle that arrives early has to wait until the beginning of the time windows. Any vehicle that arrives late will incur a penalty, and the delay costs of damaged goods are serious. Therefore, this study uses mixed time windows to measure fresh good distribution penalty costs. Thus, the penalty function can be shown in the formula: aRikikk ~LTi) Tie tik ,tik Pi(fik) represents penalty costs if vehicle k transgresses the time window of customer i. The lower bound Tie represents the earliest arrival time that a customer can endure when a service starts earlier than ETi. Similarly, the upper bound Ta represents the latest arrival time that the customer can endure when the service starts later than L^. Here, a represents the penalty coefficient that is within the actual time but earlier than the optimal satisfactory time and ft represents the penalty coefficient that is within the actual time but later than the optimal satisfactory time. a > 1, ft > 1. The corresponding penalty function is shown in Fig.1. Penalty^ k la aq¡k (ET - sik ) q,t - LT Y T„ ET: LTt T, Time Fig. 1 Relationship between arrival time, time-windows and penalty cost • Average freshness The products average freshness is related to the quality and freshness factor. Thus, the average freshness is defined as: t\ ¿2 ^ qikr(tik)yik/^qik ÍEN kEK i = 0 (5) Then, the mathematical model MO-VRPMTW-P can be given as follows. Objective functions: maxZ1= ^ Ydfx0]k+ ^ c0dijYjXijk+wYJYJqik(l~e r0, ViEN (14) tjk = (tik + di]/v0)xi]k, Vi,j EN+,kEK (15) xijk E {0,1}, i,jEN+,kEK (16) In the first objective function Eq. 6, total costs are minimized, namely: the fixed costs, transportation costs, damage costs, and penalty costs. In the second objective function Eq. 7, the average remaining freshness of products to be delivered is maximized. Eq. 8 is the number of vehicles constraint. Eq. 9 states that each vehicle should leave and return to distribution centre. Eq. 10 is a vehicle capacity constraint. Eqs. 11 and 12 represent that each customer should be serviced, and each can only be serviced once time. Specially, Eq.13 defines freshness function of the perishable products, and Eq. 14 ensures the lowest level of freshness that the customer can accept. Eq.15 defines the time that vehicle k takes from leaving customer i to arrival at customer j. Eq. 16 represents that vehicle k serves customer i before customer j. 3. Used methods The MO-VRPMTW-P problem is an NP-hard problem. Multiple objectives need to be optimized at the same time. In the combination optimization problem, GA is an efficient global optimal algorithm and VNS is an efficient local search algorithm. In this work, we combine the improved GA with the VNS algorithm. Traditional algorithms mainly consider the customer spatial location relationship; however, they rarely take orders the time and space characteristic constraints of orders into consideration. Since the orders have obvious ST characteristics, we use the k-means method to cluster the nodes to obtain the initial solution considering the ST strategy in the first stage. Then, in the second stage, we adopt VNSGA to optimize the distribution route. 3.1 Generate initial solution Calculating spatio-temporal distance In the process of delivery, each perishable order has a corresponding demand. Considering orders with spatio-temporal distance may solve the problem more effectively than just considering the distance. So, we use the definition of ST distance from the literature [26] to cluster orders. Use to denote the ST distance. D?, Dfj represent the Euclidean distance and temporal distance between customer i and j, respectively. The transportation time of the points is related to the Euclidean distance, which means =ttj. Here, [a, b] and [c,d\ are the time windows of customer i and j, the specific arrival time at customer j is t' E(a',b'), a' = a + t^, b' = b + ttj. Use Savij^t') to denote the saved time when vehicle arrives at customer j at the moment t'. Here, A is the maximum window, and Kt, K2, K3 are parameters related to time. k2t' + k1d- (k1 + k2)c t'd A greater Sav^^t') means a smaller spatial distance. Dfj(t') is defined as: Dfj(t') = kxA —Savijit') t' E(a',b') (18) Temporal distance Dt] is defined as: -b' Dïj=( DTj(t')dt'/(b'-a') J n' = k1A - I (k2t' + k1d — (k1 + k2)c)dt' Jmin(a',c) (19) r max(min(b',d),c) + 1 (~k1t' + k1d)dt' Vin(ma x(a',c),d) rmin(b',d) + 1 (k3f+ k3d)dt'/(b'-a') ~'min(a/,d) We take the maximum distance as the temporal distance. Dfj = ma(WjM) (20) The ST distance is related to and Dfj. a1 , a2 are the and Dfj weight coefficients, respectively. a1 + a2 = 1. The ST distance can be expressed as follow: i Dîj - min (°mn) \ ( Djj - min (D£n) | m,nEC,m±n I . I J m,nEC i ,-„„ Df/ =a11 -=-=— I + a2 I-=-=— I (21) J \ max (0mn)- min (AnJ/ \ max (Aim)" min(Aim) \m,nEC,m^n m,nEC / \m,nEC,m±n m,nEC Construct initial solution After the calculation of the ST distance, we apply k-means method to cluster the orders and construct initial solution. Here, k is the number of vehicles. The orders are divided into k clusters, and the clustering centre zi(z1,z2, ...zk) of each cluster is The cluster k is defined as follow: k Min£ ^ D^ (22) 7 = 1 iEzi/{oi] k = maxY,iEN di/Q. E^i represents the total demand of the largest distribution order, and D(i, Oj) represents the distance of the sub-order i to the cluster centre zt. 3.2 Optimization solution based on VNSGA We combine the improved GA and the VNS to solve the multi-objective problem. Owing to the different mechanisms used in population search and local search, we use two different fitness functions in the selection operations. We apply non-dominance ranking and crowding distance sorting in the NSGA-II method as a global strategy and adopt an external archive to maintain the process of evolution. Then, VNS realizes the dynamic search. The flow chart is shown in Fig. 2. Integrating orders and coding chromosome (G^G^Gs **■,GK_) Generating population, size of N Initial the NDset I Initial external archive -?- Select a individual from the current population based on non-dominance ranking and crowding distance Select individuals from the population Crossover and mutation operation Variable neighborhood search Generate N new individuals Update external archive Satisfying iteration termination ? X Y Output the end Fig. 2 Flowchart of the solution algorithm Fitness function We adopt two different fitness functions in population and local search operation selection. First, in the population selection operation, we propose the ranking and crowding degree method based on the Pareto dominance. fit1{x) is the population selection fitness. Use rank{x) to denote the relationship of Pareto dominance, rank(x) < rank(y) indicates that x dominates y. The crowding distance is Crowding — distance(x). The first keyword is ascending in a particular order and the second key sort is descending according to the crowding distance. fit1(x) = (rank(x), Crowding — distane(x)) (23) Second, in the local search selection operation, we select a better solution from the neighbourhood. S(x) is the number of solutions x dominated. W(x) is the number of solutions x dominate in storage pool. fit2(x) is the local search selection operation fitness. 1+S(x) 4. Results and discussion 4.1 Data description In this section, we use numerical experiments to demonstrate the efficiency and advantages of applying our heuristic algorithm. We adopt the instances developed by Solomon for VRPTW. Six different instance types are considered: R1, R2, C1, C2, RC1, and RC2. Suppose the distributed perishable product is milk. Set the shrinkage factor as 0 = 1/200. a = 1.5, /3 = 1.5, = 0.5, a2 = 0.5. These experiments were performed on a personal computer with Intel® CoreTM i5-4460 CPU at 2.40 GHz and 8.00 GB of RAM. The computation run time unit is seconds. The stopping criterion is set to Maxt = 300. The parameters are listed in Table3. Table 3 Parameters of the experiments. Parameter Meaning Value V The speed of the vehicle (km/h) 30 T The lifecycle of the perishable product (h) 24 f The unit fixed costs per of the vehicle (yuan) 50 Q The maximum capacity of the vehicle (kg) 300 Co The average cost of travelling (yuan/km) 2.5 r0 The minimum freshness level that the customer can accept 0.75 w The unit value of the perishable products(yuan/kg) 30 K A set of vehicles 50 4.2 Effectiveness of considering spatio-temporal distance Consider the time and space characteristics of the order, we propose a method based on the spatio-temporal metrics to verify the strategy of spatio-temporal distance. We define the method that considers the customer spatio-temporal location relationship as ST-VNSGA and the method that does not consider the customer spatial location as VNSGA. Comparison results between ST-VNSGA and VNSGA are given in Table 4. Table 4 Comparison results between ST-VNSGA and VNSGA ST-VNSGA VNSGA Gap Case Time Cost Freshness Time Cost Freshness Time Cost Freshness (s) (yuan) (%) (s) (yuan) (%) (%) (%) (%) R101_25 43 771.58 93.4 73 791.24 92.1 41.10 2.48 1.41 R101_50 79 1403.31 91.7 144 1449.73 88.6 45.14 3.20 3.50 R101_100 127 2649.43 88.1 278 2784.37 85.7 54.31 4.84 2.80 C101_25 22 185.36 94.5 26 191.11 92.9 15.38 3.00 1.72 C101_50 44 367.98 92.8 54 379.45 91.7 18.52 3.02 1.20 C101_100 79 730.91 92.3 99 758.69 90.0 20.2 3.66 2.56 RC101_25 29 799.72 93.3 49 819.43 91.6 40.82 2.41 1.86 RC101_50 52 1360.98 91.4 94 1398.21 89.0 44.68 2.66 2.70 RC101_100 91 2556.67 87.5 189 2678.85 84.9 51.85 4.56 3.06 60.00% 50.00% o 40.00% -H .p .2 30.00% 4" ■I 20.00% a ° 10.00% QJ Ë 0.00% s^v^o - cvV O.V , O" Vy Ov Qv N □ R101_25 □ R101_50 □ R101_100 □ C101_25 □ C101_50 □ C101_100 □ RC101_25 □ RC 101.50 □ RC101_100 ^v ^ ^ C Fig. 3 Time optimization rate with the spatio-temporal strategy 6.00% 5.00% 4.00% 3.00% 2.00% 1.00% 0.00% fx ■ * cost —■• freshness R101_25 R101_50 R101_100 Fig. 4 R class case objectives ratio 4.00% 3.50% 3 00% 2.50% 1 00% 1.50% 1 00% 0.50% 0.00% freshness C10125 C101_50 C101_100 Fig. 5 C class case objectives ratio ' - 4.50% ■ 4.00% ■ 3.50% -3.00% ■ , 2.50% -2.00% -1.50% ■ 1.00% -0.50% ■ 0.00% - RC101_25 RC101_50 RC101_100 Fig. 6 RC class case objectives ratio Table 4 lists the computation results of ST-VNSGA and VNSGA with regard to run time, cost and freshness, as well as the optimization gap. Fig. 3 displays the optimization rate of the run time with regard to R, C, and RC classes. As can be seen from Table 4, ST-VNSGA can obtain a better solution in less time compared to VNSGA. For example, the ST-VNSGA run times of R101_25, R101_50, and R101_100 are 43, 79, and 127, respectively. The VNSGA run times are 73, 144, and 278, and the improvement rates of run time are 41.10 %, 45.14 %, and 54.31 %, respectively. The ST strategy can add the customer to the path where the distance is as close to the customer as possible. Thus, the ST strategy can both effectively reduce the search scope and reach a better solution faster. Fig. 4 shows R class optimization rate in cost and freshness. Similarly, Fig. 5, Fig. 6 show the results for the C class and RC class, respectively. Conclude from Table 4 and Figs. 4-6 , the cost and freshness of ST-VNSGA-P are also optimized to some extent; for example, the cost of R101_100 reduces to 4.84 % and the freshness of RC101_100 increases to 3.06 %, which proves that the considered spatio-temporal approach has certain guiding significance. At the same time, Table 4 and Fig. 3 indicate that the proposed algorithm run time is reduced greatly compared to VNSGA, especially in R and RC problems, where the run time is reduced by more than 50 %. R class problem optimization rates are 41.10 %, 45.14 %, and 54.31 %, RC class problem optimization rates are 40.82 %, 44.68 %, and 51.85 %, respectively. The customer points in C class almost appear as an aggregated distribution. The improvement in spatiotemporal distance is not obvious. However, the points in R and RC class randomly spread, so the efficiency of ST-VNSGA solution significantly improved. Thus, the strategy proposed in this study is more suitable for the optimization of the disperse region. Meanwhile, we can see from Fig. 3 that considering spatio-temporal distance has good potential in solving large-scale VRPs. For example, RC class problem run time optimization rates are 40.82 %, 44.68 %, and 51.85 %, all successively increasing. Therefore, the ST strategy optimizes obviously effect on large-scale VRPs. To validate the ST strategy optimization effect on the fresh product distribution model and algorithm, we compared without ST strategy (VNSGA) with ST-VNSGA in different order environments. Six numerical examples were created for testing. The impact of the ST strategy on run time and objective values is listed in Table 5. As indicated by Table 5, comparing the run time rate R101_100 (54.31 %) with R201_100 (46.15 %), C101_100 (20.2 %) with C201_100 (15.46 %), and RC101_100 (51.85 %) with RC201_100 (42.15 %), ST-VNSGA excels VNSGA algorithm in run time optimization, especially with the narrow time windows. In VNSGA, clustering and optimization proceed at the same time. ST-VNSGA optimizes the procedure after clusters. The strict time window constraint interferences the progress clustering. The larger time windows are, the constraints are smaller. Table 5 Comparison results of different orders environment ST-VNSGA VNSGA Gap Case Time Cost Freshness Time Cost Freshness Time Cost Freshness (s) (yuan) (%) (s) (yuan) (%) (%) (%) (%) R101_100 127 2649.43 88.1 278 2784.37 85.7 54.31 4.84 2.80 R201_100 161 1756.43 85.8 299 1802.98 83.1 46.15 2.58 3.25 C101_100 79 730.91 92.3 99 758.69 90.0 20.2 3.66 2.56 C201_100 82 368.84 90.1 97 380.16 88.8 15.46 2.98 1.46 RC101_100 91 2556.67 87.5 189 2678.85 84.9 51.85 4.56 3.06 RC201_100 199 1794.49 84.1 344 1934.89 80.5 42.15 7.26 4.47 We conclude that time windows have a major impact on the delivered perishable products; moreover, the cost and freshness of ST-VNSGA increase at different rates; for example, the cost of R201_100 decreases by 7.26 % at most, and the freshness of RC201_100 increases by 4.47 %. 4.3 Effectiveness of ST-VNSGA For a better analysis of the effectiveness of the proposed heuristic algorithm combined GA and VNS, we compare the NSGA-II algorithm (ST-NSGA-II) with ST-VNSGA. Comparison results between ST-VNSGA and ST-NSGA-II are provided in Table 6. Table 6 Comparison results between ST-VNSGA and ST-NSGA-II ST-VNSGA VNSGA Gap Case Time Cost Freshness Time Cost Freshness Time Cost Freshness (s) (yuan) (%) (s) (yuan) (%) (%) (%) (%) R101_25 43 771.58 93.4 38 789.35 91.9 13.16 2.25 1.63 R101_50 79 1403.31 91.7 69 1427.97 88.9 14.49 1.7 3.15 R101_100 127 2649.43 88.1 98 2801.66 86.3 29.59 5.43 2.09 C101_25 22 185.36 94.5 19 189.26 93.8 15.79 2.06 0.75 C101_50 44 367.98 92.8 37 383.81 90.1 18.91 4.12 3.00 C101_100 79 730.91 92.3 65 762.49 88.9 21.54 4.14 3.82 RC101_25 29 799.72 93.3 23 811.18 90.1 26.09 1.41 3.55 RC101_50 52 1360.98 91.4 43 1451.15 87.4 20.93 6.21 4.58 RC101_100 91 2556.67 87.5 72 2720.21 83.1 26.39 6.01 5.29 7.00% 0 6.00% I 5.00% ^ 4.00% I 3.00% Ä 2.00% jn ° 1.00% 0.00% i ♦ 1 M / V ' '- " ' - ■ 4 i 4 freshness v ^ ^(>V>V Fig. 7 Objective ratio between ST-VNSGA and ST-NSGA-II under different order environments We can further derive the validity of ST-VNSGA in Table 6. Result shows that for the run time in R, C, and RC classes, ST-VNSGA takes on an increasing format, for example, R101_25 (13.16 %), R101_50 (14.49 %), and R101_100 (29.59 %). Because the VNS algorithm needs to search mass neighbourhood structures, which increases the run time dramatically, and the ranking based on Pareto dominance further increases the operation time. However, we can see that ST-VNSGA ameliorates the quality of satisfactory solutions; for example, the cost optimal rate of RC101_100 is 6.01 %, the freshness optimal rate of RC101_100 is 5.29 %, indicating that the local search strategy has excellent optimization in seeking the best solution. ST-VNSGA has an advantage in solving multiple objectives, which shows that the proposed algorithm is effective and efficient. The objective value improves between ST-VNSGA and ST-NSGA-II with regard to cost and freshness, as shown in Fig. 7. In terms of the algorithm solving efficiency, the ST-VNSGA has a significant advantage, especially with the growing number of customers. 4.4 Contrastive analysis of results According to the contrastive analysis of optimization results, we have found several findings: • Compared with the method using spatial clustering, the strategy that considers the spatiotemporal distance distribution can achieve better solutions in a shorter period of time, especially in the medium or large-scale distribution problem where it can reach 54.31 % (Table 4, R101_100). This also shows that the algorithm has good potential in solving large-scale VRPs. The strategy proposed in this paper, ST-VNSGA, has obvious advantages for dispersed customer distributions. The results also provide effective decision support to solve the fresh distribution practical problem. • In fresh product distribution, the heuristic algorithm calculation run time with narrow time window constraints is better than the algorithm based on spatial distance in accordance with the actual needs of customers. In addition, it can save the cost of logistics and improve the service with regard to freshness to a certain extent. • From the results (Table 6), we can see that ST-VNSGA ameliorates the quality of satisfactory solutions, has excellent optimization in seeking the best solution, and has advantages in solving the multiple objectives, especially with the growing number of customers, with regard to cost and freshness, which show that ST-VNSGA is effective and efficient. 5. Conclusion In this study, we established the MO-VRPMTW-P model to minimize the distribution costs and maximize the freshness of perishable products. We considered the time-sensitive freshness of perishable products and the high cost of delay in mixed time windows. Then, in view of the fresh product order space and time characteristics, we designed a heuristic algorithm that considers spatio-temporal distance (ST-VNSGA) to solve the fresh product distribution problem. Several numerical examples were presented to demonstrate the effectiveness and efficiency of the proposed algorithm. It was demonstrated that these algorithms can lead to a substantial decrease in run time and major improvements in solution quality, which reveals the importance of considering a spatio-temporal strategy with mixed time windows. It is worth noting that some areas require improvement; for example, we will focus on interference management in the process of perishable product distribution in the next step. Acknowledgement This work is supported by the National Natural Science Foundation of China (Grant No. 71471025, Grant No. 71531002), Science and Technology Plan Projects of Yangling Demonstration Zone (Grant No. 2016RKX-04), China Postdoctoral Science Foundation (Grant No. 2016M600209), and China Ministry of Education Social Sciences and Humanities Research Youth Fund Project titled as "Delivery Optimization of Maturity-Based Fruit B2C Ecommerce-Taking Shaanxi Kiwifruit as the Example" (Project No. 16YJC630102). References [1] Hwang, H.-S. (1999). A food distribution model for famine relief, Computers & Industrial Engineering, Vol. 37, No. 1-2, 335-338, doi: 10.1016/S0360-8352(99)00087-X. [2] Teunter, R.H., Flapper, S.D.P. (2003). Lot-sizing for a single-stage single-product production system with rework of perishable production defectives, OR Spectrum, Vol. 25, No. 1, 85-96, doi: 10.1007/s00291-002-0105-3. [3] Amorim, P., Almada-Lobo, B. (2014). The impact of food perishability issues in the vehicle routing problem, Computers & Industrial Engineering, Vol. 67, 223-233, doi: 10.1016/j.cie.2013.11.006. [4] Song, B.D., Ko, Y.D. (2016). A vehicle routing problem of both refrigerated- and general-type vehicles for perishable food products delivery, Journal of Food Engineering, Vol. 169, 61-71, doi: 10.1016/i.ifoodeng.2015.08.027. [5] Hu, H., Zhang, Y., Zhen, L. (2017). A two-stage decomposition method on fresh product distribution problem, International Journal of Production Research, Vol. 55, No. 16, 4729-4752, doi: 10.1080/00207543.2017.1292062. [6] Tarantilis, C.D., Kiranoudis, C.T. (2001). A meta-heuristic algorithm for the efficient distribution of perishable foods, Journal of Food Engineering, Vol. 50, No. 1, 1-9, doi: 10.1016/S0260-8774(00)00187-4. [7] Tarantilis, C.D., Kiranoudis, C.T. (2002). Distribution of fresh meat, Journal of Food Engineering, Vol. 51, No. 1, 8591, doi: 10.1016/S0260-8774(01)00040-1. [8] Zhang, G., Habenicht, W., Spieß, W.E.L. (2003). Improving the structure of deep frozen and chilled food chain with tabu search procedure, Journal of Food Engineering, Vol. 60, No. 1, 67-79, doi: 10.1016/S0260-8774(03) 00019-0. [9] Faulin, J. (2003). Applying MIXALG procedure in a routing problem to optimize food product delivery, Omega, Vol. 31, No. 5, 387-395, doi: org/10.1016/S0305-0483(03)00079-3. [10] Chen, H.-K., Hsueh, C.-F., Chang, M.-S. (2009). Production scheduling and vehicle routing with time windows for perishable food products, Computers & Operations Research, Vol. 36, No. 7, 2311-2319, doi: 10.1016/i.cor.2008. 09.010. [11] Osvald, A., Zadnik Stirn, L. (2008). A vehicle routing algorithm for the distribution of fresh vegetables and similar perishable food, Journal of Food Engineering, Vol. 85, No. 2, 285-295, doi: 10.1016/i.ifoodeng.2007.07.008. [12] Hsu, C.-I., Hung, S.-F., Li, H.-C. (2007). Vehicle routing problem with time-windows for perishable food delivery, Journal of Food Engineering, Vol. 80, No. 2, 465-475, doi: 10.1016/i.ifoodeng.2006.05.029. [13] Ruan, J., Shi, Y. (2016). Monitoring and assessing fruit freshness in IOT-based e-commerce delivery using scenario analysis and interval number approaches, Information Sciences, Vol. 373, 557-570, doi: 10.1016/i.ins.2016. 07.014. [14] Rong, A., Akkerman, R., Grunow, M. (2011). An optimization approach for managing fresh food quality throughout the supply chain, International Journal of Production Economics, Vol. 131, No. 1, 421-429, doi: 10.1016/i.iipe. 2009.11.026. [15] Amorim, P., Günther, H.-O., Almada-Lobo, B. (2012). Multi-obiective integrated production and distribution planning of perishable products, International Journal of Production Economics, Vol. 138, No. 1, 89-101, doi: 10.1016/i.iipe.2012.03.005. [16] Prindezis, N., Kiranoudis, C.T., Marinos-Kouris, D. (2003). A business-to-business fleet management service provider for central food market enterprises, Journal of Food Engineering, Vol. 60, No. 2, 203-210, doi: 10.1016/ S0260-8774(03)00041-4. [17] Hasani, A., Zegordi, S.H., Nikbakhsh, E. (2012). Robust closed-loop supply chain network design for perishable goods in agile manufacturing under uncertainty, International Journal of Production Research, Vol. 50, No. 16, 4649-4669, doi: 10.1080/00207543.2011.625051. [18] Govindan, K., Jafarian, A., Khodaverdi, R., Devika, K. (2014). Two-echelon multiple-vehicle location-routing problem with time windows for optimization of sustainable supply chain network of perishable food, International Journal of Production Economics, Vol. 152, 9-28, doi: 10.1016/i.iipe.2013.12.028. [19] Claassen, G.D.H., Gerdessen, J.C., Hendrix, E.M.T., van der Vorst, J.G.A.J. (2016). On production planning and scheduling in food processing industry: Modelling non-triangular setups and product decay, Computers & Operations Research, Vol. 76, 147-154, doi: 10.1016/i.cor.2016.06.017. [20] Pahl, J., Voß, S. (2014). Integrating deterioration and lifetime constraints in production and supply chain planning: A survey, European Journal of Operational Research, Vol. 238, No. 3, 654-674, doi: 10.1016/i.eior.2014. 01.060. [21] de Keizer, M., Akkerman, R., Grunow, M., Bloemhof, J.M., Ha^ema, R., van der Vorst, J.G.A.J. (2017). Logistics network design for perishable products with heterogeneous quality decay, European Journal of Operational Research, Vol. 262, No. 2, 535-549, doi: 10.1016/i.eior.2017.03.049. [22] Albrecht, W., Steinrücke, M. (2018). Coordinating continuous-time distribution and sales planning of perishable goods with quality grades, International Journal of Production Research, Vol. 56, No. 7, 2646-2665, doi: 10.1080/ 00207543.2017.1384584. [23] Devapriya, P., Ferrell, W., Geismar, N. (2017). Integrated production and distribution scheduling with a perishable product, European Journal of Operational Research, Vol. 259, No. 3, 906-916, doi: 10.1016/i.eior.2016.09.019. [24] Moon, I., Lee, S. (2000). The effects of inflation and time-value of money on an economic order quantity model with a random product life cycle, European Journal of Operational Research, Vol. 125, No. 3, 588-601, doi: 10.1016/S0377-2217(99)00270-2. [25] Chen, J.-M., Lin, C.-S. (2002). An optimal replenishment model for inventory items with normally distributed deterioration, Production Planning & Control, Vol. 13, No. 5, 470-480, doi: 10.1080/09537280210144446. [26] Qi, M., Lin, W.-H., Li, N., Miao, L. (2012). A spatiotemporal partitioning approach for large-scale vehicle routing problems with time windows, Transportation Research Part E: Logistics and Transportation Review, Vol. 48, No. 1, 248-257, doi: 10.1016/i.tre.2011.07.001. Advances in Production Engineering & Management Volume 13 | Number 3 | September 2018 | pp 333-344 https://doi.Org/10.14743/apem2018.3.294 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper Game theoretic analysis of supply chain based on mean-variance approach under cap-and-trade policy He, L.F.a, Zhang, X.a*, Wang, Q.P.b*, Hu, C.L.c* aCollege of Management and Economics, Tianjin University, Tianjin, P.R. China bSchool of Management Science and Engineering, Hebei University of Economics and Business, Shijiazhuang, P.R. China cCollege of Intelligent Manufacturing, Tianjin Sino-German University of Applied Sciences, Tianjin, P.R. China A B S T R A C T A R T I C L E I N F O In recent years, carbon emission problem occurred by carbon dioxide as one of the main greenhouse gases, has become the focus due to its great influence on human life. With the increase of consumers' low-carbon consciousness, this paper studies the supply chain which consists of a single supplier and a single manufacturer in presence of market low-carbon preference. First, we establish the mean-variance analysis model. Second, we study the optimal decisions of channel members considering the risk factor in three situations: traditional supply chain without emission reduction, individual emission reduction by manufacturer and supply chain collaborative emission reduction. Finally, the equilibrium results are demonstrated by numerical studies. The results show that chain members' profits are not only affected by their own risk-aversion level, but also by other chain member's risk aversion level. More important is that there exist the optimal carbon emission reduction level and profits in system collaborative emission reduction. The research makes operation mechanism of low-carbon supply chain clearer and provides a theoretical reference for supply chain members on pricing and investment strategy of emission reduction. © 2018 CPE, University of Maribor. All rights reserved. Keywords: Supply chain; Cap-and-trade policy; Carbon emission; Game theoretic analysis; Mean-variance model *Corresponding authors: feigemse@gmail.com (Zhang, X.) qinpeng@heuet.edu.cn (Wang, Q.P.) chenglinhu8@gmail.com (Hu, C.L.) Article history: Received 21 January 2018 Revised 25 August 2018 Accepted 7 September 2018 1. Introduction In recent years, the quick boosting of industry economy has brought many environmental problems, such as the greenhouse effect that threats human health seriously. In this regard, it is universally accepted that massive carbon emission is the first dominant factor resulting in greenhouse effect Dinan [1], so many countries attempt to use scientific approaches to control carbon emissions. In Kyoto Protocol, the idea of 'cap-and-trade' policy was therefore designed to control the emission of greenhouse gases, which is a significant step to mitigate the negative effect of carbon emission. It is noteworthy that controlling carbon emission in the production process is the key tache of managing low-carbon supply chain. Thereby, a lot of manufacturing enterprises pay attention to low-carbon production and implement effective emission reduction measures. Generally, controlling carbon emission is of great significance to promote the sustainable development of environment. A number of researches on low-carbon supply chain have focused on performance evaluation and management, such as, Yin et al. [2], Chen [3]. For in-depth research, the studies on low carbon supply chain can be divided into three categories: carbon footprint distribution, e.g. Wang et al. [4], Yang et al. [5], production and operation decision, e.g. Wang et al. [6], Nie et al. [7], Zhao et al. [8]; supply chain structure under carbon emission constraints, e.g. Cholette et al. [9], Du et al. [10]. For ones interested in low-carbon supply chain, please refer to Jharkharia [11] to learn more progress. One of main mathematical concepts to characterize our study is the use of mean variance. Chen et al. [12] studied the basic inventory model based on the newsvendor model by applying mean variance model. Wu et al. [13] researched newsvendor model with mean-variance analysis considering stockout cost, and the results reveal that the risk-averse newsvendor may even order more products than the risk-neutral. Ray et al. [14] conducted the research on purchasing strategy by mean-variance analysis considering disruption risk, which is the first time applying this approach to manage supply chain under disruption. Considering different risk attitudes, Choi et al. [15] investigated the newsvendor problem by utilizing mean-variance approach. Similarly, Yamaguchi et al. [16] developed mean-variance model, discussed the optimal order quantity following the three risk attitudes: risk-neutral attitude, risk-averse attitude and risk-prone attitude. Zhuo et al. [17] explored the two-echelon supply chain with option contracts under the mean-variance framework. More information about supply chain with mean-variance, please see Chiu et al. [18]. There is an increasing body of literature on cap-and-trade policy, which is relevant to the present study. In 1997 the cap-and-trade system was put forward in the Kyoto Protocol, the cap is proposed by government, and green organization has become a supplier of carbon emission. Different from traditional supply chain, the carbon emission permit can be traded via carbon market so that the carbon emission is utilized and controlled effectively. Considering stochastic demand, Zhang et al. [19] addressed the optimal production strategy of manufacturer in cap-and-trade system. Furthermore, they investigated purification efficiency in the case of multitime. Du et al. [20] explored the low carbon supply chain consisting of one manufacturer who results in massive carbon emission and single emission permit supplier in a single period. Du et al. [21] discussed the carbon emission policy by utilizing the Stackelberg game in cap-and-trade system, furthermore, investigated the impact of cap on manufacturer, supplier and supply chain. Du et al. [22] constructed optimal production model for manufacturer considering environmental and preference with fixed cap and price of emission permit. According to the sizes of trade price and purification cost, the manufacturer decides whether to invest purification cost or trade emission permits via market for maximizing the profit. Zhao et al. [23] researched the problem of coordination mechanism and design of supply chain consisting of a sing manufacturer and a sing retailer considering the cap of carbon emission in a low carbon environment. Furthermore, Yuan et al. [24] studied the coordination of supply chain with revenue sharing contract in cap-and-trade system. In view of quick response strategy under newsvendor setting, Lee et al. [25] investigated the optimal pricing decisions of the company and policy maker who proposed the cap and trading price. Our study is also related to consumers' preference to low-carbon product, such as Ma et al. [26] and Wang et al. [27]. Similarly, Zhang [28] assumed that demand of products is a linear function with product prices and carbon emission and explored the optimal decisions of channel members in three situations: traditional supply chain without emission reduction, manufacturer's individual emission reduction and collaborative emission reduction. Similar researches are also conducted in He et al. [29] and He et al. [30]. However, the risk factor of channel members wasn't considered in the literature. Choi et al. [15] discussed maximum profit of supply chain under risk free constraint and risk constraint, and the result showed that supply chain risk actively affects the decision-making of supply chain members, so it cannot be ignored. Differently, this paper simultaneously considers the mean and variance to investigate the problem of low carbon supply chain, which makes the research more consistent with the reality, and ensure rigor of the research problem. Comprehensively, the extant literature on low carbon supply chain mainly focus on optimization in firm and supply chain level but without incorporating risk concerns, which is exactly the gap we attempt to fill. The Mean-variance model considers both mean and variance to capture risk concerns, which is close to the reality. However, the existing literature on low carbon supply chain solely takes mean into account. In response to the enhancement of people's awareness of low carbon, countries are also actively taking action to reduce carbon emission, such as, the designing cap-and-trading policy. To echo above problems, this paper exploits the mean-variance analysis to construct the game of supply chain under cap-and-trade policy, then the optimal decisions of chain members are explored and some numerical examples are provided to illustrate the models. 2. Preliminaries and notation The necessary assumptions used in this paper are listed as follows: • The market information is symmetrical and complete, and the market can completely be cleared. Thereby, there is no excess supply or excess demand, and the Pareto optimal of equilibrium price is ensured. • Assume a unit material from the supplier can produce a unit product by the manufacturer. • According to the study of Abdallah et al. [31], we assume that the government has implemented a more market-oriented economic means, namely, cap-and-trade system. In order to implement the strategy fairer and effectiveness, the government sets the carbon cap of per unit product which is assumed as Xg. According to the international practice and combined with China's national conditions, the carbon trading price pc is allowed in the carbon emission trading market. • Assume the demand is linearly decreasing in retail price but increasing in the low-carbon preference with random interruption. We therefore give the liner function D = a — /3p + ye, where a is intrinsic demand, + Cm)■ • According to the standard hypothesis in classical model, namely, there is a quadratic function relationship between cost and R & D investment), we use C(e) = de2/2 to represent emission reduction cost, where d is the parameter of emission reduction cost. Zhang [28] and Raz et al. [32] also use the function to represent the relationship between carbon emission reduction level and emission reduction cost. The notations used in this paper are summarized in Table 1. Table 1 Notation summarization Notation Meaning P Unit retail price Pc Unit trading price of carbon emission D The demand of product Cm Unit cost of manufacturer Cs Unit cost of supplier w Unit wholesale price eo Initial carbon emission of unit product Xg Carbon emission cap of unit product e Parameter of emission reduction cost The risk-averse level of manufacturer As The risk-averse level of supplier 3. Optimal decision of supply chain considering risk attitude 3.1 Optimal decision-making in traditional supply chain without emission reduction In cap-and-trade system, if the manufacturer want to keep his traditional production without emission reduction, he must purchase carbon permits to meet the cap set by the government. In this situation, the manufacturer's profit is: nm = (p-w- 00 + f-Pp) + Pc&g -e0)(a + ^ ~Pp) The first term in the right hand side (RHS) of above function represents the sales revenue of manufacturer, and the second denotes the cost for purchasing carbon emission permit In addition, we assume that the cap proposed by the government is less than the initial carbon emission of per unit product Ag (y — ßPc)2, there are the optimal solutions of Eq. 12. Let the Eqs. 14 and 15 equal zero, then the optimal decisions of manufacturer can be expressed as follows .. [A - PcO + ffpc)]Q- Ama) — (y2 + ßpcy~ ß9)[w + cm -pc(Ag -g0)] 20(0- 2pcy) -(y- f3pc)2 = (/ + PVc)[a - Ama-p(w + cm-pc{Xg -g0))] 61 ~ 2(0.9-2pcy)-(y-ppcy ( J Substituting the Eqs. 16 and 17 into the Eq. 13, we can rewrite the utility function of supplier as below Us = 2p(d-2^JcY~)C-(y-pPc)^[p6[a ~ P(W + c™ + + PVcYUm° (18J ' ' -[20(0- 2Vcy) -(y- (3pc)2]Asa} Solving the first derivative and second derivative with wholesale price of Eq. 18, we can obtain dUs w-cs -WV- ^ - (y—^m*} - w-ïS-b-tef d2Us _ 2/329 d^ï =~ 2(3(9- 2pcy) -(y- (3pc)2 < 0 Let Eq. 19 equal zero, the optimal wholesale price is expressed as below 2* a cm-cs-Vc(Xg-e0) \fie-{y + PVc)2}Ama-[2p{e-2VcY)-(Y-pVc)2]As w =------I--(20) 2p 2 2p28 1 J Finally, substituting the Eq. 20 into Eqs. 16 and 17, we rewrite the optimal solutions as below 2* _Y(Y + PPc)-P0_ for .of . „ ^ V =- 2^2(3(9- 2pcy) -(y- 0pc)2] ^ + ^ + ~g°»l + [a + p(cm + cs-pc(2(3(9- 2pcy) - (y- ppc)2))]Ama (21J =_Y + PPc_ a-p[cm + cs-pc(2p(9- 2pcy) - (y- 0pc)2)] 61 20(0 — 2pcy) — (y — PPc)2 2 (22J [(300 ~(y + PPc)2)]°~ [2(3(9- 2pcy) - (y- 0pc)2Ra L J 2(39 ' 3.3 Optimal decision-making in channel collaborative emission reduction In order to improve the market competitiveness and get more profits, manufacturer conducts low carbon production considering consumers' preference of low carbon products. Then the market demand has been promoted with the behavior of manufacturers' emission reduction. Finally, the profit of supplier has been improved with more materials ordered by manufacturer. Thereby, the rational supplier proposes that he is willing to cooperate with manufacturer and bear a part of carbon emission reduction cost. Under the mode of cooperative emission reduction, what is the changes in decision-making between supply chain members? Does the collaborative emission reduction between the supply chain members maximize the economic benefits of each other? This section analyzes these problems. First, we get the chain members' profits as follows. 1 nm = (p-w- cm)D + pc(Ag -e0 + e)D- ^(1 - S)Qe2 1 ns = (w- Cs^D-^SOe2 Where the parameter 5 denotes the emission reduction cost sharing proportion of supplier. The expected profits of chain members are expressed as follows 1 Enm = [p-w- cm + pc(Ag — e0 + g)](a - ftp + ye) - ^(1 - S)6e2 (23] 1 Ens = (w-cs)(a -/3p+ ye)-^Sde2 (24) The first term of RHS in Eq. 23 represents the revenues consisting of both sales revenue and carbon emission trading revenue or cost and the second term the carbon emission reduction cost undertaken by manufacturer. Similarly, the first term of Eq. 24 expresses the sales revenue, and the last term the carbon emission reduction cost undertaken by supplier. Then we can get the utility functions of chain members as follows MV: Um = [p-w-cm + pc(Ag -e0 + g)](a -f3p + ye- Ama) - i(1 - 8)6e2 (25) 1 MV: Us = (w- cs)(a - pp + ye - Asa) -~8Be2 (26) Similarly, we get the first partial derivatives of Um with respect to the retail price and carbon emission reduction level as follows da m = -2ßp + ß[w + cm-pcaq —e0)] + a + yet- Ama [27) H¿ = dp dUm = y[p-w-cm+ pc(Ag -e0 + e)] + pc(a-^p + ye- Ama) - (1 - S)9e [28) Solving the Hessian matrix of the Eq. 25 in terms of the retail price and carbon emission reduction level, we get " -20 Y-PPc -Y-PPc 2pcy-(1-S)6_ When 20[(1 — 5)6 — 2pcy] >(y — fipc)2, there are the optimal solutions of Eq. 25. Let the Eqs. 27 and 28 equal zero, the optimal decisions of manufacturer are expressed as follows = [(1- Ô)6-Pc(y + /?pc)](q - Ama) - [y2 + (3pcy-(3{1-8)6][Ag -e0 + g] [2g) = (/ + PPc)[a - Ama-p(w + cm-pc{Xg -g0))] [30) 2$\{\-8)e-2pcy\-{y-$pc)2 [ ) Substituting the Eqs. 29 and 30 into the Eq. 26, we can rewrite the utility function of supplier as below w — cq 9 Us = „--;rT2 {Pd(1-8)[Xg -eo] + [£9(1-S)-(y + ppc)2]Ama 2ß(_8-2pcy)-(y-ßpcy , „„ n—Rfu, + r —n fi [31) [ßea-5)-(y + ßPc)2]Ama}-28e"r 22 ^ 13Q(Y + ßPc) [a-Amo-ß(w + cm-pc(Äg-eo))Y 2 [2ß(y- ßPc)-(y-ßPc)2f Solving the second derivative with wholesale price of Eq. 31, we can get d2us =__2/329(1-S)___p2es(y + pPc)2_ dw2 2£[(1-50)-2pcr]-(7-£pc)2 [2^((1-50)-2pcr)-(r-^Pc)2]2 Therefore, the optimal wholesale price is generated by letting the first derivative with wholesale price of Eq. 31 equal zero 2. a cm-cs-Vc(Xg —e0) [ßS - (y + ßpc)2]Ama- [2ß{9- 2pcy) -(y- ßpc)2]A W — -------r + (32) 2P 2 2p2e [(2^(1- S)0- 2pey) -(y- pPc)2)2 8)9- 2pey) -(y- jSpe)2)+ 2jS2e25(1-5)]^mff ^20[2(2^((1- 5)0- 2pcy) -(y- ^pc)2) + 5(y + pPcy] Finally, substituting the Eq. 32 into the Eqs. 29 and 30, we get 3. = [(1-S)d-pc(y + pPcMa - Ama) - [y2 + pPcy- p(1-S)9][w3* + cm-pc(Ag -e0)] P 2p[(1-5)6-2pcy]-(y-ppc)2 ( J = (/ + PVc)[a - Ama~P(w3* + cm-pc(Ag -g0))] (34) 2^[(1 — 5)6 — 2pcy] ~(y — PPc)2 It is difficult to elegantly express the p3* and e in explicit format. Thereby, we substitute w3* into Eqs. 29 and 30 to represent the equilibrium results. Further research on the problem will be carried out in the numerical analysis section. 4. Computational study and discussion In this section, we conduct a series of numerical computation to demonstrate previous models and obtain some managerial insights. Proposition 1: The carbon emission reduction level in chain members' reduction cooperation is greater than that in emission reduction by manufacturer solely. The result reveals that the reduction cooperation of chain members can not only meet the low carbon preference of consumers, but promote the environment sustainable development Additionally, the carbon emission reduction level is increasing in supplier's risk-averse level, while decreasing in manufacturer's risk-averse level, as shown in Fig. 1 and Fig. 2, respectively. the reason resulting in this phenomenon can be explained like that the risk-averse manufacturer wouldn't like to input more costs to conduct low-carbon production, while the risk-averse supplier is just opposite. 4.5'L 3.5 ■ 2.5 0.2 0.4 0.6 A 0.8 0.4 0.6 A Fig. 1 The carbon emission reduction level with versus Fig. 2 The carbon emission reduction level with versus Aq variation Am variation Proposition 2: The retail price is decreasing in risk aversion level of supply chain members. As shown in Figs. 3 and 4, we can explain the phenomenon in this way. The manufacturer reduces retailer price to attract price sensitive consumers for getting high profit. Furthermore, observing Fig. 3, we find there are different relationships among the three retail prices with the Fig. 3 The retail price with versus v4s variation Fig. 4 The retail price with versus Am variation variation of supplier's risk aversion level. While the retail price in collaborative emission reduction is the lowest price no matter how manufacturer's risk aversion level changes, as shown in Fig. 4. Proposition 3: The wholesale price is closely related to the risk aversion level of chain members. It is decreasing in supplier's risk aversion level, while it shows opposite result in manufacturer's risk aversion level, as shown in Figs. 5 and 6. The reason resulting in this phenomenon is that the risk-averse supplier reduces the wholesale price in order to promote the order of manufacturer. For risk-averse manufacturer, as we described aforesaid, he reduces retail price aim to stimulate consumption with more order occurred. Finally, the upstream supplier takes this chance to improve wholesale price. Fig. 5 The wholesale price with versus variation Fig. 6 The wholesale price with versus Am variation Proposition 4: The profits of chain members are decreasing in their own risk aversion levels. In addition, there are always the optimal profits of chain members in emission reduction cooperation when supplier keeps the certain risk aversion level. Observing Figs. 7, 8, 9, and 10, the chain members all get the optimal profits in collaborative emission reduction when the risk aversion level of supplier is less than 0.4. Additionally, risk-averse chain members' profits are increasing in their own risk aversion levels. the cause of this kind of appearance is that risk-averse chain members decreases the prices (wholesale price or retail price) to stimulate consumption for maximizing the profits, nevertheless, the profits are still decreasing. Fig. 7 The profit of manufacturer with versus variation Fig. 8 The profit of manufacturer with Am variation A= A Fig. 9 The profit of supplier with versus variation Fig. 10 The profit of supplier with Am variation Proposition 5: The profits of chain members are increasing in the cap imposed by government. The cap proposed by the government directly impact on the economic benefits of the supply chain members in the risk-neutral supply chain. With increase of the cap, the economic benefits of supply chain members have been improved, and it is more beneficial to the chain members to cooperate with each other about carbon emission. The result also reveals another information that the cap provided by the government can be directly converted into economic income for the chain members. As shown in Table 2. The parameters are summarized in Table 3. Table 2 The profits of chain members xq uL y2 um u3 um ui 5 7528.8 7803.5 8071.2 15057.5 15607.0 15660.4 6 7600.2 7877.6 8147.4 15200.4 15755.1 15809.1 7 7672.0 7952.0 8224.0 15344.0 15903.9 15958.4 8 7744.1 8026.7 8301.0 15488.3 16053.5 16108.5 Table 3 Parameter Parameter a ß Y cs e a Pc e0 S xq Value 200 0.3 0.4 10 8 80 10 3 10 0.4 5 5. Conclusion Considering consumers' low-carbon preference and the risk concerns of chain members, this study exploits mean-variance approach to explore the optimal decisions of supply chain agents. Three different cases have been discussed: traditional supply chain without emission reduction, emission reduction by manufacturer individually and chain members' collaborative emission reduction.Some meaningful conclusions are obtained through analytical analysis and computa- tional study. The results show the emission reduction level in collaborative emission reduction is greater than that by manufacturer individually. Moreover, there is an optimal profits of supply chain members in collaborative emission reduction, which indicates that in the cap-and-trade system, collaborative emission reduction of chain members can not only increases the profits of chain members, but also improves the emission reduction level. Additionally, the optimal decisions are affected by both manufacturer's risk level and supplier's risk level. The research provides a cooperative emission-reduction way for chain members constrained by low-carbon policy. The significance of this paper is reflected in three aspects. First, considering the chain member's risk, we explore the equilibrium results in three cases as described aforesaid, and figure out the operational mechanism for a risk-aversion supply chain in the cap-and-trade policy. Second, we study the influence of risk attitude on the profits of chain members, which reveals that the level of risk aversion plays an important role in the decision-making of supply chain members. The research compensates the defects of existing literature that assumes supply chain members are risk-natural. Finally, we consider the risk of chain members and construct mean-variance analysis model to let the research close to the reality. The proposed models can be extended in several ways. For instance, extending the model to information asymmetry, which is common and important in reality. Furthermore, the study may be extended to more-than-two-tier supply chain. Additionally, studying the low carbon supply chain in multiple periods may be also a potential direction in the future. Acknowledgment The authors sincerely appreciate the anonymous referees and editors for their time and patience devoted to the review of this paper as well as their constructive comments and helpful suggestions. This work is partially supported by NSFC Grants (No. 91646118, 71701144, 71602142, 71528002). It is also supported by Natural Science Foundation of Hebei Province of China under Grant No. G2017207015 and Humanity and Social Science Youth Foundation of Ministry of Education of China under Grant No. 18YJC630183. Reference [1] Dinan, T. (2008). Policy options for reducing CO2 emissions, Washington, D.C., The Congress of the United States, Congressional Budget Office. [2] Yin, Z., Zhang, X. (2014). A review of the research on low carbon supply chain in the context of open economy, Inquiry Into Economic Issues, Vol. 9, 154-159, doi: 10.3969/j.issn.1006-2912.2014.09.025. [3] Chen, J. (2012). Study on supply chain management in a low-carbon era, Journal of Systems & Management, Vol. 21, Vol. 6, 721-728, doi: 10.3969/j.issn.1005-2542.2012.06.002. [4] Wang, L., Zheng, Y., Gao, Y. (2014). Study on multi stage carbon footprint optimization of manufacturing supply chain, Journal of the Party School of the Central Committee of the C.P.C., Vol. 18, No. 4, 92-97, doi: 10.14119/j.cnki. zgxb.2014.04.057. [5] Yang, C., Li, X., Shao, J. (2013). Optimization of carbon footprint of reverse supply chain in complex uncertain environment, Low-carbon Economy, Vol. 4, doi: 10.13529/j.cnki.enterprise.economy.2013.04.007. [6] Wang, W.-B., Deng, W.-W., Bai, T., Da, Q.-L., Nie, R. (2016). Design the reward-penalty mechanism for reverse supply chains based on manufacturers' competition and carbon footprint constraints, Journal of Industrial Engineering and Engineering Management, Vol. 30, No. 2, 188-194, doi: 10.13587/j.cnki.jieem.2016.02.023. [7] Nie, J.-J., Wang, T., Zhao, Y.-X., Zahng, L.-N. (2015). Collecting strategies for the closed-loop supply chain with remanufacturing in the constraint of carbon emission, Journal of Industrial Engineering and Engineering Management, Vol. 29, No. 3, 249-256, doi: 10.13587/j.cnki.jieem.2015.03.028. [8] Zhao, D.-Z., Yuan, B.-Y., Xia, L.-J., Xie, X.-P. (2014). Dynamic game study in supply chain with manufacturers' competition under the constraint of productions' emission, Industrial Engineering & Management, Vol. 1, No. 8, 65-71, doi: 10.3969/j.issn.1007-5429.2014.01.011. [9] Cholette, S., Venkat, K. (2009). The energy and carbon intensity of wine distribution: A study of logistical options for delivering wine to consumers, Journal of Cleaner Production, Vol. 17, No. 16, 1401-1413, doi: 10.1016/j. jclepro.2009.05.011. [10] Du, S., Ma, F., Fu, Z., Zhu, L., Zhang, J. (2011). Game-theoretic analysis for an emission-dependent supply chain in a 'cap-and-trade'system, Annals of Operations Research, Vol. 228, No. 1, 135-149, doi: 10.1007/s10479-011-0964-6. [11] Das, C., Jharkharia, S. (2018). Low carbon supply chain: A state-of-the-art literature review, Journal of Manufacturing Technology Management, Vol. 29, No. 2, 398-428, doi: 10.1108/IMTM-09-2017-0188. [12] Chen, F., Federgruen, A. (2000). Mean-variance analysis of basic inventory models, Working Paper, Columbia Business School, New York, USA. [13] Wu, J., Li, J., Wang, S., Cheng, T.C.E. (2008). Mean-variance analysis of the newsvendor model with stockout cost, Omega, Vol. 37, No. 3, 724-730, doi: 10.1016/j.omega.2008.02.005. [14] Ray, P., Jenamani, M. (2016). Mean-variance analysis of sourcing decision under disruption risk, European Journal of Operational Research, Vol. 250, No. 2, 679-689, doi: 10.1016/j.ejor.2015.09.028. [15] Choi, T.-M., Li, D., Yan, H. (2008). Mean-variance analysis of a single supplier and retailer supply chain under a returns policy, European Journal of Operational Research, Vol. 184, No. 1, 356-376, doi: 10.1016/j.ejor.2006. 10.051. [16] Yamaguchi, S., Goto, H., Kusukawa, E. (2017). Mean-variance analysis for optimal operation and supply chain coordination in a green supply chain, Industrial Engineeering & Management Systems, Vol. 16, No. 1, 22-43, doi: 10.7232/iems.2017.16.1.022. [17] Zhuo, W., Shao, L., Yang, H. (2018). Mean-variance analysis of option contracts in a two-echelon supply chain, European Journal of Operational Research, Vol. 271, No. 2, 535-547, doi: 10.1016/j.ejor.2018.05.033. [18] Chiu, C.-H., Choi, T.-M. (2016). Supply chain risk analysis with mean-variance models: A technical review, Annals of Operations Research, Vol. 240, No. 2, 489-507, doi: 10.1007/s10479-013-1386-4. [19] Zhang, J.-J., Nie, T.-F., Du, S.-F. (2011). Optimal emission-dependent production policy with stochastic demand, International Journal of Society Systems Science, Vol. 3, No. 1-2, 21-39, doi: 10.1504/HsSs.2011.038931. [20] Du, S., Ma, F., Fu, Z., Zhu, L., Zhang, J. (2011). Game-theoretic analysis for an emission-dependent supply chain in a 'cap-and-trade' system, Annals of Operations Research, Vol. 228, No. 1, 135-149, doi: 10.1007/s10479-011-0964-6. [21] Du, S., Zhu, L., Liang, L., Ma, F. (2013). Emission-dependent supply chain and environment-policy-making in the 'cap-and-trade' system, Energy Policy, Vol. 57, 61-67, doi: 10.1016/j.enpol.2012.09.042. [22] Du, S., Hu, L., Song, M. (2014). Production optimization considering environmental performance and preference in the cap-and-trade system, Journal of Cleaner Production, Vol. 112, Part 2, 1600-1607, doi: 10.1016/j.jclepro. 2014.08.086. [23] Zhao, D., Wang, B., Xu, C. (2014). Study on the coordination mechanism of supply chain considering the restriction of product carbon emissions, FORECAST, Vol. 33, No. 5, 76-80. [24] Yuan, J., Ma, J., Yang, W. (2016). Revenue-sharing contract for supply chain under a cap and trade system, In: Proceedings of 2016 International Conference on Logistics, Informatics and Service Sciences (LISS), Sydney, Australia, 1-6, doi: 10.1109/LISS.2016.7854442. [25] Lee, J., Lee, M.L., Park, M. (2018). A newsboy model with quick response under sustainable carbon cap-n-trade, Sustainability, Vol. 10, No. 5, doi: 10.3390/su10051410. [26] Ma, C., Liu, X., Zhang, H., Wu, Y. (2016). A green production strategies for carbon-sensitive products with a carbon cap policy, Advances in Production Engineering & Management, Vol. 11, No. 3, 216-226, doi: 10.14743/ apem2016.3.222. [27] Wang, Q., He, L. (2018). Managing risk aversion for low-carbon supply chains with emission abatement outsourcing, International Journal of Environmental Research and Public Health, Vol. 15, No. 2, 367-387, doi: 10.3390/ijerph15020367. [28] Zhang, J. (2016). The game analysis of low carbon supply chain collaboration in emission reduction research and development, Master thesis, Hunan University, Changsha, Hunan Province, China. [29] He, L., Zhao, D., Xia, L. (2015). Game theoretic analysis of carbon emission abatement in fashion supply chains considering vertical incentives and channel structures, Sustainability, Vol. 7, No. 4, 4280-4309, doi: 10.3390/ su7044280. [30] He, L., Hu, C., Zhao, D., Lu, H., Fu, X., Li, Y. (2016). Carbon emission mitigation through regulatory policies and operations adaptation in supply chains: Theoretic developments and extensions, Natural Hazards, Vol. 84, Supplement 1, 179-207, doi: 10.1007/s11069-016-2273-5. [31] Abdallah, T., Farhat, A., Diabat, A., Kennedy, S. (2012). Green supply chains with carbon trading and environmental sourcing: Formulation and life cycle assessment, Applied Mathematical Modelling, Vol. 36, No. 9, 4271-4285, doi: 10.1016/j.apm.2011.11.056. [32] Raz, G., Druehl, C.T., Blass, V. (2013). Design for the environment: Life-cycle approach using a newsvendor model, Production and Operations Management, Vol. 22, No. 4, 940-957, doi: 10.1111/poms.12011. Advances in Production Engineering & Management Volume 13 | Number 3 | September 2018 | pp 345-357 https://doi.Org/10.14743/apem2018.3.295 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper Decision-making strategies in supply chain management with a waste-averse and stockout-averse manufacturer Jian, M.a, Wang, Y.L.a* aSchool of Transportation and Logistics, Southwest Jiaotong University, Chengdu, P.R. China A B S T R A C T A R T I C L E I N F O Behavioral preferences is an important factor that affects the decision-making strategies of enterprises. Usually, the behavioral preferences will lead to decision-making that deviates from profit maximization. In this study, we investigate the influence of a dominant manufacturer's behavioral preferences on decision-making and subsequent impact on profits. This study looks at the profits of the manufacturer, retailer and the system as a whole. We construct a two-stage supply chain involving a retailer and a manufacturer who may have risk-neutral (RN), stockout-aversion (54), waste-aversion (WA), and stockout- and waste-aversion (5W) preferences. Through a comparison and analysis of the four cases, we find that the manufacturer's wholesale price increases (decreases) with the 54 (WA) coefficient, while the retailer's order quantity is completely the opposite. The manufacturer's wholesale price is the highest in the WA model, followed by the RN, 54 and SW models, in that order. The retailer's order quantity is the largest and smallest in the 54 and WA models, respectively, while the size of the order quantity between the RN and SW models depends on the ratio m (the ratio of the 54 to the WA). Moreover, we also explore the changing trends of the decision-making and profits of the participants and the system profit with the degree of 54 and WA, comparing the profits of the four cases. © 2018 CPE, University of Maribor. All rights reserved. Keywords: Decision-making strategy; Supply chain management; Waste-averse preferences; Stockout-averse preferences *Corresponding author: wangyonglong@my.swjtu.edu.cn (Wang, Y.L.) Article history: Received 6 July 2018 Revised 20 August 2018 Accepted 22 August 2018 1. Introduction The supply chain decision-making is a very important problem in supply chain (SC) management and therefore is also a very difficult problem for enterprises. In practice, many retailers must look at overordering cost, underordering cost and predicted demand to determine their ordering decisions before a sales season [1]. In addition, the manufacturer also must decide the pricing strategy for its products according to its own business interests. Traditional research on decision-making strategies are based mainly on risk neutrality, as in [2-3]. However, we know from behavioral economics that decision-makers in the real world may have behavioral preferences, which is also an important reason for the existence of "decision bias" [1]. Therefore, it is more practical to study the decision-making strategies of decision-makers by integrating behavior preferences into the SC. In the existing literature, behavioral preferences research focuses mainly on loss aversion and risk aversion. A study of loss aversion by Niu et al. [4] considers the sustainability of the with alternative power structures, and the impact of the supplier's attitudes to loss on the 's sustainability and profitability. Choi [5] studies the newsboy problem, where the decision-maker is loss-averse under carbon emission constraints. When there is competition between manufacturers, research by Li and Li [6] show an ordering decision problem in the case of a supply disrup- tion. Xu et al. [7] study the optimal order quantity in the loss-averse newsvendor model with backordering. In the two-period game, research by Yang and Xiao [8] study the pricing, ordering and quality of retailers' service decisions when loss-averse consumers are sensitive to service quality. They also construct a SC coordination mechanism. A study of risk aversion by Giri [9] simultaneously considers the case of a retailer ordering from a manufacturer with random output and a manufacturer with a fixed output and studies the inventory decision problem of a risk-averse retailer. Zhu et al. [10] consider the global optimization problem when the manufacturer has two sales channels. Wu et al. [11] study the influence of capacity uncertainty on the order decision of risk-averse managers. Additionally, Oh et al. [12] consider the impact of uncertainty in service cost on the product pricing decision when the service provider of the downstream SC is risk-averse. Egging et al. [13] discuss the influence of risk aversion on gas investment decisions. Xiao et al. [14] consider the retailer's service and pricing decisions when both system members and consumers are risk-averse. Zhou et al. [15] study a coordination problem when chain members and service provider are risk-averse. Zheng et al. [16] study the pricing decisions of two competitive shipping companies, one of which is risk-averse. Liu [17] analyzes the influence of risk aversion on the decision-making and expected profits of dual-channel members when information asymmetry exists in the risk aversion. For more research on risk aversion, see [18-20]. However, there are few studies available on the decision-making strategy problems that arise when decision makers have S^ or WA preferences. Schweitzer and Cachon [1] study the ordering decision when the decision-maker in the downstream SC has either a S^ or WA preferences. In contrast to the previous literature cited, in this paper, we introduce the S^ and WA coefficients under the manufacturer-led Stackelberg game. We also consider that the manufacturer may have RN, SA, WA, and SW preferences. This paper mainly studies the following three problems: • How will the four different behavioral preferences of the manufacturer affect the decision-making and expected profits of the chain members, as well as the expected profit of the system? • How will the change in the manufacturer's S^ and WA preferences affect the decisions of the chain members? • How will the change in the manufacturer's S^ and WA preferences affect the expected profits of the chain members and system? To solve these problems, we build a Stackelberg game model composed of a retailer and a manufacturer, who is the leader and assume the demand is random. The retailer pursues its own expected profit maximization by making order decisions, while the manufacturer pursues its own expected utility maximization by setting wholesale pricing. We compare the decision-making behavior and expected profits of the chain members, and the expected profit of the system, under the four different behavioral preferences. At the same time, we also analyze the influence of the change in the S^ and WA coefficients on the decision-making and expected profits of the chain members, as well as the expected profit of the system. The above conclusions can be used to some extent to guide decision-makers on how to adjust their decision-making strategies according to the manufacturer's S^ and WA preferences. This paper is organized as follows. Section 2 provides the assumptions and notations of the model. Section 3 presents the model analysis under the four kinds of behavioral preferences. Section 4 offers the comparison and analysis of the four cases. Section 5 provides the numerical study, and Section 6 summarizes the full text. 2. Assumptions and notations 2.1 Model assumptions The mathematical model presented in this paper is based on the following assumptions: Assumption 1: The entire SC consists of a manufacturer and a retailer, in which the manufacturer is the leader and the retailer is the follower. Assumption 2: The cumulative distribution function and probability density function of random demand D are F(x) and f(x), respectively, and the expected demand is E(x) =p. Assumption 3: Because market demand D is random, the retailer may overorder or underorder. To simplify the model, we ignore any salvage value or shortage loss. This assumption is common in the existing literature (e.g., [1, 4, 16, 20]) and does not have any effect on the conclusions. Assumption 4: The manufacturer may have four behavioral preferences: RN, SA, WA, and SW. Assumption 5: p>w >c. 2.2 Notations The meanings of the other symbols used in this article are in Table 1. Table 1 The listing of notations (j e {RN, WA.SA.SW}) Symbol Meaning Symbol Meaning w Wholesale price c Production cost P Retail price R Order quantity ß WA coefficient X Stockout-aversion coefficient RN Risk-neutral WA Waste-averse SA Stockout-averse SW Stockout- and waste-averse SC Supply chain n1R The retailer's expected profit in case j nJs The supplier's expected profit in case j The expected profit of supply chain in case j ul The supplier's expected utility in case j * Optimal value 3. Model analysis In this section, we will study the optimal decisions of the participants under the four behavioral preferences models. In the manufacturer-led Stackelberg game, the manufacturer has priority decision-making power and first determines the wholesale price w to maximize its expected utility. Then, the retailer determines its order quantity q in response to the manufacturer's decision-making and the predicted market demand D. The expected profits of the manufacturer and the retailer are as follows: n]s = (w — c)q (1) fq nJR(q) = pEmin(q,D) — wq = (p — w)q — p I F(x)dx (2) Jo where j 6 {RN, WA, SA, SW}. Superscripts RN, WA, S^ and SW indicate that the manufacturer has risk-neutral, waste-aversion, stockout-aversion and stockout- and waste-aversion preferences, respectively. Case 1 With reference to Schweitzer and Cachon [1], when the manufacturer has SW preferences, the manufacturer's expected utility function is Uiw(w) = wq-cq- E{fi(q- D)+ -A(D- q)+] = (w-c + A)q-An + (J3+A) [ F(x)dx (3) Jo In the Eq. 3, wq is the product sales income of the manufacturer, cq is the production cost of the product, — D)+is the waste penalty loss and A(D — q)+is the shortage penalty loss. Among these variables, higher values of coefficients fi and A indicate higher degrees of WA and SA, respectively, by the manufacturer. d2nsw The inverse induction method is used to solve the model. From Eq. 2, we have „ 1 = aqi dnsw —pf(q) <0; then, n:|wis concave in q. Let J* = 0, and we can obtain the following: = f-i [4) By substituting q w*into Eq. 3, the manufacturer's expected utility optimization problem is Let = 0, and we can obtain the following: dq ! (P~w\ _ (w-c + A)(p-v) + (p- w)(A + fj) _o m P J- (p_v)2/[F-1(E^)] =0 Eq. 6 is too complex to obtain intuitionistic results, in order to gain more valuable conclusions, we follow prior literature [10, 20] and assume that the random demand variable obeys a uniform distribution, that is, x~U(0, B). Through Eqs. 4 and 6, we can now obtain the following: ww* = P(P + c + fO qSW* = B(p-c + A) f7) 2p + fi + A 2p+p+A LJ By substituting Eq. 7 into Eqs. 1 and 2, we can obtain the following: w* = B(p-c + X)\p(p + p)-c(p + p+X) w* = Bp(p-c + A)2 % (ip + p + xy ,Ur 2(2p + p + xy f8) sw* B(p — c + D[p(3p + 2P + A)-C(3p + 2/? + 21) nc -ns +nR - 2(2p + p+XY f9) Case 2 When the manufacturer has only WA preferences, that is, 2 = 0 , then the optimal decisions and expected profits of the participants are shown in the third row of Table 2. Case 3 When the manufacturer has only preferences, that is, fi = 0 , then the optimal decisions and expected profits of the participants are shown in the fourth row of Table 2. Case 4 When the manufacturer has RN preferences, that is, 2 = 0 and fi = 0, then the optimal decisions and expected profits of the participants are shown in the fifth row of Table 2. Table 2 Optimal decisions and expected profits of the manufacturer and retailer fj e {RN, WA,SA,SW}) Preferences RN WA SA SW sw qSW* ma£U;jw = (w-c + A)qsw*-An + (ß + A) í F(x)dx (5) Jo q¡* w¡* u1* ns nR B(P - c) p + c B(p- c)2 B(p- c)2 2p 2 4p 8p B(p- c) p(p + c+ß) B{p-c)2{p + ß) Bp(p — c)2 2p + ß 2p+ß (2p+ß)2 2{2p+ß)2 B(p-c + A) B{p + c) ß(p-c+l)[p2-c(p+l)] Bp(p — c+X) 2p + X 2p + Â (2p+A)2 2(2p+X)2 B(p — c + X) B(p + c + ß) B(p-c+Ä)[(p-c)(p+ß)-cÄ] Bp(p — c+X) 2p+ß + Ä 2p+ß+X (2p+ß+X)2 2{2p + ß+X) Table 3 Optimal expected profit of the SC nscA* 3B(p- c)2 8p B(p-c)2(3p + 2ß) 2(2p + ß)2 B(p-c + A)[p(3p + X)~ c(3p + 21)] 2(2p+A)2 B(p-c+X)[p(3p + 2ß+Ä)-c(3p + 2Ä + 2ß)] 2(2p + Ä + ß)2 4. Comparison and analysis of the four behavioral preferences mode This section mainly analyzes the conclusions obtained in the previous section. We analyze the influence of the change in S^ coefficient A and WA coefficient p on the decisions and expected profits of the participants, as well as the expected profit of the system. Additionally, the decisions and expected profits of the participants and the expected profit of the system are compared and analyzed under the four cases. The equilibrium results in the four cases are summarized in Tables 2 and 3. We first analyze the effect of A and p on the qsw* and wsw*. From Eq. 7, we can obtain dqsw*(A,P) _ B(p + p + c) dqsw*(A,P)_ B(p+A-c) dX _ (2p + p+A)2 >0, dp __ (2p + p + A)2 <0 (10) dwsw*gP) = p(p+ p + c) dwsw*(A,p) = p(p+A-c) , dA (2p+p+A)2 , dp (2p + p + A)2 ( ) Therefore, we can obtain Proposition 1. Proposition 1 (1) The retailer's order quantity increases with the A, while it decreases with the p. (2) The manufacturer's wholesale price decreases with the A, while it increases with the p. Proposition 1 shows that the S^ (WA) preferences motivates (weakens) the retailer's ordering enthusiasm. As A (P) increases, the order quantity will increase (decrease). On the other hand, the dominant manufacturer will adjust the price of the product according to the ordering decisions of the retailer in order to maximize its expected utility. So as A (P) increases, the wholesale price will decrease (increase). Next, we will analyze the influence of A and p on the expected profits of the participants. From Table 2, we can obtain dnsRw*^A,p) _ Bp(p-c + A)2 dnsRw\A,p) _ Bp(p-c + A)(p + c + p) dp ~~~2TpTW<0, dX _ (2p + p + A)3 >0 (12) dnjA'-(X) _ BAjp + c)2 dn^'XP) _ Bpjp- c)2 dA (2p + A)3 0 dp (2p + PY (13) dnssw\A,p) B(p — c + A)[p(A-p) + c(A + p)] dp " (2p + p + A)3 (14) dn|w*(À,p) _ B(p + c + p)[p(p~ À)-c(À + p)] dÀ (2p + p + X)3 (15) dnsw* dnsw* From Eq. 14, there are —-—>0when 2>21and —-—< 0when 0<2<21, where A1 = n dp 1 dp i' i Rirt-rdwSW* dwSW* From Eq. 15, there are >0 when P>PX and <0 when 0

A1 (P> Pt), the manufacturer's expected profit increases with the p (A); conversely, the manufacturer's expected profit decreases with the p (A). As A (P) increases, the retailer's expected profit increases (decreases), which is due to an increase (decrease) in order volume. When the manufacturer has only WA (S^) preferences, it is understandable that the manufacturer's expected profit decreases as the order quantities de- crease. However, the manufacturer's expected profit decreases as order volume increases, which seems counterintuitive; this is because the gains from the increase in order quantity are not sufficient to compensate for the loss caused by the reduction in product price. In addition, the relationship between manufacturer's expected profit and the A (0) is related to the 0 (2) in the SW model. From Table 3, we can obtain dnscw*(.A,f3) = B(p + c + 0)[p(0 + p)-c(p+1 + f3)} dA (2p + 0 + A)3 (16) dn^w*(A,p) _ B(p — c + A)[p(/3 + p)-c(0when /3 > 02and—£— < Owhen 0 <0 <02, where02 = OA OA r(n+2\-r,2 fiuSW* fin™* c(p ) v . From Eq. 17, there are2^ >0 when A>A2 and^- <0 when 0 0 ( p2-c(p + A)>0 (18) dnscA* _ B(p + c)[p2 -c(p + A)] , dn™A* _ B(p — c)2(p + 0) dA (2p + A)3 0 dp (2p+A)3 (19) From the above analysis, we can obtain Proposition 3. Proposition 3 (1) When the manufacturer has only S^ (WA) preferences, the expected profit of SC increases (decreases) with the A (ft). (2) When the manufacturer has SW preferences, the relationship between the expected profit of SC and the 0 (A) is determined by A( ft). There is a threshold A2 (fi> fi2), and when the value of A (ft) is large, that is, A> A2 (fi> fi2), the expected profit of SC increases with the 0 (A); in contrast, the expected profitof SC decreases with the 0 (A). Proposition 3 shows that the system profit is increasing (decreasing) with the A (0) in the S^ (WA) model. Because for the whole system, the order quantity increases (decreases) with the increase (decrease) in A (0), which can (cannot) satisfy the potential market demand as much as possible. In addition, if A (¡3) is larger, there will be a risk of overordering or underordering in the SW model. Therefore, reducing (increasing) the order amount is beneficial to increase the profit of the system. From Eqs. 10 and 11, we can obtain dwsw*{A,p) <0_( wsw*{A,p) 0_f wsw*(A,ß) wWA*(0) = w 0_( qSW*{A,ß) >qsw*{0,ß) = qWA*{ß) I dA 0) = qSA*{A) >qSA*{0) = q dqsw*{A,ß) <0_f qsw*{A,ß) qK'v*when m>mi and qbW* qWA* mx, the order quantity is larger in the SW model; otherwise, the order quantity is larger in the RN model. B = 100, p = 10, c = 1,0 = 5. ■o o « 45 J3 Fig. 1 The effect of ratio m on the retailer's order quantity From the Eq. 12, we can obtain the following: dnsRw* <0 dß dnsRw\X,ß) dX Let m = 2/ß, then „SW* „RN* _ nR nR — nsRw*{X,ß)0 n sw*. inR nö"" nsRw*(.0,ß) = n™A*(ß) .Rt R SW*(A, 0) = nsRA*(X) >^*(0) = nRN* or n%w* njjwhen m>m1 and ^^ when 0m1 rsw* SA* n WA* 0 m1, the expected profit is larger in the SW model; conversely, the expected profit is larger in the RN model. 70 65 C 60 re 40 35 0 0.5 1.5 m 1.5 m (a) (b) Fig. 2 The effect of ratio m on the expected profit of manufacturer (a) and retailer (b) n, WA* Proposition 6 < ■s -SA* nS m4 where m2 = -2p2-c(2p+ß)+^4p3(p+ß)+4cp2(2p+ß)+c2(4p2+ß2) 2 P{p+C) The proof for Preposition 6 is in Appendix A. m3 = pjp-c) p2+c(p+ßY m4 = (2p+ß)(p-c) p2+c(p+ß) As shown in Fig. 2(a), Proposition 6 indicates that the manufacturer's expected profit is highest in the RN model. However, the expected profit difference between the WA, SW and SA models depends on the ratio m. There are three thresholds, m2 , m3 and m4, where m2 m4, the expected profit is the highest and lowest in the WA and SA models, respectively. nWA* mt: where mr = and mfi = 3 p+c 0 2p2+pß-2c(3p+ß)+^p2(2p+ß)2+4cp(2p2+pß-ß2)+4c2(p2+ß2) 4cß . The proof for Preposition 7 is in Appendix A. 200 190 140 2 2.5 3 0 0.5 2 2.5 3 a a. f/ / / // m 5 C SA* nC -SW* 6 8 10 12 14 16 18 Fig. 3 The effect of ratio m on the expected profit of SC 80 RN > 75 n C 70 65 C 60 m 6 55 0 20 As shown in Fig. 3, Proposition 7 shows that the expected profit of SC is lowest in the WA model. However, the expected profit difference between the RN, SW and SA models depends on the ratio m. There are two thresholds, m5 and m6, where m5 < m6. When 0 m6, the expected profit is the highest and lowest in the SW and RN models, respectively. 5. Results and discussion A numerical simulation was used to intuitively show the impact of the A and the ft on the decisions and expected profits of the chain members and the expected profit of the system. At the same time, we also want to verify the correctness of the above conclusions. Therefore, this section uses a numerical simulation method for further analysis. Without the loss of generality, the parameters are set as follows: B = 100,p = 10, c = 1. From Fig. 4(a), we can see that the manufacturer's wholesale price is decreasing with the A and increasing with the ft It can be seen from Fig. 4(b) that the retailer's order quantity is increasing with the A and decreasing with the ft a o IB CO ® o ■C is a a. 3 CO c ro 3 U" V 0 0.5 1 1.5 2 2.5 3 □ □ (a) (b) Fig. 4 The effect of X and p on wholesale price (a), and order quantity (b) From Fig. 5(a), we can see that the retailer's expected profit is increasing with the A and decreasing with the ft Fig. 5(b) shows that when the manufacturer has only SA (WA) preferences, the manufacturer's expected profit is decreasing with the A (ft. From Figs. 6(a) and 7(a), we can see that when the manufacturer has SW preferences, the relationship between the expected profit of the manufacturer (SC) and the A depends on the p. If ft is large, that is, /3 = 4.0 , then the expected profit is increasing with A; in contrast, if ft is small, that is, p = 0.5 , then the expected profit is decreasing with A. Figs. 6(b) and 7(b) show that the relationship between the expected profit of the manufacturer (SC) and the /3 depends on the A. If A is large, that is, A = 3.0 , then the expected profit is increasing with ft in contrast, if A is small, that is, A = 0.4 , then the expected profit is decreasing with p. Fig. 8 shows that when the manufacturer has only the SA (WA) preferences, the expected profit of the SC is increasing (decreasing) with the A (ft. 4.4 0 2 3 4 2 120 a. ■a S o 110 0.5 1 2.5 3 . ■o . . 3 CO 0.5 1 1.5 2 2.5 3 □ □ or □ (a) (b) Fig. 5 The effect of X and p on the expected profits of retailer (a) and manufacturer (b) o o CD a X CD a a 3 to «= o o CD a X CD CO a a □ □ (a) (b) Fig. 6 The effect of X and p on the expected profit of the manufacturer c CO ■C o > a a 3 a T3 o -»-» o o a 100 105 110 115 120 405.5 405 404.5 404 403.5 403 402.5 402 401.5 401 400.5 0 □ □ (a) (b) Fig. 7 The effect of X and p on the expected profit of the SC 203 202 201 200 90 80 70 0 1.5 2 0 202.5 202.5 202 202 201.5 201.5 201 201 200.5 200.5 200 200 199.5 3 199.5 199 199 1.5 2 2.5 3 1.5 2 2.5 3 405.05 404.95 404.9 404.85 404.8 404.75 404.7 404.65 95 2 3 Fig. 8 The effect of X or /3 on the expected profit of the SC 6. Conclusion This paper studies the influence of a dominant manufacturer's behavioral preferences on the decision-making and expected profits of the chain members, as well as the expected profit of the system. The behavioral preferences investigated include RN, S^, WA and SW preferences. Through the analysis of this paper, we can draw the following conclusions: • The manufacturer's wholesale price is increasing (decreasing) with the S^ (WA) coefficient, while the retailer's order quantity (expected profit) is the opposite; • When the manufacturer has only the S^ (WA) preferences, the manufacturer's expected profit decreases with the S^ (WA) coefficient; • When the manufacturer has the SW preferences, the relationship between the expected profit of the manufacturer (SC) and the S^ coefficient depends on the WA coefficient. The relationship between the expected profit of the manufacturer (SC) and the WA coefficient depends on the S^ coefficient; • The manufacturer's wholesale price is the highest in the WA model, followed by the RN, S^ and SW models, in that order; • The retailer's order quantity (expected profit) is the largest and smallest in the S^ and WA models, respectively, while the size of the order quantity (expected profit) between the RN and SW models depends on the ratio m, and there is a ratio threshold of m^, • The manufacturer's expected profit is the largest in the RN model. However, the size of the expected profit among the WA, SW and S^ models depends on the ratio m, and there are three ratio thresholds of m2, m3 and m4, where m2 < m3 < m4; • The expected profit of the SC is the lowest in the WA model. However, the size of the expected profit among the RN, SW and S^ models depends on the ratio m, and there are two ratio thresholds of m5 and m6, where m5 < m6. Acknowledgement This work was supported by the Natural Social Science Foundation of China (18BGL104). References [1] Schweitzer, M.E., Cachon, G.P. (2000). Decision bias in the newsvendor problem with a known demand dis-tribtion: Experimental evidence, Management Science, Vol. 46, No. 3, 404-420, doi: 10.1287/mnsc.46.3.404. 12070. [2] Jian, M., Fang, X., Jin, L.-Q., Rajapov, A. (2015). The impact of lead time compression on demand forecasting risk and production cost: A newsvendor model, Transportation Research Part E: Logistics and Transportation Review, Vol. 84, 61-72, doi: 10.1016/j.tre.2015.10.006. [3] Zhu, X.D., Li, B.Y., Wang, Z. (2017). A study on the manufacturing decision-making and optimization of hybridchannel supply chain for original equipment manufacturer, Advances in Production Engineering & Management, Vol. 12, No. 2, 185-195, doi: 10.14743/apem2017.2.250. [4] Niu, B., Chen, L., Zhang, J. (2017). Sustainability analysis of supply chains with fashion products under alterntive power structures and loss-averse supplier, Sustainability, Vol. 9, No. 6, 1-19, doi: 10.3390/su9060995. [5] Choi, S. (2018). A loss-averse newsvendor with cap-and-trade carbon emissions regulation, Sustainability, Vol. 10, No. 7, 1-12, doi: 10.3390/su10072126. [6] Li, X., Li, Y. (2016). On the loss-averse dual-sourcing problem under supply disruption, Computers & Operations Research, Vol. 100, 301-313, doi: 10.1016/i.cor.2016.12.011. [7] Xu, X., Wang, H., Dang, C., Ji, P. (2017). The loss-averse newsvendor model with backordering, International Journal of Production Economics, Vol. 188, 1-10, doi: 10.1016/i.iipe.2017.03.005. [8] Yang, D., Xiao, T. (2017). Coordination of a supply chain with loss-averse consumers in service quality, International Journal of Production Research, Vol. 55, No. 12, 3411-3430, doi: 10.1080/00207543.2016.1241444. [9] Giri, B.C. (2011). Managing inventory with two suppliers under yield uncertainty and risk aversion, International Journal of Production Economics, Vol. 133, No. 1, 80-85, doi: 10.1016/i.iipe.2010.09.015. [10] Zhu, L., Ren, X., Lee, C., Zhang, Y. (2017). Coordination contracts in a dual-channel supply chain with a risk-averse retailer, Sustainability, Vol. 9, No. 11, 1-21, doi: 10.3390/su9112148. [11] Wu, M., Zhu, S.X., Teunter, R.H. (2013). The risk-averse newsvendor problem with random capacity, European Journal of Operational Research, Vol. 231, No. 2, 328-336, doi: 10.1016/i.eior.2013.05.044. [12] Oh, S., Rhodes, J., Strong, R. (2016). Impact of cost uncertainty on pricing decisions under risk aversion, European Journal of Operational Research, Vol. 253, No. 1, 144-153, doi: 10.1016/i.eior.2016.02.034. [13] Egging, R., Pichler, A., Kalv0, 0.I., Walle-Hansen, T.M. (2017). Risk aversion in imperfect natural gas markets, European Journal of Operational Research, Vol. 259, No. 1, 367-383, doi: 10.1016/i.eior.2016.10.020. [14] Xiao, T., Choi, T.M., Yang, D., Cheng, T.C.E. (2012). Service commitment strategy and pricing decisions in retail supply chains with risk-averse players, Service Science, Vol. 4, No. 3, 236-252, doi: 10.1287/serv.1120.0021. [15] Zhou, Y.-W., Li, J., Zhong, Y. (2018). Cooperative advertising and ordering policies in a two-echelon supply chain with risk-averse agents, Omega, Vol. 75, 97-117, doi: 10.1016/i.omega.2017.02.005. [16] Zheng, W., Li, B., Song, D.-P. (2017). Effects of risk-aversion on competing shipping lines' pricing strategies with uncertain demands, Transportation Research Part B: Methodological, Vol. 104, 337-356, doi: 10.1016/i.trb.2017. 08.004. [17] Liu, M., Cao, E., Salifou, C.K. (2016). Pricing strategies of a dual-channel supply chain with risk aversion, Transportation Research Part E: Logistics and Transportation Review, Vol. 90, 108-120, doi: 10.1016/i.tre.2015.11.007. [18] Downward, A., Young, D., Zakeri, G. (2016). Electricity retail contracting under risk-aversion, European Journal of Operational Research, Vol. 251, No. 3, 846-859, doi: 10.1016/i.eior.2015.11.040. [19] Ohmura, S., Matsuo, H. (2016). The effect of risk aversion on distribution channel contracts: Implications for return policies, International Journal of Production Economics, Vol. 176, 29-40, doi: 10.1016/i.iipe.2016.02.019. [20] Li, B., Chen, P., Li, Q., Wang, W. (2014). Dual-channel supply chain pricing decisions with a risk-averse retailer, International Journal of Production Research, Vol. 52, No. 23, 7132-7147, doi: 10.1080/00207543.2014.939235. Appendix A Proof of Proposition 6 Let m = A/fi; from Eq. 13 and Table 2, we can obtain dn™* (!) dA dnfA* (ß) <0 ^ n%A* (A) < n%A* (0) = n trn* ls dß <0 ^ n™A* (ß) < n™A* (0) = nRN* n sw* RN. Bß2[c(1 + m) — p(1 -m)]2 Us =---—~—-—~--< 0 _sw* ns „SA* _ ns = 4p(2p + ß + mß)2 ?A* „SW* „RN*} _ ,iis , u s } — Bß2(p — c + mß)[p[p(2m — 1) + mß2] + c[p + 2mp + mß(1 + m)]} max[nssA*, n™A*, n™*, n™*} = n™* „WA* „54* _ ns ns — (2p + mß)2(2p + mß + ß)2 B(p — c)2(p + ß) B(p — c + mß)[p2 — c(p + mß)] „sw* ns „WA* _ ns — (2p + ß)2 (2p + mß)2 B(p — c + mß)[(p — c)(p + ß) — cmß] B(p — c)2(p + ß) (2p + ß + mß)2 (2p + ß)2 (A1) (A2) (A3) (A4) From Eqs. A1 to A4, we can obtain Proposition 6. Where m2, m3, and m4 satisfy Eqs. A2, A3, and A4, respectively. A,,... „ _ p(p-c) _ -2p2-c(2p+P)+^4p3(p+P)+4cp2(2p+P)+c2(4p2+P2) Additionally, m3 — p2+c(p+p), ™2 — 2^—) , and m4 = (2p+ß)(p-c) p2+c(p+ß) . Proof of Proposition 7 Let m — A/fi; from Eqs. 19 and Table 3, we can obtain B(p - c)2(3p + B(p-c + m^)[p(3p + 20 + mfi) - c(3p + 2m0 + 20)] nWA* - nsw* _ d(n £ 2(2p + ß)2 2(2p + mß +ß)2 WA* _ nSW*) Bß(p + c + ß)[p(p + ß) - c(p + mß + ß)] dm (2p + mß + ß)3 < 0 max{n™A* - nscw*} = (n™A* - nscw*)m=0 = 0 ^ n™A* < nscw* dnSA* w dX c (ß) gnWA* n RN* dß ~sw* _ llr — > 0 ^ nscA*(X) > nscA*(0) = n£N* <0 ^ n^A*(ß) < n^A*(0) = n%N* „WA* ^ „RN* ^ „SA* nC ^ nC ^ nC n. SA* Bß[c(l + m) - p(l - m)][c[4p + 3ß(l + m)] - p[4p + ß(3 + m)]} 8p(2p + mß + ß)2 Bß(p -c + mß)[p[4p2 + pß(3 + 2m) + mß2] -c[4p2 + 2mß2(1 + m) + 3p(ß + 2mß)]} „SW* _ llr — 2(2p + mß)2(2p + mß + ß)2 From Eqs. A5 to A8, we can obtain nWA* < nSW* < nRN* < nScA*> 0 m6 Where m5and m6satisfy Eqs. A7 and A8, respectively. Additionally, m5 = and m6 = 2p2+pP-2c(3p+P)+^p2(2p+P)2+4cp(2p2+pP~P2)+4c2(p2+P2) 4c/3 . (A5) (A6) (A7) (A8) APEM jowatal Advances in Production Engineering & Management Volume 13 | Number 3 | September 2018 | pp 358-365 https://doi.Org/10.14743/apem2018.3.296 ISSN 1854-6250 Journal home: apem-journal.org Original scientific paper Genetic programming method for modelling of cup height in deep drawing process Gusel, L.a*, Boskovic, V.b, Domitner, J.b, Ficko, M.a, Brezocnik, M.a aUniversity of Maribor, Faculty of Mechanical Engineering, Maribor, Slovenia bGraz University of Technology, Institute of Materials Science, Joining and Forming, Graz, Austria A B S T R A C T A R T I C L E I N F O Genetic programming method for modelling of maximum height of deep drawn high strength sheet materials is proposed in this paper. Genetic programming (GP) is an evolutionary computation approach which uses the principles of Darwin's natural selection to develop effective solutions for different problems. The aim of the research was the modelling of cylindrical cup height in deep drawing process and analysis of the impact of process parameters on material formability. High strength steel sheet materials (DP1180HD and DP780) were formed by deep drawing using different punch speeds and blank holder forces. The heights of specimens before cracks occur were measured. Therefore, four input parameters (yield stress, tensile strength, blank holder force, punch speed) and one output parameter (cup height) were used in the research. The experimental data were the basis for obtaining various accurate prediction models for the cup heights by the genetic programming method. Results showed that proposed genetic modelling method can successfully predict fracture problems in a process of deep drawing. © 2018 CPE, University of Maribor. All rights reserved. Keywords: Metal forming; Deep drawing; Modelling; Genetic programming *Corresponding author: leo.gusel@um.si (Gusel, L.) Article history: Received 19 December 2017 Revised 22 February 2018 Accepted 7 June 2018 1. Introduction Deep-drawing processes are frequently used in the sheet metal forming industry, because of the achievable formability. In deep drawing the sheet blank is deformed by tensile and compressive loads in different directions of action. The sheet thickness should remain as constant as possible. Many parameters influence the deep drawing process and should be carefully selected for effective and economical production Among them, the influence of punch speed and blank holding force are of great importance, since the shortened process duration leads to an increase in productivity. Therefore, gaining knowledge of the influence of these parameters is of great interest when achieving the greatest possible productivity in producing fault-free parts. Accurate models, which describe the influences of different parameters in deep drawing, can be obtained by many different modelling methods, and serve for optimization and prediction purposes in the process. In widely used deterministic methods like multi regression method, the form and size of a model is determined in advance. These modelling methods have one common feature: all of them optimize a given model of a problem. That is, each genotype determines a particular combination of values of free variables of the model, and the interpretation of those variables within a model is fixed. When applied to a different problem instances, they often turn out to be infea-sible or obtain much worse fitness. Much better results are often achieved by non-deterministic modelling methods, such as genetic programming method (GP). In GP the model is not known in advance, and only the constraints of the model are given, e.g., instruction set, the maximum size of organisms etc. [1]. Hence GP creates and optimizes entire model of the problem. In other words, GP evolves an executable code that inputs a given problem instance and outputs a solution to the instance. In this sense GP's solutions are active, i.e. adapt to the given problem instance [2]. Majority of papers dealing with modelling methods in metal forming processes, describing neural network and genetic algorithms (GA) methods while only a few are dealing with GP. In paper [3] authors describe prediction and optimization of sealing cover thinning in deep drawing process by using GA as very accurate and effective method. In [4] GA have been used to develop an optimization strategy for choosing blank holder force and punch force for enable fracture free and wrinkle-free production of deep drawn cups while in [5] optimum blank shape and process parameters were defined and calculated for deep drawing process using Taguchi optimization method and GA. Pareto optimal solution search techniques in GA to reduce excessive thinning and wrinkling occurrence was described in [6] and [7]. An evolutionary structural optimization method for the sheet blank shape optimization in deep drawing process was presented in [8]. Many authors have used combination of neural network and GA approach for modelling and predicting deep drawing parameters and some properties of drawn products [9, 10]. In [11], the effect of hydro-mechanical deep drawing process parameters was investigated by FE simulations and neuro-fuzzy modelling method to predict the maximum thinning of sheet material, while in [12] a fuzzy control algorithm for optimization of loading profiles and drawing ratio was proposed. An overview of recent applications of evolutionary computing in manufacturing industry was presented in [13]. Very few papers (and especially sheet forming processes) dealing with modelling by the GP modelling method. In papers [14, 15] authors compared genetic algorithm models and GP models for the distribution of effective strain and stress and also some mechanical properties in forward extruded alloy. GP method performed well in developing more accurate genetic models. In [16] authors described the application of bi-objective GP modelling for prediction of the size of austenite grain as a function of heating time. GP method was also used very successfully for modelling of bending capability of titan-zinc sheet in [17]. Different bending parameters and their impact on bending capability were analysed and optimized and a variety of accurate models were developed by GP. This paper proposed a GP modelling of maximum height of cup deep drawn high strength steel specimens. Maximum height achieved in deep drawing process depends on many parameters and is excellent indicator for formability level of chosen material. Experimental data measured during the deep drawing processes served as an environment for evolution process in genetic programming. Blank holder force, punch speed, tensile and yield stress were used as independent input variables, while maximum height of deep drawn specimens before first cracks occur was a dependent output variable. 2. Materials and methods 2.1 Genetic programming Genetic programming is a sub brunch of evolutionary computation (EC) emulating the natural evolution of species and is probably the most general approach among EC methods. To imitate the evolutionary process in GP, certain components must be defined, such as mathematical functions, problem variables and genetic operations like reproduction, mutation or crossover [1]. The GP process starts with random generation of initial population of organisms (computer programmes). Terminal genes, such as variables and constants, and function genes (arithmetic operations) compose each programme and must be carefully selected. After that the fitness function must be determined for evaluation of programme adaptation to the environment (e.g., to the experimental data). Adaptation is a main force in natural selection. The change and improvement of programmes fitness during GP process is enabled by genetic operations, such as reproduction, mutation and crossover. The right selection of genetic operations and their probability is of vital importance for successful GP process (and often vary regarding the problem to be solved), because genetic operations provide an increasing diversity and genetic exchange among computer programmes. The mutation also introduces new code fragments into population and is used as a common workaround for loose of diversity and stagnation, especially in small populations. The last step of the process is the definition of termination criterion which is usually prescribed number of generations. If the termination criterion is fulfilled the evolution is then terminated. In general, many independent GP runs are needed for successful and accurate problem solutions. 2.2 Experimental details The main goal of the experiments was determination of the impact of punch speed and blank holder force on the maximum height of deep drawn sheet metal cups before crack occurs. The two chosen parameters are very important for efficient and quality process of deep drawing. Deep drawn cups are cylinders that are closed on one end and open on the other, with or without a flange on the open end. Two different sheet materials were used (DP1180HD and DP780) which are new advanced high strength steels with good cold formability developed for automotive industry. These two materials are also suitable for welding and show excellent characteristics by crash tests. With the help of tensile tests two important mechanical characteristics were determined: Rm (tensile strength) and yield stress Rpo.2 for DP1180HD (Rpo.2 = 1077 N/mm2, Rm = 1269 N/mm2) and for DP780 (Rpo.2 = 490 N/mm2, Rm = 840 N/mm2). Material thickness was s = 1.5 mm, diameter of round sheet specimens was D = 215 mm and punch diameter was d = 100 mm. For deep drawing process special experimental tool with a cylinder shaped punch and hydraulic press SHC-400 were used [18], Fig. 1(a). The process starts with the application of the same amount of Wisura FMO 5020 lubricating oil on both sides of the specimens. Because of cup geometry, anisotropy doesn't have any influence so position of the blank doesn't have to be considered. The measurement of the deep drawn cup's height was executed by laser and acoustic measurement devices. Laser measures drawing height and acoustic sensor device, which was attached in the inner side of the die, detects the moment when crack occurs. With the interaction of these two devices and special software it is possible to detect exact drawing height of the cup just before first crack occurs. For punch speed the experimental range was set from 50 mm/s to maximum speed of 150 mm/s. The latter is usually the highest punch speed tolerated in praxis for successful and economical deep drawing process of high strength steel. Blank holder force was set to 200 kN and 600 kN respectively. In order to provide reliable results, three experiments for the same parameters were performed (60 experiments for both materials) and then average value of the height was calculated. The shape of deep drawn cups for both materials after the experiment is shown in Fig. 1(b). Therefore, we obtained a total of 20 combinations of experimental data. The results for measured cup's height are listed in Table 1. Table 1 Experimental results of deep drawing process Experiment Yield stress Tensile strength Blank holder Punch speed Cup height No. Rp0.2 (X1) Rm (X2) force Fb (X3) v (X4) H (N/mm2) (N/mm2) (kN) (mm/s) (mm) 1 490 840 200 50 20.086 2 1077 1269 200 50 15.596 3 490 840 200 75 20.043 4 1077 1269 200 75 15.360 5 490 840 200 100 19.853 6 1077 1269 200 100 15.220 7 490 840 200 125 19.423 8 1077 1269 200 125 14.626 9 490 840 200 150 18.623 10 1077 1269 200 150 13.073 11 490 840 600 50 19.863 12 1077 1269 600 50 15.580 13 490 840 600 75 19.667 14 1077 1269 600 75 15.140 15 490 840 600 100 19.370 16 1077 1269 600 100 14.500 17 490 840 600 125 18.936 18 1077 1269 600 125 14.293 19 490 840 600 150 18.360 20 1077 1269 600 150 13.256 Fig. 2 shows the influence of punch speed v and blank holder force Fb on maximum cup height It is obvious that an increase of the punch speed leads to a decrease of the maximum drawing height, i.e. to decrease of the formability. This is valid for both examined high strength steels. The highest difference between measured height when blank holder force Fb = 200 kN was used was 7.2 % for DP780. The difference was much higher for DP1180HD, i.e. 16 % decrease of height when punch speed v = 150 mm/s was used compared to the height achieved at v = 50 mm/s. The increase of blank holder force also leads to a decreasing drawing height but to a smaller extent. A change of the blank holder force does not affect the shape or tendency of the curve. Since increasing punch speed has confirmed the tendency of the formability to decrease, but high speeds are needed to achieve a certain level of productivity, it is important to optimize the process parameters with the goal to improve the efficiency and quality of the deep drawing process. Fig. 2 The influence of punch speed v and blank holder force Fb on maximum cup height 3. Results and discussion Evolutionary parameters that were used for GP modelling processes were: population size was 400 for all GP runs, maximum number of generations varied from 500 to 1000, reproduction probability was set to 0.1, probability of crossover varied from 0.1 to 0.3, while probability of mutation varied from 0.6 to 0.8. Such a high probability of mutation has been chosen because of the relatively small allowable maximal depth of organisms. In this way, sufficient diversity of organisms in the population is preserved. Maximum allowed depth for organisms created in the initial generation was 6 and maximum depth after crossover was 10 (in some cases only 8). Experimental results from Table 1 were used as a training data set for GP process. Additional experiments were also performed for testing data purpose and were not used in training data set. Three different function sets were applied for GP modelling. Each function set contained function genes. First set contained three operations (function genes): addition, subtraction and multiplication, i.e. F = {+, -, *}, while division operation was added to the second function set, i.e. F = {+, -, *, /}. In third function set, natural exponential function was added to the second function set {+, -, *, /, ZEXP}. All function genes were protected against extreme values that can be occurred during simulated evolution process [1]. Terminal set comprised of terminal genes. In our case, the terminal consisted of four input variables, and random real numbers, i.e. T = {X1, X2, X3, X4, M}. X1 is yield strength (#p0.2), X2 is tensile strength (#m), X3 is blank holder force (Fb), X4 is punch speed (v), and M are random generated real constants from the interval -1 to 1. By applying the first function set and terminal set, the most accurate genetically developed model for a cup height H developed by genetic programming was: (+ (* -8.71014 -2.13249) (* (+ (* (+ (- (+ (- (* -9.70734 X1) (* X1 X4)) (* (- X2 X1) (+ X2 9.80336))) (+ (* (+ 9.80336 X4) (+ X1 X1)) (* X1 (+ X4X1)))) (+ (- (- (* -9.70734 X1) (* X4X3)) (+ (* X3 X4) (* X4 X4))) (- (- (* -9.70734 X1) (* X4 X4)) (* (+ 9.80336 X4) (+ X1 X3))))) (- -2.07369 -2.07544)) (- (+ (+ (- (- X2 1.62782) (* -0.0546093 X4)) (+ 9.76417 1.167)) (- (+ (- (-X2 1.62782) (* -0.0546093 X4)) (+ 9.76417 1.167)) X1)) (- (- (- (+ X1 (* -0.0546093 X4)) 9.80336) 9.80336) (- (- (- X2 1.62782) (* -0.0546093 X4)) (+ X1 (+ (* 3.08974 -6.87342) (- -7.08458 4.51225))))))) (- -2.07369 -2.07544)))) The above model is presented in prefix notation of programming language LISP. After simplifying, the above model for a cup height H can be written as the mathematical expression (please note that in Eq. 1 instead of variable names X1, X2, X3, X4 original symbols are used): 18.6958 - 3.0625 * 10~6fln2 + 0.00528002 Rm + 3.0625 * 10~6Rr 0.0000300228 Fb +ßp(—0.00545928 - 3.0625 * 10"6 R, -9.1875 * 10~6Fhv - 6.125 * 10"6 v2 0.0000153125 v) + 0.000382265 v (1) Generation Fig. 3 Deviation of best GP model using first function set {+, -, *} and training data 0 100 200 300 400 500 600 700 800 900 1000 Generation Fig. 4 Deviation of best GP model using second function set {+, -, *, /} and training data The model in Eq. 1 was generated in generation 500. It has a total of 127 genes, 63 functional genes and the depth of this model is 10. Probability of reproduction was 0.1, probability of mutation 0.6, and probability of crossover 0.3. Maximal allowed depth after crossover was 10. The average difference between experimental results and the results predicted by the genetic developed model in Eq. 1 was 5 = 1.14 % and 5 = 1.18 % for training data set and testing data set, respectively. Fig. 3 presents the average deviation (error) of best GP models when applying first function set. In first few generations there is a large deviation between the prediction results of the GP developed model and experimental data but after several generations natural selection randomly added new genes. Thus, diversity get bigger and deviation becomes better, i.e. lower. After generation 250, the deviation of best models improves very slowly and does not change much until final generation. The most accurate GP model of all experiments for modelling of cup height was developed by applying the second function set (+, -, *, /). In the prefix LISP form it can be expressed as: (+ (/ (+ (- (/ (- (* 9.40782 X2) (* X2 -5.06327)) (+ (- X3 9.69337) (* -0.346372 X3))) (+ (+ (3.49304 5.8424) -8.94781) (+ (- X1 0.583014) -8.27698))) (- (/ (- (- X2 -5.60735) (- X3 9.69337)) 4.48845) (+ (/ (* X3 3.73772) (- X2 X1)) (+ (- -9.68001 3.38485) (- X4 -0.0676581))))) (+ (+ (/ (+ (* X1 3.13824) (* X1 3.13824)) (- (+ 7.6819 5.58638) (- X3 X4))) (+ (/ (+ X2 X3) (- X1 X2)) (+ -1.18787 X2))) (/ (+ (- (* X1 -8.97815) (* X1 X4)) (* (- X4 -8.75927) (* X1 -2.32527))) (- (* (-X2 8.5169) 0.560745) (+ (/X1 X4) (/X1 X4)))))) (+ (/ 5.35335 (/ (+ (- (+ X1 -8.27698) (/ X2 X4)) (/ (+ X2 5.21486) -9.66085)) (+ X2 (/ (* X1 8.08171) (+ -6.90552 X3))))) 8.66904))) The above model can be written in the infix mathematical expression as: 8.08171 R„ \ 5.35335 IRm + - x -6.90553 + Fi 8.66904+ ---l- -8.81677 + R„- 0.103511 Rm--^ P m y í Í 14 4711 \ 3 73772 Fu \ (2) (36.5633 -Rp+Rp (0.222794 - 9.69337 J 0.653628 FJ ~ a222794 + Rp+Fb , 6.27648 fip fip(-29.3458 - 3.32527 v) v -1.18787 + Rm + tt—ît1 -F-n— + ■ RP ~Rm 13.2683 -Fb + v^ -2Rp - 4.77581 v + 0.560745 Rm v In this process we increased the highest allowed number of generations from 500 to 1000. Also the values of probability for some genetic operations were changed. The obtained model (Eq. 2) has an average error of 5 = 0.65 % (5 = 0.68 % for testing set) and was generated in generation 1000, has a total of 139 genes, 69 functional genes and the depth of this model is 8 which was also the highest allowed depth after crossover. In this case the values of main evolutionary parameters were: probability of reproduction was 0.1, probability of crossover 0.1 while mutation probability was set to 0.8. This model is more complicated compared to GP model in Eq. 1 but it is also much more accurate. Fig. 4 presents the average deviation (error) of the best GP model when second function data set was applied for GP process. It can be seen that after generation 50 all calculated average errors of best GP models are smaller than 3 % and after generation 360 they are all under 1 %. From that point and up to final generation the percentage deviation decreases very slowly. If more generations would be used in genetic programming the accuracy of the models would not increase significantly, but processing time would be much longer. Some very simple genetic models were also obtained by GP process. One of the best simple model with average error of 1.31 % was developed with first function set (+, -, *) and is presented in Eq. 3. -0.30418 - 0.0469128 Rp + 0.0532937 Rm - 0.000797616 Fb - 0.0170231 v (3) The above model is the evidence that GP process is capable of developing not only complex genetic models, but also very simple models with satisfactory accuracy. Simple models are easy to use, but when accuracy of the model is priority, such as in our research, the more complex and more accurate models should be used for prediction and simulation purposes. We have also performed a few runs of GP modelling with third function set by adding natural exponential function (ZEXP) to four basic math functions. The best obtained genetic models were complex but not so accurate. It was obvious that adding ZEXP function does not results in getting better genetic models, in some cases natural selection in GP even eliminates exponential function and so the best developed models were without ZEXP function. The reason for this happening could be in relatively simplicity of the studied problem. 4. Conclusion In the paper a GP modelling method for maximum height of cup deep drawn high strength steel specimens, which is an indicator for cold formability level of the material, was presented. The research showed that an increase of the punch speed and blank holder force leads to a decrease of the maximum drawn height. By GP modelling it was possible to develop many different and very accurate genetic models. By using and combining different values for evolution parameters the optimal and most suitable GP models were developed. In our case, a high probability of mutation was vital for good convergence of the algorithms, because it preserves high diversity of a population. With lower value of probability of mutation, the convergence becomes either very slow due to low diversity of organisms in a population or organisms do not even converge to a sufficiently good solution. In the paper only three of the best developed genetic models were presented. The most accurate GP model was also very complex, but this is not an obstacle because modern production processes are all supported by high performance computers, so it is easy to use even the most complex models for accurate prediction of chosen parameters. In mass production processes, such as deep drawing, sometimes even the tiny improvement in optimization of process parameters can lead to massive reduction of production costs and consequently highly accurate models are desired. In our future work we intend to perform experiments with much more different materials and also additional parameters with the goal to optimize the input parameters for achieving better formability in deep drawing process. Acknowledgement The authors thank the Institute of Materials Science, Joining and Forming at Graz University of Technology, Austria for possibility of using their equipment for the experimental work, and the Slovenian Research Agency ARRS for partial financing of Slovenian research team. References [1] Koza, J. (1992). Genetic programming, The MIT Press, Massachusetts, USA. [2] Pawlak, T.P. (2015). Competent algorithms for geometric semantic genetic programming, PhD thesis, University of Technology, Poznan, Poland, doi: 10.13140/RG.2.1.1240.7760. [3] Kakandikar, G.M., Nandedkar, V.M. (2016). Prediction and optimization of thinning in automotive sealing cover using genetic algorithm, Journal of Computational Design and Engineering, Vol. 3, No. 1, 63-70, doi: 10.1016/j. jcde.2015.08.001. [4] Gharib, H., Wifi, A.S., Younan, M., Nassef, A. (2006). Optimization of the blank holder force in cup drawing, Journal of Achievements in Materials Manufacturing Engineering, Vol. 18, No. 1-2, 291-294. [5] Sener, B., Kurtaran, H. (2016). Optimization of process parameters for rectangular cup deep drawing by the Taguchi method and genetic algorithm, Materials Testing, Vol. 58, No. 3, 238-245, doi: 10.3139/120.110840. [6] Wei, L., Yuying, Y. (2008). Multi-objective optimization of sheet metal forming process using Pareto-based genetic algorithm, Journal of Materials Processing Technology, Vol. 208, No. 1-3, 499-506, doi: 10.1016/j.jmatprotec. 2008.01.014. [7] Di Lorenzo, R., Ingarao, G., Marretta, L., Micari, F. (2009). Deep drawing process design: A multi objective optimization approach, Key Engineering Materials, Vol. 410-411, 601-608, doi: 10.4028/www.scientific.net/KEM.410-411.601. [8] Naceur, H., Guo, Y.Q., Batoz, J.L. (2004). Blank optimization in sheet metal forming using an evolutionary algorithm, Journal of Materials Processing Technology, Vol. 151, No. 1-3, 183-191, doi: 10.1016/i.imatprotec.2004. 04.036. [9] Singh, D., Yousefi, R., Boroushaki, M. (2011). Identification of optimum parameters of deep drawing of a cylindrical workpiece using neural network and genetic algorithm, International Journal of Mechanical, Aerospace, Industrial, Mechatronic and Manufacturing Engineering, Vol. 5, No. 6, 987-993, doi: 10.1999/1307-6892/2043. [10] Zhao, J., Wang, F. (2005). Parameter identification by neural network for intelligent deep drawing of axisymmet-ric workpieces, Journal of Materials Processing Technology, Vol. 166, No. 3, 387-391, doi: 10.1016/i.imatprotec. 2004.08.020. [11] Tinkir, M., Dilmef, M., Türköz, M., Halkaci, H.S. (2015). Investigation of the effect of hydromechanical deep drawing process parameters on formability of AA5754 sheets metals by using neuro-fuzzy forecasting approach, Journal of Intelligent & Fuzzy Systems, Vol. 28, No. 2, 647-659, doi: 10.3233/IFS-141346. [12] Öztürk, E., Türköz, M., Halkaci, H.S., K09, M. (2017). Determination of optimal loading profiles in hydro-mechanical deep drawing process using integrated adaptive finite element analysis and fuzzy control approach, The International Journal of Advanced Manufacturing Technology, Vol. 88, No. 9-12, 2443-2459, doi: 10.1007/ s00170-016-8912-x. [13] Oduguwa, V., Tiwari, A., Roy, R. (2005). Evolutionary computing in manufacturing industry: An overview of recent applications, Applied Soft Computing, Vol. 5, No. 3, 281-299, doi:10.1016/i.asoc.2004.08.003. [14] Brezocnik, M., Kovacic, M., Gusel, L. (2005). Comparison between genetic algorithm and genetic programming approach for modeling the stress distribution, Materials and Manufacturing Processes, Vol. 20, No. 3, 497-508, doi:10.1081/AMP-200053541. [15] Gusel, L., Rudolf, R., Brezocnik, M. (2015). Genetic based approach to predicting the elongation of drawn alloy, International Journal of Simulation Modelling, Vol. 14, No. 1, 39-47, doi: 10.2507/IISIMM14(1)4.277. [16] Halder, C., Madei L., Pietrzyk, M., Chakraborti, N. (2015). Optimization of cellular automata model for the heating of dual-phase steel by genetic algorithm and genetic programming, Materials and Manufacturing Processes, Vol. 30, No. 4, 552-562, doi: 10.1080/10426914.2014.994765. [17] Kovacic, M., Uratnik, P., Brezocnik, M., Turk, R. (2007). Prediction of the bending capability of rolled metal sheet by genetic programming, Materials and Manufacturing Processes, Vol. 22, No. 5, 634-640, doi: 10.1080/10426910701323326. [18] Radakovics, S.M. (2017). Einfluss der Tiefziehgeschwindigkeit auf die Kaltumformbarkeit von HSS und AHSS Stählen, Bachelor's thesis, Fakultät für Maschinenbau und Wirtschaftswissenschaften, TU Graz, Austria (in German). Calendar of events • Industrial Engineering and Manufacturing Technologies, Phuket, Thailand, September 21-23, 2018. • 26th International Conference on Materials and Technology, Portorož, Slovenia, October 3-5, 2018. • International Conference on Additive Technologies, Maribor, Slovenia, October 10-11, 2018. • European Simulation and Modelling Conference (ESM 2018), Ghent, Belgium, October 24-26, 2018. • 5th International Conference on Materials, Mechatronics, Manufacturing and Mechanical, Kuching, Malaysia, November 2-4, 2018. • International Conference on Industrial Engineering and Engineering Management (IEEM), Bangkok, Thailand, December 16-19, 2018. • 3rd International Conference on 3D Printing Technology and Innovations, Rome, Italy, March 25-26, 2019. • 9th IFAC Conference on Manufacturing Modeling, Management, and Control, Berlin, Germany, August 28-30, 2019. Notes for contributors General Articles submitted to the APEM journal should be original and unpublished contributions and should not be under consideration for any other publication at the same time. Manuscript should be written in English. Responsibility for the contents of the paper rests upon the authors and not upon the editors or the publisher. Authors of submitted papers automatically accept a copyright transfer to Chair of Production Engineering, University of Maribor. For most up-to-date information on publishing procedure please see the APEM journal homepage apem-journal.org. Submission of papers A submission must include the corresponding author's complete name, affiliation, address, phone and fax numbers, and e-mail address. All papers for consideration by Advances in Production Engineering & Management should be submitted by e-mail to the journal Editor-in-Chief: Miran Brezocnik, Editor-in-Chief UNIVERSITY OF MARIBOR Faculty of Mechanical Engineering Chair of Production Engineering Smetanova ulica 17, SI - 2000 Maribor Slovenia, European Union E-mail: editor@apem-journal.org Manuscript preparation Manuscript should be prepared in Microsoft Word 2007 (or higher version) word processor. Word .docx format is required. Papers on A4 format, single-spaced, typed in one column, using body text font size of 11 pt, should not exceed 12 pages, including abstract, keywords, body text, figures, tables, acknowledgements (if any), references, and appendices (if any). The title of the paper, authors' names, affiliations and headings of the body text should be in Calibri font. Body text, figures and tables captions have to be written in Cambria font. Mathematical equations and expressions must be set in Microsoft Word Equation Editor and written in Cambria Math font. For detail instructions on manuscript preparation please see instruction for authors in the APEM journal homepage apem-journal.org. The review process Every manuscript submitted for possible publication in the APEM journal is first briefly reviewed by the editor for general suitability for the journal. Notification of successful submission is sent. After initial screening, and checking by a special plagiarism detection tool, the manuscript is passed on to at least two referees. A double-blind peer review process ensures the content's validity and relevance. Optionally, authors are invited to suggest up to three well-respected experts in the field discussed in the article who might act as reviewers. The review process can take up to eight weeks on average. Based on the comments of the referees, the editor will take a decision about the paper. The following decisions can be made: accepting the paper, reconsidering the paper after changes, or rejecting the paper. Accepted papers may not be offered elsewhere for publication. The editor may, in some circumstances, vary this process at his discretion. Proofs Proofs will be sent to the corresponding author and should be returned within 3 days of receipt. Corrections should be restricted to typesetting errors and minor changes. Offprints An e-offprint, i.e., a PDF version of the published article, will be sent by e-mail to the corresponding author. Additionally, one complete copy of the journal will be sent free of charge to the corresponding author of the published article. APEM journal Chair of Production Engineering (CPE) University of Maribor APEM homepage: apem-journal.org Advances in Production Engineering & Management Volume 13 | Number 3 | September 2018 | pp 233-368 Contents Scope and topics 236 Sustainable manufacturing - An overview and a conceptual framework for continuous 237 transformation and competitiveness Hussain, S.; Jahanzaib, M. A hybrid grey cuckoo search algorithm for job-shop scheduling problems under fuzzy conditions 254 Yang, F.; Ye, C.M.; Shi, M.H. Design, finite element analysis (FEA), and fabrication of custom titanium alloy cranial 267 implant using electron beam melting additive manufacturing Ameen, W.; Al-Ahmari, A.; Mohammed, M.K.; Abdulhameed, O.; Umer, U.; Moiduddin, K. Dynamic integration of process planning and scheduling using a discrete particle swarm 279 optimization algorithm Yu, M.R.; Yang, B.; Chen, Y. An integral algorithm for instantaneous uncut chip thickness measuring in the milling process 297 Li, Y.; Yang, Z.J.; Chen, C.; Song, Y.X.; Zhang, J.J.; Du, D.W. Change-point estimation for repairable systems combining bootstrap control charts and 307 clustering analysis: Performance analysis and a case study Yang, Z.J.; Du, X.J.; Chen, F.; Chen, C.H.; Tian, H.L.; He, J.L. Multi-objective optimization for delivering perishable products with mixed time windows 321 Wang, X.P.; Wang, M.; Ruan, J.H.; Li, Y. Game theoretic analysis of supply chain based on mean-variance approach under cap-and-trade policy 333 He, L.; Zhang, X.; Wang, Q.P.; Hu, C.L. Decision-making strategies in supply chain management with a waste-averse and 345 stockout-averse manufacturer Jian, M.; Wang, Y.L. Genetic programming method for modelling of cup height in deep drawing process 358 Gusel, L.; Boskovic, V.; Domitner, J.; Ficko, M.; Brezocnik, M. Calendar of events 366 Notes for contributors 367 Copyright © 2018 CPE. All rights reserved. apem-journal.org 9771854625008