An efficient Lagrangian-based heuristic to solve a multi-objective sustainable supply chain problem

TL;DR

This paper introduces a novel Lagrangian-based heuristic, $AugMathLagr$, for solving complex multi-objective sustainable supply chain problems, demonstrating competitive and outstanding performance through extensive testing and case studies.

Contribution

The paper presents a new Lagrangian matheuristic method, $AugMathLagr$, specifically designed for multi-objective sustainable supply chain optimization problems.

Findings

$AugMathLagr$ outperforms existing methods in artificial instance tests.

The method shows superior results in a real-world case study.

Computational experiments validate the efficiency and competitiveness of $AugMathLagr$.

Abstract

Sustainable Supply Chain (SSC) management aims at integrating economic, environmental and social goals to assist in the long-term planning of a company and its supply chains. There is no consensus in the literature as to whether social and environmental responsibilities are profit-compatible. However, the conflicting nature of these goals is explicit when considering specific assessment measures and, in this scenario, multi-objective optimization is a way to represent problems that simultaneously optimize the goals. This paper proposes a Lagrangian matheuristic method, called $AugMathLagr$, to solve a hard and relevant multi-objective problem found in the literature. $AugMathLagr$ was extensively tested using artificial instances defined by a generator presented in this paper. The results show a competitive performance of $AugMathLagr$ when compared with an exact multi-objective method limited by time and a matheuristic recently proposed in the literature and adapted here to address the studied problem. In addition, computational results on a case study are presented and analyzed, and demonstrate the outstanding performance of $AugMathLagr$.

Tables

Table 1. Table 1: Notations for the economic, environmental and social functions.

Description	Notation
Economic function
Max NPV	NPV
Min costs	C
Max profit	P
Global after tax profit	ATP
Environmental function
ReCiPe 2008	ReCiPe
Min GHG emissions or CO2 emissions	GHG
Min energy consumption	EC
Min waste	W
Eco-indicator 99	EI99
Environmental impact indicator	IC
Water consumption	WC
Social function
Jobs created	JC
Maximize justice based employment	EJ
Workers’ damage	WD
Workers’ health	WH
Social Benefit Indicator	SB
Min proximity of employees to work	EP
Min hires and layoff to incentive stability	ES
Max the distance between undesirable facilities and customers	UF
Equity and workload in work	WE
Workers’ safety	WS
Consumer risk	CR
Maximizes humanitarian forces	HF
Min total dissatisfaction of served customer zones	CD
Maximize food security	FS
Charity donations or support	Ch
Training hours	WT
Community service hours	CS
Social responsibility	SR
Social impact to open plants	SI
Hazardous by-products	HP
Uncollected wasted cooking oil	UWCO

Table 2. Table 2: SSC optimization models and solution methods.

Ref.

Problem

CL

IQ

PQ

MP

PT

MT

MH

Eco

Env

Soc

Method

Type

Solutions

This

Eletronic components

yes

yes

yes

yes

yes

yes

yes

NPV

ReCiPe

SB

AUGMECON and Lagrangian heuristic

Matheuristic

Generative

[63]

Numerical

yes

no

no

no

no

yes

no

C

CO2

JC,WD, WU

GA and CS

Meta-heuristic

Generative

[60]

Paper industry

no

yes

yes

yes

yes

no

yes

C

GHG, EC, W

EP,ES

AUGMECON2

Exact

Generative

[57]

Electronic components

yes

yes

yes

yes

yes

yes

yes

C

ReCiPe

JC

AUGMECON and relaxation-based heuristic

Matheuristic

Generative

[51]

Bio-energy SC

no

no

yes

yes

yes

no

yes

ATP

GHG

JC

AUGMECON

Exact

Generative

[50]

Numerical

no

yes

yes

yes

yes

no

yes

P

GHG

JC,WD

$ϵ$-Constraint

Exact

Generative

[47]

Numerical

yes

no

no

no

yes

no

no

C

-

JC

Fuzzy solution approach and weighted sum

Fuzzy

A priori

[43]

Numerical

yes

no

yes

yes

yes

no

no

P

Green degree

CD

Game theory combined with NSGA-II and MOICA

Meta-heuristic

Generative

[36]

Food bank SC

no

no

no

yes

yes

no

no

C

FW, CO2

Ch

Lexicographic ordering

Exact

A priori

[27]

Glass industry

yes

no

yes

no

no

no

no

C

-

JC,WD

Several meta-heuristic, among them ICA, GA, SA and TS

Meta-heuristic

Generative

[25]

Numerical

no

no

yes

yes

no

yes

yes

C

GHG

JC,WE

Swarm intelligence algorithms and VNS

Meta-heuristic

Generative

[23]

Electricity generation

no

yes

yes

no

yes

no

no

C

GHG

JC

AUGMECON

Exact

Generative

[20]

Numerical

yes

no

no

no

no

yes

yes

C

CO2

UF

NSGA-II and SA

Meta-heuristic

Generative

[13]

Numerical

yes

no

yes

yes

no

yes

no

P

CO2

JC,WT, CS,Ch

GP and fuzzy solution approach

Fuzzy

A priori

[59]

Numerical

no

no

yes

no

no

no

yes

C

CO2

JC,WD

Fuzzy solution approach

Fuzzy

Interactive

[54]

Numerical

yes

no

yes

yes

no

no

yes

C

ReCiPe

JC,WS HP

Meta-heuristics and $ϵ$ constraint

Meta-heuristic

Generative

[52]

Numerical

yes

no

no

no

no

yes

yes

C

-

JC

CS

Meta-heuristic

Generative

[48]

Switchgrass-based bioenergy SC

no

yes

yes

yes

yes

yes

no

C

GHG

JC

AUGMECON and TOPSIS

Exact

Generative

[46]

Numerical

yes

no

no

no

yes

no

no

C

-

JC

NSGA-II

Meta-heuristic

A priori

[42]

Numerical

yes

yes

yes

yes

yes

no

no

P

GHG

JC

MOICA

Meta-heuristic

Generative

[41]

Electronic components

yes

yes

yes

yes

yes

yes

yes

NPV

ReCiPe

SB

$ϵ$-Constraint

Exact

A priori

[56]

Biomass SC

no

no

yes

yes

no

yes

yes

P

IC

JC,WH, WS

PCA and AHP

Decomposition

A priori

[21]

Switchgrass-based bioethanol SC

no

yes

yes

yes

yes

yes

yes

C

ReCiPe

SR

fuzzy solution method

Fuzzy Interactive

Generative

[16]

Blood bank SC

no

yes

no

yes

yes

yes

no

C

-

JC

$ϵ$-Constraint

Exact

Generative

[11]

Biofuel SC from coffee crop residues

no

yes

yes

yes

yes

no

no

NPV

CO2

JC,FS

$ϵ$-Constraint

Exact

A priori

[61]

Wine industry

no

no

no

no

no

yes

no

C

CO2

SI

AUGMECON

Exact

Generative

[53]

Healthcare SC for natural disaster response

no

no

yes

no

no

no

yes

C

-

HF

Neural network, GA and PSO

Neural network, Meta-heuristic and GA

A priori

[44]

Bioethanol SC

no

no

yes

yes

yes

no

yes

P

GHG

JC

AUGMECON and Benders decomposition

Matheuristic

Generative

[32]

Textile industry

yes

no

yes

yes

no

yes

no

C

WC

EJ

MOVDO

Meta-heuristic

Generative

[2]

Global manufacturer

no

yes

yes

no

yes

yes

no

C

GHG, W

JC*

GP and $ϵ$-Constraint

Exact

Generative

[64]

Biodisel SC

no

no

no

no

no

no

no

P

CO2

UWCO

Decomposition GA and weighted sum

Heuristic/ Decomposition

A priori

[39]

Bioethanol SC

no

yes

yes

yes

yes

no

yes

C

Eco-costs

JC

GP

Exact

A priori

[3]

Bioethanol SC

no

yes

yes

yes

yes

yes

yes

P

EI99

JC

Robust possibilistic programming

Exact/ Robust

Interactive

[40]

Battery industry

yes

yes

no

yes

yes

no

no

C

ReCiPe

SB

$ϵ$-Constraint

Exact

A priori

[6]

Lumber industry

no

yes

yes

yes

yes

yes

yes

C

GHG

EP

Weighted sum

Exact

A priori

[45]

Medical needle and syringe SC

yes

no

yes

yes

no

no

yes

C

ReCiPe

JC,CR

Benders decomposition and fuzzy solution approach

Matheuristic

Generative

[14]

Glass industry

yes

no

yes

no

no

no

yes

C

GHG, W

JC,WD

ICA based and VNS

Meta-heuristic

Generative

Table 3. Table 3: Decision variables of the SSC formulation introduced in [ 41 ] .

Variable	Description
Continuous decision variables
$S_{mit}$	Amount of product $m$ stocked in entity $i\in I_f\cup I_w$ in time period $t\in T$.
$X_{maijt}$	Amount of item $m\in M$ transported from entity $i\in I$ to entity $j\in I$.
	by transportation mode $a$ in time period $t\in T$.
$P_{mgit}$	Amount of product $m\in M$ produced by technology $g\in G_{prod}$ in factory $i\in I_f$ in time period $t\in T$.
$R_{mgit}$	Amount of product $m$ remanufactured by technology $g\in G_{rem}$ in factory $i\in I_f$ in time period $t\in T$.
$Y{C}_i$	Capacity of entity $i\in I_f\cup I_w$.
$YC{T}_{it}$	Effective use of capacity in entity $i\in I_f\cup I_w$ in time period $t\in T$.
$K{T}_{ait}$	Upper bound to the number of transportation modes $a$ from entity $i\in I$ to another in time period $t\in T$.
Integer decision variables
$K_{ai}$	Number of transportation modes $a\in A$ used to transport products from entity $i\in I$.
$Q_{aijt}$	Number of trips from entity $i\in I$ to $j\in I$ by transportation mode $a$ in time period $t\in T$.
Binary decision variables
$Y_i$	value $1$ indicates that entity $i\in I$ is installed and $0$, otherwise.
$Z_{gmi}$	value $1$ indicates that the technology $g$ is selected to produce product $m$ in factory $i\in I_f$
	and $0$, otherwise.

Table 4. Table 4: Data parameters for the constraints of the SSC formulation introduced in [ 41 ] .

Parameter	Description
Demands and specification of items
$dm{d}_{mit}$	Demand of customer $i\in I_c$ for product $m\in M_{fp}$ in a period $t\in T$.
$BO{M}_{mng}^{prod}$	Bill of materials $m\in M_{rm}$ to produce a product $n\in M_{fp}$
	using technology $g\in G_{prod}$.
$BO{M}_{mn}^{rem}$	Bill of recovery products $m\in M_{rp}$ to remanufacture a final product $n\in M_{fp}$.
$p{w}_m$	Weight of a raw material, final product or recovery product $m\in M_{rm}\cup M_{fp}\cup M_{rp}$ ($kg$).
$ap{u}_m$	Necessary area per unit of $m\in M_{rm}\cup M_{fp}\cup M_{rp}$ ($m^2$).
$ret{F}_m$	Return rate of a recovery product $m\in M_{rp}$.
$LF$	Lifetime in time periods of final products.
Maximum and minimum capacities of entities
$i{c}_{mi}^{max}$	Storage capacity for a product $M_{fp}$ in a entity $i\in I_f\cup I_w$.
$i{c}_{mi}^{min}$	Minimum stock for a product $m\in M_{fp}$ in a entity $i\in I_f\cup I_w$
$p{c}_g^{max}$	Maximum capacity of a technology $g\in G$.
$p{c}_g^{min}$	Minimum use required by a technology $g\in G$.
$e{a}_i^{max}$	Maximum installation area of a entity $i\in I_f\cup I_w$ $\left(m^2\right)$.
$e{a}_i^{min}$	Minimum installation area of a entity $i\in I_f\cup I_w$ $\left(m^2\right)$.
$s{c}_{mi}^{max}$	Maximum amount of a raw material $m\in M_{rm}$ that can be supplied by the supplier $i\in I_p$ in each time period.
$s{c}_{mi}^{min}$	Minimum order of a raw material $m\in M_{rm}$ from a supplier $i\in I_p$ in each time period.
$dis{t}_{ij}$	Distance between two entities $i\in I$ and $j\in I$ ($km$)

Table 5. Table 5: Data parameters for the objective functions of the SSC formulation introduced in [ 41 ] .

Parameter	Description
Costs and related parameters
$te{c}_g$	Cost to install a technology $g\in G$ (€).
$op{c}_g$	The operating cost of a technology $g\in G$ (€).
$ps{u}_m$	Selling price of a final product $m\in M_{fp}$ (€).
$s{c}_m$	Inventory price of a final product $m\in M_{fp}$ (€).
$rp{c}_m$	Cost of recovery product $m\in M_{rp}$ (€).
$rm{c}_m$	Cost of raw material $m\in M_{rm}$ (€).
$av{c}_a$	Average vehicle consumption for $a\in {𝒜}_{truck}$ ($l$ per $100km$).
Social parameters
$sqm{c}_i$	Construction costs for each entity $i\in I_f\cup I_w$ (€).
$l{c}_i$	Labor cost in an entity $i\in I$ (€).
${\mu }_i^{gdp}$	Inverse of the country GDP per capita (GDPPC) where entity $i\in I$ is located.
$pWor{k}_i$	Working condition index in an entity $i\in I$.
$w{i}_i$	Minimum number of workers in entity $i\in I_f\cup I_w$.
$w{a}_a$	Number of workers in transportation mode $a\in A$.
$w{g}_g$	Number of workers to operate technology $g\in G$.
$wps{q}_i$	Minimum number of workers per square feet in entity $i\in I_f\cup I_w$.
Environmental impact parameters
$e{i}_c$	Entity installation impact for category $c\in C$ (per $m^2$).
$e{s}_{ic}$	Impact from supplier $i\in I_p$ for category $c\in C$.
$e{g}_{mgc}$	Production impact for manufacturing $m\in M_{fp}$ using technology $g\in G$
	for category $c\in C$ (per product).
$e{t}_{ac}$	Impact of the transportation mode $k\in K$ for the category $c\in C$ (per $kg$).

Table 6. Table 6: Results for optimizing each objective function individually. Very small and very large values are expressed in E notation.

Function to	Time(s)	Values of the objective functions
be optimized		$f_{eco}$	$f_{env}$	$f_{soc}$
$\max f_{eco}$	0.306	1.096e+06	6.185e+05	11.704
$\min f_{env}$	0.256	-1.094e-07	3.885e+05	50.697
$\max f_{soc}$	0.091	-2.503e-08	4.593e+05	499.081

Table 7. Table 7: Details about the entities of the supply chain.

Entity	Label	Location	Region	GDPPC (PPP)	Construction	Labor
					cost (€)	cost (€)
Supplier	$S_1$	Verona, Italy	Europe	0.98	-	28.1
	$S_2$	Hannover, Germany	Europe	1.24	-	30.4
	$S_3$	Leeds, United Kingdom	Europe	1.06	-	15.3
Factories	$F_1$	Verona, Italy	Europe	0.98	595	28.1
	$F_2$	Hannover, Germany	Europe	1.24	661	30.4
	$F_3$	Leeds, United Kingdom	Europe	1.06	601	15.3
Warehouses	$W_1$	Verona, Italy	Europe	0.98	595	28.1
	$W_2$	Hannover, Germany	Europe	1.24	661	30.4
	$W_3$	Leeds, United Kingdom	Europe	1.06	601	15.3
	$W_4$	Zaragoza, Spain	Europe	0.95	373	21
	$W_5$	Lisbon, Portugal	Europe	0.75	318	12.2
	$W_6$	São Paulo, Brazil	Brazil	0.355	538	8.98
	$W_7$	Recife, Brazil	Brazil	0.355	538	8.98
	$W_8$	Budapest, Hungary	Europe	0.67	282	7.5
	$W_9$	Sofia, Bulgaria	Europe	0.47	270	3.7
Customers	$C_1$	Italy	Europe	0.98	-	28.1
	$C_2$	Germany	Europe	1.24	-	30.4
	$C_3$	United Kingdom	Europe	1.06	-	15.3
	$C_4$	Spain	Europe	0.95	-	21
	$C_5$	Portugal	Europe	0.75	-	12.2
	$C_6$	São Paulo, Brazil	Brazil	0.355	-	8.98
	$C_7$	Recife, Brazil	Brazil	0.355	-	8.98
Airports	$Ai{r}_1$	Zaragoza	Europe	1.19	-	21
	$Ai{r}_2$	Paris-Charles de Gaulle	Europe	1.08	-	32.4
	$Ai{r}_3$	Kortrijk-Wevelgem	Europe	0.95	-	37.2
	$Ai{r}_4$	São Paulo	Brazil	0.355	-	8.98
Seaports	$Se{a}_1$	Hamburg, Germany	Europe	1.24	-	30.4
	$Se{a}_2$	Santos, Brazil	Brazil	0.355	-	8.98

Table 8. Table 8: Optimization of the individual objective functions by an exact strategy, MathLagr and MathFix .

			Function to be optimized
			$\max f_{eco}$	$\min f_{env}$	$\max f_{soc}$
Exact method	Time (s)		154.218	15.676	3742.222
	Values of the	$f_{eco}$	1.405e+09	1.052e+09	-5.713e+09
	objective	$f_{env}$	9.766e+08	9.108e+08	9.249e+08
	functions	$f_{soc}$	895.804	2176.588	40800.432
MathLagr	Time (s)		43.035	22.469	150.603
	Values of the	$f_{eco}$	1.404e+09	8.906e+08	-5.646e+09
	objective	$f_{env}$	9.774e+08	9.108e+08	9.721e+08
	functions	$f_{soc}$	943.260	2573.369	40432.385
	$GAP$ (%)		0.041%	2.273e-3%	0.902%
	Speed up		3.584	0.698	24.848
MathFix	Time (s)		51.515	34.231	152.861
	Values of the	$f_{eco}$	1.401e+09	1.040e+09	-5.669e+09
	objective	$f_{env}$	9.820e+08	9.108e+08	9.288e+08
	functions	$f_{soc}$	828.169	2038.037	40532.255
	$GAP$ (%)		0.230	6.038e-3	0.657
	Speed up		2.994	0.458	24.481

Table 9. Table 9: Results of the multi-objective metrics for the Pareto frontier approximations obtained by AUGMECON2, AugMathLagr and AugMathFix .

	Number	$\mathrm{𝐚𝐌𝐈𝐃}$	$\mathrm{𝐚𝐒𝐍𝐒}$	$\mathrm{𝐑𝟐}$	Time (h)
AUGMECON2	31	12.177	3.010	0.000	153.359
AugMathLagr	31	13.252	2.943	5.558e-04	2.656
AugMathFix	26	11.076	3.021	-6.096e-04	4.913

Table 10. Table 10: Case study: location and installation capacities of factories and warehouses in the solutions obtained by the exact method, MathLagr and MathFix .

		Function to be optimized
		$\max f_{eco}$			$\min f_{env}$			$\max f_{soc}$
		Location	Maximum	Minimum	Location	Maximum	Minimum	Location	Maximum	Minimum
			capacity	capacity		capacity	capacity		capacity	capacity
Factories	$F_1$	$A,L,M$	-	-	$A,L,M$	-	-	$A,L,M$	-	-
	$F_2$	$A,L,M$			$A,L;M$	$A,L,M$		$A,L,M$	$A,L,M$
	$F_3$	$A,L,M$	$A,L,M$		$A,L,M$	$A,L,M$		$A,L,M$	$A,L,M$
Warehouses	$W_1$	$A,L,M$	-	-	$A,L,M$	-	-	$A,L,M$	-	-
	$W_2$				$L,M$	$L,M$		$A$	$A$
	$W_3$				$A,L,M$	$A,L,M$		$A,L,M$	$A,L,M$
	$W_4$				$A$	$A$		$A,L,M$	$A,L,M$
	$W_5$	$M$			$A$	$A$		$A,L,M$	$A,L,M$
	$W_6$				$A,L$	$A,L$		$A,L,M$	$A,L,M$
	$W_7$				$L,M$	$L,M$		$A,L,M$	$A,L,M$
	$W_8$							$A,M$	$A,M$
	$W_9$	$A,L$						$A,L,M$	$A,L,M$

Table 11. Table 11: Number of entities, items and periods in each set of instances.

Instance		$\left|I\right|$	$\left|M_{rm}\right|$	$\left|M_{fp}\right|$	$\left|M_{rp}\right|$	$\left|T\right|$
Prefixes	Suffix
STD, TECHC, RAWC, SUP, CAP	_I17_M4_T3	17	2	1	1	3
STD, TECHC, RAWC, SUP, CAP,	_I17_M4_T6	17	2	1	1	6
STD, TECHC, RAWC, SUP, CAP	_I17_M4_T12	17	2	1	1	12
STD, TECHC, RAWC, SUP, CAP	_I17_M8_T3	17	4	2	2	3
STD, TECHC, RAWC, SUP, CAP	_I17_M8_T6	17	4	2	2	6
STD, TECHC, RAWC, SUP, CAP	_I17_M8_T12	17	4	2	2	12
STD, TECHC, RAWC, SUP, CAP	_I25_M4_T3	25	2	1	1	3
STD, TECHC, RAWC, SUP, CAP	_I25_M4_T6	25	2	1	1	6
STD, TECHC, RAWC, SUP, CAP	_I25_M4_T12	25	2	1	1	12
STD, TECHC, RAWC, SUP, CAP	_I25_M8_T3	25	4	2	2	3
STD, TECHC, RAWC, SUP, CAP	_I25_M8_T6	25	4	2	2	6
STD, TECHC, RAWC, SUP, CAP	_I25_M8_T12	25	4	2	2	12

Equations

$$\sum _{\begin{array}{c}a\in A,j\in I_f\end{array}}{X}_{maijt}\le s{c}_{mi}^{max}{Y}_i,i\in I_{sup},m\in M_{rm},t\in T$$

$$\sum _{\begin{array}{c}a\in A,j\in I_f\end{array}}{X}_{maijt}\ge s{c}_{mi}^{min}{Y}_i,i\in I_{sup},m\in M_{rm},t\in T$$

$$\sum _{m\in M,a\in A,j\in I}{X}_{maijt}\le e{c}_i^{max}{Y}_i,i\in I,t\in T$$

$$\sum _{m\in M,a\in A,j\in I}{X}_{majit}\le e{c}_i^{max}{Y}_i,i\in I,t\in T$$

$$S_{mit}\le i{c}_{mi}^{max}{Y}_i,m\in M_{fp},i\in I_f\cup I_w,t\in T$$

$$S_{mit}\ge i{c}_{mi}^{min}{Y}_i,m\in M_{fp},i\in I_f\cup I_w,t\in T$$

$$Y{C}_i\le e{a}_i^{max}{Y}_i,i\in I_f\cup I_w$$

$$Y{C}_i\ge e{a}_i^{min}{Y}_i,i\in I_f\cup I_w$$

$$\sum _{m\in M,a\in A,i\in I,t\in T}{X}_{maijt}\ge Y_j,j\in I$$

$$\sum _{m\in M,a\in A,j\in I,t\in T}{X}_{maijt}\ge Y_i,i\in I$$

$$Q_{aijt}\le BigM{Y}_i,a\in A,i,j\in I,t\in T$$

$$Q_{aijt}\le BigM{Y}_j,a\in A,i,j\in I,t\in T$$

$$K_{ai}\le BigM{Y}_i,a\in A_{truck},i\in I$$

$$P_{mgit}\ge p{c}_g^{min}{Z}_{gmi},i\in I_f,(m,g)\in H_{prod},t\in T$$

$$P_{mgit}\le p{c}_g^{max}{Z}_{gmi},i\in I_f,(m,g)\in H_{prod},t\in T$$

$$R_{mgit}\ge p{c}_g^{min}{Z}_{gmi},i\in I_f,(m,g)\in H_{rem},t\in T$$

$$R_{mgit}\le p{c}_g^{max}{Z}_{gmi},i\in I_f,(m,g)\in H_{rem},t\in T$$

$$\sum _{g:(m,g)\in H_{prod}}{Z}_{gmi}\le Y_i,m\in M_{fp},i\in I_f$$

$$\sum _{g:(m,g)\in H_{rem}}{Z}_{gmi}\le Y_i,m\in M_{fp},i\in I_f$$

$$\min$$

$$s.t.$$

$$E^{{}^{\prime }}u+E^{{}^{\prime \prime }}v\le d$$

$$u\in \text{mathds}{R}_{\ge 0},v\in \text{mathds}{Z}_{\ge 0}^{n{}^{\prime \prime }},w^{\prime }\in \{0,1{\}}^{b^{\prime }},w^{\prime \prime }\in \{0,1{\}}^{b^{\prime \prime }}$$

$$\begin{array}{r}{f}_{soc}=\sum _{i\in I_f\cup I_w}{\mu }_i^{gdp}pWor{k}_i(w{i}_i{Y}_i+wps{q}_{i}Y{C}_i)+\\ \sum _{i\in I_f}{\mu }_i^{gdp}pWor{k}_i\sum _{(m,g)\in H}w{g}_g{Z}_{gmi}+\sum _{i\in I}{\mu }_i^{gdp}pWor{k}_i\sum _{a\in A_{truck}}w{a}_a{K}_{ai}+\\ \sum _{i,j\in I}{\mu }_i^{gdp}pWor{k}_{i}dis{t}_{ij}\sum _{m\in M}p{w}_m\sum _{a\in A_{plane}\cup A_{ship}}\frac{w{a}_a}{yth}\sum _{t\in T}{X}_{maijt}\end{array}$$

$$\text{mbox}(MOP):\min$$

$$s.t.$$

$$E^{{}^{\prime }}u+E^{{}^{\prime \prime }}v\le d$$

$$f_{env}(u,v,w^{\prime },w^{\prime \prime })+l_{env}={ϵ}_{env}$$

$$-f_{soc}(u,v,w^{\prime },w^{\prime \prime })+l_{soc}={ϵ}_{soc}$$

$$u\in \text{mathds}{R}_{\ge 0},v\in \text{mathds}{Z}_{\ge 0}^{n{}^{\prime \prime }},w^{\prime }\in \{0,1{\}}^{b^{\prime }},w^{\prime \prime }\in \{0,1{\}}^{b^{\prime \prime }},$$

$$l_{env},l_{soc}\in \text{mathds}{R}_{\ge 0}$$

$$\text{mbox}(PRL):\min L(u,v,w^{\prime },w^{\prime \prime },\lambda )=$$

$$+$$

$$s.t.$$

$$E^{{}^{\prime }}u+E^{{}^{\prime \prime }}v\le d$$

$$f_{env}(u,v,w^{\prime },w^{\prime \prime })+l_{env}={ϵ}_{env}$$

$$-f_{soc}(u,v,w^{\prime },w^{\prime \prime })+l_{soc}={ϵ}_{soc}$$

$$u\in \text{mathds}{R}_{\ge 0},v\in \text{mathds}{Z}_{\ge 0}^{n{}^{\prime \prime }},w^{\prime }\in \{0,1{\}}^{b^{\prime }},w^{\prime \prime }\in \{0,1{\}}^{b^{\prime \prime }}$$

$$l_{env},l_{soc}\in \text{mathds}{R}_{\ge 0}$$

$$GAP(f_i(S^{{}^{\prime }}),f_i(S^{\ast }))=\frac{|f_i(S^{{}^{\prime }})-f_i(S^{\ast })|}{\max (f_i(S^{\ast }),f_i(S^{{}^{\prime }}))}$$

$$I_{R2}(P^A,U)=\frac{1}{|U|}\left(\sum _{\mu \in U}\underset{S\in P^Z}{\max }\{\mu (S)\}-\sum _{\mu \in U}\underset{a\in P^A}{\max }\{\mu (a)\}\right)$$

$$u_{\gamma }(S)=\underset{i\in \{eco,env,soc\}}{\max }{\gamma }_i|f_i^{\ast }-f_i(S)|$$

$$GAPM(f(S^{{}^{\prime }}),f^I)=\sqrt{\sum _{i\in \{eco,env,soc\}}{\left(GAP(f_i(S^{{}^{\prime }}),f_i^{\ast })\right)}^2}$$

Full text

††journal: Paper accepted to European Journal of Operational Research

\externaldocument

table_related

\externaldocumenttable_labels

\externaldocumenttable_params

\nonumnote

Accepted to the European Journal of Operational Research. ©2021. This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/

1 Introduction

The change in mentality of organizations to go beyond profit maximization and consider the long-term economic success by taking into account social and environmental responsibilities has raised the interest of researchers in the subject of Sustainable Supply Chain (SSC) management [55, 17, 4]. In particular, multi-objective SSC management models have relied on the integration of social, environmental and economic issues – the so-called triple bottom line – to explicitly define the long-term planning of sustainable companies [10]. Multi-objective optimization plays an important role in approaching SSC management since a compromise among conflicting goals related to profit, environment and social issues can be reached.

A measure commonly used as an economic criterion in SSC management optimization models is the Net Present Value (NPV) [12]. There is a plethora of measures [1] to quantitatively assess the environmental impacts related to, for example, greenhouse gas emissions, pollution and resource usage. Additionally, there exist some methods specifically designed to quantify these impacts such as Eco-Indicator 99 [24] and ReCiPe 2008 [22], which are based on the Life Cycle Analysis (LCA) methodology.

In spite of the relevance of the social impacts in the SSC management, a minority of models in the literature optimize this indicator along with environmental and economic indicators, e.g. in [40, 6, 9, 41, 57]. Eskandarpour et al. [17] discuss various indicators to assess social responsibilities, most of which related to human rights and social justice laws. In particular, the authors suggest metrics that encourage the creation of job vacancies [40, 6, 41], employment stability [6] and working conditions. Recently, Mota et al. [40] introduced the Social Benefit Indicator (SB) whose calculation relies on the total number of jobs created in each location across the supply chain. This indicator weighs preferably the creation of entities in less developed countries. Mota et al. [41] have optimized this social indicator in the multi-objective SSC management problem under investigation. This indicator is recommended in situations that involve hiring and do not require layoffs [41]. In addition to optimizing the SB, the SSC management problem introduced in [41] aims at maximizing the NPV and minimizing the environmental impacts quantified by ReCiPe 2008.

The problem in question models a generic SSC management problem that manufactures multiple products in a planning horizon composed of multiple periods. The appeal of the model lies in its applicability to a wide variety of industries, since it integrates a number of strategic and tactical decisions related to, for example, the use of technologies to manufacture and refurbish products and the shipment of items through different transport modes, among others. The authors solved mono-objective problems to minimize the environmental objective function and to maximize the economic, i.e., the NPV, and social objective functions. Moreover, they have also considered two scenarios to maximize the social objective function with additional constraints for the NPV value. In these scenarios, they required the NPV value to be at least $85\%$ and $95\%$ of the value found when optimizing the mono-objective problem that maximizes the economic objective function. The computational difficulty to solve the resulting problem was the main reason why the authors did not consider the three objective functions simultaneously when optimizing the problems in the case studies.

Many studies employ heuristic methods instead of exact solution strategies to solve SSC management problems within a reasonable time [14, 54, 25]. However, to design a heuristic to find feasible solutions for the SSC management problem introduced in [41] is particularly challenging. This is due to the integration of several complicated constraints involving binary and integer variables with strategic planning decisions. In line with this, a promising approach to this problem is the relaxation of complicated constraints and decision variables. Lagrangian decomposition allows us to explore these aspects and, in specific, can provide efficient methods for this problem by relaxing complicated constraints while penalizing their violation.

In this context, the primary contributions of this paper are:

The development of an efficient Lagrangian matheuristic for SSC management problems, here called AugMathLagr, to tackle the SSC problem proposed in [41].

• 2.

The adaptation of the AugMathFix matheuristic introduced in [57] to the SSC management problem put forward in [41].

• 3.

A test bed of artificial instances and an instance generator loosely based on the study of real data.

The present paper also reports the computational experience with AugMathLagr on a case study and a set of medium-sized randomly generated instances. Results from of AugMathLagr were compared with results obtained with the Augmented $\epsilon$-Constraint method (AUGMECON2) [38], which is an enhancement of AUGMECON [37], and with the adapted AugMathFix. The comparison involved evaluating the solution quality and time-to-solution. The proposed matheuristic achieved very good results in significantly lesser computational times than AUGMECON2 and AugMathFix.

The remainder of this paper is organized as follows. Section 2 reviews the related methods for solving multi-objective SSC management problems; Section 3 presents the SSC model to which we propose a Lagrangian heuristic; Section 4 describes the proposed multi-objective matheuristic, AugMathLagr, and briefly discusses the adapted AugMathFix; Section 5 introduces the random instance generator; Section 6 presents the computational experiments and Section 7 sums up the paper drawing some conclusions and suggesting future works.

2 Related works

In this section, we review SSC management multi-objective optimization models and solution methods related to this paper. Due to the wide range of SSC problems, we narrowed the scope of our search to forward and closed-loop multi-objective SSC models that aim at optimizing economic, environmental and social objectives and that integrate both strategic and tactical planning. The integration of strategic and tactical planning decisions prevents optimization models from finding sub-optimal solutions [34, 20]. Tactical decisions, such as the inventory and production quantities in warehouses and factories, may affect long-term decisions on where to locate these entities. They represent the average values of these variables within the periods considered. For further details on optimization models and solution methods for SSC management problems, we refer to Barbosa-Póvoa et al. [4].

Consider the notations for the objective functions presented in Table 1. Table 2 shows SSC optimization models and solution methods published in the last 5 years. In this table, columns “CL” to “MH” indicate whether or not the study is characterized by: closed-loop or forward supply chains (CL); inventory decisions (IQ); production quantity choices (PQ); multiple products (MP); multi-period planning horizon (PT); multiple transportation modes (MT) and multiple technologies to manufacture the products (MH). Columns “Eco”, “Env” and “Soc” show respectively the economic, environmental and social indicators optimized by the model or solution method. The references marked as “-” in column “Env” did not specify which environmental impact indicators they use in the optimization. Column “Method” presents the optimization approach to solve the corresponding problem. Column “Type” mentions the type of solution method the studies consider among exact, meta-heuristic and matheuristic methods. Finally, column “Solutions” classifies the solution methods reported in the corresponding references as generative, a priori or interactive. Generative methods approximate the Pareto frontiers; a priori methods assign preferences to the objectives; and interactive methods require the intervention of a decision-maker.

One may observe that only the closed-loop models presented in [57] and [41] include all the strategic and tactical decisions listed in Table 2. The model introduced in [57] also considers the carbon market. Mota et al. [41], on the other hand, model cross-docking transportation through airports and seaports and further explore transportation decisions such as the number of vehicles and trips required to deliver products. Despite including all the strategic and tactical decisions presented in Table 2, Ghaderi et al. [21], Bairamzadeh et al. [3], Boukherroub et al. [6] do not address the reverse flow in their models. The variety of decisions each model takes into account suggests that generic models should include as many decisions as possible in order to represent a wider range of supply chains.

The minimization of costs and profit maximization are the most commonly employed economic objective functions in the reviewed SSC models. Most studies minimize the environmental indicators related to Greenhouse Gases (GHG) or specific CO2 emissions, despite the European Commission putting forward the recommendation that Lyfe Cycle Approaches (LCA) are the best available methodology to evaluate the environmental impact [19]. Other case-oriented indicators are the minimization of energy and water consumption besides the waste resulting from SC operations.

Social indicators are difficult to quantify [17], being commonly linked to the planning situation. The majority of references in Table 2 regard the socio-economic development of regions by optimizing job creation indicators. Mota et al. [40, 41] designed the Social Benefit Indicator (SB) to favor the job creation in less developed countries. SB assesses the development of a country according to its Gross Domestic Product (GDP) [7], which is a widely accepted index to measure the economic activity of a nation. In this context, countries with less industrial development have a lower GDP and, according to the defined SB, are less developed. The authors recommend, however, its use only in planning situations that preferably involve hiring and do not require lay-offs. Another concern about job creation is related to working conditions. To this end, other widely employed indicators take into account workers’ safety and health. When combined, these indicators can help to prevent the creation of jobs in poor working conditions.

A vast body of literature on SSC management multi-objective solution methods is based on scenario analysis, the $\epsilon$-Constraint or weighting aggregation methods and heuristics or matheuristics.

To approach a multi-objective minimization problem, the $\epsilon$-Constraint method [35, 26] optimizes one of its objective functions and adds constraints to assign upper bounds for the other objective functions. In particular, Mavrotas [37] introduced the Augmented $\epsilon$-Constraint method (AUGMECON) that solves a sequence of mono-objective problems. Later, Mavrotas & Florios [38] proposed AUGMECON2 as an improved version of AUGMECON. Examples of works that employ AUGMECON to approach multi-objective SSC management problems are found in [48, 44, 60].

Exact-based methods often employ mixed-integer linear solvers such as CPLEX to solve the resulting mono-objective problems. In practice, these methods are computationally costly to approach most SSC problems. In such cases, AUGMECON and AUGMECON2 are impracticable to approximate the Pareto frontier. To overcome this limitation, Mota et al. [41] optimized the SSC problem using scenario analysis.

Goal Programming (GP) is an exact-based method used in a number of studies presented in Table 2. GP methods aim at finding solutions by minimizing the deviance from pre-defined values for each objective function, also known as goals. In the references indicated in Table 2, there are both generative and a priori versions of GP methods. When dealing with uncertainty in SSC, fuzzy-based solution approaches are commonly employed as a priori or interactive methods.

Several multi-objective algorithms to heuristically solve SSCs have been introduced. Among them, one can find Non-dominated Sorting Genetic Algorithm (NSGA-II) [46, 20, 43], Genetic Algorithms (GA) [64, 53, 27, 63], Variable Neighborhood Search (VNS) [14, 25], Imperative Competitive Algorithm (ICA) [14, 27], multi-objective ICA (MOICA) [42, 43], Cuckoo Search (CS) [52, 63], Multi-Objective Vibration Damping Optimization (MOVDO) [32], Particle Swarm Optimization (PSO) [53, 25], Simulated Annealing (SA) [27, 20] and Tabu Search (TS) [27]. Other strategies such as Neural Networks [53], Principal Component Analysis (PCA) and Analytical Hierarchical Process (AHP) [56] can also be found in the literature.

Matheuristics are powerful solution methods that hybridize heuristic methods with mathematical programming techniques. A few matheuristics can be found in the literature to approach multi-objective SSC problems. The matheuristic proposed by Osmani & Zhang [44] combined AUGMECON with a Benders decomposition to find Pareto solutions for a SSC management problem. Zhang et al. [64] also introduced a matheuristic that employs a Benders decomposition to split a SSC problem into two subproblems associated with the strategic and tactical decisions. Tautenhain et al. [57] introduced a multi-objective matheuristic founded on AUGMECON2, called AugMathFix, which iteratively solves mono-objective SSC management problems with strategic constraints relaxed. In this case, for the method to find a feasible solution to the original problem in a reduced computational time, the variables related to the tactical decisions are fixed at the values of the solution of the relaxed problem. Moreover, at each iteration a local search strategy is applied to the feasible solutions to enhance the quality of the solutions found.

As presented in Table 2, with the exception of AugMathFix, the existing SSC heuristics do not account for all the decisions considered in the model introduced in [41]. However, the resulting integrated SSC problem has a significant amount of integer variables, forcing the authors to simplify its solution through scenario analysis [41]. Optimization solvers like CPLEX have a better performance in solving mixed-integer programs when they have a reduced number of integer variables. In this context, matheuristics that decompose the original problem into simpler ones or that relax binary and integer variables, for example, may speed up mixed-integer solvers in finding solutions.

Lagrangian decomposition methods, in particular, are potential candidates to solve SSC optimization problems since they allow us to simplify complicated constraints, such as those including binary variables, for example. These methods presented good results in finding solutions for related SC problems. Heidari-Fathian & Pasandideh [29] and Yousefi-Babadi et al. [62] studied the Lagrangian decomposition of multi-objective supply chain problems. Heidari-Fathian & Pasandideh [29] employed the $\epsilon$-Constraint method to transform the multi-objective SSC management problem into a sequence of mono-objective problems. The resulting mono-objective problems were then solved by a Lagrangian heuristic. The pair of constraints to be relaxed was the one that, when removed from the problems, resulted in the lowest computational running times.

In the introduced mathematical formulation, Yousefi-Babadi et al. [62] did not consider environmental issues in any of the objective functions, however, they included recycling centers. The authors employed a weighting strategy to aggregate the objective functions into a mono-objective function and then relaxed the most complex constraints of the problem to then solve it using an optimization solver.

Lagrangian-based heuristics are extensively studied to approach mono-objective supply chain problems [18, 33, 15, 65, 49]. Eskigun et al. [18] reduced the original capacitated network design problem into simpler and independent subproblems to approach the design of vehicle distribution centers of several instances constructed from industrial data by a Lagrangian heuristic. In the Lagrangian heuristic proposed by Lidestam & Rönnqvist [33], the authors decomposed the supply chain planning problem into two different subproblems associated with different stages of cellulose production. Elhedhli & Merrick [15] studied a supply chain network design problem that aims to minimize environmental costs due to carbon dioxide emissions. The authors employed a Langragian relaxation by decomposing the original problem according to entity types and warehouse site. Zhang et al. [65] approached a supply chain problem that required distribution centers between (i) suppliers and factories, and (ii) factories and customers. The authors relaxed the constraints which ensured that each customer and supplier was assigned to only one distribution center in the Lagrangian heuristic they suggested. Rafie-Majd et al. [49] studied a supply chain of perishable products that takes into consideration fuel consumption and product wastage. The Lagrangian heuristic presented by the authors considers the relaxation of the vehicle capacity constraints, a nonlinear constraint and a constraint related to the allocation of distribution centers to customers.

To the best of our knowledge, this is the first Lagrangian-based method to solve an integrated SSC problem that considers all the dimensions of sustainability, as presented in Table 2.

The next section describes the SSC management problem, the focus of this paper.

3 Problem description

The studied (SSC) is composed of suppliers, factories, warehouses, customers, airports and seaports, and is tailored to support multi-period planning. Moreover, there are three types of items in the SSC: raw materials, final products and recovered products.

Figure 1 illustrates a flow network whose nodes are the entities of the SSC, and the arcs represent the flow of goods between the sites. In this figure, the rectangles represent airports and seaports whereas the ellipses represent the remaining entities. The double-sided arrows indicate that the corresponding arcs represent the flow of both final products and recovered products between the sites. The rightwards arrow represents the arc indicating the flow of raw materials from suppliers to factories.

Factories employ production technologies to manufacture final products from a bill of raw materials and may also use remanufacturing technologies to reuse products. The problem allows the storage of final products sent from factories to warehouses. Factories and warehouses can ship final products to meet customer demands, indicated by the arc linking these entities. Customers return used products to factories and warehouses after the end of their lifetime, which is defined as a parameter measured in the number of periods of the planning horizon in the model. Warehouses can then return the recovered products to the factories.

Land transportation is responsible for transporting raw materials from suppliers to factories. The possible transportation modes between factories, warehouses and customers are by land, air or sea. Trucks carry items from and to entities through land transportation, whereas air and sea transportation are only allowed from/to airports or seaport hubs.

The goals of the studied SSC management are related to the triple bottom line: (i) the maximization of the Net Present Value (NPV) as the economic function; (ii) the minimization of the environmental impact evaluated by ReCiPe 2008 [22] as the environmental function; and (iii) the maximization of the Social Benefit Indicator (SB) [41] as the social function. This social indicator favors the creation of jobs in less developed countries. The Gross Domestic Product (GDP) [7] was used to assess and compare the country’s industrial base development.

SB can also be adjusted according to the context of the supply chain under study. In analyzing supply chains at the regional level, for example, the development of regions might be assessed through unemployment rates or population densities [40].

In this paper, in order to prevent job creation in poor working environments, we also considered an additional indicator linked to the labor conditions in entities in the social objective. Moreover, we included the impacts of suppliers in the environment function as well.

In the studied SSC, the set of periods is given by $T=\{1,2,\dots,|T|\}$ and the set of entities by $I=I_{sup}\cup I_{f}\cup I_{w}\cup I_{c}\cup I_{air}\cup I_{port}$, where $I_{sup},I_{f},I_{w},I_{c},I_{air}$ and $I_{port}$ are, respectively, the sets of suppliers, factories, warehouses, customers, airports and seaports. The set of items is given by $M=M_{rm}\cup M_{fp}\cup M_{rp}$, where $M_{rm},M_{fp}$ and $M_{rp}$ are the sets of raw materials, final products and recovered products, respectively. In this SSC, intermediate products required in the bill to manufacture other products are also defined as final products. The set of transportation modes is given by $A=A_{truck}\cup A_{plane}\cup A_{boat}$, where $A_{truck},A_{plane}$ and $A_{boat}$ are the set of trucks, airplanes and ships, respectively. The technologies set is given by $G=G_{prod}\cup G_{rem}$, where $G_{prod}$ and $G_{rem}$ are, respectively, the production and the remanufacturing technologies sets. Moreover, consider $H_{prod}=\{(m,g)\mbox{ such that }m\in M_{fp}\mbox{ can be manufactured using technology }g\in G_{prod}\}$ and $H_{rem}=\{(m,g)\mbox{ such that }m\in M_{fp}\mbox{ can be remanufactured using technology }g\in G_{rem}\}$. In Section 1 of the Supplementary Material, we present a brief discussion on how this formulation can be adapted to a more general network.

This model assumes data regarding economic, environmental and social criteria are available. We acknowledge, however, that the access of information in the supply chains is a problem, not often because it does not exist but because it is not treated. So, from our experience, a large period of time must be assigned to the identification of the required data. Also and in particular in what concerns the environmental and social data, assumptions must be often taken. For instance, in the environmental component, similar operations may be taken as the basis to feed certain supply chain operations that do not have available data. This need implies that a sensitivity analysis should be performed on the parameters more subject to uncertainty. This type of work allows addressing the lack of information by looking at which information is more important, that is, which one has a more significant impact on strategic and tactical decisions.

The decisions associated with the management of the studied SSC relate to the amount of raw materials acquired from each supplier; the opening of factories and warehouses and their capacities; production and remanufacturing technologies assigned to factories and recycling centers, respectively; amount of items produced and remanufactured in factories; storage levels in warehouses; and shipment of items between entities using different transportation modes. Table 3 presents a complete list of decision variables. Tables 4 and 5 present the parameters associated with the model constraints and objective functions, respectively.

The following section presents the multi-objective formulation of the SSC management problem studied in [41].

3.1 Multi-objective formulation

To fully understand the proposed method, we show the constraints associated with the Lagrangian relaxation we perform in our method in (1)-(19).

Constraints (1) and (2) ensure that the amount of raw materials acquired by the factories from the selected suppliers are within the interval $[sc^{min}_{mi},sc^{max}_{mi}]$.

[TABLE]

Constraints (3) and (4) define, respectively, $ec_{i}^{max}$ as the maximum in- and out-flow of products between a pair of installed entities $i,j\in I$.

[TABLE]

Constraints (5) and (6) ensure that the amount of final products $m$ stored at installed factory or warehouse $i$ is within the interval $[ic_{mi}^{min},ic_{mi}^{max}]$.

[TABLE]

Constraints (7) and (8) ensure that the installation area of each factory or warehouse $i$ is within the interval $[ea_{i}^{min},ea_{i}^{max}]$, where $ea_{i}^{min}$ and $ea_{i}^{max}$ are non-negative scalars.

[TABLE]

Constraints (9) and (10) guarantee that only entities selected to be installed can receive or send items.

[TABLE]

Constraints (11) and (12) guarantee that if the number of trips to transport items from/to an entity is higher than [math], the entity must be installed.

[TABLE]

Constraints (13) restrict the purchase of trucks $a\in A_{truck}$ only at installed entities $i\in I$.

[TABLE]

Let $pc^{max}_{g}$ and $pc^{min}_{g}$ be, respectively, the maximum and minimum amounts of products that technology $g\in{G}$ can produce. For each factory $i$ and time period $t$, constraints (14) and (15) ensure that the production levels of the final product $m$ using technology $g$, ${\forall}(m,g)\in H_{prod}$, are within the interval $[pc^{min}_{g},pc^{max}_{g}]$. Analogously, constraints (16) and (17) ensure that the remanufacturing levels of technology $g$ are within the interval $[pc^{min}_{g},pc^{max}_{g}]$.

[TABLE]

Constraints (18) and (19) define that production or remanufacturing technologies can only be selected in installed factories.

[TABLE]

Let $n^{{}^{\prime}},n^{{}^{\prime\prime}},b^{{}^{\prime}}$ and $b^{{}^{\prime\prime}}$ be natural numbers that define the dimension of the following decision variables: $u\in\mathds{R}^{n^{{}^{\prime}}}_{{\geq 0}},v\in\mathds{Z}_{{\geq 0}}^{n^{{}^{\prime\prime}}}$, $w^{\prime}\in\{0,1\}^{b^{\prime}}$ and $w^{\prime\prime}\in\{0,1\}^{b^{\prime\prime}}$. Moreover, consider ${f}_{eco}$, ${f}_{env}$, ${f}_{soc}:\mathds{R}^{n{{}^{\prime}}}\times\mathds{Z}_{{\geq 0}}^{n{{}^{\prime\prime}}}\times\{0,1\}^{b^{\prime}}\times\{0,1\}^{b^{\prime\prime}}\rightarrow\mathds{R}$ as the economic, environmental and social functions of the problem, respectively. The economic function ($f_{eco}$) and social function ($f_{soc}$) must be maximized to achieve the best values of NPV and of the SB, respectively. Thereby, without loss of generality, we minimize $-f_{eco},f_{env}$ and $-f_{soc}$ to describe the multi-objective minimization problem (20)-(23).

[TABLE]

where $p$ and $q$ are natural numbers, $\mathcal{A}^{{}^{\prime}}\in\mathds{R}^{p\times n^{{}^{\prime}}}$ and $E^{{}^{\prime}}\in\mathds{R}^{q\times n^{{}^{\prime}}}$ are parameters associated with the real variables; $\mathcal{A}^{{}^{\prime\prime}}\in\mathds{R}^{p\times n^{{}^{\prime\prime}}}$ and $E^{{}^{\prime\prime}}\in\mathds{R}^{q\times n^{{}^{\prime\prime}}}$ are parameters associated with the integer variables; $\beta^{\prime}\in\mathds{R}^{p\times b^{\prime}}$ and $\beta^{\prime\prime}\in\mathds{R}^{p\times b^{\prime\prime}}$ are parameters associated with the binary variables $Y$ and $Z$ of the problem, respectively; $d\in\mathds{R}^{q}$ are the parameters not linked to decision variables; $\mathcal{A}^{{}^{\prime}}u+\mathcal{A}^{{}^{\prime\prime}}v\leq\beta^{\prime}w^{\prime}+\beta^{\prime\prime}w^{\prime\prime}$ are the constraints associated with the binary variables and presented in (1)-(19); and $E^{{}^{\prime}}u+E^{{}^{\prime\prime}}v\leq d$ are the remaining constraints.

Constraints $E^{{}^{\prime}}u+E^{{}^{\prime\prime}}v\leq d$ are related to the tactical planning. They model (i) the material balance at entities to ensure that their in-flow is equal to the out-flow, considering item manufacturing, remanufacturing and stock; (ii) cross-docking at the airports and seaports; (iii) the maximum flow between entities; (iv) product delivery to meet customer demands; (v) the return of used products after $LF$ time periods; and (vi) transportation-related decisions such as the maximum contracted capacity of airports and seaports and maximum investment in trucks.

To calculate $f_{eco}$, i.e., the NPV, Mota et al. [41] considered the difference between the profits from selling products and the costs involved in the SSC. The authors also take into account the depreciation of capital invested in the SSC. The SSC costs are those related to raw material acquisition, product manufacturing, recovery product acquisition, remanufacturing activities, transportation, handling costs at hub airport and seaport terminals, costs associated with contracts with airline and seaport companies, inventory costs and labor costs.

Function $f_{env}$ is calculated according to ReCiPe 2008 as the normalized sum of the environmental impacts due to production manufacturing and remanufacturing, transportation and entity installation. In this paper, we also consider impacts related to raw material acquisition from suppliers, which are given by $\sum_{c\in CC}\eta_{c}(\sum_{t\in T,m\in M_{rm},a\in A,i\in I_{p},j\in I_{f}}es_{ic}X_{maijt})$ where $CC$ is the set of indicators in ReCiPe 2008 and $\eta_{c}$ is a normalization factor for indicator $c\in CC$.

The SB proposed in [41] as the social objective function rewards the creation of jobs due to entity installation, manufacturing and remanufacturing technology selection and transportation decisions entities in countries with lower GDP. Such an indicator composes the social objective function $f_{soc}$ of the model studied in this paper to be maximized, that also considers the working conditions in the entities, as presented in Equation (24).

[TABLE]

where $yth$ is the number of years in the planning horizon. For more details about the formulation, we refer to [41].

The following section discusses the solution methods proposed in this paper.

4 Proposed method

This section describes the proposed Lagrangian matheuristic – called AugMathLagr — to approach the multi-objective SSCM problem discussed in the previous section.

AugMathLagr heuristically solves the multi-objective SSCM problem by following the same core strategy as the Augmented $\epsilon$-Constraint Method (AUGMECON2) [38]. AugMathLagr introduces the Lagrangian-based heuristic for solving mono-objective problems as an innovation in relation to AUGMECON2.

Mavrotas & Florios [38] proposed AUGMECON2 as an improvement of the $\epsilon$-Constraint method. Its aim is to identify the Pareto set of a multi-objective problem by systematically solving a sequence of mono-objective problems ($\epsilon$-Constrained problems). In line with this, the multi-objective problem (20)-(23) can be approached by AUGMECON2 through the solution of the mono-objective problem (25)-(31), referred here to as MOP.

[TABLE]

where $l_{env}$ and $l_{soc}$ are the slack variables of the $\epsilon$-Constraints (28) and (29); $\epsilon_{env}$ and $\epsilon_{soc}$ are scalar values defined as thresholds of $f_{env}$ and $-f_{soc}$, respectively; $r_{env}$ and $r_{soc}$ are positive scalars which are the absolute value of the difference between the best and worst possible values of the functions $f_{env}$ and $-f_{soc}$, respectively; and $eps\in{\mathds{R}_{+}}$ is a small value to promote alternative optimal solutions for $-f_{eco}$ with the best possible values of $f_{env}$ and $-f_{soc}$. The coefficients $1$ and $0.1$ on the objective function mean that $f_{env}$ must be prioritized over $-f_{soc}$. Their values were defined according to Mavrotas & Florios [38].

In this paper, we estimate the worst possible value for functions $f_{env}$ and $-f_{soc}$ as the worst values they assume on the solutions that optimize each objective function. To define the values of $\epsilon_{env}$ and $\epsilon_{soc}$, AUGMECON2 creates a grid of evenly distributed points in the Cartesian plane limited by the best and worst possible values for $f_{env}$ and $-f_{soc}$. The Pareto frontier approximation is composed of the solutions of the MOP considering pre-defined values of $\epsilon_{env}$ and $\epsilon_{soc}$.

Let us consider AUGMECON2 to solve the MOP. As AUGMECON2 is an iterative strategy, an upper index is used on the functions and variables of the MOP to indicate the solutions of a given iteration. Moreover, to assign values to $\epsilon_{j},j\in\{env,soc\}$, considering an imposed number of $dg$ points equally distributed on the grid, the step values of $\epsilon_{j}$ are defined and referred to as $step_{j}=\frac{r_{j}}{dg},j\in\{env,soc\}$.

In the first iteration, iteration [math], an initial value for $\epsilon_{soc}$ must be considered, here denoted by $\epsilon_{soc}^{(0)}$, which can be, for example, the nadir point of $-f_{soc}$. According to constraints (25), ${-f_{soc}(u,v,w^{\prime},w^{\prime\prime})}+l^{(0)}_{soc}=\epsilon_{soc}^{(0)}$, hence, $l_{soc}^{(0)}=\epsilon^{(0)}_{soc}-{(-f_{soc}(u,v,w^{\prime},w^{\prime\prime}))}$. Then, to update $\epsilon_{soc}$, i.e., to define $\epsilon^{(1)}_{soc}$, $\epsilon^{(0)}_{soc}$ is decremented by an scalar $step_{soc}$ and in the next iteration, the values are updated according to: ${-f_{soc}(u,v,w^{\prime},w^{\prime\prime})}+l^{(1)}_{soc}=\epsilon^{(0)}_{soc}-step_{soc}$. First let $l_{soc}^{(0)}-step_{soc}\geq 0$ or, equivalently, $l_{soc}^{(0)}\geq step_{soc}$. Thereby, when $l_{soc}^{(1)}=l_{soc}^{(0)}-step_{soc}$, constraint (29) ensures ${-f_{soc}({u,v,w^{\prime},w^{\prime\prime}})}+l_{soc}^{(0)}-step_{soc}=\epsilon_{soc}^{(0)}-step_{soc}$, i.e., ${-f_{soc}(u,v,w^{\prime},w^{\prime\prime})}+l_{soc}=\epsilon_{soc}$, that is the same problem solved in the current iteration. Therefore, to avoid solving redundant problems, it is necessary to choose a step size whose $l_{soc}<step_{soc}$ holds. AUGMECON2 does that by selecting the step size where $\lfloor{\frac{l_{soc}}{step_{soc}}}\rfloor$.

Mavrotas [37] observed that if a problem is infeasible for a given value of $\epsilon_{soc}$, for smaller values of $\epsilon_{soc}$, it will also be infeasible. Therefore, in this case, AUGMECON2 halts at decrementing $\epsilon_{soc}$ to avoid solving unnecessary problems. This mechanism enables AUGMECON2 to investigate a lower number of problems than an enumerative $\epsilon$-Constraint method.

In the next section, we thoroughly explain the MathLagr matheuristic, which is the Lagrangian heuristic to solve mono-objective problems in AugMathLagr.

4.1 MathLagr

MathLagr finds heuristic solutions through the Lagrangian relaxation of each MOP. The constraints associated with binary variables $\mathcal{A}^{{}^{\prime}}u+\mathcal{A}^{{}^{\prime\prime}}v\leq\beta^{\prime}w^{\prime}+\beta^{\prime\prime}w^{\prime\prime}$ can be decomposed into equivalent constraints $\mathcal{A}^{{}^{\prime}}u+\mathcal{A}^{{}^{\prime\prime}}v\leq\beta^{\prime}+\beta^{\prime\prime}w^{\prime\prime}$ and $\mathcal{A}^{{}^{\prime}}u+\mathcal{A}^{{}^{\prime\prime}}v\leq\mathcal{M}w+\beta^{\prime\prime}w^{\prime\prime}$, $\mathcal{M}\in\mathds{R}^{p\times b^{\prime}}$ being a matrix whose elements are $BigM$ values. The method keeps constraints $\mathcal{A}^{{}^{\prime}}u+\mathcal{A}^{{}^{\prime\prime}}v\leq\beta^{\prime}+\beta^{\prime\prime}w^{\prime\prime}$ in the problem and relaxes constraints $\mathcal{A}^{{}^{\prime}}u+\mathcal{A}^{{}^{\prime\prime}}v\leq\mathcal{M}w^{\prime}+\beta^{\prime\prime}w^{\prime\prime}$.

As the elements of $\beta^{\prime}$ are a trivial estimation of the $BigM$ values of $\mathcal{M}$, we shall refer to the relaxed constraints as $\mathcal{A}^{{}^{\prime}}u+\mathcal{A}^{{}^{\prime\prime}}v\leq\beta^{\prime}w^{\prime}+\beta^{\prime\prime}w^{\prime\prime}$.

Let $\lambda\in\mathds{R}^{p},\lambda\geq 0$, be the $p$-dimensional vector of Lagrange multipliers associated with the relaxed constraints. Each Lagrange multiplier, $\lambda_{i}$, penalizes in the objective function the corresponding violation of constraints $a^{{}^{\prime}}_{i,}u+a^{{}^{\prime\prime}}_{i,}v\leq\beta^{\prime}_{i,}w^{\prime}+\beta^{\prime\prime}_{i,}w^{\prime\prime}$, where $a^{{}^{\prime}}_{i,}$, $a^{{}^{\prime\prime}}_{i,}$, $\beta^{\prime}_{i,}$ and $\beta^{\prime\prime}_{i,}$ indicate the $i$-th row of, respectively, $\mathcal{A}^{{}^{\prime}}$, $\mathcal{A}^{{}^{\prime\prime}}$, $\beta^{\prime}$ and $\beta^{\prime\prime}$.

MathLagr solves the Lagrangian mono-objective problems PRL defined by Equations (32)-(38). In Equation (32), $L$ is the Lagrangian function of the problem.

[TABLE]

[TABLE]

The resulting problem is nonlinear. In order to approach it numerically using mixed-integer program methods, MathLagr employs the subgradient method [30]. This iterative method assigns, in its first iteration, initial values to the Lagrange multipliers, which are usually null values. Then, at each iteration it relies on the resulting best lower bounds found along the iterations and on estimated upper bounds, possibly from heuristics incorporated in the method, to update the Lagrange multipliers in an attempt to refine them. In this paper, $\lambda^{(k)}$ refers to the values of the Lagrange multipliers at iteration $k$. Moreover, for ease of notation, we shall refer to the solutions, that is, the values of the decision variables $u,v,w^{\prime}$ and $w^{\prime\prime}$, as $x_{RL}^{(k)}$, $x_{f}^{(k)}$ or $x^{best}$.

We chose the subgradient method because it theoretically converges to the optimal solution. Even though such convergence is not guaranteed in AugMathLagr due to the linear relaxations used, the subgradient method allows the developed matheuristic to achieve high-quality solutions.

A solution to the Lagrangian problem is not necessarily feasible for the original problem. Therefore we apply a feasibility heuristic to the solution of the Lagrangian problem obtained at each iteration of the subgradient method. Since the feasibility heuristic guarantees that the variables are integer, we solve Lagrangian problems with the integer variables relaxed to speed CPLEX up. Algorithm 1 presents a pseudocode of the proposed Lagrangian heuristic named Mathlagr to solve PRL.

Algorithm 1 uses the following as input: (i) the maximum number of iterations of the method, $kMax$; (ii) the step size to adjust the Lagrange multipliers at each iteration, $st$; and (iii) an initial value for the Lagrange multipliers, $\lambda^{(0)}$. In lines 1 and 2, the upper (UB) and lower bounds (LB) for the objectives are initialized. Then, CPLEX [31] is employed to solve the Lagrangian problem PRL with the integer variables relaxed at iteration [math], $P_{RL}^{(0)}$. In line 6, problem $P_{F}^{(k)}$ is created. It is the MOP with the real variables $X_{maijt},m\in M,a\in A,i,j\in I,t\in T$, which are null in the relaxed solution of the iteration $k$, i.e $x_{RL}^{(k)}$, fixed at zero. In line 7, CPLEX returns a solution for $P_{F}^{(k)}$ and a feasible solution to the MOP, i.e., $x_{f}^{k}$. In line 8, the algorithm updates the upper bound ($UB$) with the best economic objective function value of a solution encountered up to that iteration. In line 9, $x^{best}$ is updated if the value of its economic objective function is better than $UB$. In the first iteration, $UB$ is exactly the value of ${-f_{eco}}(x_{f}^{(0)})$, since it is the first upper bound obtained by the method. The lower bound ($LB$), on the other hand, is the Lagrangian function value of the relaxed solution obtained in that iteration, that is, $L(x_{RL}^{(0)},\lambda^{(0)})$, defined in line 10. In the next iterations, $LB$ can be worse than those found in former iterations. Therefore, $LB$ is the largest lower bound up to that iteration. In line 11, $g(x_{f}^{(k)})$ is defined as the violation of Constraints (33) according to solution $x_{f}^{(k)}$. Then, the Lagrange multipliers are updated as indicated in 12. The process is repeated until either the upper bound is equal to the lower bound or the maximum number of iterations has been reached. Algorithm 1 returns the best feasible solution found over the iterations referred to as $x^{best}$ in line 14.

4.2 Multi-objective matheuristic

Algorithm 2 presents a pseudocode of the AugMathLagr solution method that creates a sequence of mono-objective problems using the multi-objective method AUGMECON2. This algorithm has as input the maximum number of iterations of the Lagrangian heuristic, $kMax$, the step size of the subgradient method, $st$, the initial value of the Lagrange multipliers, $\lambda^{(0)}$, and the number of points of the grid, $dg$.

Mavrotas [37] suggests the use of lexicographic optimization for estimating the lower and upper bounds of each objective function. However, this approach is computationally expensive. Therefore, due to the complexity of the studied model, to define these bounds in line 1 of Algorithm 2, each objective function is individually optimized by applying the Lagrangian heuristic described in Algorithm 1 to the problems $\min\mbox{ }f_{i}(x)+\lambda^{T}({g(x)}),s.t.:$ Constraints $\eqref{eq:monoLagr_ctr}$,$\eqref{eq:monoLagr_ctr2}$ and $\eqref{eq:monoLagr_ctr3}$, for $f_{i}\in\{-f_{eco},f_{env},-f_{soc}$ }. The lower bounds $fL_{i}$ and the upper bounds $fU_{i}$ are the lowest and highest values of each objective function $-f_{eco}$, $f_{env}$ and $-f_{soc}$ and are respectively applied to a relaxed solution and an incumbent solution found along the iterations.

In line 2, the algorithm estimates the ranges of the objective functions by calculating the difference between the lower and upper bounds estimated in line 1. The number of grid points required for each objective function is pre-defined by $dg$. The initial values of $\epsilon_{env}$ and $\epsilon_{soc}$ are set as the estimated upper bound values for the respective objective functions. In line 3, the algorithm calculates the step values for $\epsilon_{env}$ and $\epsilon_{soc}$ according to the number of grid points and the ranges of the corresponding objective functions. In line 4, the algorithm initializes the Pareto frontier approximation $\mathcal{P}$ as an empty set and in line 5 a control variable $gr_{env}$ indicating the point in the grid is initialized as zero.

In line 10, AugMathLagr employs Algorithm 1 to heuristically solve the Lagrangian mono-objective problem MOP, constructed in line 9. If solution $x_{f}$ is feasible, in line 12 the approximation of the Pareto frontier, $\mathcal{P}$, is updated with $x_{f}$. In lines 13 and 14, the algorithm avoids solving redundant problems by checking the ratio between the slack variable of the social constraint and the step size. Otherwise, if $x_{f}$ is infeasible, it means that further decreasing $\epsilon_{soc}$ will only result in infeasible problems. Therefore, the algorithm sets ${gr}_{soc}={dg}$ in line 15, preventing the method from solving these problems and making it proceed to the next value of $\epsilon_{env}$.

In line 20, the dominated solutions are deleted from $\mathcal{P}$ and the updated Pareto frontier approximation $\mathcal{P}$ is returned.

In addition to introducing AugMathLagr, we have also adapted the multi-objective matheuristic AugMathFix [57] to find solutions to the target SSC management problem. The next section briefly explains this adaptation.

4.3 Adaptation of AugMathFix

To better evaluate the proposed matheuristic, AugMathLagr, we compare its performance with an AUGMECON2-based matheuristic recently proposed to approach an SSC problem, the AugMathFix. As the SSC problem is different from the one studied in this paper, we had to adapt such matheuristic, the AugMathLagr, to the target SSC management problem. The adapted AugMathFix selects the same constraints to be relaxed as AugMathLagr. To find feasible solutions to relaxed problems, in addition to imposing which warehouses and entities will not be opened, AugMathFix employs the same strategy as AugMathLagr to fix decision variables.

Details about the implementation of AugMathFix to the supply chain studied in this paper are presented in Section 2 of the Supplementary Material. The mono-objective heuristic in AugMathFix is called MathFix and was also adapted in this paper.

The next section presents the instance generator introduced in this paper.

5 Instance generator

The optimization models found in the literature for sustainable supply chains are usually problem-specific [4] or based on randomly generated parameters. As a consequence, defining a generic instance generator is not an easy task since it should incorporate several traits from this diverse range of characteristics. In this context, this paper introduces a methodology to generate random instances to target the SSC formulation proposed in [41].

In order to define the ranges of the parameters, the data distribution, and so forth, we relied on real data, in particular, from the case study in [41]. This case study is discussed in Section 6.2. Section 3 of the Supplementary Material shows the guidelines for estimating the data parameters of the uniform random variables.

In the next section, we illustrate a small case study created by the instance generator introduced in this paper.

5.1 Illustrative example

This section presents a small example to illustrate how the instance generator works when considering a planning horizon with only two periods. In this example, the SSC has one supplier, one factory, two warehouses, two customers, two airports and two seaports. Moreover, it produces only one final product using two types of raw materials. The example also considers three production and remanufacturing technologies. In our example, we refer to the production technologies as $g_{0},g_{1}$ and $g_{2}$ and to the remanufacturing technologies as $g_{3}$, $g_{4}$ and $g_{5}$. The transportation modes are trucks, airplanes and boats; and there are two types of trucks to perform land transportation, each of them specified by $k_{0}$ and $k_{1}$.

Table 6 shows the results (time to solutions and upper bounds) achieved by the Lagrangian matheuristic whose pseudocode is in Algorithm 1 - MathLagr when optimizing each objective function of the MOP individually. CPLEX [31] was the tool that solved the Lagrangian relaxed problems with a stopping criterion of $1\%$ of optimality gap. The first column of Table 6 identifies the optimized function. The third, fourth and fifth columns report, respectively, the values of the economic, environmental and social functions of the heuristic solution to the mono-objective problem whose objective function is that indicated in the first column.

Figures 2 to 4 display the decision variable values regarding each time period of the planning horizon obtained by MathLagr when optimizing the economic, environmental and social functions, respectively111The values were rounded for readability purposes.. These figures illustrate the flow of items indicated by the values of the variables of production, remanufacturing and storage on the arcs of the network. Labels $S$, $F$, $W_{1}$, $W_{2}$, $C_{1}$, $C_{2}$, $Air_{1}$, $Air_{2}$, $Sea_{1}$ and $Sea_{2}$ represent, respectively, the supplier, the factory, the first warehouse, the second warehouse, the first customer, the second customer, the first airport, the second airport, the first seaport and the second seaport.

Moreover, to identify these entities in the indices of the variables, we assign integer numbers to them. Entity $S$ corresponds to index [math], $F$ to $1$, and so on until $Sea_{2}$, which is identified by index $9$ in the variables. For the same reason, technologies $g_{0}$, $g_{1}$, $g_{2}$, $g_{3}$, $g_{4}$ and $g_{5}$ are assigned to sequential indices [math] to $5$ and trucks $k_{0}$ and $k_{1}$ to indices [math] and $1$. The labels on the arcs correspond to the decision variables, according to Table 3. In this example, warehouse $W_{2}$, customer $C_{2}$, airport $Air_{2}$ and $Port_{2}$ are located in a continent different from the remaining entities. Thereby to transport items between entities, in this case, air or sea transportation is mandatory.

No warehouse is opened in the solutions that optimize the economic function, i.e., maximizes the NPV, and minimizes the environmental impact using the ReCiPe 2008 methodology. The solution that minimizes the environmental function indicates that avoiding to install warehouses contributes to reducing the environmental impacts. Since the decisions to open warehouses contribute positively to the social objective function, the solution to the maximization of this indicator indicates the opening of the two warehouses at maximum capacity.

Truck $k_{1}$ is selected for land transportation in all the entities in the environmental solution, whereas truck $k_{0}$ is preferable in the economic solution. In particular, overall, truck $k_{1}$ has a higher price and a lower environmental impact than truck $k_{0}$. In the environmental optimization solution, apart from land transportation, only seaports are selected due to the lower environmental impact they cause. In the economic, environmental and in the second period of the social optimization solutions, seaports are responsible for delivering products to customer $C_{2}$ and returning the used products to the factory. In the first period of the solution that optimizes the social function, in contrast, airports are responsible for delivering products to customer $C_{2}$. The stock in warehouse $W_{2}$ at the end of the first time period in the social optimization solution is used to meet the demand by customer $C_{2}$ in the second time period.

The environmental and social solutions employ the manufacturing technology $g_{0}$, whereas the economic solution uses the cheaper technology $g_{1}$. The remanufacturing levels at factories are higher in the environmental solution, suggesting that reverse logistics can reduce the environmental impacts of the SSC. One of the advantages of recycling is to reduce the negative environmental impacts due to raw material acquisition from suppliers.

The next section presents the computational experiments carried out with a real case study and artificially generated instances.

6 Computational experiments

In this section, we shall discuss two experiments carried out using a case study and artificial instances drawn from the instance generator presented in this paper. The results obtained by the proposed matheuristic, AugMathLagr, are contrasted with two multi-objective methods found in the literature: AUGMECON2 and an adaptation of AugMathFix. Note that we refer to the heuristics employed in AugMathLagr and in AugMathFix to solve mono-objective problems by, respectively, MathLagr and MathFix. The first experiment concerns the case study presented by Mota et al. [41]. The second presents an analysis of the computational results achieved by the methods for a set of benchmark instances.

The number of grid points was set at 10 for all the methods. The AugMathLagr parameters were: $kMax=10$; $\lambda^{(0)}=0$; and $st=1\times 10^{-10}$, $st=1\times 10^{-10}$ and $st=1\times 10^{-8}$ when optimizing the economic, environmental and social functions, respectively. All experiments were performed on a computer with an Intel Xeon E5-2680v2 2.8 GHz processor and 128 GB of main memory. CPLEX was limited to use only 15GB of memory in all the experiments. Moreover, CPLEX is limited to use 8 threads in both experiments. A maximum time limit of $1h$ was imposed for CPLEX to solve each mono-objective problem, MOP. Exceptionally for AUGMECON in the experiment with the case study, which was computationally more challenging, we imposed a time limit of 3h for CPLEX. Very large and small absolute values are expressed in E notation.

Before detailing the results of the experiments, we shall describe the measures used to evaluate the multi-objective optimization methods.

6.1 Measures of assessment

$GAP$, a measure of proximity, was used to evaluate the quality of solutions with respect to their optimal value. $GAP$ quantifies how close a solution $S^{{}^{\prime}}$ is from another $S^{*}$, considering an objective function $f_{i}$ to be minimized, i.e. $f_{i}\in\{-f_{eco},f_{env},-f_{soc}$} indicated in Equation (39).

[TABLE]

The stopping criterion for CPLEX was when it reached a solution whose $GAP$ between the upper and lower bounds was less than or equal to $1\%$. The remaining parameters were set as the default of the solver.

To better assess the results achieved by MathLagr, similar to Tautenhain et al. [57], we used a set of multi-objective metrics to evaluate the quality of the Pareto frontiers obtained by the methods. They were: the R2 Indicator and two variations of the Mean Ideal Distance (MID) and Spread of Non-Dominated Solutions (SNS).

The R2 Indicator [28] checks the quality of a Pareto frontier approximation, $P^{A}$, by comparing it with a representative set of Pareto frontier, $P^{Z}$, as shown by Equation (40). In the experiments reported in this paper, $P^{Z}$ was estimated by AUGMECON2.

[TABLE]

where $U$ is the set of utility functions. A utility function ${\mu}\in U$, ${\mu}:\mathds{R}^{3}\rightarrow\mathds{R}$, maps a solution of the multi-objective problem to a scalar value. Lower values of $I_{R2}(P^{A},U)$ indicate better approximations. Negative values of the R2 Indicator express that solutions from $P^{A}$ are closer to the Ideal point than solutions from $P^{Z}$. Let the Ideal point be defined by $f^{I}=\{f^{*}_{eco},f^{*}_{env},f^{*}_{soc}\}$, such that $f^{*}_{i}$ is the optimal value of the problem that optimizes objective function $f_{i}\in\{-f_{eco},f_{env},-f_{soc}$}.

Hansen & Jaszkiewicz [28] suggest using the Weighted Sum and weighted Tchebycheff utility functions. The Weighted Sum function only takes into account points inside the convex hull of a feasible region and therefore is unsuitable for cases where the solution space is not convex. Thereby, Brockhoff et al. [8] suggest using the weighted Tchebycheff function, formulated as presented in Equation (41).

[TABLE]

where $\gamma_{eco},\gamma_{env}$ and $\gamma_{soc}$ are weights associated with $f_{eco},f_{env}$ and $f_{soc}$, respectively.

For ease of notation, we refer to $I_{R2}(P^{A})$ as $I_{R2}(P^{A},U)$, with $U$ representing the Tchebycheff utility function.

The Mean Ideal Distance (MID) and Spread of Non-Dominated Solutions (SNS) [5] calculate, respectively, the mean and the standard deviation of the distance between the solutions of an approximation of the Pareto frontier, $P^{A}$, and the Ideal point. Due to the different scales of the objective function values, Tautenhain et al. [57] employed the adaptation of these measures to calculate the multidimensional $GAP$ ($GAPM$) instead of the distance between the solutions. Equation (42) presents how to calculate $GAPM$, where $GAP$ is defined according to Equation (39).

[TABLE]

where $f(S^{{}^{\prime}})=\{-f_{eco}(S^{{}^{\prime}}),f_{env}(S^{{}^{\prime}}),-f_{soc}(S^{{}^{\prime}})\}$.

Equations (43) and (44) present the adapted $MID$ and $SNS$ metrics, here called $aMID$ and $aSNS$, for $P^{A}$.

[TABLE]

[TABLE]

Lower values of $aMID(P^{A})$ and $aSNS(P^{A})$ indicate better solutions in $P^{A}$.

6.2 Experiment I: case study

In the first experiment, a case study of an electronics company is investigated (see the supplementary report on Mendeley Data [58] or [41] for additional details on this case study). The experiment was carried out with the original problem of [41], which did not consider suppliers’ environmental impacts and working conditions. It is worth pointing out that Mota et al. [41] selected the set of locations that could be used to install the facilities beforehand, which guaranteed adequate working conditions for the employees.

Table 7 only shows details regarding entities in this case study – location, country GDP per capita where an entity is located (GDPPC) [41], measured in Purchasing Power Parity (PPP)222PPP is an exchange rate that translates the purchasing power of different monetary units., and costs in euros related to construction and labor – in order to contrast such information with the results obtained by the methods. The first, second, third and fourth columns respectively indicate the entities, their labels, geographical locations and regions. The fifth, sixth and seventh columns respectively report the values of GDPPC and construction and labor costs of entities.

The supply chain entities are located in several countries in Europe and Brazil. The company has its factory $F_{1}$ and warehouse $W_{1}$ installed in Verona. Since the suppliers and factories are located in Europe, the company must use either air or sea transportation to carry the final products to Brazil.

The company produces only two types of final products from four different raw materials. The returned final products are modeled as recovered products. There are also four technologies available for factories to manufacture products from raw materials and two technologies for recovering used products into final products.

Table 8 depicts the results of the exact method employed in AUGMECON2, of MathLagr and of MathFix in the optimization of each objective function of the problem.

The third, fourth and fifth columns indicate the results obtained by optimizing the economic, environmental and social functions. For each method, the results reported are the computational running time to find the solution and the values of economic, $f_{eco}$, environmental, $f_{env}$, and social, $f_{soc}$, functions. To calculate the $GAPs$ for MathFix and MathLagr according to Equation (39), we assume the solutions achieved by the exact method to be the optimal ones. Moreover, the rows identified by “Speed up” report the ratio of the exact method execution time and the matheuristics execution time. All the methods estimate the ideal and nadir points using, respectively, the best and worst values of each objective solution.

MathLagr is approximately $3.5$ and $25$ times faster than the exact method for optimizing the economic and social objective functions, respectively. Despite being slower than MathLagr, MathFix was still $3$ and $24$ times faster than the exact method for optimizing the economic and social objective functions, respectively. AUGMECON was faster than MathLagr and MathFix to optimize the environmental objective function.

MathLagr and MathFix obtained solutions for the economic and environmental objective functions whose $GAP$s to the corresponding solutions obtained by the exact method were lower than $0.041\%$ and $0.230\%$, respectively. On the one hand, MathLagr was faster and obtained a smaller $GAP$ than MathFix regarding the economic and environmental function optimization solution. On the other, MathFix obtained a solution for the social objective function whose $GAP$ to the solution obtained by the exact method was only $0.657\%$.

Table 9 presents the numbers of non-dominated solutions found by the three multi-objective methods AUGMECON2, AugMathLagr and AugMathFix as well as the $aMID$, $aSNS$, and R2 Indicator values. It also reports the total computational running times in seconds that the methods took to approximate the Pareto frontier for the case study.

Table 9 shows that AugMathLagr and AugMathFix obtained solutions with lower $aMID$ and $aSNS$ values than AUGMECON2. In addition, the R2 indicator of the solutions found by AugMathLagr and AugMathFix were low. These results attest to the good quality of the Pareto frontier approximations found by the methods when compared to AUGMECON2.

On the one hand, AugMathLagr is approximately $57$ and $1.8$ times faster than AUGMECON2 and AugMathFix, respectively, to estimate the Pareto frontier. On the other, both AUGMECON2 and AugMathLagr find more non-dominated solutions than AugMathFix.

Figure 5 exhibits projections of the objective function values of the solutions in the Pareto frontier approximations achieved by AUGMECON2, AugMathLagr and AugMathFix into two dimensional spaces.

Despite AugMathLagr and AugMathFix finding very close objective function values, Figures 5(a) and 5(c) show that AugMathFix did not obtain some solutions achieved by AugMathLagr. In particular, AugMathLagr found a non-dominated solution whose $f_{env}$ value is 9.4e+08, highlighted in yellow circles in Figures 5(a) and 5(c). No non-dominated solution close to such a solution was obtained by AugMathFix.

In Figures 5(a) and 5(c), it is not possible to notice a clear trade-off between the environmental and the other objectives. Figure 5(b), however, more clearly illustrates the trade-off between the economic ($f_{eco}$) and social ($f_{soc}$) objectives, since the higher the economic function, the lower the social value. Thereby, for a given environmental objective, the higher the $f_{eco}$, the lower the $f_{soc}$. Besides, in Figures 5(a) and 5(c), there are some solutions with the same $f_{soc}$ whose $f_{eco}$ values, despite being very close, ensure that the solutions are non-dominant.

In this context, Figure 5 graphically demonstrates the good quality of the Pareto frontier approximations found by the methods, since it is possible to visualize that the non-dominated solutions are well distributed in the hyperplane. Multi-objective optimization can add value to the supply chain by allowing decision makers to select the most suitable solutions, according to their goals. When presenting solutions to decision makers, the improvement and worsening of each objective should be taken into account. For example, a slight decrease in the economic objective might represent a significant improvement in the environmental objective.

Table 10 shows which factories and warehouses are opened in the solutions found by MathLagr, MathFix and the exact method in AUGMECON2 to optimize the economic, environmental and social criteria. Additionally, in this table we point out the entities installed at maximum and minimum capacities. The marks $A,L$ and $M$ indicate, respectively, the exact method, MathLagr, and MathFix.

It is possible to observe in Table 10 that there is no difference in the selection of factories to be opened at maximum or minimum capacity in the solutions found by any of the methods. Both the exact method and MathLagr opened warehouse $W_{9}$ in addition to warehouse $W_{1}$, which is opened at a fixed capacity in every solution, in the solutions that optimize the economic function. MathFix, on the other hand, selected warehouse $W_{5}$ besides $W_{1}$.

The solutions obtained by the three methods for the environmental minimization problem open more warehouses than the economic objective maximization solutions.

When minimizing the environmental objective function, the most relevant decisions for the environmental impact are those related to manufacturing and remanufacturing. Such decisions correspond to approximately $99.94\%$ of the environmental impacts on the solution that minimizes the environmental objective function. Transportation and entity location decisions represent only $0.4\%$ and $0.2\%$ of the overall environmental impacts on such a solution, respectively.

By comparing the environmental impacts due to the installation of additional warehouses on the solution of the mono-objective problem that minimizes the environmental objective function, they were 43.43% greater than such impacts on the solution that maximizes the economic objective. Nonetheless, the environmental impacts due to manufacturing and remanufacturing on the environmental minimization solution are approximately $7\%$ lower than such impacts on the economic optimization solution. Note that, as mentioned earlier, manufacturing and remanufacturing represent the overwhelming majority of the environmental impacts. As a consequence, this difference in the overall solution is approximately $7\%$.

From these results, we can see that the model and the methods can properly identify the manufacturing and remanufacturing activities as the most representative activities for the environmental impact. Thereby, it is recommended for the company to perform a careful analysis and estimations on the parameters that affect both the manufacturing and remanufacturing technology costs and environmental impacts. We also recommend performing a sensitivity analysis on these parameters when dealing with different industrial segments.

Except for warehouse $W_{1}$, required by the case study to be opened at a fixed capacity in every solution, all the remaining warehouses are opened at maximum capacity in the solutions obtained by the exact method through the optimization of the social function. MathLagr and MathFix opened, respectively, 7 and 8 out of the 9 available warehouses when optimizing the social function.

The social solution benefits from more entities and transportation. We suggest, thereby, to analyze in detail choices associated with solutions that benefit the social criterion. Moreover, companies might have strategic goals or preferences in selecting certain entities. Therefore, it is advisable to take these issues into account when selecting the best trade-off solution involving social optimization.

All the solutions obtained by the methods to individually optimize the objective functions use the two seaports. In every social optimization solution, all airports are operational, because they positively contribute to creating jobs. Because air transportation is more expensive than road transportation, the airports are not employed in the solution that optimizes the economic objective function. In this case study, air transportation causes lower environmental impacts than one of the trucks used in road transportation [41]. Thereby, it is preferred when optimizing the environmental function. Airports $Air_{1}$ and $Air_{3}$ are operational in all the solutions that optimize this function. In the solution found by MathLagr, airport $Air_{2}$ operates as well.

Even though the manufacturing and remanufacturing represent the overwhelming majority of the environmental impacts, all the methods still prefer air transportation due to its lower environmental impacts.

6.3 Experiment II: artificial instances

In this experiment, we created artificial instances through the introduced instance generator. For this purpose, we considered instances with a predefined number of entities, items and technologies. We generated four types of instances, differentiated by the prefix in their labels “STD”, “TECHC”, “RAWC”, “SUP” and “CAP”. The instances prefixed by “STD” follow the same patterns of the model introduced in [41]. In the instances prefixed by “TECHC”, the technology acquisition costs were the same for all technologies. The manufacturing costs were the same for all the suppliers in the instances prefixed by “RAWC”. In the instances prefixed by “SUP”, no minimum order was imposed for the suppliers. There was no minimum use of technologies in the instances prefixed by “CAP”. The corresponding remaining parameters of the instances prefixed by “STD”, “TECHC”, “RAWC”, “SUP” and “CAP” were generated as described in Section 3 of the Supplementary Material.

The generated set consisted of $60$ artificial instances, whose primary characteristics are summarized in Table 11. The first and second columns of this table indicate the prefixes and suffixes for the names of the instances, respectively. The name of an instance is given by its prefix followed by the suffix. Columns “$|I|$”, “$|M_{rm}|$”, “$|M_{fp}|$”, “$|M_{rp}|$” and “$T$” express, respectively, the number of entities, raw materials, final products and periods of the instances.

Instances with 17 entities have 3 suppliers, 3 factories, 3 warehouses, 4 customers, 2 airports and 2 seaports. Instances with 25 entities have 3 suppliers, 3 factories, 9 warehouses, 4 customers, 4 airports and 2 seaports. The number of production and remanufacturing technologies was 3, totaling 6 technologies. For ease of identification, we included in the name of the instances information regarding the number of entities, items and periods in the planning horizon and specific patterns for the parameters, as described in the following paragraph. The supplementary report in [58] presents the values chosen for the parameters to generate such instances discussed in Section 3 of the Supplementary Material.

Let $P_{AUG2}$, $P_{MathLagr}$ and $P_{MathFix}$ be, respectively, the Pareto frontier approximations obtained by AUGMECON2, AugMathLagr and AugMathFix. Figures 6, 7 and 8 present, respectively, the $aMID$ measure, the $aSNS$ measure and the R2 Indicator for $P_{AUG2}$, $P_{MathLagr}$ and $P_{MathFix}$. Figure 9 presents the total computational running times for obtaining $P_{AUG2}$, $P_{MathLagr}$ and $P_{MathFix}$.

As can be noted in Figure 6, the $aMID$ measures for the solutions of AugMathLagr and AugMathFix are lower than or approximately the same as the solutions achieved by AUGMECON2 for $42$ and $32$, respectively, out of the $60$ instances. Note that these values correspond to $70\%$ and approximately $53\%$ respectively of the instances considered. Although the $aMID$ values of $P_{MathLagr}$ are significantly high in 4 other instances, their overall results indicate that the solutions provided by AugMathLagr present low $GAPM$s in comparison to the Ideal point for most of the instances.

The standard deviations of $aMID$, measured by the $aSNS$ measure, for solutions obtained by AugMathLagr and AugMathFix were considerably larger than those achieved by AUGMECON2, which were significantly high for 6 and 10 instances, respectively. AUGMECON2, in particular, presented unexpectedly large $aSNS$ values for 4 other instances.

The R2 indicator value for the solutions obtained by AUGMECON2 is null, since we considered $P_{AUG2}$ as the set of representative solutions for the measure. The maximum distance from a solution obtained by AugMathLagr and AugMathFix to the Ideal point is close to the maximum distance from a solution found by AUGMECON2 to the Ideal point.

The total computational running times required by AUGMECON2 were at least $11$ and $6$ times higher than those required by AugMathLagr and AugMathFix in instances with at least 6 periods in the planning horizon. Concerning the running times, the advantages of AugMathLagr and AugMathFix over AUGMECON2 are even more evident when considering the largest instances with 25 entities, i.e. $|I|=25$. AugMathLagr, in particular, was from $18$ to $550$ times faster than AUGMECON2 and from $0.8$ to $3.7$ times faster than AugMathFix in instances with $25$ entities and at least $6$ periods in the planning horizon.

The results of this experiment attest to the good quality of the Pareto frontier approximations obtained by AugMathLagr. Moreover, the adaptation of AugMathFix presented significantly inferior running times in comparison to AUGMECON2. AugMathLagr was even faster than AugMathFix in 81.667% of the tested instances. This enables us to infer that the proposed Lagrangian matheuristic is suitable for solving large instances.

In the supplementary report in [58], we present experiments with instances whose demands have increasing, decreasing and alternating trends along the time periods. The results of the experiments attest that the proposed method is robust in finding good quality solutions for these scenarios. Therefore, we can infer that the proposed method is suitable for dealing with supply chains whose demands vary in the planning horizon.

7 Final remarks

SSC management has recently drawn the attention of industrial and academic sectors to new research developments. Most studies approach SSC management by multi-objective optimization, focusing on the construction of optimization models based on case studies. This approach, however, is commonly very time-consuming and, even though it belongs to the strategic and tactical levels of planning, to refine, validate and thoroughly study the supply chain, it is necessary to employ efficient solution methods [17]. Perhaps due to the difficulty in solving these problems and because the studies are usually problem-oriented, little effort has been devoted to keeping a benchmark data repository. This paper introduces an artificial instance generator and an efficient Lagrangian matheuristic, here called AugMathLagr, for SSC management multi-objective problems. In particular, AugMathLagr is implemented to the optimization model introduced by Mota et al. [41]. In addition to introducing AugMathLagr, this paper also efficiently adapts the matheuristic AugMathFix proposed by Tautenhain et al. [57] to the target problem.

Experiments conducted with a case study found in the literature and with a test bed of artificial instances indicated that in comparison with a classical exact method for multi-objective problems, known as AUGMECON2, AugMathLagr and AugMathFix were significantly faster. In particular, AugMathLagr was faster than AugMathFix and performed even better when approximating the Pareto frontier in cases where larger instances were considered. We assessed the Pareto frontier approximations using three multi-objective measures to evaluate how close the solutions from the Pareto frontier approximation were from the Ideal point. These results indicated that the solutions of the Pareto frontier approximation achieved by AugMathLagr and AugMathFix were close to the Ideal point, in comparison with the solutions from AUGMECON2. In particular, the values of the R2 indicator expressed the high quality of the solutions found by AugMathLagr and AugMathFix.

In the case we studied in this paper, the solution methods identified the most representative decisions for environmental and social objectives. On the one hand, the social criterion benefits from more entities being opened. On the other, manufacturing and remanufacturing decisions represent most of the environmental impacts. We acknowledge, however, that this result might be biased towards the fact that our case study was based on an electronics component producer, for which manufacturing and remanufacturing activities are higher than for other industries.

Nonetheless, the results related to the environmental impacts are strongly dependent on the scope of the investigation and on the case under study. Since data associated with environmental and social objectives are particularly challenging to be accessed, companies should carefully study the parameters associated with the most representative decision in order to better meet the reality. In particular, data available for similar supply chains should be employed as a starting point. Additionally, sensitivity analysis can also be conducted on these parameters to ensure that the obtained supply chain results and conclusions remain unaltered.

This paper also contributes with a framework for generating artificial instances. Artificial instances are particularly useful for assessing the performance of optimization models since they allow statistical inferences on a wide variety of numerical scenarios. Due to the specificity of each optimization model, the implementation of the instance generator to other SSC management optimization models is a suggestion for future work.

Another direction for future work regards the proposed matheuristic. The primary goal of such study would be to further improve the quality of the solutions in order to obtain even closer values to the optimum ones achieved by AUGMECON2 using local search procedures.

Acknowledgments

Authors Camila P. S. Tautenhain and Mariá C. V. Nascimento would like to acknowledge the funding granted by São Paulo Research Foundation (FAPESP), grant numbers: 14/27334-9, 15/21660-4 and 16/02203-4; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), grant numbers: 306036/2018-5; and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brazil (CAPES) – Finance Code 001. The authors Ana Paula Barbosa-Povoa and Bruna Mota would like to acknowledge the financial support from FCT and Portugal 2020 FCT under the project PTDC/EGEOGE/28071/2017, Lisboa -01.0145-Feder-28071. The author Mariá C.V. Nascimento is also grateful to Leonardo V. Rosset for giving her a hand. Research carried out using the computational resources of the Center for Mathematical Sciences Applied to Industry (CeMEAI) funded by FAPESP (grant 2013/07375-0).