Online Optimization of Pickup and Delivery Problem Considering Feasibility

Abstract

A pickup and delivery problem by multiple agents has many applications, such as food delivery service and disaster rescue. In this problem, there are cases where fuels must be considered (e.g., the case of using drones as agents). In addition, there are cases where demand forecasting should be considered (e.g., the case where a large number of orders are carried by a small number of agents). In this paper, we consider an online pickup and delivery problem considering fuel and demand forecasting. First, the pickup and delivery problem with fuel constraints is formulated. The information on demand forecasting is included in the cost function. Based on the orders, the agents’ paths (e.g., the paths from stores to customers) are calculated. We suppose that the target area is given by an undirected graph. Using a given graph, several constraints such as the moves and fuels of the agents are introduced. This problem is reduced to a mixed integer linear programming (MILP) problem. Next, in online optimization, the MILP problem is solved depending on the acceptance of orders. Owing to new orders, the calculated future paths may be changed. Finally, by using a numerical example, we present the effectiveness of the proposed method.

1. Introduction

There has been much attention paid to control technologies for realizing a smart city [1]. In a smart city, many services using control technologies should be developed, such as transportation, energy distribution, healthcare, environmental monitoring, business, commerce, emergency response, and social activities. Moreover, in control technologies for a smart city, a cyber–physical system (CPS) plays an important role. A CPS is a system where physical and information components are deeply connected through a communication network. In a smart city, multiple physical components such as mobile robots are controlled by information systems, such as cloud servers. Hence, a plant in a smart city is modeled as a CPS. Thus, in a smart city and a CPS, it is important to develop a control method for multiple agents. In this paper, we focus on a control method for multiple agents.

Control of multiple agents has several applications and has been widely studied. In the control of multiple agents, there are several problem settings, such as the achievement of cooperative tasks and surveillance (patrol). A control specification in cooperative tasks may be given by linear temporal logic formulas [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]. Temporal logic is a logical system of rules and symbolism for representing, and reasoning about, propositions qualified in terms of time (for example, “I will eventually be hungry”). The surveillance problem is to find agents’ paths patrolling a given area as evenly as possible [8,17,18,19,20,21,22,23,24,25,26,27,28]. In problem settings on the control of multiple agents, model predictive control (MPC) is frequently used (see, e.g., [9,22,27]). MPC is a control method in which the control input is generated by solving at each time the finite-time optimal control problem (see, e.g., [29,30]). In other words, the control input is calculated based on the prediction using a mathematical model. MPC was frequently used in the control of chemical plants. In recent years, with the development of computers, MPC has been used in the control of several systems, such as automobiles. In the case where a target area is unknown and is changed in time, a model-free approach such as machine learning is effective (see, e.g., [31,32]). In this paper, we suppose that a target area is fixed and use a model-based approach.

In this paper, we focus on the pickup and delivery problem as one of the control problems of multiple agents. In recent years, the food delivery market has been expanding through the coronavirus epidemic. Efficient delivery from the restaurant to the customer is becoming more and more important. Such a problem of efficiently delivering goods received at a certain location to the desired location is called a pickup and delivery problem (see, e.g., [33,34,35,36,37,38,39,40,41]). Especially, the problem of updating routes by solving optimization problems at regular intervals is called an online pickup and delivery problem, which has been studied extensively in recent years (see, e.g., [42,43,44]).

In pickup and delivery problems, it is important to use electric vehicles and drones as agents. Electric vehicles and drones have the disadvantage of severe battery constraints. Hence, fuel constraints must be considered in pickup and delivery problems. In [45,46], a modeling method to handle the fueling of vehicles has been proposed. A fuel constraint in the proposed method is basically the same as in [45,46]. Refueling is different. In this paper, the fuel of a vehicle immediately becomes full by changing the battery. Furthermore, when electric vehicles and drones are used, demand forecasting is important (see, e.g., [47,48]). This is because more efficient moves are needed due to battery constraints. For example, many agents can be placed in advance at points where demand is concentrated. To the best of our knowledge, few results have been obtained.

In this paper, we propose an online pickup and delivery problem in which demand forecasting and fuel constraints are considered. First, we formulate the pickup and delivery problem. In this paper, the target area is given by an undirected graph. Constraints including fuel constraints are explained. Next, after the demand forecast is explained, the cost function is introduced. In a certain term in the cost function, we use the result of demand forecasting. Then, agents’ paths can be calculated according to the result of demand forecasting. That is, based on the policy of MPC, the real order and demand forecast are used. Even if the result of demand forecasting is different from real orders, the problem can be solved. This is because the result of demand forecasting is used in the cost function, not a constraint. Hence, we can guarantee the feasibility. Third, we explain the procedure of online optimization. Finally, we present a numerical example. Through a numerical example, we demonstrate that the feasibility is guaranteed by the proposed method.

The conference paper [49] is a preliminary version of this paper. In [49], the authors have considered both demand forecasting and fuel constraints. However, depending on the result of the demand forecasting, the optimization problem becomes infeasible. In this paper, this technical issue is overcome.

This paper is organized as follows. In Section 2, we formulate the pickup and delivery problem studied in this paper. After a mathematical model of the target area and the definition of the orders are explained, the constraints considered in this paper are explained. Demand forecasting in this paper is also explained. Finally, the cost function is explained. In Section 3, we explain the procedure of online optimization of the pickup and delivery problem. In Section 4, we present a numerical example to show the effectiveness of the proposed method. In Section 5, we conclude this paper.

The main contributions of this paper are highlighted as follows:

(i)

For pickup and delivery problems, fuel constraints are introduced to utilize electric vehicles and drones as agents.

(ii)

A pickup and delivery problem considering demand forecasting is newly formulated.

(iii)

The effectiveness of the proposed method is clarified through a numerical example.

2. Pickup and Delivery Problem

In this section, we formulate the pickup and delivery problem studied in this paper. First, we describe the preliminaries of the problem, which are the target area and the order notation. Next, we describe the constraints in the pickup and delivery problem. Finally, we describe the cost function in the pickup and delivery problem.

2.1. Preliminaries

We introduce some notations.

First, the target area in the pickup and delivery problem is given by the following undirected graph:

$$\mathcal{G}=(\mathcal{V},\mathcal{E}),$$

where there are the following:

• $\mathcal{V}=\mathcal{S}\cup \mathcal{C}\cup \mathcal{R}$: the set of vertices ($\mathcal{S}\cap \mathcal{C}\cap \mathcal{R}=\varnothing$);
• $\mathcal{S}=\{s_1,s_2,\dots ,s_N\}$: the set of stores;
• $\mathcal{C}=\{c_1,c_2,\dots ,c_M\}$: the set of customers;
• $\mathcal{R}=\{r_1,r_2,\dots ,r_J\}$: the set of intermediate points between vertices;
• $\mathcal{E}\subseteq \mathcal{V}\times \mathcal{V}$: the set of edges.

Without the loss of generality, the set $\mathcal{V}$ can be represented by

$$\mathcal{V}=\{\underset{\mathcal{S}}{\underbrace{1,2,\dots ,N}},\underset{\mathcal{C}}{\underbrace{N+1,N+2,\dots ,N+M}},\underset{\mathcal{R}}{\underbrace{N+M+1,\dots ,N+M+J}}\}.$$

(1)

Let $\mathcal{A}=[a_{ij}]\in {\{0,1\}}^{(N+M+J)\times (N+M+J)}$ denote the adjacency matrix of $\mathcal{G}$, where if there is the edge between $v_i$ and $v_j$, then $a_{ij}=1$; otherwise, $a_{ij}=0$. An example of an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ is shown in Figure 1. We suppose that according to a given graph, an agent moves from some vertex to other one in unit time.

Next, notations about orders are introduced. We assume that new orders are received sequentially. Let ${\tau }_i$, $i=0,1,2,\dots$ denote the time that the delivery route may be changed. The initial update time ${\tau }_0$ is set as ${\tau }_0=0$. We assume that at the update time ${\tau }_i$, there is at least one new order. Let ${\mathcal{O}}^i$ denote the set of orders that are newly added at time $t={\tau }_i$. The element $o_k$ of ${\mathcal{O}}^i$ is given by the following form:

$$o_k=(s_k,c_k,T_{o_k}^s,T_{o_k}^c),$$

(2)

where there are the following:

• $s_k\in \mathcal{S}$: the store that received the order;
• $c_k\in \mathcal{C}$: the customer who placed the order;
• $T_{o_k}^s$: the time that the goods are available for pickup;
• $T_{o_k}^c$: the time limit to complete the delivery.

2.2. Constraints

We formulate the following constraints in the pickup and delivery problem at the update time ${\tau }_i$.

• Constraints on Agent Location: an agent can be located at only one vertex.
• Constraints on Agent Movement: an agent can move according to a given graph.
• Constraints on Orders: an order is assigned to an agent.
• Constraints on the Set of Orders: an agent must receive the goods before delivery and must deliver those until a certain time.
• Constraints on Agent Fuel: the fuel remaining for an agent must be equal to or greater than a certain value.
• Constraints on Luggage: the number of goods that an agent can handle once is constrained.

We explain the details of these constraints.

2.2.1. Constraints on Agent Location

First, let $\mathcal{N}$ denote the finite set of agents. We introduce a binary variable vector $p_l\left(t\right)\in {\{0,1\}}^{\left|\mathcal{V}\right|}$ for the agent $l\in \mathcal{N}$, where $t\in \{0,1,2,\dots ,\}$ is discrete time. The j-th element $p_{l,j}\left(t\right)$ of $p_l\left(t\right)$ is 1 if the agent l is located at the vertex $v_j$; otherwise, it is 0. The constraint on the location of the agent l is given by

$$\begin{array}{c}\hfill {\displaystyle \sum _{j=1}^{\left|\mathcal{V}\right|}}{p}_{l,j}\left(t\right)=1,l\in \mathcal{N},t\in \{{\tau }_i,\dots ,{\tau }_i+T_f\},\end{array}$$

(3)

which implies that at time t, the agent l is located at only one vertex.

2.2.2. Constraints on Agent Movement

The movement of the agent l is constrained by the following inequality:

$$\begin{array}{c}\hfill p_l\left(t\right)\le \mathcal{A}{p}_l(t+1),t\in \{{\tau }_i,\dots ,{\tau }_i+T_f-1\},\end{array}$$

(4)

where the inequality sign (≤) holds element-wise. Equation (4) represents that the agent l moves only between vertices for which edges exist (see, e.g., [50,51,52]). We remark that (4) becomes the constraint on the agent movement when (3) is imposed.

For example, consider the undirected graph shown in Figure 1. Suppose that at time t, the agent l is located at the vertex $s_1$. Then, the location candidates at time $t+1$ are the vertices $s_1$, $r_1$, $r_2$, and $r_5$. Based on the notation of (1), the vertices $s_1$, $r_1$, $r_2$, and $r_5$ correspond to 1, 5, 6, and 9, respectively. Because the first row of the adjacency matrix $\mathcal{A}$ is given by

$$\left[\begin{array}{cccccccccc}1& 0& 0& 0& 1& 1& 0& 0& 1& 0\end{array}\right],$$

the above situation can be represented by the following inequality:

$$p_{l,1}\left(t\right)\le p_{l,1}(t+1)+p_{l,5}(t+1)+p_{l,6}(t+1)+p_{l,9}(t+1).$$

When at time t the agent l is located at the vertex $s_1$, $p_{l,1}\left(t\right)=1$ holds. We remark that from (3), $p_{l,j}\left(t\right)=0$, $j\ne 1$ holds. Hence, under (3), either $p_{l,1}(t+1)$, $p_{l,5}(t+1)$, $p_{l,6}(t+1)$, or $p_{l,9}(t+1)$ must be 1. This implies that the location candidates at time $t+1$ are the vertices $s_1$, $r_1$, $r_2$, and $r_5$.

2.2.3. Constraints on Orders

For the agent l and the order $o_k$, we introduce a binary variable $y_{l,k}$ as follows:

$$y_{l,k}=\left\{\begin{array}{cc}1\hfill & \mathrm{if} \mathrm{the} \mathrm{agent} l \mathrm{is} \mathrm{responsible} \mathrm{for} \mathrm{the} \mathrm{order} o_k,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

As variables related to the agent l’s order receipt and delivery at time t, we also introduce binary variables $z_{l,k}^{\mathrm{pick}}\left(t\right)$ and $z_{l,k}^{\mathrm{deli}}\left(t\right)\in \{0,1\}$ as follows:

$$\begin{array}{cc}{z}_{l,k}^{\mathrm{pick}}\left(t\right)& =\left\{\begin{array}{cc}1\hfill & \mathrm{if} \mathrm{the} \mathrm{agent} l \mathrm{receives} \mathrm{the} \mathrm{goods} \mathrm{of} \mathrm{the} \mathrm{order} o_k \mathrm{at} \mathrm{time} t,\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.\hfill \\ z_{l,k}^{\mathrm{deli}}\left(t\right)& =\left\{\begin{array}{cc}1\hfill & \mathrm{if} \mathrm{the} \mathrm{agent} l \mathrm{delivers} \mathrm{the} \mathrm{goods} \mathrm{of} \mathrm{the} \mathrm{order} o_k \mathrm{at} \mathrm{time} t,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

Assume that the vertex $s_k$ in the order $o_k$ corresponds to the vertex $v_i$ (i.e., the i-th element $p_{l,i}\left(t\right)$ of $p_l\left(t\right)$). Then, receiving the goods of the order $o_k$ at time t can be represented by the following expression:

$$\begin{array}{c}\hfill [z_{l,k}^{\mathrm{pick}}\left(t\right)=1]\Rightarrow [y_{l,k}=1]\wedge [p_{l,i}\left(t\right)=1],\forall l,\forall t.\end{array}$$

(5)

Similarly, if $c_k$ corresponds to the i-th element $p_{l,i}\left(t\right)$ of $p_l\left(t\right)$, then delivering the goods of order $o_k$ at time t can be represented by the following expression:

$$\begin{array}{c}\hfill [z_{l,k}^{\mathrm{deli}}\left(t\right)=1]\Rightarrow [y_{l,k}=1]\wedge [p_{l,i}\left(t\right)=1],\forall l,\forall t.\end{array}$$

(6)

The relationship between $y_{l,k}$, $z_{l,k}^{\mathrm{pick}}\left(t\right)$, and $z_{l,k}^{\mathrm{deli}}\left(t\right)$ is represented by the following expressions:

$$\begin{array}{cc}\hfill [y_{l,k}=1]& \iff {\displaystyle \sum _{t=T_{o_k}^s}^{{\tau }_i+T_f}}{z}_{l,k}^{\mathrm{pick}}\left(t\right),\forall k\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill [y_{l,k}=1]& \iff {\displaystyle \sum _{t=T_{o_k}^s}^{{\tau }_i+T_f}}{z}_{l,k}^{\mathrm{deli}}\left(t\right),\forall k\hfill \end{array}$$

(8)

We define two sets of orders $\mathcal{P}$ and $\mathcal{D}$ to represent the status of each order. The set $\mathcal{P}$ is the set of orders where the agent has not received the goods at the store. The set $\mathcal{D}$ is the set of orders where the agent has received the goods at the store and is in the process of delivery. Here, because there is only one agent responsible for each order, the following constraints must be imposed:

$$\begin{array}{c}\hfill {\displaystyle \sum _l}{y}_{l,k}=1,\forall o_k\in \mathcal{P}\cup \mathcal{D}.\end{array}$$

(9)

2.2.4. Constraints on the Set of Orders

The constraint on $\mathcal{P}$ is that the agent receives the goods before delivery and that the delivery time of the goods is equal to or earlier than the time limit ${\tau }_i+T_f$. This constraint is represented by the following expressions:

$$\begin{array}{c}{\displaystyle \sum _l}{\displaystyle \sum _{t={\tau }_i}^{t_1}}{z}_{l,k}^{\mathrm{pick}}\left(t\right)-{\displaystyle \sum _l}{\displaystyle \sum _{t={\tau }_i}^{t_1}}{z}_{l,k}^{\mathrm{deli}}\left(t\right)\ge 0,t_1\in \{{\tau }_i,\dots ,{\tau }_i+T_f\},\forall o_k\in \mathcal{P},\hfill \end{array}$$

(10)

$$\begin{array}{c}{\displaystyle \sum _l}{\displaystyle \sum _{t={\tau }_i}^{{\tau }_i+T_f}}{z}_{l,k}^{\mathrm{deli}}\left(t\right)=1,\forall o_k\in \mathcal{P}.\hfill \end{array}$$

(11)

On the other hand, the constraint on $\mathcal{D}$ is that the time to deliver the goods must be equal to or earlier than the time limit ${\tau }_i+T_f$. Hence, the constraint on $\mathcal{D}$ is expressed by the following expression:

$$\begin{array}{c}\hfill {\displaystyle \sum _l}{\displaystyle \sum _{t={\tau }_i}^{{\tau }_i+T_f}}{z}_{l,k}^{\mathrm{deli}}\left(t\right)=1,\forall o_k\in \mathcal{D}.\end{array}$$

(12)

2.2.5. Constraints on Agent Fuel

When agents are implemented by using electric vehicles and drones, we need to consider fuel constraints.

We introduce the continuous variable $q_l\left(t\right)\in \mathcal{R}$ as a variable representing the fuel of the agent l at time t. Let $q_0$ denote the agent’s maximum fuel. Suppose that the initial fuel is given by $q_l\left(0\right)=q_0$. The fuel $q_l\left(t\right)$ at time t is defined by

$$\begin{array}{c}\hfill q_l\left(t\right)=\left\{\begin{array}{cc}{q}_0\hfill & \mathrm{if} \mathrm{the} \mathrm{agent} l \mathrm{is} \mathrm{located} \mathrm{at} \mathrm{the} \mathrm{vertex} \mathrm{in} \mathcal{S},\hfill \\ q_l(t-1)-1\hfill & \mathrm{otherwise},\hfill \end{array}\right.\end{array}$$

(13)

which implies that the fuel is fully charged at the store vertex and decreases by one at the vertex, except for the store vertex.

The constraint on the fuel of the agent l is given by

$$\begin{array}{cc}\hfill & q_l\left(t\right)\ge c,t\in \{{\tau }_i,\dots ,{\tau }_i+T_f-1\}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & q_l({\tau }_i+T_f)\ge c_f,\hfill \end{array}$$

(15)

where $c,c_f\ge 0$ are given constants. According to the policy of model predictive control, the terminal constraint may be different with other constraints. This difference can be realized via a method for choosing c and $c_f$.

2.2.6. Constraints on Luggage

We introduce a variable $X_l\left(t\right)$ that represents the number of pieces of luggage carried by agent l at time t. Here, we assume that the weights of the packages are the same. Using $z_{l,k}^{\mathrm{pick}}\left(t\right)$ and $z_{l,k}^{\mathrm{deli}}\left(t\right)$, the variable $X_l\left({\tau }_i\right)$ can then be defined as

$$\begin{array}{c}\hfill X_l\left({\tau }_i\right)=\left\{\begin{array}{cc}{\displaystyle \sum _k(z_{l,k}^{\mathrm{pick}}\left(0\right)-z_{l,k}^{\mathrm{deli}}\left(0\right))}\hfill & \mathrm{if}i=0,\hfill \\ {\left(X_l\left({\tau }_i\right)\right)}_{i-1}\hfill & \mathrm{otherwise} (i\ge 1),\hfill \end{array}\right.\end{array}$$

(16)

where ${\left(X_l\left({\tau }_i\right)\right)}_{i-1}$ is the value obtained via optimization at $t={\tau }_{i-1}$. In addition, $X_l\left(t\right)$ can be expressed as

$$\begin{array}{c}\hfill X_l\left(t\right)=X_l(t-1)+{\displaystyle \sum _k}(z_{l,k}^{\mathrm{pick}}\left(t\right)-z_{l,k}^{\mathrm{deli}}\left(t\right)),\end{array}$$

(17)

when $t={\tau }_i+1,\dots ,{\tau }_i+T_f$. Assume that the agent l can carry up to $c_w$ pieces of luggage, where $c_w$ is a given constant. The constraint on the number of pieces of luggage that the agent l carries at time t can be expressed as

$$\begin{array}{c}\hfill X_l\left(t\right)\le c_w,\forall l,\forall t.\end{array}$$

(18)

The variable $X_l\left(t\right)$ implies the number of pieces of luggage. In the case where the weights of the packages are different, the variable $X_l\left(t\right)$ is regarded as a total weight for the agent l via a simple modification.

2.3. Demand Forecast

In the conventional pickup and delivery problem, the paths of the agents are calculated at ${\tau }_i,{\tau }_{i+1},\dots ,{max}_k{T}_{o_k}^c$ (${max}_k{T}_{o_k}^c$ represents the slowest time limit in the set of orders). In this paper, we consider using the demand forecast. It is expected that better paths are obtained by considering not only the time interval where there are orders but also the future where there are no orders.

In this paper, we assume that the number of agents (denoted by $d_j\left(t\right)$) required in the store $s_j$ at time t is given as the result of the demand forecast. In the next subsection, we explain how to use $d_j\left(t\right)$.

2.4. Problem Formulation

Under the above preparations, we formulate the pickup and delivery problem at the update time ${\tau }_i$. The cost function to be minimized is defined as follows:

$$\begin{array}{c}\hfill J_i=c_{i,1}+c_{i,2}+c_{i,3}.\end{array}$$

(19)

We describe $c_{i,1}$, $c_{i,2}$, and $c_{i,3}$ in detail. At first, $c_{i,1}$ is defined as follows:

$$\begin{array}{c}\hfill c_{i,1}=w_1{\displaystyle \sum _{t={\tau }_i}^{{\tau }_i+T_f}\sum _{l\in \mathcal{N}}\sum _{j\in \mathcal{C}\cup \mathcal{R}}}{p}_{l,j}\left(t\right).\end{array}$$

(20)

$p_l\left(t\right)\in {\{0,1\}}^{\left|\mathcal{V}\right|}$ is a binary variable vector for the agent $l\in \mathcal{N}$, where $t\in \{0,1,2,\dots ,\}$ is discrete time. The j-th element $p_{l,j}\left(t\right)$ of $p_l\left(t\right)$ is 1 if the agent l is located at the vertex $v_j$; otherwise, it is 0. The scalar $w_1$ is a constant representing the weights and is determined through trial and error. In other words, Equation (20) implies the sum of the time spent at the vertices, except for the store for all agents. By reducing the value of this term, the running cost of the agents can be reduced.

Next, $c_{i,2}$ is defined as follows:

$$\begin{array}{c}\hfill c_{i,2}=w_2{\displaystyle \sum _k}{\alpha }_k(T_k-T_{o_k}^c),\end{array}$$

(21)

where $T_k$ is a variable that represents the time when the agent completes the delivery of order $o_k$ to the customer. Time $T_k$ can be expressed as

$$\begin{array}{c}\hfill T_k={\displaystyle \sum _l}[T_{o_k}^c,...,{\tau }_i+T_f]\left[\begin{array}{c}{z}_{l,k}^{\mathrm{deli}}\left(T_{o_k}^c\right)\\ ⋮\\ z_{l,k}^{\mathrm{deli}}({\tau }_i+T_f)\end{array}\right],\end{array}$$

(22)

where $z_{l,k}^{\mathrm{deli}}\left(t\right)$ is a binary variable that is 1 if agent l delivers order $o_k$ at time t and 0 if it does not. In (22), the column vectors are all zero, except for one element. Therefore, from this expression, we can find $T_k$ when the value of $z_{l,k}^{\mathrm{deli}}\left(t\right)$ is equal to 1. The scalar ${\alpha }_k$ is a variable that determines whether the delivery time of the order $o_k$ to be delivered meets the delivery deadline. The ${\alpha }_k$ is defined as follows:

$$\begin{array}{c}\hfill {\alpha }_k=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{T}_k\ge T_{o_k}^c,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

(23)

where ${\alpha }_k=1$ implies that the time $T_k$ when the order is delivered to the customer is later than the time limit $T_{o_k}^c$ to complete the delivery. The scalar $w_2$ is a constant representing the weights and is determined through trial and error.

In other words, Equation (21) is a summation of the delays in the delivery time. By reducing the value of this equation, we can obtain agent trajectories that satisfy the delivery time constraints of the order.

Finally, $c_{i,3}$ is defined as follows:

$$\begin{array}{c}\hfill c_{i,3}=w_3{\displaystyle \sum _{t={\tau }_i}^{{\tau }_i+T_f}\sum _{j\le \left|\mathcal{S}\right|}}{\alpha }_j\left(t\right)({\sigma }_j\left(t\right)-d_j\left(t\right)),\end{array}$$

(24)

where ${\sigma }_j\left(t\right)$ is a decision variable that represents the number of agents in store $s_j$ at time t. The constant $d_j\left(t\right)$ is given in advance (see Section 2.3). We assume that the demand forecast is a prediction of when, where, and how many orders will be placed. The ${\alpha }_j\left(t\right)$ is a binary variable that determines whether the number of agents in store $s_j$ at time t satisfies the required number of agents in store $s_j$ obtained from the demand forecast. The ${\alpha }_j\left(t\right)$ is defined as follows:

$$\begin{array}{c}\hfill {\alpha }_j\left(t\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{d}_j\left(t\right)\ge {\sigma }_j\left(t\right),\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.\end{array}$$

(25)

where ${\alpha }_j\left(t\right)=1$ implies that the number of agents in store $s_j$ is less than the number of agents required for store $s_j$ obtained from the demand forecast. The scalar $w_3$ is a constant representing the weights and is determined through trial and error. In other words, Equation (24) is the sum of the number of agents that are less than the number of agents required for each store. By reducing the value of this equation, the agent trajectory can be obtained with more consideration given to demand forecasting.

The pickup and delivery problem at the update time ${\tau }_i$ is given as follows.

Problem 1.

$$\begin{array}{ccc}\hfill \mathrm{minimize}& & \left(19\right)\hfill \\ \hfill \mathrm{subject}\mathrm{to}& & \left(3\right)\text{–}\left(18\right),\left(22\right),\left(23\right),\left(25\right).\hfill \end{array}$$

Via a simple calculation (see, e.g., [52]), the constraint conditions (3)–(18), (22), (23), and (25) can be transformed into linear constraint conditions. In addition, the cost function is also transformed into a linear cost function. Hence, this problem is reduced to a mixed integer linear programming (MILP) problem.

3. Online Optimization

For ${\mathcal{O}}^i$, $\mathcal{P}$, and $\mathcal{D}$, the online optimization of the pickup and delivery problem is performed as follows.

Procedure for online optimization of the pickup and delivery problem:

• Step 1: If $i=0$, set $\mathcal{P},\mathcal{D}=\varnothing$, and $t=0$.
• Step 2: For each $o_k\in \mathcal{P}\cup \mathcal{D}$, calculate the values of ${\displaystyle \sum _l\sum _{t={\tau }_{i-1}}^{{\tau }_i}{z}_{l,k}^{\mathrm{pick}}\left(t\right)}$ and ${\displaystyle \sum _l\sum _{t={\tau }_{i-1}}^{{\tau }_i}{z}_{l,k}^{\mathrm{deli}}\left(t\right)}$.
• Step 3: For $o_k\in \mathcal{P}\cup \mathcal{D}$, perform the following procedure.

  $$\begin{array}{c}\hfill \left\{\begin{array}{cc}\mathcal{P}:=\mathcal{P}\cup \left\{o_k\right\}\hfill & {\displaystyle \mathrm{if}\sum _l\sum _{t={\tau }_{i-1}}^{{\tau }_i}{z}_{l,k}^{\mathrm{pick}}\left(t\right)=0,}\hfill \\ \mathcal{P}:=\mathcal{P}\setminus \left\{o_k\right\},\mathcal{D}=\mathcal{D}\cup \left\{o_k\right\}\hfill & {\displaystyle \mathrm{if}[\sum _l\sum _{t={\tau }_{i-1}}^{{\tau }_i}{z}_{l,k}^{\mathrm{pick}}\left(t\right)=1]\wedge [\sum _l\sum _{t={\tau }_{i-1}}^{{\tau }_i}{z}_{l,k}^{\mathrm{deli}}\left(t\right)=0],}\hfill \\ \mathcal{P}:=\mathcal{P}\setminus \left\{o_k\right\},\mathcal{D}=\mathcal{D}\setminus \left\{o_k\right\}\hfill & {\displaystyle \mathrm{if}\sum _l\sum _{t={\tau }_{i-1}}^{{\tau }_i}{z}_{l,k}^{\mathrm{deli}}\left(t\right)=1.}\hfill \end{array}\right.\end{array}$$
• Step 4: Update $\mathcal{P}:=\mathcal{P}\cup {\mathcal{O}}^i$.
• Step 5: For each order, if $o_k\in \mathcal{D}$, then update $y_{l,k},l\in \mathcal{N}$ to ${\left(y_{l,k}\right)}_{i-1}$ obtained via optimization at $t={\tau }_{i-1}$, that is,

  $$\begin{array}{c}\hfill y_{l,k}:={\left(y_{l,k}\right)}_{i-1}.\end{array}$$

  (26)
• Step 6: For $o_k\in \mathcal{P}$, impose Equations (10) and (11) as a constraint. For $o_k\in \mathcal{D}$, impose (12) as a constraint.
• Step 7: Solve the pickup and delivery problem at time $t={\tau }_i$, adding the constraints in Equation (26).
• Step 8: Based on the calculation results, move the agent. Update $t:=t+1$.
• Step 9: If there are no orders at time t, return to Step 8. If there is an order, return to Step 2.

The outline of this procedure is shown in Figure 2.

In this procedure, the agents’ predicted paths may be changed by solving the optimization problem that is generated by receiving a new order. Because the agents’ predicted paths are changed online (real time), we call this procedure an online optimization.

4. Numerical Example

In this section, we present a numerical example where the optimization problem is solvable even when the demand forecast deviates significantly from the actual order by considering the delivery constraints in the cost function. Here, we compare the proposed method with the existing method [49].

4.1. Preliminaries

Suppose that the target area is given by Figure 3. The number of vertices is set to $\left|S\right|=10, \left|C\right|=10, \mathrm{and} \left|R\right|=16$, and the graph is randomly arranged. The number of agents is set to 10, and the initial position is assumed to be one agent in each store. The maximum fuel of an agent is set to $q_0=20$. All $w_1$, $w_2$, and $w_3$ in (20), (21), and (24) are set to 1. In addition, we set $T_f=20$. We used a computer with an AMD Ryzen 5 5600X 6-Core Processor 3.70 GHz CPU and 32 GB of memory. The optimization problem at each update time was written using PuLP, which is a modeler written in Python. The solver used to compute the optimization problem was IBM ILOG CPLEX 20.1.0, which can be called from Python. The simulation result is outputted as the computation result of the Python program.

The order contents are given in

Table 1. Order details.

Order	${\mathit{s}}_{\mathit{k}}$	${\mathit{c}}_{\mathit{k}}$	${\mathit{T}}_{{\mathit{o}}_{\mathit{k}}}^{\mathit{s}}$	${\mathit{T}}_{{\mathit{o}}_{\mathit{k}}}^{\mathit{c}}$
$o_1$	$s_3$	$c_2$	0	7
$o_2$	$s_2$	$c_3$	5	12
$o_3$	$s_1$	$c_4$	5	12
$o_4$	$s_4$	$c_5$	8	13
$o_5$	$s_4$	$c_6$	8	13
$o_6$	$s_4$	$c_8$	8	13
$o_7$	$s_4$	$c_9$	8	13

. In the case of this order, the optimization problem is solved at times 0, 5, and 8. However, in the existing method, ${\tau }_i+T_f$ in the Equations (11) and (12) representing the delivery time constraints is $T_{o_k}^c$. In other words, the existing method constrains the order delivery to be completed by the order delivery date, not within the time interval. In [49], the cost function $J_i$ is defined as follows, and the delivery time constraint is not considered in the cost function.

$$\begin{array}{c}\hfill J_i={\displaystyle \sum _{t={\tau }_i}^{{\tau }_i+T_f}\sum _l\sum _{j>\left|\mathcal{S}\right|}}{p}_{l,j}\left(t\right)\end{array}$$

(27)

The existing methods for demand forecasting consider the following inequalities.

$$\begin{array}{c}\hfill {\displaystyle \sum _l}{p}_{l,j}\left(t\right)\ge d_j\left(t\right),(j=1,2,\dots ,N)\end{array}$$

(28)

The $d_j\left(t\right)$ used in Equations (24) and (28) are set to $d_5\left(8\right)=4$ and 0 in all other cases. This means that at time 8, store $s_5$ needs four agents, and there are no constraints for the other times or stores. However, the actual orders are not for store $s_5$ but for store $s_4$, which means that the predictions of the stores that will receive the orders in the demand forecast are all wrong.

4.2. Calculation Results

Table 2. Results with the existing method.

	Results at Time 0			Results at Time 5			Results at Time 8
Order	Agent	Pickup Time	Delivery Time	Agent	Pickup Time	Delivery Time	Agent	Pickup Time	Delivery Time
$o_1$	3	0	5	3	0	5	–	–	–
$o_2$	–	–	–	2	5	9	–	–	–
$o_3$	–	–	–	1	5	10	–	–	–
$o_4$	–	–	–	–	–	–	–	–	–
$o_5$	–	–	–	–	–	–	–	–	–
$o_6$	–	–	–	–	–	–	–	–	–
$o_7$	–	–	–	–	–	–	–	–	–

shows the results obtained for the case solved by the existing method.

Table 3. Values of the cost function at each time in the existing method.

Time	Cost Function ${\mathit{J}}_{\mathit{i}}$
0	5
5	12
8	Infeasible

shows the cost function values at each time for the existing method.

Table 4. Results with the proposed method.

	Results at Time 0			Results at Time 5			Results at Time 8
Order	Agent	Pickup Time	Delivery Time	Agent	Pickup Time	Delivery Time	Agent	Pickup Time	Delivery Time
$o_1$	3	0	3	3	0	3	3	0	3
$o_2$	–	–	–	2	5	12	2	5	14
$o_3$	–	–	–	1	5	9	1	5	9
$o_4$	–	–	–	–	–	–	3	10	13
$o_5$	–	–	–	–	–	–	8	9	13
$o_6$	–	–	–	–	–	–	2	10	13
$o_7$	–	–	–	–	–	–	6	9	14

shows the results obtained for the case solved by the proposed method.

Table 5. Values of the cost function at each time in the proposed method.

Time	${\mathit{c}}_{\mathit{i},1}$	${\mathit{c}}_{\mathit{i},2}$	${\mathit{c}}_{\mathit{i},3}$	Cost Function ${\mathit{J}}_{\mathit{i}}$
0	5	0	0	5
5	10	0	0	10
8	17	3	0	20

shows the value of the cost function at each time for the proposed method. Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 shows the obtained trajectories of the agents. The trajectories of the agents obtained at times 0 and 5 for the existing methods are shown in Figure 4 and Figure 6. The trajectories of the agents obtained at times 0, 5, and 7 for the proposed method are shown in Figure 5, Figure 7, and Figure 8. The solid line in Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 shows that the agent has luggage, while the dotted line shows that the agent does not have luggage. The thickness of the arrows in the figure indicates the number of pieces of luggage the agent is carrying.

From Table 2 and Table 4, only the proposed method can solve the problem at time 8 for the orders given in Figure 1. This is because in the case of the existing method, Equations (11) and (12) are not satisfied for orders $o_2$ and $o_7$. On the other hand, in the case of the proposed method, we add to the cost function the amount of time past the delivery time constraint $T_{o_k}^c$ for the orders $o_2$ and $o_7$ from Table 4 and Table 5, respectively. Thus, the result at time 8 is solvable. Therefore, even when the demand forecast is significantly off and the delivery time constraints for all the orders are not met, an optimal solution that satisfies the conditions as much as possible can be obtained.

5. Conclusions

In this paper, we proposed an online pickup and delivery problem in which demand forecasting and fuel constraints were considered. First, we formulated the pickup and delivery problem with fuel constraints. Next, we introduced the cost function considering demand forecasting. After that, we explained the procedure of online optimization. Finally, through a numerical example, we presented the effectiveness of the proposed method.

In future work, we will consider practical applications. The problem formulation in this paper is relatively simpler than real food delivery services. It is important to extend the proposed method to the practical setting. We will consider many real scenarios and validate the proposed method. For large-scale scenarios, the computation time to solve an MILP problem becomes large (in this paper, as the first step, an MILP problem is solved by a conventional solver). It is also important to develop a distributed solution method, such that the problem is decomposed into small subproblems, and a parallel algorithm based on swarm optimization.