The BAP is a rectangle packing problem where a set of rectangles, denoted as O , are optimally placed in a larger rectangle, denoted as O , as shown in Figure 3. The rectangle packing problem is an NP-hard problem and can be used to describe many real-life problems (de Bruin, 2013; Murata et al., 1995). In some of these problems, the dimensions of O are held constant, such as in the problem of packing modules on a chip, where the widths and heights of the rectangles correspond to the physical widths and heights of the modules (Murata et al., 1995). Other problems, such as the one presented in this work, allow either the horizontal or vertical edge of each rectangle in O to vary. For example, suppose the vessel lengths are predefined (i.e., vertical edges are static), but the service time is allowed to vary (i.e., horizontal edges are dynamic) (Buhrkal et al., 2011).

The BAP solves the problem of optimally assigning incoming vessels to berthing positions in order to be serviced, as shown in Figure 1. To relate to the rectangle packing problem, the width and height of O represent the time horizon T seconds and the berth length L meters, respectively. Similarly, the widths and heights of each element in O represent the time spent to service vessel i and the space taken by docking vessel i , respectively. In the BAP, the vessel characteristics (length of the vessel, arrival time, handling time, and desired departure time) are assumed to be known for all vessels. A representation of a BAP solution is shown in Figure 4. The x- and y-axes represent the time horizon and berthing space, respectively. The gray squares, labeled A, B, C, and D, represent scheduled vessels. The width of the boxes represents the time spent being serviced, and the height of the boxes represents the amount of space the vessel requires on the berth. The vertical line adjacent to the “arrival time” represents the time at which the vessel arrives and is available to be berthed. The “berthing time” is the time the vessel is berthed and begins to be serviced. The “completion time” represents the time at which the berthing space becomes available again.

The BAP forms the basis of the PAP; however, there are some differences in the way the variables are interpreted. The relevant variables are given in Table 1. The starting service time, u i s, is viewed as the initial charge time, and the service time is the total elapsed time spent on the charger. Similarly, for the spatial term, v i ∈ [ 0 , L ] , the berth location is interpreted as the initial position on the charger. There are also a few clarifying concepts about how the system is modeled. The PAP models the set of chargers as one continuous line; that is, the natural behavior of the PAP model is to allow vehicles to be queued anywhere along [ 0 , L ] . Similarly, the charge times are continuous and can be placed anywhere on the time horizon, [ 0 , T ] , as long as the allocated times do not interfere with other scheduled charge times.

The PAP formulation’s parameters can be divided into two categories: input parameters and decision variables. The following parameters are assumed to be known inputs for the MILP. L defines the length of the charger in meters. As stated previously, it is modeled as a continuous bar, meaning that a vehicle can be placed anywhere in the range [ 0 , L ] . It is assumed that the time horizon, T s, is known so that vehicles may be placed temporarily in the range [ 0 , T ] . The total number of visits to the station over the time horizon is represented by n V . The arrival time for each visit is represented by a i s, and the required charge time is represented by s i s. The width of vehicle i is represented by l i m.

The decision variables provide the means by which the solver may optimize the problem. The initial and final charge times for vehicle i are u i and d i s, respectively. The starting position on the charger is denoted as v i ∈ [ 0 , L ] m. The temporal ordering of vehicles i and j is determined by σ i j ∈ { 0,1 } , where σ i j = 1   ⇒   i arrives before j for all 1 ≤ i , j ≤ n V . Similarly, ψ i j ∈ { 0,1 } determines the relative position of vehicles i and j on the charger: ψ i j = 1   ⇒   v i < v j for all 1 ≤ i , j ≤ n V .

To determine the values for each of these decision variables, a MILP was formulated by Qarebagh et al. (2019). The formulation is shown in its entirety for completeness. The problem to be solved is min   ∑ i = 1 n v d i − a i , (1)

subject to u j − u i − s i − σ i j − 1 T ≥ 0 , (2a) v j − v i − l i − ψ i j − 1 L ≥ 0 , (2b) σ i j + σ j i + ψ i j + ψ j i ≥ 1 , (2c) σ i j + σ j i ≤ 1 , (2d) ψ i j + ψ j i ≤ 1 , (2e) s i + u i = d i , (2f) a i ≤ u i ≤ T − s i , (2g) σ i j ∈ 0,1 ,   ψ i j ∈ 0,1 , (2h) v i ∈ 0 , L , (2i) i , j = 1 … n V ;   i ≠ j . (2j) The objective function, Equation 1, minimizes the idle and service time by summing over the differences between the departure time, d i , and arrival time, a i , for all visits. In other words, the objective function is to search for the schedule that removes each vehicle from the service queue as quickly as possible.

Equations 2a– e are used to ensure that individual rectangles do not overlap. In terms of the PAP, this implies that there are no conflicts in the schedule spatially or temporally. Equation 2a establishes temporal ordering when active ( σ i j = 1 ) in the manner described previously by utilizing big-M notation. Similarly, Equation 2b establishes spatial ordering when active ( ψ i j = 1 ) . Equations 2c, e enforce spatial and temporal ordering between each queue/vehicle pair. Equation 2d and Equation 2e enforce the validity of the assignments. For example, if Equation 2d resulted in a value of 2, that would imply that both vehicles i and j are scheduled before and after each other temporally, which is impossible. In the case of Equation 2e being equal to 2, it would mean that vehicles i and j are scheduled both before and after one another on the charging strip, which is again impossible.

The last constraints force relationships between arrival time, initial charge time, and departure time. Equation 2f states that the initial charge time, u i , plus the total charge time, s i , must equal the departure time, d i . Equation 2g enforces the arrival time, a i , to be less than or equal to the service start time, u i , which, in turn, must be less than or equal to the latest time the vehicle may begin charging and stay within the time horizon. Equation 2h simply states that σ i j and ψ i j are binary terms. Equations 2i, 2j ensures that the assigned value of v i is within the range [ 0 , L ] .

The queuing constraints ensure that the buses entering the charging queues are assigned feasibly. There are three sets to differentiate between different entities. B = { 1 , … , n B } is the set of bus indexes, with index b used to denote an individual bus, Q = { 1 , … , n Q } is the set of queues, with index q used to denote an individual queue, and V = { 1 , … , n V } is a set of visits to the station, with i and j used to refer to individual visits. The mapping Γ : V → B is used to map a visit index, i , to a bus index, b . The notation Γ i is used as a shorthand to refer to bus index b for visit i .

The actual physical dimensions of the BEB are ignored, and it is assumed that each BEB will be assigned to charge at a particular charger. Because of this assumption, the PAP spatial variable, l i , may be removed, and v i is made to be an integer corresponding to which queue visit i will be used, v i ∈ Q . That is, the queue position is now discretized over n Q queues, where a BEB occupies a single charge queue. Thus, when ψ i j = 1 , vehicle j is placed in a charging queue with a larger index than that of vehicle i , i.e., v j > v i . The charger length L is likewise replaced with n Q . Note that n Q = n B + n C , where n B is the number of buses and n C is the number of chargers. The rationale for adding additional idle queues is to allow BEBs to be “set aside” if no additional charge is required. Adding one idle queue for each BEB ensures that the constraints will be satisfied if multiple buses sharing overlapping times at the station are placed in idle queues. This method will be applied when defining the parameters in Section 4. The modified queuing constraints can be written as shown in Equation 3: v i − v j − ψ i j − 1 n Q ≥ 1 (3a) d i ≤ τ i , (3b) s i ≥ 0 , (3c) v i ∈ Q . (3d)

The constraint in Equation 3a is nearly identical to Equation 2b, but rather than viewing the charger as a continuous strip of length L , it is discretized into n Q queues, each with a width of unit length 1. A BEB is also assigned a unit length of 1, which is reflected in Equation 3a by ⋅ ≥ 1 . Equation 3b ensures that the time the BEB is detached from the charger, d i , is before its departure time, τ i s. Note that the introduction of the new variable τ i exists to allow the final charge time to be independent, similar to the manner in which to the initial charge time is independent of the arrival time, a i ≤ u i ≤ d i ≤ τ i . Equation 3c ensures that the change time is non-negative. Equation 3d defines v i to be an element from the set of queues.

This section introduces the battery dynamic constraints. Two constraints are enforced on the SOC for each BEB: the SOC must always remain above a specified percentage to guarantee sufficient charge to execute their respective routes, and each bus must end the day with a SOC above a specified threshold, in preparation for the next day.

The SOC upon arrival for visit i is denoted as η i kWh. Because the SOC for visit i is dependent on its previous visits, the mapping ϒ : V → V ⋃ { ∅ } is used to determine the next visit that corresponds to the same bus, with ϒ i being the shorthand notation. Thus, Γ j and Γ ϒ i , for ϒ i = j , would both map to the same bus index as long as ϒ i is not the null element, ∅ . The null element is reserved for BEBs that have no future visits.

To define the time spent on the charger, s i , as well as the initial, final, and intermediate bus charges for each visit i , the sets for the initial and final visits must be defined. Let the mapping of the first visit by each bus be denoted as Γ 0 : B → V . The resulting value of the mapping Γ 0 represents the index for the first visit of bus b . Similarly, let Γ f : B → V map the indexes for the final visits for each bus b ∈ B . Let the shorthand for each mapping be denoted as Γ b 0 and Γ b f . The initial and final bus charge percentages, α and β , can then be represented by the constraint equations η Γ b 0 = α b κ b and η Γ b f = β b κ b , respectively. The intermediate charges must be determined during runtime.

It is assumed that the charge received is proportional to the time spent charging. The rate for charger q is denoted as r q kW. Note that a value of r q = 0 corresponds to a queue where no charging occurs. A bus in such a queue simply waits at the station for the departure time. The queue indexes are ordered such that the first n B queues have r q = 0 to allow an arbitrary number of buses to sit idle at any given moment in time. The next n C queues are reserved for the slow and fast chargers. The amount of discharge between visits i and ϒ i , the next visit of the same bus, is denoted as Δ i kWh. If visit i occurred at charger q , the SOC of the BEB’s next arrival, ϒ i , would be η ϒ i = η i + s i r q − Δ i .

The binary decision variable w i q ∈ { 0,1 } is introduced to indicate the active charger for visit i in vector form. The form of the SOC for the next visit, ϒ i , can be written using the following constraints: η ϒ i = η i + ∑ q = 1 n Q s i w i q r q − Δ i , (4a) ∑ q = 1 n Q w i q = 1 , (4b) w i q ∈ 0,1 . (4c)

Where Equation 4a represents the SOC for the next visit of BEB b , Equation 4b ensures that the BEB for visit i is assigned to a single queue, and Equation 4c specifies that w i q is a binary value.

The choice of queue for visit i becomes a slack variable and is defined in terms of w i q as v i = ∑ q = 1 n Q q w i q . (5)

Maximum and minimum values for the charges are included to ensure that the battery is not overcharged and guarantee sufficient charge for subsequent visits. The upper and lower battery charge bounds for bus b are κ b and ν b κ b , respectively, where κ b is the battery capacity and ν b is a percent value. The upper and lower bounds for the current SOC are written as follows: η i + ∑ q = 1 n Q s i w i q r q ≤ κ Γ i , (6a) η i ≥ ν Γ i κ Γ i . (6b)

Equation 6a ensures that the BEB SOC does not exceed the battery capacity, and Equation 6b ensures that the initial SOC for each visit is above the threshold of ν Γ i κ Γ i . Note that the term s i w i q is a bilinear term. A standard method for linearizing a bilinear term that contains an integer variable is by introducing a slack variable with an either/or constraint (Chen et al., 2010; Rodriguez and Vecchietti, 2013). Allowing the slack variable g i q s to be equal to s i w i q , g i q can be defined as g i q = s i w i q = 1 0 w i q = 0 . (7)

Equation 7 can be expressed as a mixed-integer constraint using big-M notation with the following four constraints: s i − 1 − w i q M ≤ g i q , (8a) s i ≥ g i q , (8b) M w i q ≥ g i q , (8c) 0 ≤ g i q , (8d) where M is a large unitless value. If w i q = 1 , then Equations 8a, b become s i ≤ g i q and s i ≥ g i q , forcing s i = g i q , with Equation 8c being inactive. If w i q = 0 , Equation 8a is inactive, and Equations 8c, d force g i q = 0 .

The goal of the MILP is to utilize chargers as little as possible to reduce energy costs, with fast charging being penalized more to avoid its adverse effects on battery health and to account for the higher usage costs. Thus, assignment cost m q and usage cost ϵ q are associated with each charger, q . These unitless weights can be adjusted based on the charger type or time of day when the visit occurs. The assignment term takes the form w i q m q , and the usage term takes the form g i q ϵ q . The resulting BEB charging problem is defined in Equation 9: min ∑ i = 1 N ∑ q = 1 n Q w i q m q + g i q ϵ q , (9)

subject to the constraints u j − u i − s i − σ i j − 1 T ≥ 0 , (10a) v j − v i − ψ i j − 1 n Q ≥ 1 , (10b) σ i j + σ j i + ψ i j + ψ j i ≥ 1 , (10c) σ i j + σ j i ≤ 1 , (10d) ψ i j + ψ j i ≤ 1 , (10e) s i + u i = d i , (10f) η Γ b 0 = α Γ i κ Γ i , (10g) a i ≤ u i ≤ T − s i , (10h) d i ≤ τ i , (10i) η i + ∑ q = 1 n Q g i q r q − Δ i = η γ i , (10j) η i + ∑ q = 1 n Q g i q r q − Δ i ≥ ν Γ i κ Γ i , (10k) η i + ∑ q = 1 n Q g i q r q ≤ κ Γ i , (10l) η Γ b f ≥ β Γ f κ Γ f , (10m) s i − 1 − w i q M ≤ g i q , (10n) s i ≥ g i q , (10o) M w i q ≥ g i q , (10p) 0 ≤ g i q , (10q) v i = ∑ q = 1 n Q q w i q , (10r) ∑ q = 1 n Q w i q = 1 , (10s) w i q , σ i j , ψ i j ∈ 0,1 , (10t) v i , q i ∈ Q , (10u) i ∈ V . (10v)

Equations 10a– i are reiterations of the queuing constraints from Equation 2; Equation 3. Equations 10g–m provide the battery charge constraints. Equations 10n–q define the charge duration of every visit/queue pairing. Equations 10r, 10s are reiterations of Equation 4b and Equation 5, respectively. Equations 10t–v define the sets of valid values for each variable.

To display the capabilities of the model, an example scenario is presented. The scenario was run over a time horizon of T = 24 h, utilizing n B = 35 buses with n V = 338 visits divided between the n B buses. As stated before, the route times are sampled from a set of routes provided by the Utah Transit Authority (UTA). Each bus has a battery capacity of κ b = 388 kWh, which is required to stay above an SOC of ν b = 20% (77.6 kWh) ∀ b ∈ B to ensure that each BEB can complete its next route while maintaining battery health. Each bus is assumed to begin the working day with an SOC of a l p h a b = 90% ∀ i ∈ V (349.2 kWh). Additionally, each bus is required to end the day with a minimum SOC of β b = 70% (271.6 kWh) ∀ b ∈ B . This assumes that overnight charging can account for the 20% SOC deficit. Furthermore, the value selected for β b   ∀ b ∈ B is primarily for illustrative purposes and has no particular significance other than to demonstrate the user’s ability to specify the desired end-of-day SOCs. Each bus is assumed to discharge at a rate of ζ b = 30 kW. Note that many factors play a role in the rate of discharge; however, for the sake of simplicity, since the discharge calculation is beyond the scope of this work, an average rate is used. A total of n C = 30 chargers are utilized, where 15 of the chargers are slow charging (30 kW) and 15 are fast charging (911 kW). The technique to minimize the total charger count is used.

To encourage the MILP PAP to utilize the fewest number of chargers, the value of m q in the objective function, Equation 9, is set as follows: ∀ q ∈ { 1,2 , … , n B } ; m q = 0 and ∀ q ∈ { n B + 1 , n B + 2 , … , n B + n C } ; m q = 1000 q . The charge duration scalar, ϵ q , is defined as ϵ q = r q to create a consumption cost term, g i q ϵ q kWh. By consumption cost, it is meant that the total energy consumed by the charge schedule will be accounted for in the objective function. This method encourages the model to minimize active charger times, particularly for the fast chargers.

Another heuristic-based optimization strategy, referred to as the Qin-modified method, is also used to compare the results of the MILP PAP. The Qin-modified strategy is based on the threshold strategy proposed by Qin et al. (2016). The strategy has been modified slightly to accommodate the case of multiple charger types without an exhaustive search for the best charger type. The heuristic is based on a set of rules that revolve around the initial SOC of the bus visit i . There are three different thresholds, low (85%), medium (90%), and high (95%). Buses below the low threshold ( SOC ≤ 85 %) are prioritized to fast chargers and are then allowed to utilize slow chargers if no fast chargers are available. Buses between the low and medium threshold (85% < SOC ≤ 90%) prioritize slow chargers first and utilize fast chargers only if no slow chargers are available. Buses above the medium threshold and below the high threshold (90% < SOC ≤ 95%) will only be assigned to slow chargers. Buses above the high threshold ( SOC > 95%) will not be placed in a charging queue. Once a bus has been assigned to a charger, it remains on the charger for the duration of the time it is at the station or reaches an SOC of 95% charge, whichever comes first.

The total number of constraints resulted in 7,511 continuous and 328,282 integer/binary constraints. The optimization was performed using the Gurobi MIP solver (Gurobi Optimization, LLC, 2021) on a machine equipped with an AMD Ryzen 9 5900 × 12 Processor (24 core) running at 4.95 GHz. The solver was allowed to run for 7.200 s.

The schedule generated by the Qin-modified strategy and the MILP PAP is shown in Figures 5A, B, respectively. The x-axis represents the time in hours. The y-axis represents the assigned charging queue. Rows between 0 and 14 are active times for slow chargers, and rows in the range of 15 and 29 are active times for fast chargers. The unique color/symbol-styled vertices represent the starting charge time for a bus b , with the line to the vertical tick signifying the region of time the charger is active. The lines connecting the points represent the charge sequence for each BEB.

The first observation is in the choice of preferred chargers between the Qin-modified and MILP scheduler. Figures 6A, B show that the Qin-modified schedule uses at most four fast chargers and three slow chargers at the same time, whereas the MILP schedule uses at most one fast charger and six slow chargers at the same time. Both the Qin-modified and MILP schedule used the fast chargers in short bursts (˜0.2–0.5 h). The main difference lies in the utilization strategy of the slow chargers. The Qin-modified method, for the most part, opted for shorter bursts for the slow chargers (˜0.3–0.7 h), most heavily placed on the first slow charger. The MILP utilized the slow chargers in short bursts; however, the solver recognized moments where a BEB being placed in a slow-charging queue for a longer duration was more cost-effective (with respect to the objective function) than the BEB placed in a fast-charging queue. Although one of the objectives of the MILP was to minimize the amount of chargers used, the Qin-modified schedule ended up using fewer chargers than the MILP. Note that the MILP schedule packed the first queue for the fast and slow chargers more effectively than the Qin-modified schedule. Although both schedules generated are valid, no comparison of the quality of the schedule can be made directly from Figures 5A, B.

Figures 7A, B depict the charge for every bus over the time horizon for the Qin-modified and MILP schedules, respectively. Every vehicle begins with a SOC of α b = 90%, finishes with a SOC of β b = 70% in the MILP PAP schedule, and never goes below 20% in the intermediate arrivals, as stated in Equation 10. There is no guarantee for this in the Qin-modified strategy, which can be observed by some intermediate charges reaching a SOC of 0% and the distribution of final charges, the minimum being 0% and the maximum being 94.845%. It is of note that the SOC plummeting to 0% is due to the reactive nature of the Qin-modified method. Specifically, the schedule utilized in this example contains some BEBs that begin with routes that deplete the SOC by small amounts but do not fall below the threshold. These routes are followed by small charge windows, which are then preceded by large routes depleting the SOC by a significant amount, resulting in the SOC reaching 0%. The only guarantee that the Qin-modified method supplies is its predictability for the intermediate visits due to its heuristic nature (i.e., if the BEB charge is within the low threshold, a fast charger will be prioritized), whereas MILP places a bus in the queue that “makes sense” in the context of the larger picture. The MILP PAP does not have an obvious decision-making process because its weighted objective function is affected by the accumulation of previous decisions.

Another important measure for the chargers is to compare the amount of power and energy consumed. Figure 8 depicts the power consumption throughout the time horizon. It can be seen that the Qin-modified power consumption is steadily less or the same as that of the MILP schedule. This can be accounted for by the MILP’s constraints to keep the bus SOC above 20% and to reach a final SOC of 70% at the end of the working day. Similarly, the accumulated energy consumed is shown in Figure 9. The MILP schedule is more efficient up until about the 11th hour. Again, this can be accounted for by the fact that the MILP accommodates the extra constraints. Due to these constraints, the MILP PAP consumes approximately 0.1 ⋅ 1 0 4 kWh more than the Qin-modified method. The overlap of the MILP PAP can be accounted for by referencing Figures 6A, B. Between the 5th and 10th hours, the MILP schedule mainly uses slow chargers, increasing the rate at which power is being consumed. Afterward, the MILP schedule continues to use at least the same number of chargers as the Qin-modified schedule. However, due to the added constraints, the MILP schedule must utilize more resources to remain within the specified bounds.