The third party logistics provider freight management problem: a framework and deep reinforcement learning approach

Abstract

In many large manufacturing companies, freight management is handled by a third-party logistics (3PL) provider, thus allowing manufacturers and their suppliers to focus on the production of goods rather than managing their delivery. Provided their pivotal supply chain role, in this work we propose a general framework for what we term as “the 3PL freight management problem” (3PLFMP). Our framework identifies three primary activities involved in 3PL freight management: the assignment of orders to a fleet of vehicles, efficient routing of the fleet, and packing the assigned orders in vehicles. Furthermore, we provide a specific instantiation of the 3PLFMP that considers direct vs. consolidated shipping strategies, one dimensional packing constraints, and a fixed vehicle routing schedule. We solve this instantiated problem using several Reinforcement Learning (RL) methods, including Q-learning, Double Q-learning, SARSA, Deep Q-learning, and Double Deep Q-learning, comparing against two benchmark methods, a simulated annealing heuristic and a variable neighborhood descent algorithm. We evaluate the performance of these methods on two datasets. One is fully simulated and based on past work, while another is semi-simulated using real-world automobile manufacturers and part supplier locations, and is of our own design. We find that RL methods vastly outperform the benchmark heuristic methods on both datasets, thus establishing the superiority of RL methods in solving this highly complicated and stochastic problem.

Appendices

Appendix A Nomenclature

Notation	Definition
C	Set of customers (indexed by i)
$n_c, n_s, n_p, n_v, n_t,n_{ipj}$	Number of customers, suppliers, parts, vehicles, time periods, and flows
S	Set of suppliers (indexed by j)
P	Set of products (indexed by p)
V	Set of vehicles (indexed by v)
T	Set of time periods (indexed by t)
D	Demand domain (indexed by d)
$o_{ijp}^{(t)}$	Order of product p from the supplier j by customer i at time t
$t_{o^{(t)}_{ijp}} \text {(or} \, t_o)$	Due date of order $o_{ijp}^{(t)}$
$O^{(t)}$	Set of unassigned orders at time t
$\Sigma$	Set of consolidation strategies (indexed by $\sigma$)
$\Omega$	Set of solution strategies (indexed by $\omega$)
$\Psi$	Set of problem considerations
$u_i (u'_i)$	Weight (volume) of customer i’s order
$U_v (U'_v)$	Weight (volume) capacity of vehicle v
${\tilde{O}}^{(t)}$	Set of assigned orders to vehicles
$(l_p, w_p, h_p)$	Length, width, and height of product p
$(L_v, W_v, H_v)$	Length, width, and height of vehicle v’s container
${\hat{O}}$	Data structure containing the arrangements of orders in containers
$f(\cdot )$	Function capturing the overall transportation costs
$TC^{(t)}$	Transportation costs at time t
R	Reward function
${\mathcal {X}}$	The set of decision variables
${\mathbb {G}}$	The constraint set
$x_{ov}$	A binary variable that takes the value of 1 if order o is assigned to vehicle v, and 0 otherwise
$x_{iv} (x_{jv})$	A binary variable that takes the value of 1 if vehicle v visits customer i (supplier j) on it’s route, and 0 otherwise
$t_{vi} (t_{vj})$	The time that vehicle v next visits customer i (supplier j)
$V'$	Feasible vehicle/action set ($V'\subset V$) created by removing invalid actions/vehicles
$\theta$	RL method parameters
$\gamma$	Bellman equation discount parameter
$\pi$	Penalty function
$\beta$	Barrier function
$s \in {\mathcal {S}}$	State s and state space ${\mathcal {S}}$
$a \in {\mathcal {A}}$	Action a and action space ${\mathcal {A}}$
$\vartheta$	RL policy
${\tilde{v}}$	Direct shipment vehicle

Appendix B An expanded state definition

In real-world applications, there are challenges in defining the elements of an MDP because of infinite state space and computational and data efficiency concerns (Sutton & Barto, 2018). Therefore, the state representation needs to be defined in a way that contains the most salient interaction aspects between the agent and the environment. Therefore, the state representation should give the agent the proper domain knowledge and also make generalization possible, but must also necessarily be limited in scope (Sutton & Barto, 2018).

In light of this, we experimented with using more complicated state representations by including each vehicle’s current weight and/or volume, the day of the week, chronological time step (expressed as both a scalar and a binary subvector), and various combinations of these. Denote the state representation with additional features as $s'^{(t)}$ and the set of extended features as $y^{(t)}$. With this notation, we have $s'^{(t)} = s^{(t)} \oplus y^{(t)} = (ipjd)^{(t)} \oplus y^{(t)}$ where $\oplus$ is the concatenation operator, and $y^{(t)}$ is an arbitrary set of extended features, e.g., vehicles’ weight (at time t), vehicles volume (at time t), etc.

We experimented with $y^{(t)} \in Y= \{ \{\overrightarrow{W}^{(t)}\}, \{\overrightarrow{V}^{(t)}\}, \{t\}, \{\overrightarrow{W}^{(t)} \oplus \overrightarrow{V}^{(t)}\}, \{\overrightarrow{W}^{(t)} \oplus \overrightarrow{V}^{(t)} \oplus t \}, \{t_c\}, \{\overrightarrow{W}^{(t)} \oplus \overrightarrow{V}^{(t)} \oplus t_c \} \}$ where $\overrightarrow{W}$ and $\overrightarrow{V}$ are normalized vectors of vehicle weights and volumes, respectively, t is the chronological time step, and $t_c$ is the day of the week. We did not observe performance or convergence improvements by including any $y^{(t)} \in Y$ and thus elected to keep the state representation simple by excluding them.

Appendix C Discrete demand scenario details

One could alternatively adopt a discrete state space representation by discretizing the demand since customer i, part p, and supplier j are already discretely represented. Such a representation would permit tabular RL formulations and produce a less complex state space. While we did not experiment with such representations we wished to provide details on these should one choose to adopt our framework.

There are two discrete demand representations that one may wish to consider:

• Binary demand satisfaction: The demand of each ipj can only have two levels ($n_d = 2$), which is interpreted as either not satisfied or satisfied, i.e. $d\in \{0,1\}$. The maximum state space size would therefore be $\vert {\mathcal {S}}\vert = 2*n_{ipj} * n_t$ in this scenario.

• Ordinal demand bins: The demand of each ipj is generated with a countably finite number of levels, i.e., $1< n_d < \infty$ and such that the $n_d$ demand levels coincide with binned continuous demand levels, each having a governing lower and upper limit—e.g., $d^{(1)}_l < d^{(1)}_u \le d^{(2)}_l,\dots ,d^{(n_d)}_u$. This representation allows the incoming demands, and coinciding stochasticity, to be modeled in a discretized fashion. The maximum state space size is $\vert {\mathcal {S}} \vert = n_d *n_{ipj} * n_t$ in this scenario.

Appendix D CDF1 3PLFMP algorithms

In this section, we provide the specific details of the ASSIGN procedure for each of our adopted methods, Q-learning, Double Q-learning, SARSA, Simulated Annealing (SA) heuristic, and Variable Neighborhood Descent (VND) algorithm.

1.1 Q-learning algorithm

Algorithm 2 provides the ASSIGN Q-learning pseudocode for the CDF1 3PLFMP setting.

Algorithm 2

ASSIGN Q-learning Algorithm for the CDF1 3PLFMP.

1.2 Deep Q-learning algorithm

Algorithm 3 provides the pseudocode for ASSIGN using Deep Q-learning with the CDF1 3PLFMP setting.

Algorithm 3

ASSIGN Deep Q-learning Algorithm for the CDF1 3PLFMP.

1.3 Double Q-learning algorithm

The full ASSIGN pseudocode of the CDF1 3PLFMP Double Q-learning (DQL) method is shown in Algorithm 4.

Algorithm 4

ASSIGN Double Q-learning Algorithm for the CDF1 3PLFMP.

1.4 Double deep Q-learning algorithm for the CDF1 3PLFMP

Algorithm 5 shows the ASSIGN Double Deep Q-learning pseudocode for the CDF1 3PLFMP setting.

Algorithm 5

ASSIGN Double Deep Q-learning Algorithm for the CDF1 3PLFMP.

dummy

1.5 SARSA algorithm

The full pseudocode of SARSA adapted to the CDF1 3PLFMP setting is shown in Algorithm 6.

Algorithm 6

ASSIGN SARSA Algorithm for the CDF1 3PLFMP.

1.6 Benchmark simulated annealing and variable neighborhood descent heuristic details and algorithms

Heuristic methods are capable of finding reasonable solutions to large, stochastic optimization problems and have been extensively employed to handle NP-Hard problems, such as TSP, BPP, etc. Indeed, as shown in our literature review in Section 2, they are exclusively employed to handle pairings (i.e., two) of the three overarching problems that constitute the 3PLFMP. Therefore, such methods are suitable benchmarks against which we can compare our proposed RL methods. Furthermore, we can also assess how well heuristics perform on this specific 3PLFMP.

A popular and well-known heuristic method is Simulated Annealing (SA), introduced by Kirkpatrick et al. (1983). The method is directly motivated and derived from the metallurgy annealing process. The annealing process involves initially heating a substance to a high temperature, which is then gradually decreased. In metallurgy, this process allows the molecules in the metal to reach a crystalline (optimal) form, which takes the system from a starting state with a high level of internal energy to a state with minimum levels of energy. The SA heuristic applies this idea to optimization problems by moving from one solution state to another through an arbitrary perturbation algorithm, defined based on the problem at hand, with the goal of reaching an optimal objective function value. To reflect this process, SA has a temperature parameter that is initially set to a starting value that is decreased over the course of algorithmic execution using a “cooling schedule”. The temperature parameter affects the probability of accepting non-improving solutions and thus acts as a control on the exploration-exploitation trade-off. The temperature is initially high to encourage exploration early on. Like the RL methods, the SA algorithm is shown to converge to optimality under certain conditions (Bertsimas & Tsitsiklis, 1993).

The pseudocode of the SA algorithm for the CDF1 3PLFMP is shown in Algorithm 7. First, the parameters are initialized and the temperature is set to an initial value (line 2). Subsequently, the algorithm iterates over orders (line 3). To obtain an initial solution (state of the system), orders are assigned to vehicles randomly (line 6). Then a swap heuristic (the perturbation mechanism; lines 10 and 11) is used to create a neighborhood of the current solution at each inner iteration. If the new solution does not improve the reward, it may be set as the incumbent solution probabilistically relative to the current temperature (calculated in line 16). In the outer loop, the temperature is reduced using the cooling schedule (line 21; $T = T_0 \zeta ^ {it_{outer}}, (0< \zeta < 1)$). This procedure repeats until all orders are assigned.

Algorithm 7

Simulated Annealing for the 3PLFMP.

Another popular optimization algorithm is Variable Neighborhood Descent (VND), first proposed by Mladenović and Hansen (1997). This method works by iteratively exploring different neighborhoods around the current solution, where each neighborhood corresponds to a specific perturbation of an inputted solution. VND’s flexibility in exploring diverse solution spaces and its iterative refinement process make it a powerful method for solving optimization problems. The pseudocode of the VND algorithm for the CDF1 3PLFMP is shown in Algorithm 8 (Note that the method is technically Variable Neighborhood Ascent since we are maximizing the reward. This is however equivalent to minimizing the negative reward and we thus opt to maintain the original “descent“ nomenclature.). The algorithm begins by first initializing the parameters and neighborhood structures (lines 1–2). Subsequently, the algorithm iterates over orders (line 3). To obtain an initial solution, orders are assigned to vehicles randomly (line 6). Subsequently, swap and shift heuristics (the neighborhood structure function in our problem; line 11) are used to create alternative neighborhood solutions at each inner iteration. If there is an improvement in the reward (lines 15–16), the solution with the largest improvement (line 17) will be set as the incumbent solution (line 18). This procedure repeats until all orders are assigned.

Algorithm 8

Variable Neighborhood Descent for the 3PLFMP.

Appendix E Experiment details

In this section we provide some additional experiment details. In particular, we provide visual illustrations of the customers, suppliers, and depots of both datasets, describe how the flows used in our experiments are generated in greater detail, how orders are generated, and how vehicle routing schedule, weight, and volume capacities are determined. We also provide additional details on the parameterization the 3PLFMP framework, MDP formulation, and RL methods.

1.1 Dataset visualizations

A map showing the geolocated customers, suppliers, and depot of the fully-simulated dataset is shown in Figure 9.

Fig. 9

The locations of the customers, suppliers, and depot of the fully-simulated data

A map showing the geolocated customers, suppliers, and depot of the semi-simulated dataset is shown in Figure 10.

Fig. 10

The locations of the customers, suppliers, and depot of the semi-simulated data

1.2 Flow generation

The ipj flows of an instantiated problem are randomly determined. A user first specifies the number of customers ($n_c^{\prime } \le n_c$), parts ($n_p^{\prime } \le n_p$), and suppliers ($n_s^{\prime } \le n_s$) that should be in a given problem’s logistics network. To create the flows, the following procedure is carried out:

• 1. Create set $J^{'}$: for each supplier $j=1,\dots ,n_s$, select a random customer $i^{\prime } \sim {\mathcal {U}}\left( \{1,\dots ,n_c\}\right)$ and a random part $p^{\prime } \sim {\mathcal {U}}\left( \{1,\dots ,n_p\}\right)$; $J^{'} \leftarrow J^{'} \cup \{i^{\prime }p^{\prime }j\}$. Set $J^{'}$ includes all $n_s$ suppliers paired with random customers and products.

• 2. Create set $P^{'}$: for each part $p=1,\dots n_p$, select a random customer $i^{\prime } \sim {\mathcal {U}}\left( \{1,\dots ,n_c\}\right)$and a random supplier $j^{\prime } \sim {\mathcal {U}}\left( \{1,\dots ,n_s\}\right)$; $P^{'} \leftarrow P^{'} \cup \{i^{\prime }pj^{\prime }\}$. Set $P^{'}$ includes all $n_p$ parts paired with random customers and suppliers.

• 3. Create set $I^{'}$: for each customer $i=1,\dots ,n_c$, select a random part $p^{\prime } \sim {\mathcal {U}}\left( \{1,\dots ,n_p\}\right)$ and a random supplier $j^{\prime } \sim {\mathcal {U}}\left( \{1,\dots ,n_s\}\right)$; $I^{'} \leftarrow I^{'} \cup \{ip^{\prime }j^{\prime }\}$. Set $I^{'}$ includes all $n_c$ customers paired with random parts and suppliers.

• 4. $flows \leftarrow J^{'} \cup P^{'} \cup I^{'}$. The number of elements in this set is $n_{ipj} = \vert flows \vert$, representing the total number of flows used in the problem instance.

Depending on whether the flows are being generated from the semi- or fully-simulated data, the customer and supplier coordinates are selected from the dataset, starting from the first customer and supplier. This ensures that when the problem size increases, the customers and suppliers in previous (i.e., smaller) instances are also in the larger-sized instance.

The following table provides the number of customers $n_c^{'}$, suppliers $n_s^{'}$, and parts $n_p^{'}$ that produce the total number of flows $n_{ipj}$ for each problem instance:

1.3 Order generation

To generate the orders for each episode, we assign each $ipj \in flows$, corresponding to an eligible order $o_{ipj}$, a random arrival time $t_o^{'} \sim {\mathcal {U}}(\{1,\dots ,\left\lfloor .5*n_t \right\rfloor \})$. Based on this arrival time we then randomly generate an order due date $t_o = \sim {\mathcal {U}}(\{t_o^{'}+1,\dots ,t_o^{'}+\left\lfloor .4*n_t \right\rfloor \})$.

We also randomly generate a demand for each ipj at each episode. This demand represents the amount of part p from supplier j that is in demand by customer i at the time $t_o^{'}$ the ipj order $o_{ipj}$ is placed. Demand $d_{ipj}$ is drawn randomly from a normal distribution as follows:

$$\begin{aligned} d_{ipj} \sim {\mathcal {N}}\left( \tau _{ipj},\rho _{ipj}\right) , \end{aligned}$$

(E1)

where $\tau _{d}$ and $\rho _{d}$ are determined by $\tau _{d} \sim {\mathcal {U}}(\{600,\dots ,4000\})$ and $\rho _{d} = \frac{\tau _{d}}{6}$. This procedure produces a fair amount of stochasticity in the demand-generating process since the average of the demand-producing distribution has a large and random range of mean values, as well as a large and varying standard deviation, which is dependent upon the random mean selected.

1.4 Vehicle parameterization

For each vehicle $v \in V$ we need to generate (simulate) weight capacity, volume capacity, days available (in a week), fixed routing schedule, starting location, transportation cost, and transportation time. These properties, however, must be generated such that the final problem instance is feasible. For example, if we set very low capacities the vehicles may not be able to support the generated ipj demands.

Therefore, we set the vehicle volume and weight capacities by summing over the mean ipj demands, denoted $\tau$, plus one standard deviation, denoted $\rho$, and dividing by the number of vehicles. Furthermore, we multiply by an expansion factor, $\eta (\eta \ge 1)$ to alter capacity constraint tightness. Thus, assuming each product has one unit of weight (volume), the weight (volume) capacity for each vehicle is $U_v (U^{\prime }_v) = \frac{\sum _{ipj \in flows}{(\tau _{ipj} + \rho _{ipj})}*\eta }{n_v}$. The weight and volume capacities of the direct shipping vehicle are set to $\infty$.

The fixed routing schedule for each vehicle is set by selecting a number of suppliers $n_s^{(v)} \sim {\mathcal {U}}\left( \{2,\dots ,\frac{n_s^{\prime }}{2}\}\right)$ and then randomly selecting this number of suppliers for inclusion on v’s route – i.e., $\{Route_v \leftarrow Route_v \cup \{j\}: j \sim {\mathcal {U}}\left( \{1,\dots ,n_s^{\prime }\}_{\setminus Route_v}\right) ,1,\dots ,n_s^{(v)}\}$. All customers i associated with any supplier j through an $ipj \in flows$ are then added to the vehicle’s route to ensure there are no spurious supplier-only visits (i.e., a vehicle would not visit a supplier without there being a corresponding customer to deliver to). Appending the customers to the end of the route also ensures that suppliers are visited before their customers. The direct shipment vehicle is assumed to visit all customers and suppliers. The starting location of all vehicles is initialized to a depot, located at the centroid of the geocoded customers and suppliers for the semi-simulated dataset; the depot is already given in the fully-simulated dataset.

The transportation distance for each leg of vehicle’s tour is calculated using the Euclidean metric. By setting the average speed for each vehicle, the transportation time can also be calculated. It is worth noting that when orders are delivered using direct shipment, the transportation time is smaller, but more vehicles are being used so the total transportation distance (and therefore cost) is higher.

1.5 Parameterization

In this subsection, we provide the details on how our 3PLFMP framework is parameterized, as well as the parameterization of our adopted RL methods.

1.5.1 Framework and MDP parameterizations

We specify a variety of framework- and MDP-specific parameters in this subsection. These specifications encompass both the broader 3PLFMP framework, disclosed in Sect. 3, as well as the specific implementation we adopted for the CDF1 3PLFMP.

First, we note that we incorporate the constraints ${\mathbb {G}}({S}^{(t)})$ into our reward function, through both a barrier term and a penalty term, and also incorporate these into $\Psi$, although we do not explicitly enforce them. The $\Psi$ constraints are, however, used to inform vehicle eligibility when assigning orders. We set the exponent of the penalty term (Eq. (13)) $\xi =10$. The following are specific framework component parameterizations:

• $\mathbf {ASSIGN_\sigma ^\omega (\cdot )}$: As we mentioned, our specific implementation of the 3PLFMP explores the use of an order consolidation strategy, which is what $\sigma$ is set equal to. On the other hand, the solution strategy $\omega \in \{\text {Q-learning}, \text {Double Q-learning}, \text {SARSA}, \text {SA}\}$.

• $\mathbf {ROUTE_\sigma ^\omega (\cdot )}$: We set $\sigma =fixed$, $\omega = \emptyset$, because the vehicle routes are fixed. The procedure therefore returns the fixed routing schedule. This pre-determined routing schedule is one of the elements included in $\Psi$.

• $\mathbf {PACK_\sigma ^\omega (\cdot )}$: We set $\sigma =1D$, $\omega = \emptyset$ since we only consider one-dimensional weight and volume constraints. That is, an order is considered packed at the time of assignment so long as the addition of the order does not cause the vehicle’s weight and volume capacities to be exceed. These weight and volume constraints are included in $\Psi$.

1.5.2 Method parameterization

There are a few parameters associated with each of our methods that must be specified. Therefore, we perform a variety of parameter tuning experiments and report the results in Table 8, below. Table 8 shows the average of average rewards over problem instances 01–07 (see Table 7) for different parameter combinations for each RL method on the semi-simulated data. We observe that parameter setting 0 gives the best reward for all RL algorithms (we thus also adopt this parameterization for the fully simulated data as well). Accordingly, we parameterize the three RL algorithms with a future discount reward factor of $\gamma = 0.8$, a learning parameter $\alpha = 0.001$, a batch size of 512, and $\epsilon =0.90$ for $\epsilon$-greedy selection.

Table 7 The $n^{\prime }_c$, $n^{\prime }_p$, and $n^{\prime }_s$ values used to generate the $n_{ipj}$ flows of each problem instance

Table 8 Parameter tuning of RL methods on the semi-simulated data

Further, note that to improve the stability of the RL methods, we employ reward clipping (see, e.g., van Hasselt et al., 2016; Sutton & Barto, 2018). We clipped the rewards between $[-8000,0]$ and then normalized the value to a range $[-1,0]$ when computing the loss and updating the parameters. In addition to the above-mentioned parameters for all RL algorithms, when training the two deep RL algorithms, we used stochastic gradient descent due to its efficiency and a hard update scheme with a fixed target network, updated every 50 iterations.

Fig. 11

3PLFMP subproblem benchmark results. Semi-simulated dataset. $n_{ipj} = 53$. Equation 12 reward function moving average (window = 50). a Baseline, b No VRP c No BPP, and d No VRP + BPP

Appendix F 3PLFMP subproblem benchmark experiments

In this section we conduct additional experiments that showcase the importance of each component of the 3PLFMP, thus highlighting the inter-related nature of the three combinatorial subproblems we identify. Doing so creates a benchmark between what we propose and past works that consider only a single subproblem or pairings of these subproblems (we review these works in Section 2.2).

To conduct these benchmark experiments we remove certain elements from Eq. 12 to simulate the removal of a certain combinatorial subproblem component of the 3PLFMP. We consider the following cases:

• Baseline: Original 3PLFMP (no alteration to Eq. 12).

• No VRP: Remove the routing penalties.

• No BPP: Remove the packing penalties.

• No VRP + BPP: Removing both the routing and packing penalties.

For each of the modified reward functions above we then train the various RL methods and plot the results using the original reward function to assess whether the method can still account for the intricacies of the full 3PLFMP problem setting. We opt to use a problem size of $n_{ipj}=53$ and the semi-simulated dataset.

The results of these benchmark experiments are provided in Fig. 11.

We observe that removal of any of the subproblem components (VRP, BPP, VRP+BPP) prevents all RL methods from appropriately accounting for the removed component; all averaged rewards are flat and unable to improve their solution quality. This result showcases the importance of considering all three subproblem elements that make up the 3PLFMP framework.

Appendix G Robustness to demand distribution

Prior researches used various continuous and discrete distributions to model uncertainty in demand (see, e.g., the survey by Oyola et al., 2018). As described in Section E, our obtained results are generated using a normal distribution, which is widely used in the literature of stochastic VRP. However, to test the robustness of our framework, in this section we explore Poisson as a discrete distribution, another common VRP distribution - e.g., in Mendoza et al. (2016) and Sluijk et al. (2023). The experiments will assess the effects of this change on our results, and observe the robustness of our proposed framework. Note that we set the mean values of the Poisson distribution similar to the mean values mentioned in Section E for normal distribution. Table 9 summarizes the percentage of times each algorithm wins, and Table 10 shows the results of semi- and fully-simulated test instances along with the results of the two comparator algorithms. Looking at both tables, we observe that in all of the larger instances, RL algorithms perform the best compared to SA and VND algorithms. VND algorithm could provide the best solution in only one small instance (instance #2) but it is dominated in all other instances by RL algorithms. However, SA algorithm never outperformed other algorithms in the experiments. Thus, we find that the results obtained using a Poisson distribution are comparable to those obtained used a normal distribution.

Table 9 Percentage of semi- and fully-simulated problem instances that each algorithm performs the best based on average rewards (demand is Poisson-distributed)

Table 10 Semi-simulated (first 12 rows) and fully-simulated (second 12 rows) test results for Poisson demand