Simulation optimisation of a sustainable copper mining closed-loop supply chain

Abstract

This paper addresses the problem of designing a supply chain (SC) according to sustainability aspects. It identifies a research gap where an optimisation model to address the location, inventory and transportation decision in a sustainable SC applied to a copper mining industry by complementarily using a simulation model to assess SC performance in different scenarios is needed. For this purpose, we propose a simulation model based on system dynamics (SD) to assess SC performance and to support decision making. The basis of the simulation model structure is a multi-objective mixed-integer linear programming model that aims to minimise total economic, emissions and social costs, and to also contemplate social impacts. We consider a real copper mining industry problem to illustrate this. We obtain a solution using a mathematical programming model and a simulation model. The optimisation results show improved SC performance in cost and emission reduction terms, and an improved social impact. The simulation model presents a near-optimal result and allows the possibility of anticipating adverse scenarios. Future research is oriented to other real-world applications, and to: consider alternative inventory policies; contemplate a stochastic approach; add new production and routing decisions; develop a hybrid multi-agent SD model.

Keywords:

• Closed-loop supply chain
• sustainability
• optimisation
• simulation
• copper mining industry

Introduction

Economic globalised systems have complex supply chains (SCs) with environmental and social impacts that need to be managed in line with different stakeholders’ expectations and for mitigating sustainability-related risks (Rebs et al., 2019). In a complex system such as the copper industry, with its multiple inputs and system changes behave in a non-linear way, the system contains extensive feedback. The analysis of this behaviour is time-consuming and procedurally sensitive (Mohammadi et al., 2022). Both inter- and intra-organisational SCs are part of economic and social systems, which are consumers of natural resources. Sustainable SC management has become a growing concern given limited available resources and an increasing population. Production and logistics activities consume these resources, and generate waste and greenhouse gas emissions. In addition, labour conditions and other social factors often place pressure on local communities, especially if benefits are inequitably distributed (Rebs et al., 2019). Therefore, a framework is needed to analyse multiple scenarios, support decision making, select appropriate solutions and discover how different variables relate to one another over time. By considering all this, as an experimental method SC simulation can explain how SC performance indicators react to changes in both controllable internal factors and uncontrollable external factors (Campuzano & Mula, 2011). Systems dynamics (SD) is a technique that can address these complex systems, and it allows the structural causes of the studied system's behaviour to be understood (Sterman, 2000). According to practitioners and researchers, simulation models can contribute to: (i) study system changes in the model; (ii) verify analytical solutions; (iii) provide a view about key variables and how they interact; (iv) experiment with new situations that involve risk or uncertainty; (v) test new policies and decision rules (Buschiazzo et al., 2020; Campuzano-Bolarín et al., 2013, 2020; Esteso et al., 2019; Freile et al., 2020).

In their literature review, Rebs et al. (2019) found that most sustainable SC models deal with macroscopic analysis levels, while models for inter-organisational SCs are less prominent. They also observed the need for hybrid models that integrate different simulation, optimisation or multicriteria decision-making models.

This paper presents a simulation optimisation model of a sustainable SC in the copper mining industry. The main objective is to validate the optimisation model to design a sustainable SC that considers the three sustainable (3S) aspects (economic, environmental and social) in location, inventory and transport (LIT) decisions in an integrated way. Here the optimisation model that was proposed originally by Becerra et al. (2023) is adopted as the basis of this research, and is dubbed as 3S-LIT. Besides, the simulation model aims to support decision making in different simulation scenarios in a real copper mining industry context. Here the main contribution is to create a simulation model based on SD to study and optimise the SCs behaviour in different scenarios and its resilience based on the aforementioned optimisation model. This model also shows major improvement in computing times. Thus SC simulation with SD allows the copper mining SC's behaviour to be evaluated in complex scenarios, such as supplier or market disruption, which are problems that industry constantly faces.

The remainder of this paper is organised as follows. The second section presents the literature review on simulation models for sustainable SCs. The third section indicates the problem description of designing an SC according to sustainability aspects. The fourth section describes the mathematical programming (MP) model formulation, namely 3S-LIT, and compares it to the current procedure in the real mining company under study. The fifth section details the SD model formulation and presents the simulation results. The sixth section provides the conclusions.

Literature review

The selected literature corresponds to optimisation-based simulations, which allow SC behaviour to be evaluated in different situations and scenarios based on an optimal result, which is the main proposal of this paper. The integration of a direct and reverse logistics system, including the use of traditional SC structures of the direct movement of goods to consumers and the application of specialised operations for reverse SC activities, builds a closed-loop SC (CLSC) structure (Kuo, 2011). Thus incorporating reverse logistics into SC modelling includes a number of activities, such as collection, cleaning, disassembly, testing and sorting, storage, transportation and recovery operations, among others (Bostel et al., 2005). For example, Calmon and Graves (2017) optimise inventory in a reverse logistics system that supports warranty returns and replacements for a consumer electronic device. Chen et al. (2015) develop a decision support tool for a product recovery strategy under uncertainties. Deng et al. (2014) build a model by considering reverse logistics cost and customer time satisfaction to solve the location-inventory-routing problem. Jindal et al. (2015) propose a model to design and optimise a CLSC in an uncertain environment. Rezapour et al. (2015) evaluate the CLSC's performance behaviour by focusing on the impacts of strategic facility location decisions and transport and inventory decisions. Tighazoui et al. (2019) determine the optimal capacities of manufacturing and remanufacturing stocks, purchasing warehouses and transport vehicles, and define the optimal percentage of end-of-life returned products. Specifically, the goal of incorporating reverse and closed-loop logistics into the SC design is to achieve a certain level of sustainability (Becerra et al., 2021).

Very few papers that appear in the recent literature have developed hybrid models based on optimisation and simulation for sustainable SCs, and we highlight the following: Azadeh and Vafa Arani (2016) develop an SD MP approach to design and plan a biodiesel SC from biomass fields to consumption markets. The system combines certain social and macro-economic causal relations and public policies to provide the main input parameters of the MP model. Sudarto et al. (2017) develop an efficient and flexible long-term capacity planning policy for reverse logistics (RL) social responsibility by using an SD approach and an MP model, which consist of an efficient flexible capacity planning policy delivered by solving a multi-objective mixed-integer non-linear programming model (MO-MINLP) in conjunction with the SD model and an optimum seeking grid procedure. Liu et al. (2018) build a hybrid of multi-objective optimisation and SD simulation to optimise the structure of straw-to-electricity SC and to design motivational mechanisms to enhance its sustainability. Specifically, after obtaining the optimal SC, the SD model introduces governance motivation measures and identifies key stakeholder relationships. Motevalli-Taher et al. (2020) propose a multi-objective mathematical model to minimise network costs and water use, and maximise job opportunities, where the simulation model is first applied to forecast each customer's product demand. Ekren et al. (2021) propose an inventory model to optimise the reorder and up-to (s; S) inventory levels of e-groceries for predefined sharing policies by using a simulation optimisation approach, which is based on the scatter search (SS) metaheuristic.

Furthermore, a two-stage stochastic mixed integer programming model is proposed by Aloui et al. (2021) to integrate key location-allocation, inventory and routing planning decisions, and Monte Carlo simulation to consider epidemic disruptions and demand uncertainty. All this is applied to construct a smaller sample average problem instead of solving the problem with all the possible scenarios. Table 1 summarises these hybrid simulation optimisation models for sustainable SCs. Thus we analyse (i) the problems faced, namely LIT decisions; (ii) the optimisation modelling approach, namely MP or heuristic algorithms; (iii) the simulation modelling approach, namely SD or discrete-event simulation; (iv) the sustainability approach, namely the economic, environmental or social aspect.

Table 1. Literature review of sustainable SCs’ hybrid simulation-optimisation models.

	Problems			Optimisation approach		Simulation approach		Sustainability approach
References	L	I	T	MP	H	SD	DES	EC	EN	S
Azadeh and Vafa Arani (2016)	✓	✓		S-MIP		✓		✓	✓
Sudarto et al. (2017)	✓			MO-MINLP		✓		✓	✓	✓
Liu et al. (2018)		✓	✓	MO-MINLP		✓		✓	✓
Motevalli-Taher et al. (2020)	✓	✓	✓	MO-MILP			✓	✓	✓	✓
Ekren et al. (2021)		✓			SS		✓	✓	✓
Aloui et al. (2021)	✓	✓	✓	S-MILP			✓	✓	✓	✓
Our model	✓	✓	✓	MO-MILP		✓		✓	✓	✓

L: location; I: inventory; T: transport; MP: mathematical programming; H: heuristic; SD: system dynamics; DES: discrete-event simulation; EC: economic; EN: environmental; S: social; S-MIP: stochastic mixed-integer programming; MO-MINLP: multi-objective mixed-integer non-linear programming; MO-MILP: multi-objective mixed-integer linear programming; SS: scatter search.

Our proposal is aligned with the research methodology proposed by Liu et al. (2018) which, after finding the optimal SC design, applies an SD simulation model to evaluate the impact of government motivation measures. From the literature review, a research gap is identified where, to the best of our knowledge, no research proposes an optimisation model to address the LIT decision in a sustainable SC applied to a copper mining industry by complementarily using an SD simulation model to assess SC performance in different scenarios. According to Rebs et al. (2019), an SD model is chosen for its predominant strategic character and its relevance for modelling sustainable intra- and interorganisational SCs at the operational level. Hence the flexibility and simplicity with which business systems can be modelled by SD tools are highlighted by Sharif (2005). Nevertheless, Tako and Robinson (2012) do not find any evidence to support the notion that SD is used more for strategic problems than for operational/tactical issues.

Our proposal also considers the 3S for each decision, which is the main novelty of the MO-MILP model proposed by Becerra et al. (2023). Hence the proposed SD model can be considered an initial model, called a nominal model, of a digital twin, which is usually a physical-based model of the system that has been verified, validated and calibrated (Chakraborty et al., 2021).

Problem description

The MP model adopted as the basis of this work was proposed originally by Becerra et al. (2023), and is dubbed 3S-LIT. It corresponds to an MO-MINLP model whose objectives consider minimising economic, carbon emissions and social costs, and maximising the social impact of SC operations. It also incorporates LIT decisions in an integrated manner, along with the 3S aspects, into each named decision.

The model contemplates a single production plant and a single type of finished product: copper cathodes. A single waste type is considered, as is the recycling of only copper scrap from that waste. The disposal of non-recycled waste is not covered by the model. To measure route hazardousness, the number of people potentially affected by a possible accident on the route is considered. For inventory decisions, a fixed size of shipments defined by contracts with distributors is envisaged. Finally, shortages are taken as lost demand and are penalised due to non-compliance with contracts.

In this paper, copper ore suppliers and the supply of ore from the mine itself are both considered. Suppliers correspond to the small mining sites located in the north-central area of Chile, which supply a copper cathode production plant located in the municipality of Salamanca, Chile. After processing ore, the product, which corresponds to high-quality copper cathodes (LME grade A), is sent to distributors, which are responsible for globally marketing them to different customers. Our proposal incorporates elements into the global chain that are related to the circularity and recycling of copper waste, and which refer to any discarded items containing copper, to foster a more sustainable mining SC through the 3S-LIT model. To do so, we propose installing collection and repair centres that receive damaged products, which receive industrial products that have deteriorated during transport or old products discarded by end consumers and nearby communities to be either repaired or sent for treatment in recycling centres. Copper scrap refers to old, obsolete post-consumer or externally sourced scrap from copper waste, and is obtained at these recycling centres to be sent to warehouses for storage and subsequent shipment to the processing plant to be incorporated as raw material into a new process. Figure 1 presents the relation among damaged copper products, waste and scrap.

Figure 1. Damaged copper products, waste and scrap relations.

Figure 2 presents the structure of the modelled copper mining industry closed-loop SC (CLSC). As a novelty, here the simulation model aims to study how the behaviour of suppliers, distributors and RL facilities in ore supply, production lead times and demand terms affects the environmental, economic and social impacts of the CLSC. It also aims to observe CLSC performance when faced with the various disruptions and scenarios to which it may be exposed.

Figure 2. A sustainable copper mining CLSC structure.

The procedure followed in this research work is presented in Figure 3.

Figure 3. Research procedure from the MP model to the SD model.

Comparing current procedures in a copper mining industry to the 3S-LIT model

In this section, the current procedures of the SC in the copper mining industry are compared to those proposed by the 3S-LIT model, which is adopted as the basis of our SD model, through both MP models. We used data based on a real-world Chilean copper mining company. Table 2 presents each problem size considered in the present SC under study and in the CLSC by Becerra et al. (2023).

Table 2. Problem sizes.

Indices	Current SC	CLSC
Suppliers (I)	12	12
Processing centres (J)	1	2
Distributors (K)	1	2
Customers (R)	6	6
Collection and repair centres (C)	0	2
Recycling centres (D)	0	2
Scrap warehouses (S)	0	2
Raw materials (M)	1	1
Time periods (T)	12	12

Becerra et al. (2023) address a CLSC structure that incorporates all the 3S-LIT model entities by considering RL through collection and repair centres, recycling centres and scrap warehouses, where the material that flows backwards in the chain is incorporated into the production process as raw material. Suppliers are the copper ore producers with whom the company has contracts. The company's own mine is also considered an ore supplier. The shortest routes among locations are determined using the shortest route between the points algorithm in the QGIS^® geographic information system. The production plant has a processing centre with an operating capacity below its maximum capacity. This plant's idle capacity is due mainly to poor own ore production because the company's underground mine is not operating. Copper cathodes are distributed by a single company, which has a contract for a fixed order volume at a set price. The distribution company collects cathodes and then distributes them to its customers. Customer demand is presented in Table 3.

Table 3. Customer demand (tonnes of copper cathodes) per time period.

t	r = 1	r = 2	r = 3	r = 4	r = 5	r = 6
1	265.71	83.93	303.98	81.45	67.57	37.52
2	284.82	89.97	325.85	87.31	72.43	40.22
3	269.47	85.12	308.28	82.6	68.53	38.05
4	336.62	106.33	385.11	103.18	85.61	47.54
5	300.58	94.94	343.87	92.13	76.44	42.45
6	198.98	62.85	227.64	60.99	50.6	28.1
7	328.91	103.89	376.29	100.82	83.65	46.45
8	298.72	94.36	341.75	91.56	75.97	42.18
9	239.08	75.52	273.51	73.28	60.8	33.76
10	285.45	90.17	326.57	87.5	72.59	40.31
11	284.13	89.75	325.06	87.09	72.26	40.12
12	272.53	86.08	311.78	83.54	69.31	38.49

Moreover, the current SC of the industry under study corresponds to a cojoined structure with four steps: suppliers, production plant, distributors and consumers. These members of the current SC follow the same procedures as the previously described CLSC. To adopt the 3S-LIT to the present SC under study, it is modified by adjusting it to the current SC characteristics and, specifically, collection and repair centres, recycling centres and scrap warehouses are eliminated (Table 2).

Both the SCs, the current one based on Becerra et al. (2023) and the original 3S-LIT CLSC, are modelled using Pyomo in the Python language. Then they are solved by commercial solver Gurobi on a computational server with 256 Gb RAM and two AMD EPYC 7402 24-Core Processor processors, 2.80 GHz frequency and a Windows Server 2022 standard operating system. Using the same input data for both models, the following results are obtained and presented in Table 4.

Table 4. Comparison of the present SC to the CLSC when applying the 3S-LIT model.

	Current SC	CLSC	Cost savings	Increase in impacted people	% variation
Total costs (TEUR)	126056.144	47576.768	78479.38		−62.26%
Configuration cost (TEUR)	36136.068	43435.650	−7299.58		20.20%
Inventory costs (TEUR)	86548.833	267.639	86281.19		−99.69%
Order cost (TEUR)	10.768	64.653	−53.88		500.42%
Holding cost (TEUR)	78.065	184.376	−106.31		136.18%
Shortage cost (TEUR)	86460.000	18.610	86441.39		−99.98%
Transport cost (TEUR)	3371.244	3873.479	−502.24		14.90%
Total emissions costs (TEUR)	260646.119	260224.413	421.71		−0.16%
Facilities emissions costs (TEUR)	260481.123	259968.037	513.09		−0.20%
Inventory emissions costs (TEUR)	0.211	0.873	−0.66		313.59%
Transport emission costs (TEUR)	164.786	255.503	−90.72		55.05%
Social cost (TEUR)	92.897	247.370	−154.47		166.28%
Social Impact (people)	9190	12,732		3541.80	38.54%
Economic development (people)	9696	14,807		5111.00	52.71%
Route hazardousness (people)	506	2122		1616	319.36%

The results obtained when comparing both models show lower total costs due to the drop in inventory costs, specifically shortage costs. Although the inventory holding and ordering costs increase, it is not enough to compensate for the drop in inventory costs. Here copper demand is met by the cathodes produced in the plant and the cathodes repaired at the collection and repair centres (see Figure 4). There is a slight increase in configuration costs because new facilities are contemplated in the CLSC. Transport costs rise due to new routes. An increase in transport and inventory emissions is observed, but it is not enough to offset the reduction in emissions caused by installations, which implies slightly lower total emissions costs. Finally, social costs increase because, as the amount of stored product grows, the number of occupational accidents that results from these activities increases. However, the positive social impact is stronger because there are more direct and indirect jobs in the localities where facilities are located than the people who may be potentially affected on transport routes.

Figure 4. Copper demand, the current SC, and the CLSC produced and repaired cathodes production per time period.

Next a simulation model based on SD is used to test other alternative scenarios to validate the 3S-LIT model in SC performance terms in different circumstances. For this purpose, the optimal CLSC design is applied, obtained by applying the above-described MP model. The selected suppliers, operational facilities and customers served during the 12 time periods and represent the optimal CLSC design are presented in Table 5.

Table 5. Optimal CLSC design per time period.

	Time periods
	1	2	3	4	5	6	7	8	9	10	11	12
Suppliers (i)
supplier1	0	0	0	0	0	0	0	0	0	0	0	0
supplier2	1	1	1	1	1	0	1	1	1	1	1	1
supplier3	1	1	1	1	1	0	1	1	1	1	1	1
supplier4	1	1	1	1	1	0	1	1	1	1	1	1
supplier5	1	1	1	1	1	0	1	1	0	1	1	1
supplier6	1	1	1	1	1	0	1	1	1	1	1	1
supplier7	0	0	0	0	0	0	0	0	0	0	0	0
supplier8	1	1	1	1	1	0	1	1	1	1	1	1
supplier9	1	1	1	1	1	0	1	1	1	1	1	1
supplier10	1	1	1	1	1	0	1	1	1	1	1	1
supplier11	1	1	1	1	1	0	1	1	1	1	1	1
supplier12	1	1	1	1	1	0	1	1	0	1	1	1
Process Centres (j)
processCentre1	1	1	1	1	1	0	1	1	1	1	1	1
processCentre2	1	1	1	1	1	1	1	1	1	1	1	1
Distributors (k)
distributor1	1	1	1	1	1	1	1	1	1	1	1	1
distributor2	0	0	0	1	0	0	1	0	0	0	0	0
Customers (r)
customer1	1	1	1	1	1	1	1	1	1	1	1	1
customer2	1	1	1	1	1	1	1	1	1	1	1	1
customer3	1	1	1	1	1	1	1	1	1	1	1	1
customer4	1	1	1	1	1	1	1	1	1	1	1	1
customer5	1	1	1	1	1	1	1	1	1	1	1	1
customer6	1	1	1	1	1	1	1	1	1	1	1	1
Collect and repair centres (c)
collect1	1	1	1	1	1	1	1	1	1	1	1	1
collect2	1	1	1	1	1	1	1	1	1	1	1	1
collect3	1	1	1	1	1	1	1	1	1	1	1	1
Recycling centres (d)
recycling1	1	1	1	1	1	1	1	1	1	1	1	1
recycling2	0	0	0	0	0	0	0	0	0	0	0	0
Scrap Warehouses (s)
scrapwarehouse1	0	0	0	0	0	0	0	0	0	0	0	0
scrapwarehouse2	1	1	1	1	1	1	1	1	1	1	1	1

When simulating the CLSC, the same data are used as in the application of the 3S-LIT model to initialise the SD model, which allows a profounder analysis of potential SC improvements.

SD model formulation

The SD model is used to measure the performance of the optimal sustainable CLSC for the copper mining industry by applying the 3S-LIT model (Table 5). The simulation model is built with the following methodology: (i) construct the casual loop diagram; (ii) propose the flow diagram to represent the process; (iii) obtain the equations that define the SD model's behaviour; (iv) validate and perform SD.

The causal loop diagram (Figure 5) shows the cause–effect relations among the different system variables, which support their understanding to subsequently draw the flow diagram of the simulation model. As the causal loop diagram depicts, the quantity of produced copper cathodes remains constant when a push inventory strategy is adopted. The produced quantities are stored with distributors as on-hand inventory, which demonstrates their positive relation. Every time simulation begins, a quantity of available product remains in the distributor's warehouse, which constitutes the initial stock. The inventory on-hand with distributors decreases as consumer orders increase, which demonstrates their negative relation. The same occurs in scrap warehouses, where the orders to supply scrap to the plant and the scrap inventory are negatively related.

Figure 5. Causal loop of the 3S-LIT CLSC material and information flows.

When collecting damaged cathodes, which can be repaired or sent for recycling, the bigger the number of cathodes sent to customers, the more damaged cathodes are sent to collection and repair centres and, thus, the larger the number of repaired cathodes, which shows a positive relation. The same is also true for the copper ore supply from suppliers because, if production increases (if production capacity is higher), the ore required from suppliers will increase.

Finally, the larger loop, which includes the production of cathodes (i.e. the storage of cathodes with distributors, the cathodes sent to customers, the scrap sent to the collection centres, recycled scrap and the scrap stored in warehouses and, finally, the scrap sent to the plant to be processed) corresponds to a reinforcement loop. This means a directly proportional relation among the named elements.

From the causal loop diagram (Figure 5), the relations among the variables to be studied are identified in Table 6. This allows the next step in the process of modelling the CLSC to be done, namely the construction of the Forrester diagram (Figure 6) by the SD methodology (Campuzano & Mula, 2011), which represents the system under study and allows the simulation of the 3S-LIT model. To do this, we first identify the model's level, flow and auxiliary variables. Hence from the 3S-LIT model, the SD model incorporates the same objectives to evaluate CLSC performance as the level variables. Moreover, some of the 3S-LIT decision variables that define the CLSC design are taken as auxiliary variables in the SD model, but constraints are not considered because the optimal 3S-LIT results are used as initialisation values (Table 5). Finally, equations are introduced into the SD variables according to the original 3S-LIT model by Becerra et al. (2023).

Figure 6. Flow chart of the copper mining industry CLSC design by the 3S-LIT model.

Table 6. Notation.

Level variables
AIWk	Available inventory of finished product in distributors’ k warehouse (tonne)
AISWs	Available inventory of scrap in warehouse s (tonne)
Flow variables
Qjk	Quantity of finished product shipped from processing centre j in the plant to distributors k (tonne/month)
Qkr	Quantity of finished product shipped from distributors k to customers r (tonne/month)
Qds	Quantity of scrap shipped from recycling centres d to scrap warehouses s (tonne/month)
Qsj	Quantity of scrap shipped from scrap warehouses s to processing centre j in the plant (tonne/month)
Auxiliary variables
avgINk	Average inventory of finished product with distributors k per time period (tonne)
avgINs	Average inventory of scrap in warehouses s per time period (tonne)
capj	Production capacity of each processing centre j (tonne)
demr	Demand of finished product per customer r (tonne/month)
DPAc	It identifies if the product discarded by customers is assigned to collection and repair centre c (Dmnl)
OAkr	Finished product orders to distributor k from customer r (tonne)
ogi	Average ore grade of each supplier i (%)
OSAs	Orders to scrap warehouses s from the production plant (tonne)
PADk	It identifies if the produced quantity in the plant is allocated to distributor k (Dmnl)
Prod	Total produced quantity of finished product shipped to distributors (tonne)
Qd	Waste shipped to recycling centres d (tonne/month)
QDr	Products discarded by customer r (tonne/month)
Qi	Quantity of ore supplied by suppliers i to the plant (tonne/month)
QLCc	Local community waste collected by collection/repair centres c (tonne/month)
Qr	Repaired products shipped to customers r (tonne/month)
Qc	Waste from products discarded by customers shipped to collection/repair centres c (tonne/month)
QRc	Repaired products at each collection/repair centre c (tonne/month)
Qs	Quantity of scrap supplied by scrap warehouses s to the plant (tonne/month)
QURc	Unrepaired products at each collection/repair centre c (tonne/month)
recj	Metallurgical recovery of processing centre j = 1 (%)
RPAr	It identifies if the repaired product by collection and repair centres is assigned to customer r (Dmnl)
SAs	It identifies if the scrap generated by recycling centres is assigned to scrap warehouse s (Dmnl)
shr	Unmet demand of each customer r (tonne/month)
WAd	It identifies if the waste generated by collection/repair centres is assigned to recycling centre d (Dmnl)
αij	Conversion rate of the raw material from suppliers i in the finished product at processing centre j = 1 in the production plant (%)
γr	Percentage of the finished product discarded by customer r (%)
δc	Conversion rate of scrap into finished product at collection/repair centre c = 2 (%)
ϵc	Repair rate at collection/repair centre c of waste from customers (%)
θd	Percentage of waste recycled as scrap at recycling centre d (%)
FC	Total fixed configuration cost of facilities (€)
FCj	Fixed configuration costs of process centre j (€)
FCk	Fixed configuration costs of distributor k (€)
FCc	Fixed configuration costs of collection/repair centre c (€)
FCd	Fixed configuration costs of recycling centre d (€)
FCs	Fixed configuration costs of scrap warehouse s (€)
IC	Total inventory cost of distributors and scrap warehouses (€)
SHcost	Shortage cost of distributors for unmet demand (€/tonne)
hk	Holding costs of distributor k (€/tonne)
hs	Holding costs of scrap warehouse s (€/tonne)
ork	Ordering costs of distributor k (€/tonne)
ors	Ordering costs of scrap warehouse s (€/tonne)
TC	Total transport cost (€)
ctij	Unit transportation cost from supplier i to processing centres j in the production facility (€/tonne)
ctjk	Unit transportation cost from processing centres j in the production plant to distributor k (€/tonne)
ctkr	Unit transportation cost from distributor k to customer r (€/tonne)
ctrc	Unit transportation cost from customer r to collection/repair centre c (€/tonne)
ctcr	Unit transportation cost from collection/repair centre c to customer r (€/tonne)
ctcd	Unit transportation cost from collection/repair centre c to recycling centre d (€/tonne)
ctds	Unit transportation cost from recycling centre d to scrap warehouse s (€/tonne)
ctsj	Unit transportation cost from scrap warehouse s to processing centres j in the production facility (€/tonne)
FEN	Facilities environment cost (€)
ecc	Cost of the CO2 equivalent of GHG emissions per unit of waste processed at collection and repair centre c (€/tonne)
ecd	Cost of the CO2 equivalent of GHG emissions per unit of unrepaired waste processed at recycling centre d (€/tonne)
eck	Cost of the CO2 equivalent of GHG emissions per unit of finished product inventory held by distributor k (€/tonne)
ecmi	Cost of the CO2 equivalent of GHG emissions per unit of raw material m produced by supplier i (€/tonne)
ecmj	Cost of the CO2 equivalent of GHG emissions per unit of raw material m processed at processing centre j in the production facility (€/tonne)
ecs	Cost of CO2 equivalent of GHG emissions per unit of scrap inventory held by scrap warehouse s (€/tonne)
TEN	Total transport cost of the CO2 equivalent of GHG emissions (€)
etmij	Cost of the CO2 equivalent of GHG emissions per transported unit of raw material m from supplier i to processing centre j in the production facility (€/tonne)
etcr	Cost of the CO2 equivalent of GHG emissions per transported unit of finished product from collection/repair centre c to customer r (€/tonne)
etjk	Cost of the CO2 equivalent of GHG emissions per transported unit of finished product from processing centres j in the production facility to distributor k (€/tonne)
etkr	Cost of the CO2 equivalent of GHG emissions per transported unit of finished product from distributor k to customer r (€/tonne)
etcd	Cost of the CO2 equivalent of GHG emissions per transported unit of finished product from collection/repair centre c to recycling centre d (€/tonne)
etds	Cost of the CO2 equivalent of GHG emissions per transported unit of finished product from recycling centre d to scrap warehouse s (€/tonne)
etrc	Cost of the CO2 equivalent of GHG emissions per transported unit of finished product from customer r to collection/repair centre c (€/tonne)
etsj	Cost of the CO2 equivalent of GHG emissions per transported unit of finished product from scrap warehouse s to processing centre j in the production facility (€/tonne)
etmij	Cost of the CO2 equivalent of GHG emissions per transported unit of raw material m from supplier i to processing centre j in the production facility (€/tonne)
etcr	Cost of the CO2 equivalent of GHG emissions per transported unit of finished product from collection/repair centre c to customer r (€/tonne)
IEN	Total inventory emissions cost (€/)
eck	Cost of the CO2 equivalent of GHG emissions per unit of the finished product inventory held by distributor k (€/tonne)
ecs	Cost of the CO2 equivalent of GHG emissions per unit of scrap inventory held by scrap warehouse s (€/tonne)
SC	Total social cost (€)
irk	Injury cost per unit of the finished product in the inventory at distributor s (€/tonne)
irs	Injury cost per unit of scrap in the inventory in scrap warehouse k (€/tonne)
ED	Positive social impact (people)
edc	Economic development measured as the number of direct and indirect jobs generated by collection/repair centre c (people)
edd	Economic development measured as the number of direct and indirect jobs generated by recycling centre d (people)
edi	Economic development measured as the number of direct and indirect jobs generated by supplier i (people)
edj	Economic development measured as the number of direct and indirect jobs generated by centre processing j in the production facility (people)
edk	Economic development measured as the number of direct and indirect jobs generated by distributor k (people)
eds	Economic development measured as the number of direct and indirect jobs generated by scrap warehouse s (people)
HC	Negative social impact (people)
hccd	Route hazard factor from collection/repair centre c to recycling centre d (people)
hccr	Route hazard factor from collection/repair centre c to customer r (people)
hcds	Route hazard factor from recycling centre d to scrap warehouse s (people)
hcjk	Route hazard factor from processing centres j in the production facility to distributor k (people)
hckr	Route hazard factor from distributor k to customer r (people)
hcmij	Route hazard factor from supplier i to processing centre j in the production facility that transports raw material m (people)
hcsj	Route hazard factor from scrap warehouse s to processing centre j in the production facility (people)

Figure 6 depicts the flow diagram of the CLSC designed with the 3S-LIT model. This model allows experiments to study the CLSC's performance and behaviour in different scenarios. The Vensim® simulation software is used to implement the model on a computational server with 256 Gb RAM and two AMD EPYC 7402 24-Core Processor processors, 2.80 GHz frequency and a Windows Server 2022 standard operating system. Table 6 presents the notations and descriptions of the variables employed in the model; that is, the variables that represent the SC structure and the SC economic, emissions and social costs. The relations of the variables that allow the models’ objectives to be calculated (i.e. total costs, total emissions, social impacts) are presented in the equations describing the model.

When simulation starts, an initial inventory of zero is considered. For the following time periods, suppliers provide the production plant with ore, the final product is transported to distributors, which manage the inventory, and then customers are served the final product.

Products discarded as scrap are repaired and sent to customers or are recycled in recycling centres. Recycled waste is transformed into copper scrap and then stored in scrap warehouses, where its inventory is managed and sent to the processing plant.

The product orders per customer to each distributor, as well as the scrap orders to each warehouse, are known when simulation starts. In the same way, the product allocations among the different facilities are previously determined from the MP model.

In the first instance, the model replicates the optimal CLSC design defined by the 3S-LIT model, the allocation of orders and the selection of the facilities to be operated. By simulating the chain's material and information flow, at the same time CLSC performance is calculated from the economic, environmental and social objectives set out in the 3S-LIT model.

To gain a better understanding of the simulated CLSC, the equations comprising the model are presented.

Equation (1) defines the quantity of finished product to be produced in the plant that is to be shipped to distributors during each time period t. It is calculated as the sum of the production of each processing centre j, multiplied by its allocation in each distributor k. In this case, we consider two processing centres. (1) $\begin{array}{r}{Q}_{\mathrm{jk}}(t)=\text{PA}{\text{D}}_k(t)\ast {\displaystyle \sum _{j=1}^2\text{pro}{\text{d}}_j(t);\mathrm{\forall k}}\end{array}$(1) Equations (2) and (3) define the calculation of the production at each processing centre j. In both cases, production is conditioned by the production capacity of each j. Here two processing centres and 12 ore suppliers are considered. (2) $\begin{array}{rl}\text{Pro}{\text{d}}_{j=1}(t)& =\{\begin{array}{c}\text{ca}{\text{p}}_1,\text{ if }{\displaystyle \sum _{i=1}^{12}{\alpha }_i\ast Q_i(t)>\text{ca}{\text{p}}_1}\\ {\displaystyle \sum _{i=1}^{12}{\alpha }_i\ast Q_i(t),\text{ otherwise }}\end{array}\end{array}$(2) (3) $\begin{array}{rl}\text{Pro}{\text{d}}_{j=2}(t)& =\{\begin{array}{c}\text{ca}{\text{p}}_2,\mathrm{if}{\displaystyle \sum _{s=1}^2{\delta }_s\ast Q_s(t)>\text{ca}{\text{p}}_2}\\ {\displaystyle \sum _{s=1}^2{\delta }_s\ast Q_s(t),\text{ otherwise }}\end{array}\end{array}$(3) Equation (4) defines the level variable AIW_k as the amount of finished product inventory on hand at the beginning of time period t. It is defined as the quantity of product that enters in each distributor k, minus the quantity served to each customer r. Here the initial inventory on hand is considered null. (4) $\begin{array}{rl}& \text{AI}{\text{W}}_k(t)={\int }_{t_0}^t\left[{\displaystyle \sum _j{Q}_{\mathrm{jk}}(t)-{\displaystyle \sum _r{Q}_{\mathrm{kr}}(t)}}\right]\mathrm{dt};\\ & \text{AI}{\text{W}}_k(t_0)=0\mathrm{\forall k}\end{array}$(4)

The quantity of finished product shipped from each distributor k to customer r is defined in Equation (5) as a finished product replenishment order. (5) $\begin{array}{r}{Q}_{\mathrm{kr}}(t)=\{\begin{array}{c}{\displaystyle \sum _{r}O{A}_{\mathrm{kr}}(t),\text{ if }{\displaystyle \sum _{r}O{A}_{\mathrm{kr}}(t)}}\\ \le \text{AI}{\text{W}}_k(t)+{\displaystyle \sum _j{Q}_{\mathrm{jk}}(t)}\\ \text{AI}{\text{W}}_k(t)+{\displaystyle \sum _j{Q}_{\mathrm{jk}}(t),\text{ otherwise }}\end{array};\mathrm{\forall k}\end{array}$(5) The average inventory is also defined in Equation (6) as the final inventory in distributor k, minus the initial inventory on hand in k, which is the quantity of finished product received from processing centres j, divided by two because we consider the simplest average inventory formula. (6) $\begin{array}{r}\text{avgI}{\text{N}}_k(t)=\frac{{\sum }_j{Q}_{\mathrm{jk}}(t)-\text{AI}{\text{W}}_k(t)}{2};\mathrm{\forall k}\end{array}$(6)

Equation (7) defines the level variable AISW_s as the amount of scrap inventory on-hand in at the beginning of time period t. It is defined as the quantity of scrap that enters each warehouse s, minus the quantity served to the production plant. Here the initial scrap inventory on-hand is considered null. (7) $\begin{array}{rl}& \text{AIS}{\text{W}}_s(t)={\int }_{t_0}^t\left[{\displaystyle \sum _d{Q}_{\mathrm{ds}}(t)-{\displaystyle \sum _j{Q}_{\mathrm{sj}}(t)}}\right]\mathrm{dt};\\ & \text{AIS}{\text{W}}_s(t_0)=0,\mathrm{\forall s}\end{array}$(7) The quantity of scrap shipped from recycling centres d to scrap warehouses s is defined in Equation (8) as the product of the total waste processed as scrap at recycling centres d and its allocation in each scrap warehouse s. (8) $\begin{array}{r}{Q}_{\mathrm{ds}}(t)={\displaystyle \sum _d{Q}_d(t)\ast S{A}_s(t);\mathrm{\forall s}}\end{array}$(8)

Equation (9) defines the amount of scrap shipped from scrap warehouses s to the processing plant as a scrap product replenishment order. (9) $\begin{array}{r}{Q}_{\mathrm{sj}}(t)=\{\begin{array}{c}\text{OS}{\text{A}}_s(t),\text{if OS}{\text{A}}_s(t)\le \text{AIS}{\text{W}}_s(t)\\ +{\displaystyle \sum _d{Q}_{\mathrm{ds}}(t)}\\ \text{AIS}{\text{W}}_s(t)+{\displaystyle \sum _d{Q}_{\mathrm{ds}}(t),\text{otherwise }}\end{array};\mathrm{\forall s}\end{array}$(9) The average scrap inventory is also defined in Equation (10) as the final inventory of scrap, minus the initial scrap inventory on-hand, divided by two, because we consider the simplest average inventory formula. (10) $\begin{array}{r}\text{avgI}{\text{N}}_s(t)=\frac{\text{I}{\text{N}}_s(t)-\text{AIS}{\text{W}}_s(t)}{2};\mathrm{\forall s}\end{array}$(10)

To determine the amount of waste shipped to collection and repair centres c, we use auxiliary variable Q_rc, which is defined in Equation (11) as the product of the total amount of product discarded by customers and its allocation to each collection and repair centre c. (11) $\begin{array}{r}{Q}_c(t)=\text{DP}{\text{A}}_c(t)\ast {\displaystyle \sum _{r}Q{D}_r(t);\mathrm{\forall c}}\end{array}$(11) The quantity of discarded product is defined in Equation (12) as the portion of finished product discarded by customer r, multiplied by the total quantity of the finished product served to customer r. (12) $\begin{array}{r}\text{Q}{\text{D}}_r(t)={\gamma }_r\ast {\displaystyle \sum _k{Q}_{\mathrm{kr}}(t);\mathrm{\forall r}}\end{array}$(12)

Equation (13) defines the amount of repaired products. It is calculated as the product of the repair rate and the amount of waste processed by collection and repair centres c. (13) $\begin{array}{r}\text{Q}{\text{R}}_c(t)={\epsilon }_c\ast \left({\displaystyle \sum _r{Q}_{\mathrm{rc}}(t)+\mathrm{Ql}{c}_c(t)}\right);\mathrm{\forall c}\end{array}$(13) The repaired products shipped to customers r are defined in Equation (14) as the product of the total amount of repaired products and their allocation to each customer r. (14) $\begin{array}{r}{Q}_r(t)=\text{RP}{\text{A}}_r(t)\ast {\displaystyle \sum _{c}Q{R}_c(t);\mathrm{\forall c}}\end{array}$(14)

Equation (15) defines the amount of product sent to be recycled. It is calculated as the multiplication of the unrepair rate and the amount of waste processed by collection and repair centres c. (15) $\begin{array}{r}\text{QU}{\text{R}}_c(t)=(1-{\epsilon }_c)\ast \left({\displaystyle \sum _r{Q}_{\mathrm{rc}}(t)+\text{Ql}{\text{c}}_c(t)}\right);\mathrm{\forall c}\end{array}$(15) An unrepaired product is shipped to recycling centres d as waste and is defined in Equation (16) as the product of the total amount of unrepaired product and its allocation to each recycling centre d. (16) $\begin{array}{r}{Q}_d(t)=\text{W}{\text{A}}_d(t)\ast {\displaystyle \sum _c\text{QU}{\text{R}}_c(t);\mathrm{\forall d}}\end{array}$(16)

Equations (17–30) provide the economic results obtained during simulation. Equation (17) defines the total economic costs. This is calculated as the sum of the fixed configuration, inventory and transportation costs. (17) $\begin{array}{r}\text{EC}(t)=\text{FC}(t)+\text{ IC}(t)+\text{TC}(t)\end{array}$(17) Fixed configuration costs are defined in Equation (18) as the sum of the fixed costs of facilities j, k, c, d and s. (18) $\begin{array}{rl}\text{FC}(t)& ={\int }_{t_0}^t[{\displaystyle \sum _j\text{FC}{\text{J}}_j(t)+{\displaystyle \sum _k\text{FC}{\text{K}}_k(t)}}\\ & +{\displaystyle \sum _c\text{FC}{\text{C}}_c(t)+{\displaystyle \sum _d\text{FC}{\text{D}}_d(t)}}\\ & +{\displaystyle \sum _s\text{FC}{\text{S}}_s(t)}]\mathrm{dt};\text{FC}(t_0)=0\end{array}$(18)

Here each fixed cost per time period is taken into account if facilities j, k, c, d and s are selected. Equations (19–23) determine how each fixed cost is considered if facilities j, k, c, d and s are operating. Here fixed costs are considered null if there is no production or there are no quantities of finished products or scrap to be shipped through the different facilities. (19) $\begin{array}{rl}\text{FC}{\text{J}}_j(t)& =\{\begin{array}{c}\text{F}{\text{C}}_j,\text{if pro}{\text{d}}_j>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall j}\end{array}$(19) (20) $\begin{array}{rl}\text{FC}{\text{K}}_k(t)& =\{\begin{array}{ll}\text{F}{\text{C}}_k,& \text{if }{\displaystyle \sum _j{Q}_{\mathrm{jk}}>0}\\ 0,& \text{otherwise}\end{array};\mathrm{\forall k}\end{array}$(20) (21) $\begin{array}{rl}\text{FC}{\text{C}}_c(t)& =\{\begin{array}{c}\text{F}{\text{C}}_c,\text{ if }{Q}_c>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall c}\end{array}$(21) (22) $\begin{array}{rl}\text{FC}{\text{D}}_d(t)& =\{\begin{array}{c}\text{F}{\text{C}}_d,\text{ if }{Q}_d>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall d}\end{array}$(22) (23) $\begin{array}{rl}\text{FC}{\text{S}}_s(t)& =\{\begin{array}{c}\text{F}{\text{C}}_s,\text{ if }{\displaystyle \sum _d{Q}_{\mathrm{ds}}>0}\\ 0,\text{ otherwise }\end{array};\mathrm{\forall s}\end{array}$(23)

To define inventory costs, holding inventory, ordering and shortage costs are taken into account. Equation (24) calculates the total inventory costs for distributors k and scrap warehouses s per time period. $\begin{array}{rl}\text{IC}(t)& ={\int }_{t_0}^t[\text{SH }(t)+{\displaystyle \sum _k{H}_k(t)+\text{O}{\text{R}}_k(t)}\\ & +{\displaystyle \sum _s{H}_s(t)+\text{O}{\text{R}}_s(t)}]\mathrm{dt};\end{array}$ (24) $\begin{array}{rl}& \text{IC}(t_0)=0\end{array}$(24)

Equations (25–29) define the calculation of each inventory cost element, the cost of holding inventory from the average inventory, the cost of ordering from the number of orders and the shortage costs from unmet demand. (25) $\begin{array}{rl}\text{SH}(t)& =\mathrm{sh}(t)\ast \text{SHcost}\end{array}$(25) (26) $\begin{array}{rl}{H}_k(t)& =\text{avgI}{\text{N}}_k(t)\ast h_k,\mathrm{\forall k}\end{array}$(26) (27) $\begin{array}{rl}\text{O}{\text{R}}_k(t)& ={\displaystyle \sum _r{Q}_{\mathrm{kr}}(t)\ast \text{o}{\text{r}}_k,\mathrm{\forall k}}\end{array}$(27) (28) $\begin{array}{rl}{H}_s(t)& =\text{avgI}{\text{N}}_s(t)\ast h_s,\mathrm{\forall s}\end{array}$(28) (29) $\begin{array}{rl}\text{O}{\text{R}}_s(t)& ={\displaystyle \sum _d{Q}_{\mathrm{ds}}(t)\ast \text{o}{\text{r}}_s,\mathrm{\forall s}}\end{array}$(29)

Equation (30) describes the sum of all the transport costs, defined as the product of the unit cost of transporting material on the defined route and the quantity of material transported on that route. $\begin{array}{rl}\text{TC}(t)& ={\int }_{t_0}^t[{\displaystyle \sum _i{\displaystyle \sum _{j}c{t}_{\mathrm{ij}}\ast Q_i(t)+{\displaystyle \sum _j{\displaystyle \sum _{k}c{t}_{\mathrm{jk}}\ast Q_{\mathrm{jk}}(t)}}}}\\ & +{\displaystyle \sum _k{\displaystyle \sum _{r}c{t}_{\mathrm{kr}}\ast Q_{\mathrm{kr}}(t)}}\\ & +{\displaystyle \sum _r{\displaystyle \sum _{c}c{t}_{\mathrm{rc}}\ast Q_c(t)}}\\ & +{\displaystyle \sum _c{\displaystyle \sum _{r}c{t}_{\mathrm{cr}}\ast Q_r(t)+{\displaystyle \sum _c{\displaystyle \sum _{d}c{t}_{\mathrm{cd}}\ast Q_d(t)}}}}\\ & +{\displaystyle \sum _d{\displaystyle \sum _{s}c{t}_{\mathrm{ds}}\ast Q_{\mathrm{ds}}(t)}}\\ & +{\displaystyle \sum _s{\displaystyle \sum _{j}c{t}_{\mathrm{sj}}\ast Q_s(t)}}]\mathrm{dt};\end{array}$ (30) $\begin{array}{rl}& \mathrm{TC}(t_0)=0\end{array}$(30)

Equations (31–34) provide the environmental results obtained during simulation. Equation (31) defines the total environmental costs, calculated as the sum of the facility, inventory and transportation emissions costs. (31) $\begin{array}{r}\text{EN}(t)=\text{FEN}(t)+\text{TEN}(t)+\text{IEN}(t)\end{array}$(31) Emissions from facilities i, j, c and s are defined in Equation (32) as the unit cost of the emissions generated by the material processed at these facilities. (32) $\begin{array}{rl}\text{FEN}(t)& ={\int }_{t_0}^t[{\displaystyle \sum _{i}e{c}_i\ast Q_i(t)+{\displaystyle \sum _{j}e{c}_j\ast \text{pro}{\text{d}}_j(t)}}\\ & +{\displaystyle \sum _{c}e{c}_c\ast Q_c(t)}\\ & +{\displaystyle \sum _{s}e{c}_s\ast Q_s(t)}]\mathrm{dt};\text{FEN}(t_0)=0\end{array}$(32)

Equation (33) describes the transport emission costs, defined as the sum of the product of the unit emission costs of transporting material on the defined route and the quantity of material transported on that route. (33) $\begin{array}{rl}\text{TEN}(t)& ={\int }_{t_0}^t[{\displaystyle \sum _i{\displaystyle \sum _{j}e{t}_{\mathrm{ij}}\ast Q_i(t)}}\\ & +{\displaystyle \sum _j{\displaystyle \sum _{k}e{t}_{\mathrm{jk}}\ast Q_{\mathrm{jk}}(t)}}\\ & +{\displaystyle \sum _k{\displaystyle \sum _{r}e{t}_{\mathrm{kr}}\ast Q_{\mathrm{kr}}(t)+{\displaystyle \sum _r{\displaystyle \sum _{c}e{t}_{\mathrm{rc}}\ast Q_c(t)}}}}\\ & +{\displaystyle \sum _c{\displaystyle \sum _{r}e{t}_{\mathrm{cr}}\ast Q_r(t)}}\\ & +{\displaystyle \sum _c{\displaystyle \sum _{d}e{t}_{\mathrm{cd}}\ast Q_d(t)+{\displaystyle \sum _d{\displaystyle \sum _{s}e{t}_{\mathrm{ds}}\ast Q_{\mathrm{ds}}(t)}}}}\\ & +{\displaystyle \sum _s{\displaystyle \sum _{j}e{t}_{\mathrm{sj}}\ast Q_s(t)}}]\mathrm{dt};\text{TEN}(t_0)=0\end{array}$(33) Equation (34) describes the emissions generated by inventory activities at distributors k and scrap warehouses s, defined as the product of the average inventory and the unit emission costs, plus the inventory obsolescence cost, with their respective obsolescence rates. (34) $\begin{array}{rl}\text{ IEN}(t)& ={\int }_{t_0}^t[{\displaystyle \sum _k\text{avgI}{\text{N}}_k(t)\ast [e{c}_k+\beta \ast \text{Cob}{\text{s}}_k]}\\ & +{\displaystyle \sum _s\text{avgI}{\text{N}}_s(t)\ast [e{c}_s+\beta \ast \text{Cob}{\text{s}}_s]}]\mathrm{dt};\\ & \text{IEN}(t_0)=\mathrm{\forall k}\end{array}$(34)

Regarding the social impact objectives, the social cost is expressed in Equation (35) and defined as the injury co per unit of held inventory, multiplied by the average inventory for both distributors k and scrap warehouses s. (35) $\begin{array}{rl}\text{SC}(t)& ={\int }_{t_0}^t[{\displaystyle \sum _k\text{avgI}{\text{N}}_k(t)\ast \text{I}{\text{R}}_k}\\ & +{\displaystyle \sum _s\text{avgI}{\text{N}}_s(t)\ast I{R}_s}]\mathrm{dt};\text{SC}(t_0)=0\end{array}$(35) Another proposal is to measure the positive social impact defined in Equation (36) as the sum of the amount of the direct and indirect jobs generated by operational facilities i, j, k, c, d and s. The negative social impact is defined in Equation (37) as the sum of the number of people potentially affected by transporting material on populated routes. (36) $\begin{array}{rl}\text{ED}(t)& ={\int }_{t_0}^t[{\displaystyle \sum _i\text{ED}{\text{I}}_i(t)+{\displaystyle \sum _j\text{ED}{\text{J}}_j(t)+{\displaystyle \sum _k\text{ED}{\text{K}}_k(t)}}}\\ & +{\displaystyle \sum _c\text{ED}{\text{C}}_c(t)+{\displaystyle \sum _d\text{ED}{\text{D}}_d(t)}}\\ & +{\displaystyle \sum _s\text{ED}{\text{S}}_s(t)}]\mathrm{dt};\mathrm{ED}(t_0)=0\end{array}$(36) (37) $\begin{array}{rl}\text{HC}(t)& ={\int }_{t_0}^t[\text{H}{\text{C}}_{\mathrm{ij}}(t)+\text{H}{\text{C}}_{\mathrm{jk}}(t)+\text{H}{\text{C}}_{\mathrm{kr}}(t)+\text{H}{\text{C}}_{\mathrm{rc}}(t)\\ & +\text{H}{\text{C}}_{\mathrm{cr}}(t)+\text{H}{\text{C}}_{\mathrm{cd}}(t)+\text{H}{\text{C}}_{\mathrm{ds}}(t)\\ & +\text{H}{\text{C}}_{\mathrm{sj}}(t)]\mathrm{dt};\text{HC}(t_0)=0\end{array}$(37) Each positive social impact is considered if facilities i, j, k, c, d and s are operating. Equations (38–43) describe the decision of how employment development is considered. Here social impacts are considered null if there is no production or if there are no quantities of finished products or scrap to be shipped through the different facilities. (38) $\begin{array}{rl}\text{ED}{\text{I}}_j(t)& =\{\begin{array}{c}e{d}_i,\text{ if }{Q}_i>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall i}\end{array}$(38) (39) $\begin{array}{rl}\text{ED}{\text{J}}_j(t)& =\{\begin{array}{c}e{d}_j,\text{ if }\text{pro}{\text{d}}_j>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall j}\end{array}$(39) (40) $\begin{array}{rl}\text{ED}{\text{K}}_k(t)& =\{\begin{array}{c}e{d}_k,\text{ if}{\displaystyle \sum _j{Q}_{\mathrm{jk}}>0}\\ 0,\text{ otherwise }\end{array};\mathrm{\forall k}\end{array}$(40) (41) $\begin{array}{rl}\text{ED}{\text{C}}_c(t)& =\{\begin{array}{c}e{d}_c,\text{ if }{Q}_c>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall c}\end{array}$(41) (42) $\begin{array}{rl}\text{ED}{\text{D}}_d(t)& =\{\begin{array}{c}e{d}_d,\text{ if }{Q}_d>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall d}\end{array}$(42) (43) $\begin{array}{rl}\text{ED}{\text{S}}_s(t)& =\{\begin{array}{ll}e{d}_s,& \text{if }{\displaystyle \sum _d{Q}_{\mathrm{ds}}>0}\\ 0,& \text{otherwise}\end{array};\mathrm{\forall s}\end{array}$(43)

Equations (44–51) determine how each hazardousness is considered if the route is selected. The number of potentially affected people is taken into account if material is transported among SC nodes i, j, k, c, d, s and r. Here social impacts are considered null if there is no ore (Q_i), quantities of finished products (Q_jk, Q_kr), damaged products (QD_r), repaired products (QR_c), waste (QUR_c) or scrap (Q_ds, Q_sj) to be shipped through the different facilities. (44) $\begin{array}{rl}\text{H}{\text{C}}_{\mathrm{ij}}(t)& =\{\begin{array}{c}h{c}_{\mathrm{ij}},\text{ if }{Q}_i>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall i}\end{array}$(44) (45) $\begin{array}{rl}\text{H}{\text{C}}_{\mathrm{jk}}(t)& =\{\begin{array}{c}h{c}_{\mathrm{jk}},\text{ if }{Q}_{\mathrm{jk}}>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall j},k\end{array}$(45) (46) $\begin{array}{rl}\text{H}{\text{C}}_{\mathrm{kr}}(t)& =\{\begin{array}{c}h{c}_{\mathrm{kr}},\text{ i}f{Q}_{\mathrm{kr}}>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall k},r\end{array}$(46) (47) $\begin{array}{rl}\text{H}{\text{C}}_{\mathrm{rc}}(t)& =\{\begin{array}{ll}h{c}_{\mathrm{rc}},& \text{if }Q{D}_r\ast \mathrm{DP}{A}_c>0\\ 0,& \text{otherwise}\end{array};\mathrm{\forall r},c\end{array}$(47) (48) $\begin{array}{rl}\text{H}{\text{C}}_{\mathrm{cr}}(t)& =\{\begin{array}{ll}h{c}_{\mathrm{cr}},& \text{if }\text{Q}{\text{R}}_c\ast \text{RP}{\text{A}}_r>0\\ 0,& \text{otherwise }\end{array};\mathrm{\forall c},r\end{array}$(48) (49) $\begin{array}{rl}\text{H}{\text{C}}_{\mathrm{cd}}(t)& =\{\begin{array}{ll}h{c}_{\mathrm{cd}},& \text{if }\text{QU}{\text{R}}_c\ast \text{W}{\text{A}}_d>0\\ 0,& \text{otherwise }\end{array};\mathrm{\forall c},d\end{array}$(49) (50) $\begin{array}{rl}\text{H}{\text{C}}_{\mathrm{ds}}(t)& =\{\begin{array}{c}h{c}_{\mathrm{ds}},\text{ if }{Q}_{\mathrm{ds}}>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall d},s\end{array}$(50) (51) $\begin{array}{rl}\text{H}{\text{C}}_{\mathrm{sj}}(t)& =\{\begin{array}{c}h{c}_{\mathrm{sj}},\text{ if}{Q}_{\mathrm{sj}}>0\\ 0,\text{ otherwise }\end{array};\mathrm{\forall s},j\end{array}$(51)

Variables DPA_c, OAk_r, PAD_k, OSA_s, RPA_r, SAS_s, WA_d are also identified as key for determining the system's performance because they define which entities of SC products or waste are allocated and, thus, denote the route by which these materials are transported. This decision conditions the system's performance for several reasons: (i) to define the amount of material that is allocated to each facility to, thus, determine whether fixed costs are incurred, emissions are generated and the social impact of that facility; (ii) to define the amount to be transported between SC nodes to, thus, establish transport costs, transport emissions and the routes along which these materials are transported; (iii) the amount of material that is allocated to each distributor or scrap warehouse defines the results associated with the inventory management at these facilities.

Applying the SD model

The simulation proposal is applied to the same copper mining SC using the same values for the auxiliary variables and the initial values for the level variables as the MP model to compare them both. In addition, the following assumptions are considered when simulating the model:

1. − The simulation run length is fixed to 12 monthly time periods

2. − The simulation objectives are to minimise the total costs, minimise any negative environmental impacts and maximise the positive social impacts

3. − Selection of locations and allocation of material to different facilities are previously defined in the MP model

4. − Copper cathodes demand varies between 28 and 330 tonnes per month, depending on the customer

5. − Push system inventory management is the considered strategy. This means constant production to generate inventory

6. − At the start of simulation, the initial inventories with distributors and in scrap warehouses are considered null

Validation

In order to validate the simulation model, three of the tests proposed by Sterman (2000) are run. The first corresponds to a dimensional consistency test to check that the units used in models are consistent and correct. Then reproduction of known behaviours is tested, which consists of comparing the results obtained with the MP model. Finally, an extreme-conditions test is applied to analyse the model's robustness in two situations: with production equalling zero and with demand equalling zero. In these extreme situations, the inventory on-hand with distributors and inventory shortages behave as expected (Figure 7). Unlike the MP model, the simulation model defines location selection if the location is operational and if there is a material flow to or from the facility. This generates an adjustment in the results because, while seeking to optimise the positive social impact, the MP model considers the operation of some facilities, even if they do not process any material. The comparison results are shown in Table 7 for 12 and 24 running time periods.

Figure 7. (a) Inventory shortage when production equals zero in the extreme-condition test. (b) The total available inventory in distributor warehouses when demand equals zero in the extreme-condition test.

Table 7. Comparison of the CLSC MP and SD results.

	12 time periods		24 time periods
Variable	MP model (MO-MILP)	Simulation model (SD)	MP model (MO-MILP)	Simulation model (SD)
EC (TEUR)	47576.768	47540.8	94431.785	95124.6
EN (TEUR)	260224.413	260,320	521235.465	520,639
SC (TEUR)	247.370	206.602	590.397	464.444
Social Impact (People)	12,685	12273.78	25,383.43	24,547.55
ED (People)	14,807	14,300	29,416	28,600
HC (People)	2122	2026.22	4032.57	4052.45
Resolution time (sec)	26.584	2	4738.018	3.41

EC: total economic costs; EN: total emissions costs; SC: social costs; ED: employment development; HC: route hazardousness.

As previously mentioned, the simulation model presents similar results to the MP model by validating its results. It is also observed how the total costs lower because the operation of some facilities is not considered, which reduces the positive social impact. This spells an advantage for the simulation model because, by relying on flows in different facilities, the real cost generated in them is measured. Moreover, the simulation model has an advantage for the computational time, which is much shorter when more time periods are considered (Table 6).

Sensitive analysis

Sensitivity analyses are performed with the Vensim^® ‘Monte Carlo simulation’ tool. The parameter to be analysed is the percentage of waste transformed into scrap and its influence on the inventory cost per time period. For them, a minimum and a maximum value (0–100%) are selected, and uniform random distribution is chosen to generate the values to be used in 200 simulations. There are several reasons why the amount of waste that can be transformed into scrap is affected. Some have to do with technological changes, waste composition, origin of waste, etc. When analysing the results obtained by the variation in theta (θ), the impact on inventory costs is related to changes in the inventory costs per time period, as shown in Figure 8. A sensitivity analysis is also performed of the average emissions costs by varying the emissions generated by suppliers, EC_i, which may be possible due to technological improvements in extraction processes. For them, a minimum and a maximum value (100–374.4) are selected, and uniform random distribution is chosen to generate the values to be used in 200 simulations. As shown in Figure 9, the average emissions costs very much depend on this parameter.

Figure 8. Sensitivity analysis of the inventory costs per time period with variation in the percentage of waste converted into scrap.

Figure 9. Sensitivity analysis of the average total emissions costs with variation in the emissions generated by suppliers.

In addition, further sensitivity tests were run with several parameters (production capacity, repair rate, percentage of customer waste, etc.). It is concluded that these parameters do not have a significant effect on the model and suppliers’ emissions cost, EC_i,, has the strongest impact on the total emissions costs. Due to space requirements, these sensitivity analyses are not provided here.

Simulating scenarios

We propose a series of scenarios that represent a risk of SC disruption by considering the risk presented by Llaguno et al. (2022). Scenario 1 contemplates the closure of the company's own mine and supply only by ore suppliers, i.e. Q_i is null for the mine itself. Scenario 2 considers a sudden increase in demand, dem_r, of 20% from month 6 onwards. Scenario 3 recreates a situation in which a lead time (LT) of 1 month in production is taken, unlike the MP model where fixing SC lead times are more complex because there are several SC levels. So the model should be modified by incorporating scheduled arrivals at the different levels to take into account materials’ real arrival times. Scenario 4 recreates a situation that contemplates an uncertain lead time with uniform integer random distribution between 1 and 3 months. In this case, it becomes even more complicated when considering a rolling horizon to apply uncertain delay times in an MP model (Díaz-Madroñero et al., 2017). Table 8 summarises the different conditions in each scenario in relation to the initial scenario.

Table 8. Conditions of scenarios.

	Scenario 1	Scenario 2	Scenario 3	Scenario 4
demr	–	+20%	–	–
Qi	Q1 = 0	–	–	–
LT	–	–	1	1–3

For each scenario, the inventory shortage, the total economic costs, the total emission costs, social costs and social impacts are analysed. Table 9 presents the results obtained from running simulations. This set of scenarios is designed to assess the effect of potential disruptions on the SC.

Table 9. Results of scenarios.

Variable	Scenario 1	Scenario 2	Scenario 3	Scenario 4
EC (TEUR)	91,176	80961.8	65236.1	78632.3
EN (TEUR)	110,568	260,321	260,319	260,319
SC (TEUR)	121.906	184.411	192.127	258.079
ED (People)	14,288	14,300	14,290	14,270
HC (People)	2026.22	2026.22	2023.57	2018.25
Inventory shortage (Tonne)	2675.05	2026.87	1079.54	1911.12

EC: total economic costs; EN: total emissions costs; SC: social costs; ED: employment development; HC: route hazardousness.

From the results, it is concluded that a disruption in the SC, regardless of it being due to the own mine closing, a sudden increase in demand or delay in production, would increase the total economic costs, caused by increased inventory shortages, which would be critical output to be managed with possible disruptions. In Scenario 1, by closing down the own mine, the environmental impact would reduce because this activity is the largest source of the CO₂ emissions in this SC. In Scenario 4, social costs are higher because the amount of inventory handled by distributors is bigger.

Two potential disruptive risks are characterised by exceptional events, whose effect is reflected on SC performance, and are considered to observe the possible ripple effect (Ghadge et al., 2022). Here the effects of disruptions are taken into account due to risks to the own mine's ore supply and the market risk of increasing demand at some point during the simulation period. These scenarios are simulated over a longer time horizon, 36 months, to identify whether the SC is resilient.

Disruption due to supply risk

This disruption type, which represents the closure of the own mine due to either a strike or a natural disaster (e.g. earthquake), is simulated by considering zero production from the own mine for a 3-month period on the simulation time horizon. The behaviour of production in the processing plant, the available inventory and unmet demand through inventory shortage all appear in Figure 10.

Figure 10. SC disruption propagation due to supply risk.

A disruption to the ore supply by the mine's own mine shows the SC's resilience because it has other suppliers. Hence the SC compensates for this reduction in production with recycled scrap and with the suppliers not considered when its own mine operates. So available inventory, plus the quantity of production and inventory shortages, tend to stabilise over time. Figure 11 shows the impact on the performance results with regular supply and when this supply is interrupted. A reduction in emission costs is observed due to the mine not operating, which with a major source of CO₂ emissions. However, an increase in economic costs occurs, which tends to stabilise over time periods.

Figure 11. Economic and environmental costs with disruption due to supply risk.

Disruption due to the market risk

Scenario 2 represents a market risk due to increased demand from month 6 onwards. However, to observe the SC's long-term behaviour, two situations are considered: one with a 10% increase and the other with a 20% increase in demand from year 2 onwards. A change in demand implies a risk because the production rate is constant. This impacts both the inventory available held at distributors and the inventory shortage generated by increased demand, as seen in Figure 12.

Figure 12. SC disruption propagation due to the market risk.

Furthermore, the behaviour of the SC performance results on economic costs is seen in Figure 13. Increased demand results in higher economic costs due to inventory shortages because the production plant operates at maximum capacity, but cannot adapt to this situation.

Figure 13. Economic costs with disruption due to the market risk.

Both studied risks correspond to disruptions that can cause a ripple effect on the SC because they are considered to occur in an exceptional event, with medium- to long-term recovery and an impact on the total costs (Dolgui et al., 2018; Ivanov et al., 2014; Llaguno et al., 2022). The studied SC has the capability to recover from disruption to its supply because many suppliers are available. Moreover, with a market risk like increased demand, the SC does not have the capacity to face this scenario, mainly because it produces at maximum capacity without being able to cover this increase.

Additionally, readers are referred to the Supplementary File provided to be opened with the Vensim DSS® simulation model as a published version in which all the input data can be accessed. These input data have been mainly provided by the copper mining company under study, and are also representative data from the copper mining sector, which have been collected from public reports.

Discussion

The importance of developing an SD model is that it allows the shortcomings of the MP model to be addressed (Azadeh & Vafa Arani, 2016). This means that it allows copper mining industry dynamics to be faced. Therefore, considering sustainability aspects in an interrelated way in the SD model can be a way to achieve sustainability at the corporate level (Sudarto et al., 2017). In relation to the economic dimension of sustainability, the SD model enables SC behaviour to be studied in different situations, such as disruptions in its normal operation. This has relevant impacts on the CLSC's economic performance. So the most relevant variable is inventory shortage as a breach of contract results, and not only as a loss of sales but also as a penalty cost.

The copper mining industry is one of the most polluting and resource-consuming industries, e.g. water. Although progress has been made in environmental issues, the impacts of this industry need to be further studied to contribute to sustainable development. Here the main contribution of the SD model is that, with a sensitivity analysis, it identifies how the variable with the strongest impact on environmental performance is the production of suppliers, which are the largest source of carbon emissions. Given this situation, it is relevant to study the impact on the environment and nearby communities generated by extractive and mineral processing industries (Mohammadi et al., 2022).

Regarding the social dimension of sustainability, the proposal considers the job opportunities generated by operating facilities (Motevalli-Taher et al., 2020). It does not only contemplate the generated direct jobs, but also the indirect jobs produced by the operation of certain facilities like processing centres, collection and repair centres, recycling centres, to name but a few. It is also important to consider the negative impact that CLSC generates by incorporating the accident rate on routes that can be caused by transport activities (Aloui et al., 2021).

The usefulness of the modelling and the results of this proposal are to serve as a basis for a balanced scorecard tool to measure and improve sustainability aspects related to economic (location, inventory and transport), environmental (CO₂ emissions) and social (route accidents) costs in other copper mining SCs, or even in different industrial SCs.

Conclusions and future research

This study presents a simulation optimisation model with very efficient computational times to support decision making in sustainable SC designs to recreate different scenarios and SC disruptions to help decision makers. First of all, the 3S-LIT MP model proposed by Becerra et al. (2023) is applied and compared to the current situation of the studied mining SC. Then the SD simulation model is presented by applying the same dataset as in the MP model, and is compared to, thus, recreate the optimisation model that validates the SD model.

After applying the MP and SD model in the copper mining industry SC, and comparing it to the current SC situation, the main findings are described as follows:

1. Lowering total economic costs by reducing stockouts and, thus, complying with customer contracts;

2. Despite the SC size growing by incorporating RL, CO₂ emissions slightly reduce;

3. Although social costs increase, the positive social impact is stronger because more direct and indirect jobs are created by new facilities;

4. Moreover by following the SD model validation, several what-if scenarios are simulated to evaluate the system's performance under such conditions, including different time periods, fixed and uncertain lead times, supply disruptions and increased demand. The sustainability results are assessed for this reason. The proposed scenarios consist of possible disruptions in the SC; for example, the disruption of the ore supply from the company's own mine, a sudden increase in demand and delayed production. From these scenarios, it is concluded that a disruption in the SC caused by supply, demand or production directly impacts economic costs because unmet demand increases due to stockouts. In addition, the SC's resilience to a possible ripple effect when supply from its own mine is interrupted can be observed, and the inability to recover from increased demand is identified, which results in bigger stockouts and, thus, higher economic costs;

5. The proposed SD model can be considered a based platform that allows the incorporation of new parameters and sustainability variables for future research purposes.

According to these findings, general managerial implications are oriented:

1. To use a combination of optimisation and simulation models as a powerful tool to make strategic LIT decisions. Thus starting with an optimal solution during a static time period allows different scenarios to be recreated and also SC disruption during a dynamic time period when the optimal parameter values cannot be adapted to our current SC constraints in a reasonable computational time;

2. To focus on employing these simulation models as a basis prototype to develop digital twins in similar SCs;

3. To integrate the SD model into the company's information system can be promoted. It is also possible to use the tool to carry out hypothetical analyses depending on the company's needs;

4. For management purposes, and once the optimisation model is built, to assess different industrial scenarios, it is easier and faster to create by allowing for management uses such as the balanced scorecard;

5. Specifically for the SC in the copper mining industry, the main managerial implication of the model is its ability to identify the SC's inability to withstand market disruption, but the ability to handle it with supply disruption by maintaining multiple suppliers. It also allows the identification of suppliers as the main source of carbon emissions, which can lead to improvements in the technology used for ore extraction processes by opting for machinery with a milder environmental impact, such as hydrogen-powered trucks (Ahluwalia et al., 2022).

The main limitation identified in this research work is oriented to the use of a specific case study to validate the model. Hence new real-world applications would be desirable to provide more generalisable conclusions. Additionally, as a representative domain or part of the studied SC has been considered, an extension of the model to the entire SC, or even interSCs modelling, is a challenge to tackle.

Future research could consider applying the MP model and the simulation model in a different industry, e.g. iron ore mining, lithium, the agri-food industry. Additionally, different inventory policies could be considered, e.g. instead of a push inventory system strategy. For example, a pull strategy could be considered in which demand drives production. Furthermore, due to space constraints, stochastic parameters like uncertain demand or waste collection rates are no considered in the different simulated scenarios, but could be taken into account by researchers and practitioners. Finally, forthcoming works will focus on incorporating new decisions, such as production and routing, and automated and autonomous decisions, by taking a hybrid simulation approach that incorporates advanced multi-agent SD models (Barbosa et al., 2023).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The authors confirm that the data supporting the findings of this study are available in the article and in the provided Supplementary Files.

Additional information

Funding

The research leading to these results received funding from the project ‘Industrial Production and Logistics Optimization in Industry 4.0’ (i4OPT) (Ref. PROMETEO/2021/065) granted by the Valencian Regional Government; and grant PDC2022-133957-I00 funded by MCIN/AEI /10.13039/501100011033 and by European Union Next GenerationEU/PRTR; and from the Horizon Europe Framework Programme and EC-funded project DiMAT, EU under grant agreement No 101091496. It was also funded by the National Agency for Research and Development (ANID)/Scholarship Programme/Doctorado Becas en el Extranjero/2020 – 72210174.