Robust Optimization Model for Sustainable Supply Chain Design Integrating LCA

Abstract

Supply chain management is the basis for the operations in an organization. The development of realistic supply chain designs that work effectively in the presence of disturbances in a stochastic environment and incorporate sustainability factors, is a complex challenge being investigated in recent years. However, the inclusion of a methodological structured framework to evaluate environmental impacts constitutes a knowledge gap in the literature on supply chain design. This study developed a model for sustainable supply chain design, integrating Life Cycle Assessment and based on a robust optimization approach. The study follows a 4-stage methodology beginning with data collection and the execution of a Life Cycle Assessment. Then, the deterministic modeling is proposed, concluding with a robust model. A bi-objective model is proposed to maximize utility and minimize environmental impact based on demand scenarios. The model was validated with real data from a medium-sized enterprise that produces antibacterial gel, generating as a result, different configuration alternatives for the supply chain to transport the products and raw materials between its elements. The conclusions of this work highlight the importance of including sustainability factors during supply chain design, the consequences and costs of its inclusion, as well as the priority actions that promote sustainable designs.

1. Introduction

The design of a supply chain (SC) is an area of decision-making that considers parameters such as planning, demand, time and costs [1]. These designs have traditionally and mostly been based on economic factors [2]; however, increasing globalization and decreasing profit margins, among other causes, prompted modern SC to focus on providing quality products to consumers consistently and quickly. This, in turn, drove a research trend that considers non-cost metrics in network design [3]. Currently, several studies focus on developing realistic SC designs that can operate effectively in the presence of external and/or internal disruptions by incorporating deterministic and, stochastic variables [4].

On the other hand, nowadays, factors such as the decrease in natural resources, the search for competitive advantages, the pressure to reduce environmental impacts of products and processes, and global laws and agreements have generated a greater interest in sustainable development. According to the 1987 Brundtland report, sustainable development “meets the needs of the present without compromising the ability of future generations to meet their own needs” [5]. Sustainability is a broad concept, encompassing environmental, economic and social aspects [6], requiring methods and tools to quantify and compare the environmental impacts of supplying products and services.

For achieving sustainable development, industries need to design, plan and operate their SC considering a sustainability path [7], which needs to be considered as a business opportunity rather than a constraint [8]. A widely used tool is the Life Cycle Assessment (LCA), structured as a methodological framework to identify, quantify, interpret and evaluate the environmental impact of a product, process, or service in a systematic manner [9].

Quantitative mathematical programming models constitute a convenient tool to deal with strategic decisions of sustainable supply chain (SSC) design and the deterministic and stochastic variables required [10]. Deterministic models are used when the parameters are known with certainty. If parameters are unknown, stochastic models are required, where uncertainty can be encompassed in a robust optimization approach [11]. Robust optimization represents uncertainty by setting up different scenarios aiming to find a robust solution that ensures that all specified scenarios are “close” to the optimum in response to changes in the input data [12].

In this context, the objective of this study is to develop a model for SSC design, integrating LCA and based on a robust optimization approach. For the validation of the model, the antibacterial gel product has been selected and data from a Mexican company is used. The model will serve to support decision-making and recommend SC configuration alternatives to maximize economic performance and minimize the organizations’ environmental impacts.

A large increase in the global production of antibacterial gel, and its current relevance, make the product a case study of academic and industrial interest. The production of antibacterial gel in Mexico has grown exponentially since 2020 due to the COVID-19 pandemic [13]. Therefore, determining the environmental impact throughout the SC of this product will allow companies to establish immediate strategies to reduce these impacts, generate competitive advantages, improve their sustainability and, in general, contribute to goal 12 of the Sustainable Development Goals (responsible production and consumption).

For a better understanding, the document is structured as follows: Section 2 summarizes the state of the art of the study. Section 3 explains the methodology and some assumptions for the development of the model. Section 4 describes the mathematical formulation of the model. Section 5 shows the main results. Finally, Section 6 summarizes the conclusions.

2. Literature Review

Realistic SC designs incorporate critical parameters modeled as random variables with known probability distributions. One classic approach to deal with the uncertainty of these variables is the concept of scenario-based stochastic robust optimization proposed by Mulvey et al. [14], which aims to generate a series of solutions that are progressively less sensitive to the realization of the model data from a scenario set.

In recent years, many authors applied robust optimization to firm-specific problems, such as inventory management, manufacturing, production planning, facility location and network design problems. For example, Pan and Nagi [15] formulated a robust optimization model for an agile manufacturing SC under uncertain demand. Likewise, Baghalian et al. [16] presented a stochastic formulation for designing a multi-product SC network based on a robust optimization concept. Mahdi et al. [17] propose a robust multi-objective framework to improve the robustness of the water supply system in accordance with the uncertainty of rainfall and water demand. Li et al. [18] develop a robust optimization approach for prismatic lithium-ion cells, aiming to limit thermal performance uncertainties to the optimal manufacturing cost.

In addition to optimizing SC, companies develop strategies to allow adapting its dynamics and controlling disruptions. In the study of Ivanov and Sokolov [19], a control approach is presented to model SCs as multi-structural dynamic systems. Sawik [20] developed a stochastic programming model for integrated supplier selection, order quantity allocation, and customer order scheduling in the presence of SC disruption risks. Ali et al. [21] model SC disruption analysis with insufficient data in order to identify supplier disruption risks that can severely impact SC performance. Ghavamifar et al. [22] redesign a competitive SC network under operational risks and disruptions, using a multi-objective programming approach.

Regarding sustainability in SC, the last decade presented a remarkable increase in research [23], including topics such as supplier selection [24], green network designs [25] and low carbon production [26]. Considering optimization models for the design of an SSC, Mota et al. [27] present a mixed-integer linear programming model with multiple objectives that integrates LCA with economic and social decisions. Wang et al. [28] design the network of a green SC by simultaneously considering carbon emissions and inter-competitor pricing. Jaber et al. [29] studied the problem of joint emission reduction between manufacturers and retailers in a two-stage SC. Ghelichi et al. [30] use a stochastic programming approach towards the optimal design and planning of an integrated green biodiesel SC network. Mele et al. [31] employ a multi-objective model for SSC fuel, carrying out a parallel LCA study.

When coupled with the robust optimization approach and SSC network design, the literature review yields some novel results. For example, Yousefi-Babadi et al. [32] use a robust optimization approach in the redesign of a three-tier wheat, flour and bread SC to meet sustainable development considerations. Geon Kim et al. [33] develop a robust optimization model for a closed-loop SC with uncertain demand and uncertain carbon tax rates. Krishnan et al. [34] present a robust integrated multi-objective optimization model to design a food SC network considering three dimensions (economic, social and environmental).

It can be evidenced then that there is extensive literature on SC management and SSC design, highlighting the use of mathematical, linear, integer, mixed, stochastic, fuzzy and robust programming. In order to understand how quantitative models have been used in the design of SSCs and what are the main sustainability pillars addressed within the decision process, the systematic literature review (SLR) performed by Flores et al. [35] is used as a base reference.

The main SLR results show that deterministic studies prevail over stochastic studies. Regarding the solution approach, mathematical programming models with uncertainty mostly use stochastic programming, and to a lesser extent robust optimization. Concerning sustainability strategies, studies on carbon and greenhouse gas emissions prevail over LCA, which are based on deterministic data and are developed as parallel studies to the SC design.

Based on the SLR [35], current gaps such as the lack of full inclusion of LCA in the design of a SSC, were identified. This means that indicators resulting from LCA are not included in the mathematical formulation of robust optimization, which would allow estimating and evaluating the environmental impacts attributable to the entire life cycle of a product and its interaction and influence on the network design.

The relevance of LCA is recognized by several authors, especially for the environmental dimension. LCA is acknowledged as a comprehensive technique to assess the environmental impacts of products and processes [36], which still needs more in-depth research to become a sound decision-making tool for each stage of a product life-cycle [37]. Until now, the two main sectors where LCA has been most frequently applied are construction [38] and renewable energies [39].

3. Methodology

In order to develop the model for SSC design, 4 stages are proposed, as shown in Figure 1. Stages 2 to 4 are detailed in Section 3.1 and Section 3.2.

Input Data refer to the necessary information that must be acquired to prepare both the LCA and the mathematical formulation of the model. This information corresponds to the characteristics of the raw materials, energy and waste flows present in the production process. Also, data on the characteristics of the SC, logistics management, production and demand are incorporated.

3.1. Life Cycle Assessment

With the information acquired in ‘stage 1, the LCA is prepared based on the ISO 14044 methodology, which consists of four steps. An additional step was included for integrating the resulting environmental impacts in a mathematical function:

(1) Goal and Scope: The objective of the proposed LCA is to quantify the environmental impacts generated by the antibacterial gel product throughout its SC (raw material extraction, production and distribution). A 120 mL bottle is defined as a functional unit. The final disposal of the product is not considered since no data was available. Likewise, a closed-loop SC is not considered since no recovery processes are implemented in the assessed company.

(2) Lifecycle Inventory: Information is collected on the raw materials used in the production of antibacterial gel (quantity and origin), energy used in the production stage and characteristics of the transportation used for product distribution.

(3) Lifecycle Impact Assessment: This step translates the emissions and resource extractions into environmental impact scores. Gabi software 10.5.1 is used, where the entire SC of the antibacterial gel and previously collected information is modeled. The modeling and the Recipe 2016 method allow for estimating 17 midpoint environmental impact categories and three final impact categories including (a) human health damage, (b) ecosystem diversity and (c) resource availability.

(4) Interpretation: performs an analysis of the midpoint environmental impact categories with the highest values in the extraction, production and transport stages. In the first two stages (extraction and production) the end-point environmental impacts are calculated following the methodology of Hauschild and Huijbregts [40]. The resulting end-point environmental impacts for the extraction and production stages according to the impact area are human health damage (IAFH) = 1.63 × 10⁻⁷, ecosystem diversity (IAFE) = 6.06 × 10⁻⁹ and resource availability (IAFR) = 4.60 × 10⁻¹.

(5) Integration: The coefficients IAFH, IAFE and IAFR are included in the objective function 2 of the mathematical model (Min EI) (Equations (6) and (8)). For the transport stage, a linear regression model (Equation (1)) is used to express the relationship between the quantity transported, the distance traveled and the environmental impact generated (Y = climate change category).

$$Y=-1.60529+0.00704\times \mathrm{K}\mathrm{g}+0.02967\times \mathrm{M}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{a}\mathrm{s}$$

(1)

The coefficients of Equation (1): b0 (IAPIN) = −1.60529, b1 (IAPKG) = 0.00704, and b2 (IAPMI) = 0.02967, are also included in the objective function 2 of the mathematical model (Min EI) (Equations (7), (9) and (10)).

A complete description of the LCA can be found in [41].

3.2. Conceptual Framework for the Mathematical Formulation

The mathematical formulation of the model is divided into two parts, deterministic and robust.

Deterministic modeling: considers the data acquired in stage 1 and the results obtained from stage 2, and uses only deterministic data to optimize two objective functions: (1) Utility maximization (Max UTI), and (2) minimization of environmental impact (Min EI). This stage also includes the definition of capacity restrictions and flow balance between the network arcs.

Robust modeling: The demand uncertainty variable is added to the deterministic formulation, considering the approach described by Jabbarzadeh et al. [42] and Mulvey et al. [14]. A set of uncertain scenarios with their respective probabilities is built and the model of stage 3 is transformed into a robust one, which will generate global results and results per scenario.

In addition, the mathematical formulation is based on a four-layered SC encompassing multiple raw material suppliers, production plants, distribution centers and customers (see Figure 2).

Finally, the following assumptions were made during model formulation.

• A single product is produced and distributed throughout the network;
• Only one transport mean is considered in the entire network. Only a one-way trip is considered (the return trip of trucks is omitted).
• Products use diversified raw materials.
• No direct shipments can be made from the plant to the customer.
• Supplier locations and customer zones are fixed.
• Production plant and distribution center locations are fixed.
• Suppliers have a maximum availability of raw materials.
• Production plants and distribution centers have a maximum capacity.
• The cost of raw materials depends on the supplier.
• The environmental impact generated in each production plant depends on the type of raw material used.

4. Mathematical Formulation of the Model

The indexes, parameters and variables used during mathematical formulation (Stage 3 and 4 of the Methodology) are illustrated below.

Indices:

s	supplier locations, s = 1, 2 …. S
m	raw materials, m = 1, 2 …. M
i	manufacturing plant locations, i = 1, 2 …. I
j	distribution center locations, j = 1, 2 …. J
n	customer locations, n = 1, 2 …. N

Parameters:

$k_s$	fixed cost of selecting a supplier s
$r_{ms}$	maximum capacity of raw material m provided by supplier s
${prp}_{ms}$	price per unit of raw material m provided by supplier s
$t_m$	rate of conversion of raw material m to finished product
${ctp}_m$	transportation cost per mile of one unit of raw material m
${\rho p}_{si}$	distance from supplier s to plant i
$IAPIN$	coefficient b0 of the EI generated from transporting the raw material from the supplier to the plant (LCA)
$IAPKG$	coefficient b1 of the EI generated from transporting the raw material from the supplier to the plant (LCA)
$IAPMI$	coefficient b2 of the EI generated from transporting the raw material from the supplier to the plant (LCA)
$f_i$	fixed cost of producing in plant i
$l_i$	maximum production capacity of plant i
$p_i$	production cost per unit of product in plant i
$ctw$	transportation cost per mile of one unit of finished product from the plant to the distribution center
$gelkg$	weight in kilograms of one unit of finished product
${\rho w}_{ij}$	distance from plant i to distribution center j
${IAFH}_i$	EI on human health of raw material extraction and production processes of plant i (LCA)
${IAFE}_i$	EI on the ecosystem of raw material extraction and production processes of plant i (LCA)
${IAFR}_i$	EI on resource availability of raw material extraction and production processes of plant i (LCA)
$IAPCIN$	coefficient b0 of the EI generated by transporting the finished product from the plant to the distribution center (LCA)
$IAPCKG$	coefficient b1 of the EI generated by transporting the finished product from the plant to the distribution center (LCA)
$IAPCMI$	coefficient b2 of the EI generated by transporting the finished product from the plant to the distribution center (LCA)
$g_j$	fixed cost of using a distribution center j
${\phi }_j$	storage capacity of distribution center j
$h_j$	maintenance cost per product unit in distribution center j
$ctc$	transportation cost per mile of one unit of product from the distribution center to the customer
${\rho c}_{jn}$	distance from distribution center j to customer n
$IACCIN$	coefficient b0 of the EI generated from transporting the product from the distribution center to the customer (LCA)
$IACCKG$	coefficient b1 of the EI generated from transporting the product from the distribution center to the customer (LCA)
$IACCMI$	coefficient b2 of the EI generated from transporting the product from the distribution center to the customer (LCA)
$d_n$	demand of customer n
V	selling price per product unit

Decision variables:

$X_s=\left\{\begin{array}{l}1, \mathrm{if} \mathrm{supplier} s \mathrm{is} \mathrm{selected}\\ 0, \mathrm{otherwise}\end{array}\right.$

$Z_i=\left\{\begin{array}{l}1, \mathrm{if} \mathrm{a} \mathrm{plant} \mathrm{is} \mathrm{selected} \mathrm{at} \mathrm{location} i\\ 0, \mathrm{otherwise}\end{array}\right.$

$Y_j=\left\{\begin{array}{l}1, \mathit{if} a distribution center is selected at location j\\ 0, \mathrm{otherwise}\end{array}\right.$

${QS}_{msi}$	quantity of raw material m provided by supplier s and shipped to plant i
${QPP}_{si}$	total quantity of raw materials provided by supplier s and shipped to plant i
${QP}_{ij}$	quantity of products produced at plant i and shipped to distribution center j
${QW}_{jn}$	quantity of products shipped from distribution center j to customer n

The three sets of variables (X, Y, Z) are binary, while the rest of variables take positive real values only.

4.1. Deterministic Mathematical Model

Objective functions of the deterministic model:

(1)

Max UTI = sales revenue—total cost

$$Max UTI={\sum }_{j=1}^J{\sum }_{n=1}^N{QW}_{jn}·V-(CP+CM+CT)$$

(2)

The first OF refers to profit maximization, composed of sales revenue minus a total cost. The total cost is divided into three parts. The first (Equation (3)) refers to the cost of acquiring raw materials (CP), the second part (Equation (4)) covers production and storage costs (CM), and the third (Equation (5)) refers to transportation costs throughout the network (CT).

$$CP={\sum }_{m=1}^M{\sum }_{s=1}^S{\sum }_{i=1}^I{prp}_{ms}·{QS}_{msi}+{\sum }_{s=1}^S{k}_s·X_s$$

(3)

$$CM={\sum }_{i=1}^I{f}_i·Z_i+{\sum }_{i=1}^I{\sum }_{j=1}^J{p}_i{·QP}_{ij}+{\sum }_{j=1}^J{g}_j· Y_j+{\sum }_{j=1}^J{\sum }_{n=1}^N{h}_j{·QW}_{jn}$$

(4)

$$\begin{gathered} CT={\sum }_{m=1}^M{\sum }_{s=1}^S{\sum }_{i=1}^I{\rho p}_{si}·{ctp}_m·{QS}_{msi}+{\sum }_{i=1}^I{\sum }_{j=1}^J{\rho w}_{ij}· ctw·{QP}_{ij}+ \\ {\sum }_{j=1}^J{\sum }_{n=1}^N{\rho c}_{jn}· ctc·{QW}_{jn} \end{gathered}$$

(5)

(2)

Min EI = sum of the environmental impact generated in each stage of the SC.

$$Min EI=IAP+IAF+IAPC+IACC$$

(6)

The second OF refers to the minimization of the environmental impact generated throughout the chain, which is composed of four sections. The first (Equation (7)) calculates the environmental impact with respect to the climate change caused during the transport of raw material from the supplier to the plant (IAP). The second part of this OF (Equation (8)) refers to the environmental impact of the processes of extraction of raw material and production of the product (IAF), considering all the impact categories of the Recipe LCA method. The third part (Equation (9)) represents the environmental impact of the climate change category generated during transporting of finished product from the plant to the distribution center (IAPC). Finally, the fourth section (IACC) encompasses the climate change environmental impact from transporting the product from the distribution center to the customer (Equation (10)).

$$IAP={\sum }_{s=1}^S{\sum }_{i=1}^I\left[ \left(IAPKG·{QPP}_{si}\right)+({\rho p}_{si}· IAPMI)+IAPIN \right]$$

(7)

$$IAF={\sum }_{i=1}^I{\sum }_{j=1}^J{QP}_{ij}·{IAFH}_i + {\sum }_{i=1}^I{\sum }_{j=1}^J{QP}_{ij}·{IAFE}_i+ {\sum }_{i=1}^I{\sum }_{j=1}^J{QP}_{ij}·{IAFR}_i$$

(8)

$$IAPC={\sum }_{i=1}^I{\sum }_{j=1}^J\left[ \left(IAPCKG·gelkg· {QP}_{ij}\right)+({\rho w}_{ij}· IAPCMI)+IAPCIN \right]$$

(9)

$$IACC={\sum }_{j=1}^J{\sum }_{n=1}^N\left[ \left(IACCKG·gelkg· {QW}_{jn}\right)+({\rho c}_{jn}· IACCMI)+IACCIN \right]$$

(10)

Restrictions of the deterministic model:

The first three constraints (Equations (11)–(13)) represent the capacity limitations of suppliers, plants and distribution centers. Equation (11) ensures that the amount of raw material m provided by each supplier s does not exceed its capacity. Equation (12) allows the production of each plant to be less than the designed capacity. Equation (13) limits the shipment of products to the distribution centers, which must be less than their capacity.

$${\sum }_{i=1}^I{QS}_{msi}\le r_{ms}· X_s \forall m\in M, s\in S$$

(11)

$${\sum }_{j=1}^J{QP}_{ij}\le l_i· Z_i \forall i\in I$$

(12)

$${\sum }_{n=1}^N{QW}_{jn}\le {\phi }_j· Y_j \forall j\in J$$

(13)

The equations for material flow balance are represented by 5 constraints (Equations (14)–(18)). Equation (14) ensures that the amount of raw material needed to manufacture the amount of product ordered is provided. Equation (15) calculates the total amount of raw material shipped by each supplier to the plant. Equation (16) formulates the balance of incoming and outgoing products at each storage facility. Equation (17) establishes that product output from warehouse j satisfies the demand of customer n. Equation (18) guarantees that the total product output from all plants must satisfy the total demand.

$${\sum }_{s=1}^S{QS}_{msi}={\sum }_{j=1}^J{QP}_{ij}· t_m \forall m\in M, i\in I$$

(14)

$${\sum }_{m=1}^M{QS}_{msi}={QPP}_{si} \forall s\in S, i\in I$$

(15)

$${\sum }_{i=1}^I{QP}_{ij}={\sum }_{n=1}^N{QW}_{jn} \forall j\in J$$

(16)

$${\sum }_{j=1}^J{QW}_{jn}=d_n \forall n\in N$$

(17)

$${\sum }_{i=1}^I{\sum }_{j=1}^J{QP}_{ij} ={\sum }_{n=1}^N{d}_n \forall n\in N$$

(18)

The formulation was developed to allow the model to make decisions regarding five general categories: (1) quantity of raw material to be used, (2) selection of suppliers, (3) selection of production plants, (4) selection of distribution centers and (5) calculation of the flow for each arc. The model seeks to maximize profitability, minimize environmental impact and satisfy 100% of the demand.

4.2. Robust Mathematical Model

In robust optimization, there are two types of constraints: structural and control, The latter differ from the former in that they are influenced by noisy data. Similarly, two types of variables are established: design and control. The value of the design variables can be determined prior to the realization of the scenario, while the value of the control variables can be changed according to various simulations of the uncertain parameters [10].

The scenario analysis method has been widely adopted to deal with SC uncertainties, and the wait-and-see approach is the most discussed one for assuming the realization of different scenarios and making decisions in each scenario [43]. Mulvey et al. [14] were the pioneers in introducing a robust linearization method when the scenario is adopted to describe uncertainty. Over time, several authors have improved Mulvey’s proposal, where Jabbarzadeh et al. [42] stands out and serves as a reference for the proposed model.

An example of a linear programming model to separate the design and control variables is explained as follows:

$$\begin{gathered} MIN { c}^{T}x+ d^{T}y, \\ Ax=b, \\ Bx+Cy=e, \\ x, y\ge 0, \end{gathered}$$

(19)

where x denotes the vector of design variables and y denotes the vector of control variables; B, C and e represent the uncertain parameters; Ax = b denotes the constraints whose coefficients are free of uncertainties; Bx + Cy = e denotes the constraints whose coefficients are subject to uncertainties.

The robustness approach includes two aspects: the robustness of the solution and the robustness of the model. In this case, solution robustness means that the solution remains “close” to the optimum for any scenario, while model robustness implies that the solution is “almost” feasible [44]. Therefore, solution robustness can be measured by evaluating how close a solution is to the optimal value in each scenario, while model robustness can be measured by evaluating constraint violations.

We define Ω = {1, 2, …, ξ} as a set of uncertain scenarios to describe the uncertainty of certain variables and assume that each scenario occurs with probability $p_{\xi }$; where ${\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}\mathrm{p}\mathsf{\xi }=1$. A control variable ${\mathrm{y}}_{\xi }$ is introduced for each scenario ξ ϵ Ω, this means that the objective function of a model (Equation (19)) would become a random variable taking the value ${\mathsf{\psi }}_{\xi }=c^{T}x+d^T{\mathrm{y}}_{\xi }$ with probability $p_{\xi }$ under scenario ξ ϵ Ω. Similarly, an error vector ${\mathsf{\delta }}_{\xi }$ is introduced to measure the infeasibility allowed in the control constraints under scenario ξ ϵ Ω. In addition, $B_{\xi }$, $C_{\xi }$, $e_{\xi }$ represent random variables in scenario ξ. With these changes, the robust optimization model for the mathematical programming problem (Equation (19)) can be formulated as follows (Equation (20)).

$$\begin{gathered} MIN \sigma (x,y_1, y_2, \dots , y_{\xi })+ \omega p{(\delta }_1, {\delta }_2, \dots , {\delta }_{\xi }), \\ Ax=b, \\ {B}_{\xi }x+C_{\xi }{\mathrm{y}}_{\xi }+{\mathsf{\delta }}_{\xi } = e_{\xi }, \forall \xi \in \Omega , \\ x \ge 0, {\mathrm{y}}_{\xi }\ge 0, {\mathsf{\delta }}_{\xi }\ge 0, \forall \xi \in \Omega \end{gathered}$$

(20)

The first term $\sigma (x,y_1, y_2,\dots , y_{\xi })$ is a measure of solution robustness, which evaluates the closeness of a solution to optimality for any realization of the scenario ξ ϵ Ω. The second term $p{(\delta }_1, {\delta }_2,\dots ,{\delta }_{\xi })$ is a measure of model robustness, which penalizes violations of control constraints in some scenarios. The model (Equation (20)) adopts a multi-criteria objective form, where the weighting penalty Ω (called risk aversion weight) is used to express the trade-off between the robustness of the solution and the robustness of the model.

The choice of suitable functions for $\sigma (x,y_1, y_2,\dots , y_{\xi })$ and for $p{(\delta }_1, {\delta }_2,\dots {\delta }_{\xi })$, is not straightforward, so typically the mean values $\sum _{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }} {\mathsf{\psi }}_{\xi }$ and $\sum _{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }} {\mathsf{\delta }}_{\xi }$ are chosen in each case. Mulvey et al. [14] suggest that a better selection for $\sigma (x,y_1, y_2,\dots , y_{\xi })$ would be the mean plus a constant (λ) multiplied by the variance, as shown in Equation (21).

$$\sigma (x,y_1, y_2, \dots , y_{\xi })={\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }} {\mathsf{\psi }}_{\xi }+ \lambda {\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }} {\left({\mathsf{\psi }}_{\xi }-{\sum }_{\mathsf{\xi }\mathrm{\prime }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }\mathrm{\prime }} {\mathsf{\psi }}_{\xi \prime }\right)}^2$$

(21)

The objective value is less sensitive to the change of input data in all scenarios as λ (the variability weight) increases [14]. Due to the complexity represented by the squared term, Yu & Li [45] proposed a technique to transform Equation (20) into a linear programming model. Using the Yu and Li approach, the robust optimization model (Equation (20)) can be transformed into the following linear model (Equation (22)).

$$\begin{gathered} MIN {\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }} {\mathsf{\psi }}_{\xi }+ \lambda {\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }} \left[\left({\mathsf{\psi }}_{\xi }-{\sum }_{\mathsf{\xi }\mathrm{\prime }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }\mathrm{\prime }} {\mathsf{\psi }}_{\xi \prime }\right)+{2\theta }_{\xi }\right]+\omega {\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }} {\mathsf{\delta }}_{\xi }, \\ {\mathsf{\psi }}_{\xi }-{\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }} {\mathsf{\psi }}_{\xi }+{\theta }_{\xi } \ge 0, \forall \xi \in \mathsf{\Omega }, \\ {\theta }_{\xi } \ge 0, \forall \xi \in \mathsf{\Omega } \\ Ax=b, \\ {B}_{\xi }x+C_{\xi }{\mathrm{y}}_{\xi }+{\mathsf{\delta }}_{\xi } = e_{\xi }, \forall \xi \in \mathsf{\Omega }, \\ x \ge 0, {\mathrm{y}}_{\xi }\ge 0, {\mathsf{\delta }}_{\xi }\ge 0, {\theta }_{\xi } \ge 0 \forall \xi \in \mathsf{\Omega } \end{gathered}$$

(22)

Following these steps, the robust formulation of the model proposed in this study is as follows.

Objective functions of the robust model:

$$\begin{gathered} Max UTI={\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }} . \left({\sum }_{j=1}^J{\sum }_{n=1}^N{QW}_{jn}^{\xi }·V -{CP}^{\xi } - {CM}^{\xi }- {CT}^{\xi } \right)+\lambda {\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }}\left[{\sum }_{j=1}^J{\sum }_{n=1}^N{QW}_{jn}^{\xi }V-{CP}^{\xi }-{CM}^{\xi }- \\ {CT}^{\xi }-{\sum }_{\mathsf{\xi }\mathrm{\prime }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }\mathrm{\prime }}\left({\sum }_{j=1}^J{\sum }_{n=1}^N{QW}_{jn}^{\xi \prime }V-{CP}^{\xi \prime }-{CM}^{\xi \prime }-{CT}^{\xi \prime } \right)+{2\theta }_{\xi } \right]+\omega {\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }}{\mathsf{\delta }}_{\xi } \end{gathered}$$

(23)

$$\begin{gathered} Min EI={\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }} . \left({IAP}^{\xi } + {IAF}^{\xi }+ {IAPC}^{\xi }+ {IACC}^{\xi }\right)+\lambda {\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }}\left[{IAP}^{\xi }+{IAF}^{\xi }+{IAPC}^{\xi }+{IACC}^{\xi }- \\ {\sum }_{\mathsf{\xi }\mathrm{\prime }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }\mathrm{\prime }}\left( {IAP}^{\xi \prime }+{IAF}^{\xi \prime }+{IAPC}^{\xi \prime }+{IACC}^{\xi \prime } \right)+{2\tau }_{\xi } \right]+\omega {\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }}{\mathsf{\delta }}_{\xi } \end{gathered}$$

(24)

Restrictions of the robust model:

$${\sum }_{j=1}^J{\sum }_{n=1}^N{QW}_{jn}^{\xi }V-{CP}^{\xi }-{CM}^{\xi }-{CT}^{\xi }-{\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }}\left({\sum }_{j=1}^J{\sum }_{n=1}^N{QW}_{jn}^{\xi }V-{CP}^{\xi }-{CM}^{\xi }-{CT}^{\xi } \right)+{\theta }_{\xi }>0 \forall \xi ϵ \Omega$$

(25)

$${IAP}^{\xi }+{IAF}^{\xi }+{IAPC}^{\xi }+{IACC}^{\xi }-{\sum }_{\mathsf{\xi }\in \mathsf{\Omega }}{p}_{\mathsf{\xi }}\left( {IAP}^{\xi }+{IAF}^{\xi }+{IAPC}^{\xi }+{IACC}^{\xi } \right)+{\tau }_{\xi }>0 \forall \xi ϵ \Omega$$

(26)

$${\sum }_{i=1}^I{QS}_{msi}^{\xi }\le r_{ms}· X_s^{\xi } \forall m\in M, s\in S, \xi ϵ \mathsf{\Omega }$$

(27)

$${\sum }_{j=1}^J{QP}_{ij}^{\xi }\le l_i·Z_i^{\xi } \forall i\in I, \xi ϵ \mathsf{\Omega }$$

(28)

$${\sum }_{n=1}^N{QW}_{jn}^{\xi }\le {\phi }_j· Y_j^{\xi } \forall j\in J, \xi ϵ \mathsf{\Omega }$$

(29)

$${\sum }_{s=1}^S{QS}_{msi}^{\xi }={\sum }_{j=1}^J{QP}_{ij}^{\xi }· t_m \forall m\in M, i\in I, \xi ϵ \mathsf{\Omega }$$

(30)

$${\sum }_{m=1}^M{QS}_{msi}^{\xi }={QPP}_{si}^{\xi } \forall s\in S, i\in I, \xi ϵ \mathsf{\Omega }$$

(31)

$${\sum }_{i=1}^I{QP}_{ij}^{\xi }={\sum }_{n=1}^N{QW}_{jn}^{\xi } \forall j\in J, \xi ϵ \mathsf{\Omega }$$

(32)

$${\sum }_{j=1}^J{QW}_{jn}^{\xi }+{\delta }_{\xi }=d_n^{\xi } \forall n\in N, \xi ϵ \mathsf{\Omega }$$

(33)

$${\sum }_{i=1}^I{\sum }_{j=1}^J{QP}_{ij}^{\xi }+{\delta }_{\xi } ={\sum }_{n=1}^N{d}_n^{\xi } \forall n\in N, \xi ϵ\mathsf{\Omega }$$

(34)

5. Validation Results and Discussion

5.1. Description of the Case Study

The case study company is a medium-sized enterprise with 100 workers, located in Guadalajara, Mexico. The company manufactures and sells personal care products and one of its best-selling items is the 120 mL antibacterial gel. Its production process involves three mixing machines and a bottling plant, while a fleet of trucks with a capacity of 22.6 tons is used for transportation.

The antibacterial gel is made up of seven raw materials (carbomer, distilled water, ethanol, triethanolamine, glycerin, tocopherol and perfume) while its packaging consists of three parts (PET bottle, cap and label). The company currently has eight approved suppliers, which have availability of one or more raw materials. Demand is considered a control variable and three demand scenarios are generated (pessimistic, normal and optimistic).

The SC case study of the antibacterial gel used to validate the model consists of several suppliers, a single production plant, two distribution centers and five wholesale customers (see Figure 3).

5.2. Results and Discussion

For the generation of results, the model was programmed in the GAMS 40 software and, in order to obtain results of the individual functions, the MIP function of the Cplex solver is applied, which uses the Branch and Bound method. Likewise, for the bi-objective problem, Cplex is applied in the same manner through the combination of multiple objectives and the ε constraint generation method.

For the execution of the algorithm, three runs of the robust model are carried out. For Run 1, only the Max UTI OF is considered. Run 2 considers only the Min EI OF. Finally, Run 3 takes into account both OF simultaneously.

Table 1. Comparison of the results obtained with the robust model development, per run and per demand scenario.

Category	OF Max UTI (Run 1)	OF Min EI (Run 2)	Max and Min Bi-Objective (Run 3)
Robustness solution (UTI) (MXN)	$1,838,225	$1,761,742	$1,830,066
Pessimistic scenario	$1,686,786	$1,458,043	$1,546,165
Normal scenario	$1,779,928	$1,657,086	$1,705,841
Optimistic scenario	$1,873,071	$1,850,010	$1,865,517
Robustness solution (EI) (Eco points)	197,668	197,662	197,666
Pessimistic scenario	168,792	168,782	168,790
Normal scenario	187,542	187,535	187,539
Optimistic scenario	206,291	206,289	206,292
Units shipped to customers
Pessimistic scenario	364,950	364,950	364,950
Normal scenario	405,500	405,500	405,500
Optimistic scenario	446,050	446,050	446,050

shows the differences for each execution of the algorithm, detailing the overall results as well as the results for each demand scenario.

The data in Table 1 show that the model reflects robustness in the solution since its values remain “close” to the optimum for any demand scenario. Furthermore, in terms of model robustness, the solution is feasible for all demand scenarios. As for the values obtained for both utility and environmental impact, based on the results in Table 1, it can be generalized that the values of Run 3 are in the middle of Run 1 and Run 2. This indicates that the model tries to find a balance between these two opposing objectives.

As for the general environmental impact, there is a slight change between the values of the three columns, due to the fact that there is only one production plant with only one method to manufacture the product. This means that the only way to reduce this environmental impact is by selecting the shortest routes between the SC arcs. This slight decrease in environmental impact generates an average decrease of 7.2% in profit. Finally, for all demand scenarios, 100% of the demand is met for all the algorithm executions, implying that all the products requested by customers are delivered. Regarding product flows through each arc,

Table 2. Flow between distribution center j and customer n—Robust model—Pessimistic demand.

Category	OF Max UTI (Run 1)	OF Min EI (Run 2)	Max and Min Bi-Objective (Run 3)
Pessimistic Scenario (Units Shipped)
D. Center_1—Customer_1	51,750	51,750
D. Center_1—Customer_2	7650	76,500
D. Center_1—Customer_3			134,300
D. Center_1—Customer_4	40,500	40,500	40,500
D. Center_1—Customer_5		25,200	25,200
D. Center_2—Customer_1			51,750
D. Center_2—Customer_2			76,500
D. Center_2—Customer_3	171,000	171,000	36,700
D. Center_2—Customer_4
D. Center_2—Customer_5	25,200
Utilization D. Center_1	84.38%	96.98%	100.00%
Utilization D. Center_2	78.48%	68.40%	65.98%

,

Table 3. Flow between distribution center j and customer n—Robust model—Normal demand.

Category	OF Max UTI (Run 1)	OF Min EI (Run 2)	Max and Min Bi-Objective (Run 3)
Normal Scenario (Units Shipped)
D. Center_1—Customer_1	57,500	57,500
D. Center_1—Customer_2		69,500
D. Center_1—Customer_3	114,500		190,000
D. Center_1—Customer_4		45,000
D. Center_1—Customer_5	28,000	28,000
D. Center_2—Customer_1			57,500
D. Center_2—Customer_2	85,000	15,500	85,000
D. Center_2—Customer_3	75,500	190,000
D. Center_2—Customer_4	45,000		45,000
D. Center_2—Customer_5			28,000
Utilization D. Center_1	100.00%	100.00%	95.00%
Utilization D. Center_2	82.20%	82.20%	86.20%

and

Table 4. Flow between distribution center j and customer n—Robust model—Optimistic demand.

Category	OF Max UTI (Run 1)	OF Min EI (Run 2)	Max and Min Bi-Objective (Run 3)
Optimistic Scenario (Units Shipped)
D. Center_1—Customer_1	57,000	63,250
D. Center_1—Customer_2	93,500	56,450	93,500
D. Center_1—Customer_3			57,000
D. Center_1—Customer_4	49,500	49,500	49,500
D. Center_1—Customer_5		30,800
D. Center_2—Customer_1	6250		63,250
D. Center_2—Customer_2		37,050
D. Center_2—Customer_3	209,000	209,000	152,000
D. Center_2—Customer_4
D. Center_2—Customer_5	30,800		30,800
Utilization D. Center_1	100.00%	100.00%	100.00%
Utilization D. Center_2	98.42%	98.42%	98.42%

show the quantity of products sent from the distribution centers to the customers and the variations that exist for each demand scenario and model run.

According to the data shown in Table 2, Table 3 and Table 4, in the three demand scenarios, during Run 1, the distribution centers are selected, giving preference to those with lower costs. In Run 2, the focus is on selecting the distribution centers closest to plants as well as to the customers, in order to reduce the environmental impact resulting from transportation. Finally, in the bi-objective run (Run 3), the aim is to balance costs and proximity. Another noteworthy aspect of the results obtained refers to the percentages of utilization of the distribution centers, where center 1 is the most used in most cases. Therefore, consideration should be given to expanding the capacity of this center as it would increase the value of the utility and reduce environmental impact.

The results obtained also reflect the importance of including an OF that minimizes environmental impact in SSC management, agreeing with similar studies that show that the application of tools to analyze LCA, carbon emissions, carbon footprint or greenhouse gases, allows industries to guide its SC towards a path of sustainability [8,10,27,46]. A priori, these efforts do not generate economic benefits and increase SC complexity but, it is the path to follow, considering the current government pressures, regulations, global standards and consumer demands [6].

According to the data obtained, the decisions that contribute the most to building an SSC design are:

• During the selection and qualification of suppliers, it is recommended to consider aspects such as location with respect to the production plants, the means of transportation and the shipment capacity of the raw material to be used;
• Investigate the use of alternative raw materials, especially raw materials that generate less environmental impact throughout the SC;
• Investigate the application of new technology, focused on reducing the environmental impact of the production stage;
• Be willing to sacrifice a percentage of the current profit for a higher profit in the near future.

The validation results of the model allow the case study to analyze different situations occurring in the SSC design of its antibacterial gel product considering demand scenarios. Therefore, the model becomes a decision support tool that provides organizations with different alternatives for configuring their SC to move their finished products or raw materials according to their economic or environmental interests.

6. Conclusions

An SSC has been designed through the execution of a robust optimization model that considers LCA indicators, which has allowed generating and analyzing the results of different demand scenarios. The design of a dynamic network is also evaluated since the model allows the selection of different suppliers, production plants and distribution centers. In the event of internal or external disruptions that force the partial closure of one of the SC elements, the model will select the second-best option, avoiding the paralyzation of the entire SC.

The validation of the model has been carried out with data from a particular company and one of its main products, the antibacterial gel. The model has been structured in a generic way, as it considers an unlimited number of suppliers, production plants, distribution centers and customers in its mathematical formulation, so it could be adapted to other conflicting industries conflicting industries such as wood, paper, food processing and textiles [34,47], which also require tools to promote the implementation of sustainable strategies in their SCs. To adapt the model, it is necessary to start by performing the LCA of the new product of interest, gathering information on its operational processes and SC and integrating them into the model proposed.

In Section 5, special attention was paid to the comparison between the classical models that consider only costs versus the developed model that considers LCA-based environmental impact indicators. Results showed that the average percentage of profit reduction is 7.2%; however, this value will depend on the managers, who should study and analyze how much they are willing to reduce their profit, what are the future benefits and what sustainability measures are feasible to implement in their respective SCs.

The decision of whether or not to implement these alternatives should consider in parallel the economic benefits generated by SSC such as market positioning, competitive advantage, the attraction of new customers, access to new markets, demand increase, acquisition of certifications, strengthening of the corporate image and customer loyalty [48]. In this manner, the initial reduction in profit is compensated [49].

In addition, the actions simulated in this work, which are oriented towards an SSC, are just examples of the wide variety of options available. Companies must take a proactive approach to reduce sources of waste or pollution and implement measures such as ecological design, green purchases, environmental collaboration of suppliers, saving resources, reducing harmful materials, and recycling or reusing products. If these initiatives are carried out jointly, they can have a significant relationship with environmental and economic performance [48,49]. These authors demonstrated that a SSC can generate both environmental and economic benefits, without sacrificing its efficiency.

However, there are also limitations of the study. The designed SC is not a closed-loop, since it does not allow identifying the environmental impact in the post-use phase of the product. Also, only the uncertainty caused by the demand is considered and other variables such as costs are omitted. These limitations show the direction to follow in future research. In addition, the model can be further enhanced by considering a closed-loop SC with a more stochastic environment.