Environmental aspects in supplier networks-a bi-objective just-in-time truck routing problem

Freight transportation, including just-in-time (JIT) supplier networks, accounts for a substantial part of the global carbon dioxide (CO $ _2 $ 2) emissions. The JIT truck routing problem (TRP-JIT) presented in the recent literature consists of several suppliers serving a single original equipment manufacturer (OEM). A logistics provider organises the milk-run routes. The shipments are available after their release dates at the suppliers and should be delivered on their due dates at the OEM with minimal total earliness-tardiness penalties (first objective). Unlike previous research on the TRP-JIT, we focus on its environmental impact: (1) We include the weight-distance (second objective), depending on the truck's curb weight, the load, and the transportation distance. (2) We adapt a state-of-the-art large neighbourhood search (LNS) from the literature considering both objectives. (3) The LNS is embedded in bi-criterial frameworks, i.e. ε-constraint and weighted sum methods. Thereby, we estimate Pareto frontiers with at least 60 solutions in less than 25 min for instances with 99 shipments. From a managerial perspective, increasing the difference between the release and due dates for a better JIT performance may worsen the environmental impact. Lighter trucks can reduce the environmental costs without affecting the JIT performance, whereas a smaller fleet negatively affects both objectives.


Introduction
The reduction of carbon dioxide (CO 2 ) emissions is a necessary step towards tackling climate change.Consequently, measures for this purpose are the most effective if they lead to a strong decline.According to Rüdiger, Schön, and Dobers (2016), logistics and especially freight transportation cause around 5% of the global CO 2 emissions.Therefore, it is crucial to consider logistics, especially transportation, regarding the environmental impact.In this article, we focus on just-in-time (JIT) routing in supplier networks from an environmental perspective.Solutions which are favourable for JIT production may not be preferable with regard to environmental aspects, especially the CO 2 emissions.In JIT production parts are required at the assembly line at a certain time (Boysen et al. 2015).The focus of JIT routing is to organise the transportation from suppliers to an original equipment manufacturer (OEM) to feed the production accordingly.Too early or late delivery causes the highest costs in JIT production, whereas the environmental impact mainly depends on factors like the CONTACT Julian Baals julian.baals@uni-jena.deChair of Business Informatics, esp.Business Intelligence, Faculty of Economics, Friedrich Schiller University Jena, Carl-Zeiss-Str.3, Jena D-07743, Germany Supplemental data for this article can be accessed online at https://doi.org/10.1080/00207543.2023.2258237.routes' length and the weight transported.Consequently, these two objectives may not be compatible, especially if orders from geographically dispersed suppliers are required simultaneously.
The environmental impact of vehicle routing problems (VRP), initially defined by Dantzig and Ramser (1959), is discussed in recent literature, and often called Green VRP or Pollution Routing Problem (PRP), e.g. by Bektaş and Laporte (2011) and Bektaş, Demir, and Laporte (2016).However, Baals, Emde, and Turkensteen (2023) have shown that classic VRP models do not accurately capture the reality of JIT routing in supplier networks.Nonetheless, we can use the green VRP literature as a starting point for investigating the environmental impact of JIT routing.We estimate the CO 2 emissions of the transportation based on the weight-distance criterion: The weight-distance depends on the weight that is transported along the routes, including the trucks' curb weight, and the routes' length.Following the literature, we assume that (1) the weight-distance is proportional to the fuel consumption, and (2) the fuel consumption is directly related to the CO 2 emissions (Kara, Kara, and Yetis 2007;Kopfer, Schönberger, and Kopfer 2014;Y. Xiao et al. 2012).
The JIT truck routing problem (TRP-JIT), introduced by Baals, Emde, and Turkensteen (2023), is evaluated from an environmental perspective in this article.We call this the environmental TRP-JIT (envTRP-JIT).The TRP-JIT consists of a set of shipments located at different suppliers that must be delivered by third-party logistics (3PL) provider with a homogeneous truck fleet to a single OEM.Hence, it follows a many-to-one structure, which makes it different from problems in classical routing literature.The pickup at the suppliers must take place after a specific release date and the delivery at the OEM at a certain due date.The release dates impose strict constraints in the problem because the shipments cannot be picked up before their production is finished.Contrarily, we can deviate from the due date, which causes earliness or tardiness costs.These costs refer to either storage costs or additional costs for production scheduling at the OEM.Allowing deviations from the due dates enables the 3PL to plan routes, where several shipments with different due dates are consolidated on a so-called milk-run in a single truck.The objective of the TRP-JIT is to minimise the earliness-tardiness costs caused by all shipments (Baals, Emde, and Turkensteen 2023).
In this article, we consider the case, where a 3PL aims to evaluate and eventually improve its environmental impact.The 3PL requires knowledge about the JIT performance, when shifting its focus towards environmental goals.This arises two important questions for the 3PL: (1) What is the relationship between the JIT performance and the environmental impact in general?In which sense do these objectives conflict or comply?(2) How can we achieve a certain improvement of the environmental impact?What are the costs regarding the JIT performance?
Hence, we evaluate the interplay and trade-off between the earliness-tardiness and the weight-distance objective in the envTRP-JIT by applying multi-objective optimisation methods, i.e.ε-constraint methods and weighted sum methods (Laumanns, Thiele, and Zitzler 2006;Przybylski, Gandibleux, and Ehrgott 2010).Therefore, we extend the model and solution procedure of Baals, Emde, and Turkensteen (2023): First, we adapt the objective and constraints in the mixed-integer linear programming model (MIP) of the TRP-JIT for solving the envTRP-JIT optimally, including the weight-distance objective.Exact methods are only practical for small instances of the envTRP-JIT.Consequently, we extend the metaheuristic by Baals, Emde, and Turkensteen (2023), which combines a Large Neighbourhood Search (LNS) metaheuristic with a local search, by a weight-distance objective.We finally embed the MIP model using a commercial solver and the LNS solution procedure in ε-constraint and weighted sum method frameworks.In essence, the contributions of this article are: • We extend the TRP-JIT to include both JIT and environmental objectives.• We extend the LNS metaheuristic framework of Baals, Emde, and Turkensteen (2023) without increasing the computational complexity, despite the bicriterial nature of the problem.• We present an extensive numerical study to quantify the trade-off between the JIT perspective and environmental impact.We derive managerial insights to evaluate the design of a JIT logistics system according to environmental goals in terms of the planned time span between release and due dates and fleet composition.
The outline of this article is: We present the relevant literature to evaluate environmental aspects of the TRP-JIT in Section 2.Then, we extend the problem formulation of the TRP-JIT, including the MIP in Section 3. The required adaptions to the LNS metaheuristic and the bi-objective frameworks, i.e.the ε-constraint method and weighted sum method frameworks, are presented in Section 4. The computational study, including its design and results, is presented in Section 5.It focuses on the implementation of the LNS within the bi-objective frameworks and managerial insights.We summarise the outcomes and conclusions of the article in Section 6.

Literature review
The original TRP-JIT can be expressed and modelled by adapting well-studied variants of the VRP.Baals, Emde, and Turkensteen (2023) follow this approach and propose an LNS-based metaheuristic to solve the TRP-JIT.Toth and Vigo (2014) provide an overview of different variants of the VRP, i.e.VRPs with pickup and deliveries and (soft) time windows.Bräysy andGendreau (2005a, 2005b) focus on heuristic solution procedures, including local search for the VRP with time windows.Shaw (1998) provides the fundamental framework for the LNS metaheuristic.For further references related to the TRP-JIT and where it differs from common routing problems, we refer to the introduction and literature review by Baals, Emde, and Turkensteen (2023).
Additional to the economic goals, we need to define objectives if we aim to minimise the environmental impact of transportation or, more generally, supply chains.Approaches to quantify environmental aspects of transportation have been especially studied by Blanco and Sheffi (2017), Demir, Bektaş, andLaporte (2011), Demir, Bektaş, andLaporte (2014b), Turkensteen (2017), andHeni et al. (2023).Blanco and Sheffi (2017) introduce different measurements for the emissions of greenhouse gases, such as CO 2 , and propose different metrics dependent on the mode of transportation.Furthermore, Blanco and Sheffi (2017) distinguish between emissions, having an impact globally on the climate by global warming, and pollution, which affects the environment surrounding the source of emission.Demir, Bektaş, and Laporte (2011) present emission models focusing on road freight transportation.A review of recent literature on environmental aspects of road freight transportation, including different metrics for fuel consumption and emissions, is provided by Demir, Bektaş, and Laporte (2014b).We can state that the emission metrics in literature are mainly based on the vehicle's weight, speed, and distance travelled.Dekker, Bloemhof, and Mallidis (2012) provide a general literature review on green logistics as an aspect of operations research without explicitly mentioning emission metrics.Besides road freight transportation, Bauer, Bektaş, and Crainic (2010) consider emissions in rail service design as an application in intermodal freight transport.
Environmental aspects have also been considered certainly in vehicle routing.The most relevant literature is summarised in Table 1.Kara, Kara, and Yetis (2007) introduce the energy minimising VRP, following a weight-distance objective to relate to the fuel consumption.Formally, this means minimising the product of the vehicles' total weights, and the distance of the route segments this weight applies to.According to Kirby et al. (2000), the fuel consumption directly relates to CO 2 emissions.Additionally, Kopfer, Schönberger, and Kopfer (2014) elaborate on greenhouse gas reduction based on a weight-distance objective and assuming a heterogeneous vehicle fleet.
The vehicle's speed is neglected in the models by Kara, Kara, and Yetis (2007) and Kopfer, Schönberger, and Kopfer (2014).Contrarily, the impact of the vehicle's speed is studied by Fagerholt, Laporte, and Norstad (2010) for maritime shipping, including speed as a decision variable.Maden, Eglese, and Black (2010) and Jabali, Van Woensel, and de Kok (2012) evaluate the effect on the CO 2 emissions, and consider time-varying VRPs, where the travel speed and time vary due to temporary congestion on the road.Liu et al. (2023) consider the effects of velocity in a time-dependent green VRP with time windows.They apply an adaptive LNS.Figliozzi (2010) includes speed as a decision variable in the so-called emissions VRP, which also includes time windows.Figliozzi (2010) refers to Barth and Boriboonsomsin (2008) to show the link between the average travel speed and the CO 2 emissions.Similar to Turkensteen (2017), they state that the CO 2 emissions are rather stable for a broad range of speed levels.Y. Xiao et al. (2012) study models to optimise fuel consumption in VRPs.They find a linear relationship between the fuel consumption rate and the vehicle's weight.Heni et al. (2023) apply supervised learning to estimate the fuel consumption in vehicle routing.Bektaş and Laporte (2011) introduce the so-called pollution routing problem (PRP), which aims for the minimisation of the drivers' wages depending on the total travel time along the routes and the costs for emissions using a weight-distance criterion and a speed level decision.Demir, Bektaş, and Laporte (2012) present an adaptive LNS metaheuristic for the PRP.Franceschetti et al. (2013) consider the time-dependent PRP to account for traffic congestion.An overview of different models based on the PRP by Bektaş and Laporte (2011) is provided by Bektaş, Demir, and Laporte (2016).Y. Y. Xiao and Konak (2015) define the CO 2 emissions, total arrival time, and total travel time as hierarchical objectives, and the proposed solution procedure is a simulated annealing algorithm.Zhu and Hu (2019) focus on congestion and emission factors in a time-dependent VRP.
The proposed solution procedure is a combination of a genetic algorithm and particle swarm optimisation.Gupta et al. (2022) apply a genetic algorithm to solve a green VRP with fuzzy time-distances and demands split into bags.
Besides single-objective approaches or monetisation, i.e. defining a price for an emission unit like Bektaş and Laporte (2011), multi-objective methods can be applied to evaluate greenhouse gas emissions in VRPs.These methods are ε-constraint, and weighted sum methods.ε-constraint methods are based on the principle of choosing only one of the multiple objectives for optimisation and using the other objectives as constraints.Efficient implementations of these methods are studied by Ehrgott (2005a), Laumanns, Thiele, and Zitzler (2006), Mavrotas (2009), and Mavrotas and Florios (2013).In weighted sum methods, variations of the weights for the different objectives are applied to find non-dominated solutions.For more details, we refer to Cohon (1978), Przybylski, Gandibleux, and Ehrgott (2010), Kim and de Weck (2005), and Ehrgott (2005b).Additionally, Zitzler et al. (2003), Zitzler and Thiele (1998a), Zitzler andThiele (1998b), andAudet et al. (2021) study indicators to evaluate the performance of different multi-objective methods.
Demir, Bektaş, and Laporte (2014a) consider the PRP introduced by Bektaş and Laporte (2011) applying multiobjective approaches to estimate the Pareto front, i.e. different weighted sum methods, an ε-constraint method, and a combination of both, called hybrid method.The adaptive LNS metaheuristic proposed by Demir, Bektaş, and Laporte (2012) is applied to solve larger instances.As objectives, Demir, Bektaş, and Laporte (2014a) consider the travel time determining the drivers' wages on one hand, and the fuel costs directly related to the CO 2 emissions on the other hand.Costa et al. (2018) consider the same problem setup.Their main contribution is the so-called two-phase Pareto local search heuristic.Erdoğdu and Karabulut (2022) elaborate on the green vehicle routing problem by following a distance objective for the travel costs and a weight-distance objective for CO 2 emissions.Although considering different objectives, they refer to Demir, Bektaş, and Laporte (2014a) for the solution approach, including the hybrid multiobjective approach, a local search scheme, and adaptive LNS.Zarouk et al. (2022) minimise the energy consumption and maximise the customer satisfaction in a stochastic vehicle routing and scheduling problem.Their metaheuristic is based on a genetic algorithm and simulated annealing.Camacho-Vallejo et al. (2022) and Shi et al. (2023) focus on the maximisation of profits and the minimisation of CO 2 -emissions in bi-objective VRP variants.They apply tabu search metaheuristics.
Environmental aspects are discussed in literature for JIT applications, too: Karakostas, Sifaleras, and Georgiadis (2020) describe the fleet size and mix pollution location-inventory routing problem.Their objective is to minimise facilities' opening costs, inventory holding costs, general routing costs, including driver wages and vehicle usage costs, and fuel and CO 2 emissions costs.They propose a variable neighbourhood search metaheuristic as a solution procedure.Borumand and Beheshtinia (2018) consider the supply chain scheduling problem, where they minimise the delivery tardiness, the production costs, the emissions of suppliers and vehicles, and maximise the order quantity.A genetic algorithm is proposed as a solution procedure.Lee and Prabhu (2016) apply a feedback control algorithm for the JIT delivery problem, including CO 2 emissions costs, which is similar to a VRP with soft time windows.Fatemi-Anaraki et al. (2021) introduce the green delivery and pickup problem, where they evaluate the performance based on the fleet size against fuel consumption and greenhouse gas emissions by applying a K-means and genetic algorithm.Ganji et al. (2020) focus on an integrated supply chain scheduling problem aiming for the minimisation of transportation costs, tardiness costs, and environmental costs.They apply different multi-objective metaheuristics.Time windows and a heterogeneous fleet are included in the problem formulation, too.Wang et al. (2019) consider the so-called multi-depot green vehicle routing problem, where the interplay between minimal CO 2 emissions and minimal total operating costs, including costs for transportation, penalties for violating time windows, and vehicle maintenance, is surveyed.A particle swarm optimisation algorithm is used in this context.
Finally, the most relevant literature for the envTRP-JIT is summarised in Table 1.First, we refer to the VRP literature to extend the model and solution procedure of the TRP-JIT from an environmental perspective.The TRP-JIT can be seen as a routing problem with a manyto-one structure for JIT supplier networks minimising earliness and tardiness penalties.Second, we derive the definition of an environmental objective for the TRP-JIT, which aims for minimal CO 2 -emissions based on common approaches in green logistics and routing.Third, weighted sum and ε-constraint methods are proposed and applied for bi-objective studies on the economic and environmental perspective of JIT supplier networks.Consequently, this article contributes to the green routing literature by evaluating the environmental impact of JIT supplier networks.

Problem description
In this section, we formalise the envTRP-JIT based on the TRP-JIT by Baals, Emde, and Turkensteen (2023), which does not only consider the earliness-tardiness costs, but also a weight-distance objective to account for the environmental impact.We start with the formal definition of the envTRP-JIT in Section 3.1.An example is introduced in Section 3.2, focusing on the interplay and trade-off between the JIT perspective (earlinesstardiness) and environmental aspects (weight-distance).Finally, we present a MIP model for the envTRP-JIT based on its formal definition in Section 3.3.
For the envTRP-JIT, the following assumptions apply: • The fleet is homogeneous, i.e. the 3PL uses a fixed number of vehicles with the same capacity.• The vehicles' capacities are one-dimensional.
• The parameters can all be expressed by integer values, e.g. by scaling parameters if necessary.• The penalties for earliness and tardiness of a shipment can be represented by piecewise linear convex semicontinuous functions.• The triangle inequality holds for the travel time and distance matrix.• The fleet of the 3PL is also homogeneous in the curb weight of the vehicles.• The environmental impact of a route depends on the total weight of the truck, i.e. the curb weight and load, and the distance it travels.• The speed of the vehicles is predetermined by an average speed, i.e.it is not a decision variable.

Formal description
We define the following sets to describe the envTRP-JIT: • M for the m trucks, i.e.M = {1, . . ., m}, • N for the n shipments, N = {1, . . ., n}, • V for a set of all locations, which also includes the depot (node 0), n shipments, and OEM (node η = n + 1), V = N ∪ {0, η}, and • A for the arcs between the depot, shipments and the OEM, The envTRP-JIT is then defined on a graph G = [V, A, t, c], where V represents the set of nodes, and A the set of arcs, t : A → N as the function for the travel times, and additionally c : A → N as the function for the travel distances.For both, the travel times and distances, we apply the index notation c ij and t ij to express the travel distance and time along an arc [i, j] ∈ A. The parameters for the capacity, and the JIT properties are determined by: A solution for the envTRP-JIT, consists of the following components to specify m routes, which start at the depot (node 0), pickup at least one shipment i ∈ N, and end at the OEM (node η): sequence of the shipments on the routes, and • arrival times τ k to specify at which time a truck k ∈ M arrives at the OEM.
A due date in the envTRP-JIT is the date when a shipment should arrive at the OEM.Deviations from this date are possible, but cause either earliness or tardiness costs.For a single shipment this is modelled by a lower semi-continuous, piecewise-linear and convex function ρ i : N → N + 0 .More specifically, for any shipment i ∈ N, we define the penalty at arrival time τ by (1) The pickup times ti of a shipment i ∈ N are not explicitly part of a envTRP-JIT solution.Only the sequence of pickups on a route is determined.Based on the travel times t ij , ∀[i, j] ∈ A, and the release dates r i , ∀i ∈ N, the earliest pickup times t i ∀i ∈ N are recursively defined by (2) Then, the earliest arrival time τ e k for each route k ∈ M results in A solution for the envTRP-JIT is feasible if, and only if, the arrival time of each truck is at least as late as the earliest arrival time, and the total load on a route does not exceed the trucks' capacity, i.e.τ e k ≤ τ k and i∈R k q i ≤ Q, ∀k ∈ M. The consideration of environmental aspects does not affect the feasibility of a solution.
Given a solution for the envTRP-JIT, we can define the route segment distance c i,j , which is the distance between (4) The total length of a route ck , k ∈ M is either determined by the sum of the travel distances between all subsequent nodes (depot, shipments, and OEM) on a route or, alternatively, the travel distance from the depot to the first shipment on the route and the route segment distance from this shipment to the OEM.Formally, this is expressed by ( 5 ) Finally, we can define the objectives for the envTRP-JIT.The JIT objective f ET aims for minimising the total earliness-tardiness costs based on the penalty functions in (1).It is formally expressed by From an environmental perspective, we aim to minimise the CO 2 emissions represented by the weight-distance, which sums up the curb weight and shipments' weight multiplied by the distance they are respectively transported on.The weight-distance objective f WD is formalised by The expression i∈N q i c i,η refers to the total loaddistance of a certain envTRP-JIT solution, and q C k∈M ck to the total travel distance weighted by the trucks' curb weight q C .For each shipment, we consider the travel distance needed for delivery on a given route multiplied by the shipment's weight.The curb weight is included for the whole length of the routes.This is equivalent to first determining the total weight (load's weight and the curb weight q C ) on a route segment (an arc between two subsequent nodes on a route), and multiplying this weight by the segment's travel distance.We avoid combining the objectives f ET and f WD explicitly at this point, because our goal is to study the interplay between both objectives rather than assuming a certain preference.The search schemes in Section 4.1 operate differently on the objectives.

Example
We consider an example to illustrate the motivation and formal definition of the envTRP-JIT.It is depicted in Figure 1, showing the optimal solutions according to the earliness-tardiness and the weight-distance objective.The relevant parameters are mentioned in Table 2, except the travel times and distances, which are included in the appendix in Tables A1 and A2.Additionally, the trucks' curb weight is set to q C = 6400 and the load capacity to Q = 3600.The variable values for both solutions can be found in the appendix in Table A3, too.The example is based on that by Baals, Emde, and Turkensteen (2023).They modify the R107 instance of Solomon (1987) for the VRP with n = 25 suppliers.
Compared to the original TRP-JIT example, we scale the capacity demands q i , ∀i ∈ N, and the trucks' capacity Q by the factor 36.After this modification, the ratio between curb weight and the load capacity is approximately the same as in the instance introduced by Bektaş and Laporte (2011) for the PRP: They assume a curb  Note: The travel times and distances are specified in Tables A1 and A2.
weight of 6350 kg and a load capacity of 3650 kg, which gives a total weight of 10 metric tons.We assume that the average speed on each arc is set to 60 km/h.Then, the travel times are given in minutes [min] and the travel distances in kilometres [km].The travel times and distances are based on the euclidean distances.A service time of 10 minutes is included in the travel time if an edge leaves a shipment node.The travel times are rounded to the next integer, whereas two decimals' precision applies for the travel distances.When looking at the graphs of both solutions in Figure 1, it can be seen that the routes do not cross over each other for the weight-distance objective in Figure 1(b), whereas for the earliness-tardiness objective in Figure 1(a) they do.These crossovers seem to lead to rather long travel distances, i.e. non-preferred routes, from an environmental perspective.This is also reflected in the quality of both solutions according to the different objectives.The minimal earliness-tardiness objective is f (ET)  ET = 563.65 and leads to a weightdistance of f (ET) WD = 2,274,853.56,whereas the optimal weight-distance is given by f (WD)  WD = 1,677,336.44implying earliness-tardiness costs of f (WD) ET = 1357.56.Consequently, in the example, there is a certain trade-off between both objectives, which we survey in this paper.For the sake of completeness, the optimal solution for the earliness-tardiness objective is given by

Mixed-integer programming model
We adapt the MIP model by Baals, Emde, and Turkensteen (2023) to solve the envTRP-JIT.The model applies the previously introduced sets and parameters.Following Baals, Emde, and Turkensteen (2023), we set the big-M constant to P = max i∈N {max{r i , p i , t 0i }} + i∈N max j∈N∪{η} {t ij }.The MIP model is based on several variables to encounter for the different aspects of an envTRP-JIT.We include the following two binary variables to define the routes • x ijk , which is set to 1 if truck k travels directly from depot or shipment i ∈ N ∪ {0} to shipment or OEM j ∈ N ∪ {η} on its route, and 0 otherwise, and • y ik , which is set to 1 if shipment i ∈ N is assigned to the route of truck k ∈ M, and 0 otherwise.
To express the JIT objective in ( 6) by a linear equation, we include continuous variables for each shipment i ∈ N, i.e. e i for the earliness of its arrival at the OEM, d i for the tardiness of its arrival at the OEM, and ti for its pickup time at the supplier.For the arrival time at the OEM of each truck k ∈ M, we include a variable τ k .To convert the environmental objective in (7) into a linear equation, it is required to consider the continuous flow variables u ijk for each arc [i, j] ∈ A, truck k ∈ M, and i = 0 similar to Kara, Kara, and Yetis (2007).
Hence, the following objectives and constraint set provide a MIP model for the envTRP-JIT: Minimise: subject to: The objective in ( 8) is the linear equivalent to the earliness-tardiness objective in ( 6).The same applies to the weight-distance objective introduced in (7), and the linear objective in (9).The requirements that each shipment must be assigned to exactly one route, each route must start at the depot, and each shipment must be picked up by a truck exactly once are defined by the constraints in (10a)-(10c).The constraints in (10d) imply that if a truck enters a shipment node, the respective shipment is assigned to the route.The flow constraints in (10e) are required to optimise the load-distance as a part of the weight-distance.The constraints also avoid subcycles.Based on the flow constraints, the inequalities in (10f) account for the capacity limit.The constraints in (10g)-(10l) refer to the JIT perspective.We ensure by (10g) that the pickup at the supplier takes place not earlier than the release date.The constraints in (10h)-(10j) determine the pickup times at the suppliers and arrival times at the OEM such that the travel times are respected.The constraints in (10i) avoid subcycles, similar to the commodity flow constraints.The expressions in (10m)-(10q) determine the domains of the variables.We consider the valid inequalities defined by Baals, Emde, and Turkensteen (2023) for the envTRP-JIT, too, when solving it by a commercial solver.Additionally, we include valid inequalities for the commodity flow variables mentioned in Appendix B.

Solution method
In this section, we define the solution method for the bi-objective envTRP-JIT.The general framework of the procedure is outlined by two steps: (1) A bi-objective search scheme to find several nondominated (i.e.Pareto efficient) solutions by varying the parameters for the objectives.(2) A solution procedure to find a good solution for these parameters.
The bi-objective search schemes in the first step are different in the way they consider the objectives and how they vary on the parameters for the objectives.Section 4.1 focuses on these aspects.For the second step, we extend the LNS algorithm defined by Baals, Emde, and Turkensteen (2023) for the single objective TRP-JIT.We discuss this in Section 4.2.Alternatively, a commercial solver implementation of the MIP in Section 3.3 can be applied.

Bi-objective search schemes
We consider bi-objective search schemes, which either include all objectives in the objective function, but vary on the factors of the single objectives (weighted sum methods), or only use one objective in the objective function, and are constrained on the other ones (ε-constraint method).Formally, this means that we vary over one or several of the following parameters: • λ ET -weight factor of the earliness-tardiness objective in the objective function, • λ WD -weight factor of the weight-distance objective in the objective function, Moreover, we use the following definitions to formalise the algorithms of the bi-objective schemes for a certain problem: The bi-objective search schemes are generally independent of the solution method we apply to calculate solutions for a single set of parameters.Hence, in Algorithms 1-4, we use the following function to describe the procedure in an abstract manner, The output consists of the following elements: The tuple L, including the objective value tuple F of the best found solution according to the parameters λ ET and λ WD for the objective function, and the tuple S of its variable values.The set L can include several solutions depending on the solution procedure.Additional to the best solution according to the parameters λ ET and λ WD , a solution procedure may find better solutions for either the earlinesstardiness or the weight-distance objective.This becomes especially relevant when solving the problem with the bi-objective LNS algorithm.In general, the subscript of the objective function indicates the considered objective and the superscript the optimisation goal, e.g.f (WD)   ET represents the value of the earliness-tardiness objective when minimising the weight-distance.We either apply a MIP model introduced in Section 3.3 and specified in Appendix C, or a bi-objective version of the LNS algorithm by Baals, Emde, and Turkensteen (2023) presented in Section 4.2 to solve the bi-objective TRP-JIT.

Weighted sum methods
Algorithm 1 Bi-objective weighted sum method with constant intervals and normalisation. Input In the weighted sum methods, the factors λ ET and λ WD of the earliness-tardiness objective and the weightdistance objective vary.It does not include constraints on any of both objectives, meaning that ε ET = ε WD = ∞.The weighted sum method is also applied by Demir, Bektaş, and Laporte (2014a) to the bi-objective PRP.We present two schemes based on the weighted sum method.
The first scheme is formalised by Algorithm 1.The weights of the objectives consist of a normalisation factor, and a factor that we vary linearly.This means, we increase the factor for one objective, while it is decreased for the other objective with a constant stepsize.The normalisation factor is derived from the objective values when only considering one of both objectives.Based on κ max , the number of steps κ max + 1 and the step size 1 κ max are defined.
In Algorithm 2, we can find the second scheme.It is based on the theoretical work of Cohon (1978); Kim and de Weck (2005), and Przybylski, Gandibleux, and Ehrgott (2010) presenting a so-called dichotomic scheme.Here, the intention is to find all non-dominated extreme points by repeatedly searching for new solutions between already known ones.This is done by defining the weights of the objectives, λ ET = λ = Algorithm 2 Recursive bi-objective weighted sum method.
Input κ searchDepth,max Output L, Lpareto WD ) from the permutation of known Pareto efficient solutions Fpareto .We avoid reevaluating specific weighted sum combinations several times by storing all evaluated values of λ in a set .Especially for larger problems with many shipments, the number solutions to be found by the second scheme may become large.This implies a high computational time, which we aim to avoid by limiting the search depth κ searchDepth,max .It means the repetitions of searching for new non-dominant solutions between the already found ones.Contrarily, in the first scheme, the number of solutions to be calculated remains constant for problems with a larger amount of shipments.

ε-constraint methods
In the ε-constraint methods, we focus on minimizing one single objective and consider the other objective in the constraint set by varying its upper bound.Bektaş and Laporte (2011) consider the ε-constraint method for the bi-objective PRP, too.Ehrgott (2005a) and Mavrotas (2009) provide theoretical work about ε-constraint methods and propose different schemes.We define εconstraint method implementations for the envTRP-JIT in Algorithms 3 and 4. Regarding the objectives, this either means that we consider earliness-tardiness in the objective function (objective = 'ET'), and weightdistance in the constraint set (constraint = 'WD'), i.e. λ ET = 1, λ WD = 0, ε ET → ∞, and ε WD < ∞, or the opposite (objective = 'WD', constraint = 'ET'), which means, λ ET = 0, λ WD = 1, ε ET < ∞, and ε WD → ∞.In both schemes, we first calculate the solutions for each objective respectively.
In the first scheme in Algorithm 3, we then define the range to vary the upper bound ε of the constrained objective with a fixed stepsize.Depending on the variable κ starting with κ = 1 until κ = κ max , we derive the upper bound from the weighted sum of both objective values, i.e.ε = κ κ max f constraint .The number of solutions to be calculated remains constant for an increasing amount of shipments in the problem.But it may be the case that we get the same solution for consecutively calculated values of ε because we do not consider the results of previous iterations where the lower value for ε is already satisfied.We start with assuming the earliness-tardiness costs as objective and constrain the weight-distance.Afterwards, we repeat the procedure, and vice versa, assuming the weight-distance to be the objective and constraining the earliness-tardiness costs.
In the second scheme in Algorithm 4, we aim to find the complete Pareto frontier.Therefore, we recursively set the upper bound by where δ f is chosen sufficiently large, such that ε constraint < f (κ−1) constraint holds.Instead of iterating over the variable κ, we repeat the solution calculation while ε constraint > f (constraint) constraint holds.

Bi-objective large neighbourhood search metaheuristic with local search scheme
To define the bi-objective LNS metaheuristic, we extend the LNS metaheuristic of Baals, Emde, and Turkensteen (2023) by the weight-distance objective.We refer to Baals, Emde, and Turkensteen (2023) for the details about the algorithm for the single (earliness-tardiness) objective case.The bi-objective variant requires including additional slacks in the solution procedure and smaller changes in the general structure of the algorithm and its components.The procedure itself is outlined in Algorithm 5.The structure of the Algorithm 5 Large Neighbourhood Search for the envTRP-JIT ,0) is feasible then L (best) ← L (LS,0) end if // stopping criterion while s max > s and s imp,max > s imp and computation-Time < timeLimit do L (D,s) ← destroy L (LS,s) L (DR,s) ← repair L (D,s)  // acceptance criterion if eval L (DR,s) < eval L (DR,ref) or s DR ≥ s DR,max then L (DR,ref) ← L (DR,s) L (LS,s) ← localSearch L (DR,s) L (ET) , L (WD) ← updateBestSolutions L (LS,s) , L (ET) , L (WD)   if eval L (LS,s) < eval L (best) and S (LS,s)  best) , L (ET) , L (WD)  metaheuristic is defined by the constructive heuristic to get the initial solution, the destroy and repair operator, and the local search, including the acceptance criterion about whether to perform the local search.We consider f = λ ET f ET + λ WD f WD as primary objective to steer the search.But we also include the sum of the earliest arrival times k∈M τ e k as secondary auxiliary objective in a hierarchical manner as a tiebreaker.The comparison of both objectives is implied by the function eval(L) for a solution L. In the following, we focus on summarizing the main components of the algorithm, and show the differences to the single objective LNS.
The initial solution of the bi-objective TRP-JIT is generated based on the parallel insertion heuristic introduced by Potvin and Rousseau (1993).It inserts the shipments in the routes one after another by evaluating the 2-regret criterion.This means given a solution L that we determine for each shipment i ∈ U not inserted yet the difference of the evaluation function between the best and second best insertion position.Formally, this is expressed by regret(L, i) = eval(L * ,2 i ) − eval(L * ,1 i ) with L * ,1 i as the best and L * ,2 i as the second best solution for inserting i.Then we choose the shipment i * with the highest regret value to be inserted, i.e. i * = arg max i∈U {regret(L, i)}.We repeat this procedure until every shipment is assigned to a route.This results in solution L (DR,0) in Algorithm 5. Computing the initial solution requires a complexity of O(n 3 log n).
The local search is characterised by the neighbourhood operators and the evaluation criteria and slacks to evaluate neighbours in a computationally efficient manner.In Algorithm 5, performing the local search results in solution L (LS,s) .We apply the following neighbourhood operators in the local search based on a given solution: • the exchange operator to mutually replace a single shipment or sequence of shipments with another single shipment or sequence of shipments in the same or a different route, • the relocation operator to move a single shipment or a sequence of shipments to another route or position of the same route, • and the 2-opt * operator to replace two arcs between a pair of routes such that their tails, i.e. the sequences at the end of the routes, are exchanged.
In the exchange and relocation operator, the sequences can stand for a single shipment up to a sequence of n seq shipments.We consider each possible combination of two sequences in the exchange operator.The size of the neighbourhoods scales with the number of shipments n quadratically.The neighbourhood operators are similar to those applied in literature by Baals, Emde, and Turkensteen (2023), Bräysy andGendreau (2005a), Gendreau, Hertz, andLaporte (1992), Irnich, Funke, and Grünert (2006), and Toth andVigo (2002, 2014).
We compute different evaluation criteria and slacks for the earliness-tardiness and weight-distance objective to evaluate the neighbourhoods efficiently.Vidal et al. (2014) provide an overview of different slacks applied in vehicle routing problems.In general, defining slacks for routing problems has a long history, e.g. by Savelsbergh (1992) and Kindervater and Savelsbergh (2018).We use i − to describe the predecessor of a shipment i on its route R k , k ∈ M in a certain solution L. The set R i includes all shipments on the route until shipment i ∈ R k .We at first define the slacks to be stored.Then we show how to apply them on any subsequence of shipments σ = i, . . ., j in a given solution, and finally on a route σ 1 ⊕ σ 2 = i 1 , . . ., j 1 , i 2 , . . ., j 2 , which is a combination of two sequences σ 1 = i 1 , . . ., j 1 and σ 2 = i 2 , . . ., j 2 .A more detailed description of the slacks relevant for the earliness-tardiness costs, including an example, is provided by Baals, Emde, and Turkensteen (2023).
For the capacity demands, we store the slack q i = i∈ R i q i ∀i ∈ N, which stands for the load on the route until shipment i ∈ N. The required capacity of a sequence σ is then calculated by qσ = q j − q i + q i , and for the route σ 1 ⊕ σ 2 by qσ 1 ⊕σ 2 = qσ 1 + qσ 2 .
The partition-optimal arrival time τ opt k of a route k is only feasible if it is not before the earliest arrival time τ e k .To determine the earliest arrival time, we calculate the waiting time w i = r i − ( t i − + t i − i ) for each shipment i ∈ N, where we assume t i − = 0 for i − = 0. Based on the waiting times, we calculate the total forward time buffer w ij = i ∈ R j \ R i max{0, w i }, and the maximum backward time buffer ← − w ij = max i ∈ R j \ R i min{0, w i } for each possible subsequence i, . . ., j with i = j of the routes from a solution S. We use these slacks to evaluate the effect of a change of the earliest pickup time t i at the beginning of sequence σ on its end, i.e. the change of the earliest pickup time t j .This is given by the following expression When concatenating the sequences σ 1 and σ 2 , this implies the following steps to calculate the earliest arrival time τ e σ 1 ⊕σ 2 : (1) Determine the new earliest pickup time at the start of σ 1 by t new,i 1 = max{r i 1 , t 0i 1 } and its change t i 1 = t new,i 1 − t i 1 , (2) Calculate t j 1 at the end of sequence σ 1 by using the slacks, (3) Compute t new,i 2 = max{r i 2 , t j 1 + t j 1 + t j 1 i 2 } and t i 2 = t new,i 2 − t i 2 , (4) Calculate t j 2 at the end of sequence σ 2 by using the slacks, (5) Compute τ e σ 1 ⊕σ 2 = t j 2 + t j 2 + t j 2 η .
If σ 1 or σ 2 only consists of a single shipment, then t i 1 = t j 1 or t i 2 = t j 2 holds and we can skip step 2 or 4. If the partition-optimal arrival time is feasible, i.e.τ e k ≤ τ opt k holds, it is also the optimal arrival time for the route k ∈ M and we set τ * k = τ opt k .The earlinesstardiness costs correspond to the stored penalty at the certain due date.Otherwise, we state that the earliest arrival time is optimal, i.e. τ * k = τ e k .Then a binary search is required to find the closest due date before the earliest arrival time, the earliness-tardiness costs at this due date, and the slope of penalty function ρ σ 1 ⊕σ 2 at τ e k .Baals, Emde, and Turkensteen (2023) explain in detail how to calculate the earliness-tardiness costs at the earliest arrival time based on these values.
We calculate the travel distance ck of a route k ∈ M based on the route segment distance c ij for each subsequence of any route in a given solution S as introduced in Section 3.1.For the route σ 1 ⊕ σ 2 composed of subsequences σ 1 and σ 2 , this results in cσ The load distance vk of a route k ∈ M is evaluated based on the load added on a route segment v ij , which is defined by In the ε-constraint method, either the earlinesstardiness or the weight-distance objective is considered as constraints.In order not to strongly limit the search space, we allow exceeding the constraints with a penalty higher than any feasible objective value.If required, we consider the earliness-tardiness or weight-distance excess, i.e.
The slacks and evaluating criteria can be calculated in constant computational time O(1) except the definition of the partition-optimal arrival time.The binary search for the partition-optimal arrival time and possibly for the earliness-tardiness penalty at the earliest arrival time enforces a computational complexity of O(log n) to evaluate a single neighbour.Initially calculating all neighbours requires a complexity of O(n 2 log n).The computational complexity can be reduced to O(n 2 ) for any further iteration.If none of the objectives is constrained (weighted sum method), it is sufficient to store the best move for each pair of routes (∼ 0.5m(m − 1) moves) or route (∼ m moves) in each of the neighbourhoods because the slacks do not need to be reevaluated for routes not affected by the last move.Updating all neighbours affected by the manipulation of a route or pair of routes corresponds to O(n log n).If we constrain one of the objectives (ε-constraint method), it is necessary to store the complete neighbourhoods and to additionally update the excess value ( f ET,C or f WD,C ) for the constrained objective for all neighbours.Then we can avoid to reevaluate the optimal arrival times with a complexity of O(log n) for each neighbour and keep an overall complexity of O(n 2 ).Hence, when including the excess penalty (ε-constraint method), moves in a route or pair of routes affect all neighbours of other routes.
We terminate the algorithm by a stopping criterion.It is defined based on s max , i.e. the maximum number of total iterations, s imp,max , i.e. the maximum number of iterations without finding a new best solution L (best) , and a time limit.We count the total number of iterations by s, and without improvement by s imp .Within an iteration, the destroy and repair operation, and, if accepted, the local search is performed.In case of a constraint on one of the objectives (ε-constraint method), the stopping criterion for not finding an improving solution only applies if at least one feasible solution is found.A tight constraint may require more than s imp,max iterations to even find a feasible solution.
The destroy operator removes a certain amount of shipments from a given solution.The number of shipments to remove is at first defined based on a uniform distribution In the second step, this amount of shipments is randomly removed from the solution.We get the incomplete solution L (DR,s) and a set of unassigned items U from the destroy operator.
The repair operator applies the same principle like the parallel insertion heuristic for the initial solution to reassign all removed shipments to routes again.This means that we at first determine the value of the 2-regret criterion for each shipment.We then randomly choose one of the shipments to be inserted among v max shipments with the highest value of the 2-regret criterion.This results in solution L (DR,s) .Repairing a solution requires a complexity of O(n 2 log n) if the number of removed shipments is independent of n, otherwise it is O(n 3 log n).Destroy and repair operators are the core principles of the LNS introduced by Shaw (1998), and also of the adaptive LNS by Ropke and Pisinger (2006).
We apply an acceptance criterion to decide if the local search should be performed on the solution L (DR,s) after the destroy and repair operation.In case that L (DR,s) performs better than the reference solution L (DR,ref) , it gets accepted, otherwise we reject it and repeat the destroy and repair operation based on L (DR,ref) .We also accept a solution after rejecting for more than s DR,max iterations, which we count by s DR .
Instead of only tracking the solution L (best) , which is the best according to the objective function f, we also save the solutions, L (ET) and L (WD) , best according to the earliness-tardiness objective f ET and weight-distance objective f WD disregarding constraints on any objective function.Hence, the output set of solutions L is defined by L = {L (best) , L (ET) , L (WD) }.This means that in a single run of the LNS algorithm, we can find several Pareto efficient solutions, which especially affect the behaviour of the recursive weighted sum method and the output of all bi-objective search schemes.We also repeat the LNS algorithm for the same set of parameters λ ET , λ WD , ε ET , and ε WD several times defined by parameter κ rep,max .Due to variations in the search between the repetitions, we can find additional solutions.

Example
We consider the example introduced in Section 3.2 again.The bi-objective search procedures defined in Section 4.1 are applied on the example either using the MIP model implemented in a commercial solver (Gurobi) from Section 3.3 and Appendix C, or the LNS in Section 4.2.The results are illustrated in Figure 2, where the earlinesstardiness costs are plotted against the weight-distance objective.Figure 2(a) shows the Pareto frontiers estimated by the weighted sum methods and Figure 2(b) by the ε-constraint methods.We compare the Pareto frontier resulting from a particular procedure against the set of Pareto efficient solutions derived by aggregating over all procedures.Figure 2 shows that all of the bi-objective search procedures are close to or on the aggregated Pareto frontier.It can also be seen that in the example, there is a certain trade-off between the earliness-tardiness costs and the weight-distance objective representing the CO 2emissions.In Figure 2(a), it can be seen for the weighted sum methods that the LNS compared to the solver finds additional solutions close to or on the aggregated Pareto frontier because it includes several solutions in the set L. The weighted sum method performs well on finding a convex hull for the Pareto frontier, but neglects the non-convex parts of the curve (Przybylski, Gandibleux, and Ehrgott 2010).The differences between the solution procedures and the interplay between the objectives are surveyed in more detail in the computational study.

Computational study
We first define the instances and computational environment in Section 5.1.Then, we compare the different methods for the bi-objective search, and the LNS against the MIP model in Section 5.2.We discuss the relationship and trade-off between the just-in-time and environmental perspective in Section 5.3.In Section 5.4, we investigate the impact of the time window tightness, i.e. the time span between release and due date, and the fleets' composition on the environment and punctuality.

Benchmark instances, definition of the algorithmic parameters and computational environment
The computations are based on a modification of the instances for the PRP by Demir, Bektaş, and Laporte (2012).There, the locations represent different cities in Great Britain and real-world road distances.Several sets of instances are provided, comprising 20 instances per set.All instances in a set include the same number of shipments.Generally, the instances include information about the distance between the locations, demands to be delivered, curb weight, load capacity, and time windows.We derive the instances for the envTRP-JIT based on these information.Then, we must make additional assumptions, especially about the JIT perspective.The parameters for the instance generation are described in the following paragraphs.
We consider respectively 20 envTRP-JIT instances with either n = 12 or n = 99 shipments.These are modifications of the PRP instances with 15 or 100 pickup locations.The OEM is assumed to be at the last location in the PRP instances.The instances with n = 12 shipments are derived from those of the PRP with 15 pickup locations by leaving away the second and third last location.This ensures that the commercial solver can find optimal solutions in reasonable computational time.To define the weight-distance objective, we use the distances c ij , ∀[i, j] ∈ A (in meters) from the PRP instances accordingly.The same applies to the curb weight q C = 6350[kg], the load capacity of Q = 3650[kg], and the respective capacity demands q i , ∀i ∈ N.
Additionally, it is required to define the parameters to compute the earliness-tardiness costs.Therefore, we assume a constant speed of ν = 60[km/h] = 1000[m/min].The service times s i , ∀i ∈ N are distributed according to the PRP instances.Then, the values are scaled and rounded from seconds to minutes.For the depot (node 0), the service time s 0 = 0 is set.Hence, we derive the travel times by t ij := We use the same approach as Baals, Emde, and Turkensteen (2023) to define the release dates r i , due dates d i , and the penalties for earliness α i and tardiness β i , ∀i ∈ N: We set a planning horizon of t plan = 480[min].The release dates r i are uniformly distributed on the interval [0, t plan − t iη ], and the due dates d i on [max{r i , t 0i } + t iη , t plan ] for all shipments i ∈ N. The earliness penalties α i are randomly drawn from a uniform distribution on the interval [0.70, 0.99], and the tardiness penalties β i on [1.00, 1.29], rounded to two decimals for each shipment i ∈ N.
The bi-objective search schemes are performed assuming κ max = 21 for the weighted sum method with constant intervals and normalisation (WM constant) in Algorithm 1, and κ max = 11 for the εconstraint method with constant intervals (EC constant) in Algorithm 3.For the recursive weighted sum method (WM recursive) in Algorithm 2, we set the search depth κ searchDepth,max = 5 for the instances with n = 12 shipments, and κ searchDepth,max = 3 for the instances with n = 99 shipments.When applying the LNS within the biobjective search schemes, we repeat the computations five times, i.e. we set κ rep,max = 5.Within the LNS, we set the limit for the total number of iterations s max = 50,000, and the number of iterations without any improvement to s imp,max = 2500 for the stopping criterion.The parameter for the acceptance criterion is defined by s DR,max = 5.The upper and lower bound for the destroy operator are n − rem = 5 and n + rem = max{6, 0.4n }.The number of shipments to be considered in the repair operator is set to ν max = 5.In the local search, we consider sequences of up to n seq = 3 shipments.
We use C# 8.0, targeting the.NET Core 3.1 framework to implement the algorithms.The computational study was conducted on an x64 PC with Intel Core i9-10900K, 128GB RAM, and Microsoft Windows 10.We compare the performance of the LNS to an exact commercial solver (Gurobi 9.1.1).The commercial solver uses all 10 cores of the CPU, whereas the LNS is limited to one core.We repeat the LNS for the same set of parameters five times (κ rep,max = 5) on parallel threads.

Comparison of the solution methods
We compare the different bi-objective search schemes and solution procedures introduced in Section 4. We use the computation time of the different procedures to compare their performance.We count the number of different solutions (# different) found by a search scheme and derive its Pareto frontier.We determine the number of solutions that are Pareto efficient (# PE) within the solution procedure.From the different procedures, we can derive a best known Pareto frontier and the number of solutions it consists of (# PE (all)).Then, we compare the solutions of the Pareto frontier of a specific procedure against the best known one.We count the number of Pareto efficient solutions found by a certain procedure and stated to be best known (# PE (bk)).We also evaluate the hypervolume and the spacing metric of the Pareto frontiers and compare it between the different solutions and the best known Pareto frontier.The definitions of the hypervolume and the spacing metric are illustrated in Figure 3.The hypervolume V comprises the volume or area of the rectangles between the data points of the Pareto frontier and a reference point (also called the Nadir point).We define the reference point by (f ET,ref , f WD,ref ) = (1.01fET,max , 1.01f WD,max ), i.e. the point represents the worst case for both objectives, f (ET) max and f (WD) max , from the Pareto frontiers of the different methods additionally increased by 1%.Formally, this means , where we assume that the solutions are sorted in decreasing order of the earliness-tardiness penalty f ET for increasing i ∈ Fpareto .Instead of considering the absolute value of the hypervolume, we evaluate the hypervolume relative to the hypervolume of the Pareto frontier aggregated over all solution procedures V all , i.e.V = V V all .For more details about the hypervolume, see Demir, Bektaş, and Laporte (2014a) for the bi-objective PRP.Audet et al. (2021) and Zitzler et al. (2003) survey the performance of multi-objective procedures in general.
The spacing metric measures the Manhattan distance δ (i) ∀i ∈ Fpareto between each point of the normalised Pareto frontier and its closest neighbouring solution and considers the estimated standard deviation σδ of these values.In the results, we additionally include the mean μδ .Formally, these metrics are defined by δ (i) δ (i) , and Smaller values of μδ and σδ are better.
The results of the computations are summarised in Table 3 for the instances with n = 12 (DB13) and n = 99 (DB100) shipments.The detailed results are included in Appendix D. For the smaller instances, we used both the solver and the LNS to calculate the solutions for specific weights of the objectives or the upper bound for the ε-constraint.If we use the solver in the search scheme, we hierarchically optimise on both objectives.Hence, the EC recursive procedure should find all Pareto efficient solutions when using the exact solver.We refer to Baals, Emde, and Turkensteen (2023) for the evaluation of LNS metaheuristic's performance on a single objective, i.e. the earliness-tardiness penalty.The authors demonstrate an average gap to the best known solution of 1%.
We can see that the LNS clearly outperforms the solver concerning the computation time while both solution procedures achieve a similar hypervolumes.The solver only provides the best solution according to the objective function and constraints defined by the parameters λ ET , λ WD , ε ET , and ε WD .Additionally, we also track the best solution according to the earliness-tardiness costs and the weight-distance objective in the LNS.This explains the higher number of solutions found when using the LNS compared to the solver.The rate of Pareto efficient  solutions among the different output solutions #PE #different is 100% the solver and only between 26.5% and 47.65% for the LNS.More Pareto efficient solutions are still found in absolute numbers by the LNS than by the solver both within the procedure and compared to the best known Pareto frontier, except for the EC recursive scheme.This is also indicated by the rates #PE (bk)  #PE(all) in Table 3.The EC recursive scheme aims to find all Pareto efficient solutions.Both the solver and the LNS find the same amount of Pareto efficient solutions in this scheme.The rates #PE(bk) #PE show that for the solver, the share of Pareto efficient solutions of a procedure also included in the best known Pareto frontier tends to be higher.It means that for the LNS it is more likely to get solutions close to but not on the best known Pareto frontier.This effect is negligible as the number of Pareto efficient solutions found by the LNS is higher anyway -also for those being best known, and there is no significant negative effect on the hypervolume.Consequently, when using the LNS in a bi-objective scheme, a Pareto frontier can be computed substantially faster, but it may include solutions inferior but very close to the best known.In the EC recursive scheme, the LNS is faster than the solver while providing equivalent solution quality.For the larger instances, we only apply the LNS because the solver cannot provide solutions within reasonable computational time.
For the smaller instances (DB-13), the ε-constraint methods outperform the weighted sum methods regarding the solution quality, whereas for the larger instances (DB-100), we state the opposite.When looking at the computation time, especially for the larger instances, the ε-constraint method takes significantly longer than the weighted sum methods.We can state that the computation time of the LNS is notably higher if one of the objectives is considered as a constraint.For this reason, we did not perform the computations on the recursive εconstraint method (EC recursive) for the larger instances because the number of solutions on the Pareto frontier grows and affects the computational time, too.According to our performance metrics, the weighted sum methods provide Pareto efficient solutions on a better quality level within a shorter computation time for the larger instances.We base the computations for the managerial insights on this procedure.Finally, we choose the WM constant for the following considerations because it provides more than 40 Pareto efficient solutions compared to the best known frontier on average, and at least 60 Pareto efficient solutions on the estimated frontier for the larger instances.The computations were completed in less than 25 min for each instance.
The performance of the ε-constraint method seems to be weak for the larger instances when looking at the number of solutions of the procedure also included in the best known Pareto frontier but acceptable regarding the hypervolume, the mean distance, and the spacing metric.The results for the hypervolume and spacing metric lead to a similar conclusion about the performance of the procedures.When looking at the EC constant solver and the WM constant LNS method for the DB13 instances, we draw different conclusions based on the hypervolume and spacing metric for example.The Pareto frontiers of the first instance with n = 12 shipments (DB13-1) and n = 99 shipments (DB100-1) are illustrated in Figure 4 for the different bi-objective search schemes.The Pareto frontiers of the smaller instance look similar to those of the example instance in Figure 2.For the larger instance, we state that the ε-constraint method finds solutions that are close but not on the Pareto frontier found by the weighted sum methods, which explains the low value of # PE (bk).In the ε-constraint method, only either the earliness-tardiness or the weight-distance define the objective function when using the LNS.Hence, there may be better solutions according to the other objective, which is only considered as a constraint.

Managerial insights: a given scenario
We can observe in Figure 4 that the graph seems to resemble a hyperbole for the larger instances.From a managerial perspective, this means that extreme punctuality must be paid for by relatively high weight-distance, while compromising on earliness-tardiness has very little effect on weight-distance past a certain level.
We generally assume the best solution according to the earliness-tardiness costs to be the starting point to gain managerial insights because most supply chains in practice still prioritise JIT objectives over sustainability goals.We aim to evaluate the potential to reduce the environmental impact.we call the total relative increase of either the earlinesstardiness or weight-distance objective.This is formally described by . It means, on one hand, that when focusing on the JIT perspective, our costs from the environmental perspective are between 12.2% and 61.3% higher than in the best known case for the smaller instances, and between 45.3% and 58.6% for the larger instances.On the other hand, when aiming for the weight-distance objective, the earlinesstardiness costs are between 65.9% and 703.7% higher for the smaller instances, and between 231.0% and 622.7% higher for the larger instances than in the best known solution.Consequently, when trading off between both objectives, the risk of increasing the earliness-tardiness costs is relatively higher than the potential for improving the weight-distance objective.Instead of making a binary decision for one of both objectives, we consider the known Pareto efficient solutions between these two extreme points.As the starting point for decision making is the best known solution for the earliness-tardiness costs, we focus on the relative deviation from this solution.It means that we normalise the Pareto frontier by the best known solution for the earliness-tardiness costs and shift the graph such that this solution is located at the point of origin.
Formally, this is expressed by fET = for the earliness-tardiness costs f ET and by fWD = Figure 5 illustrates the normalised and shifted known Pareto frontier of instance DB100-1 from the weighted method with constant intervals and normalisation.Figure 5(a) shows the complete range of the plot, whereas Figure 5(b) focuses on the range for reducing the weightdistance up to 15%.The plot confirms that the potential of improving the weight-distance is relatively low compared to the implied increase of earliness-tardiness costs, which we already concluded in the previous from Table 4. Additionally, we state based on Figure 5 that we can improve the weight-distance objective to a certain degree without excessively increasing the earlinesstardiness costs.For example, reducing the environmental impact measured by the weight-distance by 2.0% only causes 1.2% higher earliness-tardiness costs, or a reduction by 8.6% causes around 8.3% higher costs from the JIT perspective.Any further reduction of the weightdistance objective is linked to relatively high earlinesstardiness costs.In Section 5.4, we observe a similar interrelation between both objectives for the other instances, too.
We can state that aiming for lower CO 2 emissions increases the earliness-tardiness penalties.The costs for minimizing the environmental impact are relatively high for a JIT supplier network but reasonable for a reduction to a certain degree.Hence, it is possible to lower the environmental impact of transportation in a JIT supplier network to a limited extent.

Managerial insights: potential of negotiations between suppliers, OEM and logistics provider
When negotiations between the 3PL, the suppliers, and the OEM take place, decisions on the following aspects can be discussed to improve the network's performance: (1) The suppliers adapt the release dates such that the shipments are available to be picked up earlier.
(2) The OEM delays the due dates to increase the time span available for transportation without tardiness.(3) The 3PL decides on the fleet composition defined by the number of trucks, their load capacity, and important for the environmental perspective also on their curb weight.Baals, Emde, and Turkensteen (2023) show that enabling earlier pickups by the suppliers and later due dates at the OEM have a positive effect on the earlinesstardiness costs.A smaller fleet of larger vehicles increases the earliness-tardiness costs, whereas a large fleet of small vehicles leads to the opposite.The effects of the negotiations on the environment, i.e. the weight-distance objective and the interplay with the earliness-tardiness costs, are evaluated by the following parameter variations: • We vary on the release dates r new,i := r scale • r i ∀i ∈ N by factor r scale := 0, 0.25, 0.5, 0.75, 1, such that the shipments are released at latest at time 0, 120, 240, 360, 480 [TU] during the planning horizon.Instead of only aiming on finding the best solution according the earliness-tardiness costs, our goal is to provide an estimate of the Pareto frontier for the instances with n = 99 shipments, and to derive the effect of shifting the objective towards the environmental impact.
For the sake of clarity, the results of instance DB100-1 are illustrated in Figure 6 and of instance DB100-2 in Figure 7.The plots of the remaining instances with n = 99 shipments can be found in Appendix D. We can state that earlier release dates generally imply lower earliness-tardiness costs, even if we do not solely focus on this objective.The distance between the Pareto frontiers is relatively high, when scaling the release dates from factor r new,i = 1 (initial parameter) to r new,i = 0.75 compared to the distance between r new,i = 0.75 and r new,i = 0.5.This means that guaranteeing earlier release dates by the suppliers can only significantly reduce the earliness-tardiness costs up to a certain level.The weightdistance objective remains on a similar level for earlier release dates in instance DB100-1.Contrarily, the weightdistance objective increases by up to 8% with earlier release dates in instance DB100-2, when minimizing the earliness-tardiness costs.
Delaying the due dates also improves the earlinesstardiness costs, but if we only aim for this objective, the weight-distance increases by up to 60%.The negative impact on the weight-distance is stronger than stated for earlier release dates.The reduction of the earlinesstardiness costs becomes smaller with larger offsets, similar to earlier releases.The 3PL can still reduce the earliness-tardiness costs even if a solution is chosen that does not increase the weight-distance.When delaying the due dates, we also increase the planning horizon.This enables longer routes without causing expensive delays, and may reduce the earliness-tardiness costs.
Changes in the fleet composition include three different dimensions, the number of trucks m, the curb weight q C , and the load capacity Q.In Figures 6(c  we can see the results for four different fleet compositions.The tuple fleet A represents the initial fleet.In this case, the curb weight q C is higher than the load capacity, whereas this is the opposite in the other alternative fleet compositions.Additionally, the curb weight of all three alternatives (fleet B, C, D) tested is lower than in the initial case (fleet A).This explains why the weight-distance objective of the initial fleet composition is higher than in any of the alternatives.The total curb weight m • q C is the lowest for fleet A, the highest for fleet B, and fleets C and D are in between.For the total load capacity m • Q the relation is the opposite.Fleets C and D have the same total curb weight and load capacity.We state from Figures 6(c) and 7(c) that although the total curb weight is the lowest for fleet B, its Pareto frontier is dominated by fleet D. The Pareto frontier of fleet D also dominates fleet C. The most significant difference between both fleets is the number of trucks m = 12 for fleet C and m = 16 for fleet D. Hence, the results indicate that reducing the fleet size and aggregating the shipments in a smaller number of trucks is not a viable option for a JIT supplier network.Contrarily, lighter trucks with a smaller curb weight q C reduce the environmental impact of the transportation.
Consequently, the 3PL should prefer a larger fleet of smaller trucks from the environmental and JIT perspective.Earlier release dates at the suppliers and later due dates at the OEM have a direct and positive impact on the JIT perspective.On the other hand, a negative impact on the environment is observed if we only focus on JIT aspects.Alternatively, earlier release dates or later due dates can be considered to reduce CO 2 emissions while keeping the earliness-tardiness costs constant.

Conclusion
In this article, we study the TRP-JIT from an environmental perspective, i.e. the envTRP-JIT.We derive a bi-objective problem formulation and suitable solution procedures.The objectives are defined by the earlinesstardiness costs from the JIT perspective, and the weightdistance objective from the environmental point of view.We extend the MIP model and the LNS with local search introduced by Baals, Emde, and Turkensteen (2023) accordingly.Weighted sum and ε-constraint methods are used to estimate the Pareto frontier of the problem.The computational study aims to compare different variants of the proposed solution procedures, and gain managerial insights.The outcomes of this article are summarised by the following: • For small instances with up to n = 12 shipments, a commercial solver can be used within the bi-objective search schemes to derive solutions in reasonable computation time.While achieving an equivalent solution quality, the LNS requires a shorter computation time for these instances.• The ε-constraint methods find better estimates for the Pareto frontier than the weighted sum methods for the smaller instances but require a longer computation time.For the larger instances with n = 99 shipments, the weighted sum methods are preferable because the ε-constraint methods require a longer computational time, and they are inferior in the solution quality.• To a certain degree, we can trade-off between the earliness-tardiness costs from the JIT perspective, and the weight-distance objective from the environmental point of view without strongly increasing one of both objectives.The potential for improving the weightdistance objective is relatively small compared to the increase in earliness-tardiness costs.• Increasing the time span between release and due dates reduces the earliness-tardiness costs.Especially, when only considering JIT aspects and delaying the due dates, a negative effect on the weight-distance can be stated.Earlier release dates seem to have a smaller impact on the weight-distance.Alternatively, the savings from the JIT perspective can be spent on reducing the environmental costs.• The fleet composition affects both objectives.A smaller fleet size is not preferable for any objective, whereas trucks with a smaller curb weight improve the weight-distance objective.
In future research, the following aspects should be considered.First, by applying our LNS metaheuristic, we prioritise computational time over optimality.Hence, focusing on exact methods in this bi-objective scenario is a promising option.Second, we made the assumption that the speed of the vehicles does not significantly affect the objectives.In this context, the effects of congestion and other stochastic factors, e.g.urgent unplanned deliveries for the OEM, can be the scope of further studies.Third, the assumptions of a homogeneous fleet may not hold if many parameters of the vehicles are considered, e.g. the curb weight, speed, engine options, the drivers' behaviour, or fixed costs per vehicle.After evaluating the relevance of the parameters, a heterogeneous fleet may be considered in the problem formulation and solution procedure.

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

Notes on contributor
Julian Baals received a B.Sc. and M.Sc.degree in Engineering Management/ Industrial Engineering from the University of Technology in Darmstadt, Germany in 2016 and 2019 respectively.Between 2020 and 2022, he was enrolled as a Ph.D. fellow at Aarhus University, Denmark.Since 2022, he is continuing his Ph.D. studies at Friedrich Schiller University Jena, Germany.His focus is on just-in-time logistics especially the optimisation in transportation networks by developing metaheuristic solution procedures.
for the curb weight, i.e. the weight of any empty truck k ∈ M.

Figure 1 .
Figure1.The optimal solution for both objectives of an envTRP-JIT, including the relevant travel times and distances based on the example byBaals, Emde, and Turkensteen (2023).(a) Earliness-tardiness objective and (b) Weight-distance objective.
tuple containing the objective values for the solution variables in S, • L = (F, S) -a tuple containing both, the objective values and variables of a certain solution, • L = {L (1) , . . ., L (|L|) } -a set of objective values and variables of different solutions, • Lpareto = L (1) , . . ., L (| Lpareto |) -a permutation of objective values and variables of Pareto efficient solutions sorted by increasing earliness-tardiness and decreasing weight-distance derived from a certain set of solutions L.

Figure 2 .
Figure 2. Pareto frontiers derived from the different bi-objective search schemes (WM: weighted sum method, EC: ε-constraint method) and solution procedures.(a) Weighted sum methods and (b) ε-constraint methods.

Figure 5 .
Figure 5. Normalised and shifted Pareto frontier of instance DB100-1 generated by the weighted sum method with constant intervals and normalisation.(a) Complete range and (b) Improvement of f WD up to 15%.

Figure 6 .
Figure 6.Pareto frontiers derived from the different bi-objective search schemes (WM: weighted sum method, EC: ε-constraint method) and solution procedures of instance DB100-1.(a) Release dates r i .(b) Due dates p i and (c) Fleet composition (m, q C , Q).

Figure 7 .
Figure 7. Pareto frontiers derived from the different bi-objective search schemes (WM: weighted sum method, EC: ε-constraint method) and solution procedures of instance DB100-2.(a) Release dates r i .(b) Due dates p i and (c) Fleet composition (m, q C , Q).

Table 1 .
References considering routing problems closely related to the envTRP-JIT.

Table 2 .
Parameter values to define the example setup in Figure1.

Table 4 .
Change of the objective values between the best known solutions for both objectives relative to the best known solution for the respective objective.

•
We delay the due dates p new,i := p i + p off ∀i ∈ N by the offset p off := 0, 60, 120, 180, 240, such that the shipments have to be delivered later and the planning horizon gets longer.• We evaluate the effect of the fleet composition by setting the tuple of the fleet size, curb weight, and load capacity, (m, q C , Q), for fleet A to (16, 6350, 3650), fleet B to (8, 5200, 8300), fleet C to (12, 4800, 5200), and fleet D to (16, 3600, 3900).The fleets are defined such that the total weight and load capacity of fleet B corresponds to class 7, fleet C to class 6, and fleet D to class 5 according to the official classification in the United States (Transportation Research Board and National Research Council 2010).