Battery Dimensioning for Hybrid Vehicles in a Routing Application Generalised Duality and Logic-Based Benders Decomposition Master’s thesis in Engineering Mathematics and Computational Science JONAS KINDSTRAND & LINUS NORDGREN Department of Mathematical Sciences CHALMERS UNIVERSITY OF TECHNOLOGY UNIVERSITY OF GOTHENBURG Gothenburg, Sweden 2018 Master’s Thesis 2018 Battery Dimensioning for Hybrid Vehicles in a Routing Application Generalised Duality and Logic-Based Benders Decomposition Authors: Supervisor: Jonas Kindstrand Ann-Brith Strömberg Linus Nordgren Department of Mathematical Sciences Chalmers University of Technology University of Gothenburg Gothenburg, Sweden 2018 Battery Dimensioning for Hybrid Vehicles in a Routing Application An Application of Logic-Based Benders Decomposition Jonas Kindstrand & Linus Nordgren © 2018 Jonas Kindstrand & Linus Nordgren Supervisor & examiner: Ann-Brith Strömberg Master’s thesis 2018 Department of Mathematical Sciences Chalmers University of Technology and University of Gothenburg SE-412 96 Gothenburg Sweden Telephone +46 (0)31-772 1000 Cover: An illustration of a solution to a hybrid vehicle routing problem with two vehicles. The depot node is orange and the dashed arrow indicates the freedom to not take the detour through the blue recharge station node. Typeset in LATEX Gothenburg, Sweden 2018 ii Battery Dimensioning for Hybrid Vehicles in a Routing Application An Application of Logic-Based Benders Decomposition Jonas Kindstrand & Linus Nordgren Department of Mathematical Sciences Chalmers University of Technology and University of Gothenburg Abstract The Vehicle Routing Problem (VRP), which is defined as to find optimal routes for a fleet of delivery vehicles to various customers, constitute an important class of combinatorial optimisation problems of both practical and theoretical interest. Among the various flavours of VRP, this report specifically focuses on a case with hybrid vehicles with two fuel types, with the goal of finding the optimal battery sizes which minimises the total cost. We present an exact solution method using a generalised Benders decomposition method, known as logic-based Benders decomposition. In this method, the subproblems are generalised to mixed integer linear optimisation problems. The master problem is a simple routing problem, while the subproblems concern resource constraints and battery types. The mixed integer master problem is solved by branch-and-bound, and lower bounds are generated from the solution tree. Only small instances of up to 14 customers are solved to optimality, and the performance of our algorithm is compared with more direct solution methods. As it is, the method is slower than solving the full problem directly, and further work is needed to make it competitive. Keywords: Vehicle routing problem (VRP), hybrid vehicles, battery capacity, logic-based Benders decomposition (LBBD), branch-and-bound iii Contents Abstract . iii Contents . v 1 Introduction . 1 1.1 Problem Formulation & Project Aim . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Outline of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Literature study . 5 2.1 Overview of Vehicle Routing Problems . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Benders Decomposition and Applications to VRP . . . . . . . . . . . . . . . . . . . 6 3 Theory . 9 3.1 Linear and Integer Linear Programming Preliminaries . . . . . . . . . . . . . . . . 9 3.1.1 Linear Programming, Duality and Sensitivity Analysis . . . . . . . . . . . . 9 3.1.2 Mixed Integer Linear Programming and Branch-and-Bound . . . . . . . . . 10 3.1.3 A Branch-and-Bound Example . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.1.4 The Benders Decomposition Algorithm . . . . . . . . . . . . . . . . . . . . 14 3.2 Inference Duality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2.1 The Inference Dual of a Linear Program . . . . . . . . . . . . . . . . . . . . 17 3.3 The Logic-Based Benders Decomposition Algorithm . . . . . . . . . . . . . . . . . 18 3.4 Logical Clauses and Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.4.1 Resolution Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.5 An LBBD Method for Mixed binary Linear Programming . . . . . . . . . . . . . . 21 3.5.1 Surrogate Inequalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.5.2 Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4 Mathematical Formulation of the Problem . 29 4.1 A Mixed Binary Linear Optimisation Model . . . . . . . . . . . . . . . . . . . . . . 30 4.2 A Reformulation of the Mixed Binary Linear Model . . . . . . . . . . . . . . . . . 31 5 Method . 33 5.1 Linear Programming Subproblems . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.1.1 The Time Window Subproblem . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.1.2 The Cargo Subproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.2 The Battery Charge Distribution Subproblem . . . . . . . . . . . . . . . . . . . . . 36 5.3 Default Lower Bounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.4 The Master Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 v CONTENTS 5.5 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 6 Numerical Tests and Results . 43 6.1 Implementation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 6.2 Test Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 6.3 Test Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 6.4 Test Results and Algorithm Performance . . . . . . . . . . . . . . . . . . . . . . . . 44 6.4.1 On the Strength of the Logic-Based Benders Cuts . . . . . . . . . . . . . . 47 7 Discussion . 51 7.1 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 7.2 Further Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 7.2.1 Minor Improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 7.2.2 Structural Improvements to the Benders Algorithm . . . . . . . . . . . . . . 53 7.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Bibliography . 55 A Problem Data for Test Instance P6 . 59 B Computation Times . 61 C Alternative Proof of Lemma 3.5.2 . 65 vi 1 Introduction Vehicle Routing Problems, abbreviated throughout this text as VRPs, constitute a well-studied class of problems in combinatorial optimisation and operational research. Generally, the problem consists of a road network, a set of customers with demands interspersed in the road network, and a vehicle fleet to service those demands. The goal is typically to minimise the total distance that the vehicles travel. Multiple variations of this problem have been investigated; usually some additional constraints are added to the problem to more accurately model real-life routing applications. Some well-studied variations add limited load capacities to vehicles [32 . ], time windows for deliveries [46 . ], synchronisation constraints [23 . ], or the need to visit some customers in a particular order [6 . ]. Mathematically, a VRP is represented as a weighted graph; customers, the depot, and charging stations are represented as nodes, and roads are represented as arcs. A route in the road network corresponds to a path in the graph, and road lengths are represented by the arc weights. The problem is then to find routes emerging from the depot node such that all customer nodes are included in at least one route, and to minimise the total length of the combined routes while doing so. As alluded to earlier, additional constraints may be put on this minimisation problem to represent the amount of goods to be picked up at each customer along with vehicle load capacities, time constrained nodes that must be visited during certain time windows, and so on. From an environmental point of view, minimising the total fuel use of the fleet is of great interest. Normally this is not meaningfully different mathematically from minimising total distance travelled, since each arc is assigned a constant cost in both cases. But if instead the problem contains two different fuel types and the total combined fuel cost is minimised, this is qualitatively a different problem compared to the single fuel case. For example, this is the case with hybrid vehicles which can use both conventional fuel and electricity. Hybrid vehicles are unable to drive long distances on electricity before recharging. Any VRP model containing hybrid vehicles needs to consider how this changes the problem. Normally fuel stations can be left out of the model since the time to refuel is negligible, but this is not the case for electric vehicles. Driving on electricity is typically cheaper, which means that it might be beneficial for a vehicle to take a longer route if this allows a stop at a charging station and less use of conventional fuel. All of these complications must be reflected in a solution method. Solution methods are an important aspect of VRP type problems, since the problems are often computationally heavy due to the large number of constraints and variables needed to model them. In practice, both exact and heuristic methods are used. Branch-and-bound and variants thereof are some of the most common exact solution methods. A search of literature revealed few studies which solve VRPs with hybrid vehicles, and especially few which employed exact methods. For more complex VRP models, branch-and-bound might not be enough, due to the fact that the problems become too large. Column generation (originally introduced in [28 . ]) is a type of decomposition method which has been found the be effective at solving these types of problems. Another one is Benders decomposition which was first proposed by Benders [5 . ] in 1 CHAPTER 1. INTRODUCTION 1962. It is a method for decomposing certain kinds of complicated optimisation problems. Since the original suggestion of the method by Benders, it has been generalised to work for a larger class of problems. One such generalisation is the logic-based Benders decomposition algorithm, abbreviated as LBBD, presented by Hooker and Dawande [22 . ] which is based on ideas from computational logic. One advantage of these kinds of decomposition methods is that by splitting the problem into several parts, solution methods can be tailored to each part of the problem. 1.1 Problem Formulation & Project Aim The problem studied in this project is a modification of the VRP with load capacity and time window constraints for hybrid vehicles. Multiples vehicles have to be routed out of one depot in such a way that all customers are visited by exactly one vehicle. Each customer has a time window in which they can be serviced, and requires a certain amount of cargo. All pairs of customers are connected to each other and to the depot by road links. For simplicity, the vehicles are modelled to have the same maximum cargo capacity. The vehicles use two different fuel types: Conventional fuel and electricity. It is assumed that the conventional fuel capacities are unlimited, since either refuelling takes a negligible amount of time or the distances are sufficiently small such that refuelling is not needed. Recharging batteries, however, takes a significant amount of time. Therefore maximum battery capacities are included in the model together with the location of recharge stations and recharge times. Moreover, it is assumed that the cost for each unit of distance is significantly higher when using conventional fuel than when using electricity, and that recharge stations are connected to every customer and to the depot. In addition, it is assumed that the vehicles always wait for a full recharge before leaving the recharge stations. Customers, the depot, and the recharge stations are modelled as nodes in a graph, with arcs representing routes between nodes. Combined with the constraints mentioned above, the problem can then be written as a minimisation problem. A so called three-index vehicle flow formulation [18 . ] is used as the model, and the different resource constraints are easily added to this model. The full mathematical model with all constraints is presented in Chapter 4 . . This model is solved using a Benders decomposition algorithm, which splits the problem into two layers. These are referred to as the master problem and the subproblem. The master problem is a simplified routing problem without many of the constraints from the full formulation. Cargo capacity, time windows, and battery constraints are included only in the subproblem, which turns out be decomposable into three corresponding parts. A consequence of this particular split is that the subproblem is not of the kind typically considered in the classical Benders decomposition algorithm, since it includes binary variables. This necessitates a different approach, and a so called logic-based Benders decomposition algorithm will be used instead, which has the following requirements: i. A proof scheme has to be devised for the solution to the subproblem. ii. A method for generating sensitivity information from the solution to the subproblem has to be established. The aim of this project is to formulate and implement a logic-based Benders decomposition algorithm, and in particular the two requirements above, for the problem outlined in this section. Limitations to this problem are discussed in the following section. 2 CHAPTER 1. INTRODUCTION 1.2 Limitations To simplify the model, it is assumed that vehicles are instantaneously recharged at recharge stations. The technical reason for this is so that the variables that model time and charge should not be directly dependent on each other. Further, it is possible to consider a problem with several depot nodes, but only one depot will be used in the model. The master problem from the logic-based Benders decomposition algorithm is computationally quite difficult to solve, but this is not the primary focus of this project. Instead, the focus is on the theoretical aspect of how to solve and generate sensitivity information from the subproblem, which is a difficult problem in its own right. The subproblem consists of three parts, and the focus is on the battery capacity part. The time window and cargo capacity parts are simplified variants of the battery part, and are thus not as interesting in the context of the focus chosen. Finally, there is a natural limit imposed by the size of problem instances a computer can feasibly solve in a reasonable time frame with the available hardware. For this project only instances with no more than 15 customers have been fully solved. 1.3 Outline of Thesis Following this introductory chapter, Chapter 2 . contains a literature review, focusing on the vehicle routing problem and Benders decomposition. Chapter 3 . is a presentation of the necessary theoretical concepts used later in the report: basic mixed-integer linear programming (MILP), classical Benders decomposition for linear programs (LP), and logic-based Benders decomposition for MILP. In Chapter 4 . , we present the mathematical formulation of our VRP as one large MILP. Chapter 5 . describes how this model is partitioned in order for a logic-based Benders decomposition scheme to be applied, and discusses the master problem and subproblems which arise. The chapter also contains a brief description of our computer implementation of the resulting algorithm. Chapter 6 . contains our results, along with a discussion about the efficacy of our method. Finally, in Chapter 7 . we discuss concrete directions for further research, both in terms of performance tuning and in broader changes to our method. The report ends with our conclusions in Section 7.3 . . 3 2 Literature study In the following section a brief overview of existing literature on different VRP problems is presented together with some common methods used for solving them. Section 2.2 . continues by focussing on the Benders decomposition algorithm, and extensions to it. 2.1 Overview of Vehicle Routing Problems The vehicle routing problem was originally introduced as ‘The Truck Dispatching Problem’ by Dantzig and Ramser [11 . ] as a generalisation of the so called Travelling Salesman Problem (TSP), which is the problem of finding the shortest path that visits all the nodes in a graph. Their VRP was distinct from a multiple TSP by the introduction of carrying capacities. Although the concept of an NP problem (we will not be needing concepts from complexity theory, but an overview can be found in [1 . ]) had at this point not yet been named or even formalised, the computational intractability of the travelling salesman problem for even relatively small problem instances was well recognised. As a generalisation, the VRP inherits all the difficulties of the TSP, and hence the method of Dantzig and Ramser was heuristic and could only find near-best solutions in some scenarios. An example of a VRP route is shown in Figure 2.1 . . An early exact method for the VRP appeared in Christofides and Eilon [7 . ], in the form of a so called branch-and-bound algorithm. Exact methods based on so called cutting planes had however been described for the TSP much earlier [12 . ]. These two approaches, branching and cutting planes, are the core of most exact methods. A combination of the two, branch-and-cut, forms the basis for several general purpose integer optimisation solvers, including IBM ILOG CPLEX Optimization Studio [27 . ], Gurobi [19 . ], and COIN-OR [47 . ]. A survey of different exact methods was done by Laporte and Nobert [34 . ], who present the state of the art up until the late eighties. Methods discussed included direct tree search, dynamic programming, and methods using integer linear programming (ILP) formulations. Most of the ILP was based on vehicle flow formulations, where binary variables were used to indicate if a vehicle travels between two nodes. Two different types were used: three-index formulations and two- indexed formulations. The three-index formulation has one binary variable for each combination of vehicle, origin, and destination, whereas the two-indexed formulation has integer variables counting how many vehicles travel from each origin to each destination. The three-indexed formulation was first introduced by Golden, Magnanti, and Nguyen [18 . ], and the two-indexed formulation was presented by Laporte and Nobert [32 . ] together with a commodity flow model. Letchford and Salazar-González [37 . ] give a survey of formulations for the capacitated VRP, including the three-index formulations. A set partitioning formulation of the capacitated vehicle routing problem together with an exact algorithm was presented in [2 . ]. More recent surveys of exact methods are Toth and Vigo [49 . ] and Laporte [33 . ]. Solomon [46 . ] presented a variant of the VRP with time window constraints, together with 5 CHAPTER 2. LITERATURE STUDY 0 1 3 56 7 10 11 0 2 8 9 13 12 Figure 2.1: Example of route for a VRP with two vehicles (red and blue) and one depot (0). Every node except for the depot is a customer that must be visited by at least one vehicle. All vehicles start and end in the depot. heuristic solution methods and data sets of test problem instances. Another variant of the VRP is the electrical vehicle routing problem in which the vehicles have a limited range before they have to recharge at a recharge station. Several different models for this problem have been proposed. Lin, Zhou, and Wolfson [38 . ] presented a version where battery consumption was affected by vehicle loads. They found that the effect of vehicle load on routing strategy cannot be ignored. A branch-and-price-and-cut algorithm was presented in [14 . ] as a method for solving the electrical VRP with time windows, where a labelling algorithm was used to solve the subproblems. An extension to the electrical VRP is to consider vehicles that can use both conventional fuel and electrical charge. A ‘Hybrid Vehicle Routing problem’ was introduced by Mancini [39 . ], in which a large neighbourhood search is used to solve a VRP where a vehicle switches to the more expensive conventional fuel once the battery is empty. 2.2 Benders Decomposition and Applications to VRP The (classical) Benders decomposition algorithm was introduced by Benders [5 . ] in 1962. In broad terms, the method partitions an optimisation problem into a master problem and one or more LP subproblems. The master problem is solved and yields trial values which are fed to the subproblems. The subproblems are solved, and using LP duality theory, constraints (commonly called cuts) are generated and fed back to the master problem, which is then re-solved, leading to new trial values for the subproblem. This back-and-forth continues until a solution to the original problem has been found and verified. More explicitly, classical Benders cuts are obtained by 6 CHAPTER 2. LITERATURE STUDY solving the dual of the subproblem, from which a bound on the optimal value of the subproblem is deduced. If the subproblem is an ILP or MILP, the LP dual is however not defined. An LP can still be obtained by relaxing the integrality constraints, and the bound from this relaxed problem is indeed valid for the original problem. Even so, if the duality gap is large, this bound will not be very strong [43 . ], and results in slow convergence if used in the Benders decomposition method. Geoffrion [17 . ] extended the Benders algorithm to non-linear convex subproblems. Another extension is presented by Hooker [20 . ], where an inference dual is defined for very general optimisation problems. Solving this dual amounts to examining the results of a primal method. The inference dual is used in [22 . ] to develop a generalised Benders decomposition method, resulting in a new method named logic-based Benders decomposition. This method requires a proof schema, and that a method of generating Benders cuts be devised for each class of optimisation problems. Hooker and Dawande [22 . ] and Hooker and Ottosson [24 . ] applied this to binary programming problems by interpreting the branch-and-bound tree as a proof of optimality. This was then used to perform sensitivity analysis of the problem. The analysis was extended to branch-and-bound methods which included cutting planes. Hooker [21 . ] studied a planning and scheduling problem where the problem naturally decomposes into an assignment portion and a scheduling portion. The assignment problem was solved as an MILP model, and the scheduling problem was solved as a constraint programming model. A logic-based Benders decomposition algorithm was used to link these problems, and it was found that this gave a speedup over other state of the art methods. Coban and Hooker [9 . ] showed that for pure scheduling problems, a logic-based Benders decomposition algorithm performed well when the problem size scaled up. Hooker [26 . ] contains several examples of inference rules for different classes of problems, together with a so called branching dual which obtains dual information by inspection of a solution tree formed by branching on the primal problem. Kloimüllner and Raidl [30 . ] applied a logic-based Benders decomposition algorithm to a ‘Balancing Bike Sharing Systems problem’. The problem was to balance bikes in a public bike sharing network to prevent too many rental stations running full or empty. They found that using a logic-based Benders decomposition algorithm was an improvement over previous exact methods. A branch-and-check method was introduced in Thorsteinsson [48 . ], which they describe as essentially a branch-and-cut method in which a Benders subproblem is used as what they call the separation subproblem. This combines the ideas from the logic-based Benders decomposition algorithm with a branching scheme, so that subproblems of almost any form can be considered. They present a branch-and-check method for solving a capacitated VRP with time windows. Lam and Hentenryck [31 . ] considers a branch-and-price-and-check method for solving a VRP with location congestion. Branch-and-price is used to solve the underlying VRP, and constraint programming to check the feasibility of the location congestion subproblem. It contains a nice overview of different hybrid methods for solving VRPs. Riazi, Seatzu, Wigström, and Lennartson [42 . ] proposed a logic-based Benders decomposition algorithm for solving the heterogeneous multi-vehicle routing problem, a extension of the multi- TSP. The nodes are assigned to different vehicles in a master problem and each subproblem consisted of solving a TSP for that vehicle. These subproblems were solved using a standard TSP. Fischetti, Ljubić, and Sinnl [15 . , 16 . ] describes the basic steps for designing a Benders decom- position approach to be embedded in a modern MILP solver. 7 3 Theory In this chapter, we first introduce some preliminary theory of linear programming (LP) and integer linear programming (ILP) in Section 3.1 . . This theory forms a necessary basis for understanding the methods presented in Sections 3.2 . , 3.3 . , and 3.5 . . Section 3.2 . introduces the concept of inference duals. This is then used in Section 3.3 . to generalise the classical Benders decomposition algorithm. Section 3.4 . introduces notation for logical clauses and the resolution algorithm. Finally, in Section 3.5 . , we fill in the details necessary to formulate a logic-based Benders decomposition algorithm for a problem with mixed-integer linear programming (MILP) subproblems. We adopt the convention that some vectors are in the dual (i.e. row) vector space, and denote this as v ∈ V >. This eliminates the need to write explicit transposes of vectors. 3.1 Linear and Integer Linear Programming Preliminaries The theory of (continuous) LP forms the basis for the usual treatment of exact methods for ILP and MILP. This section begins with a short review of this theory and describes its application to the branch-and-bound method for MILP. The section ends with a description of the Benders decomposition method for LP. 3.1.1 Linear Programming, Duality and Sensitivity Analysis Consider a general LP, min x ax, (3.1.1a) s. t. Ax ≥ c, (3.1.1b) x ≥ 0, (3.1.1c) where A ∈ Rm×n, x ∈ Rn, a ∈ (Rn)> and c ∈ Rm. The feasible set defined by the constraint (3.1.1b . ) together with the domain defined by constraint (3.1.1c . ) define a convex polyhedron. An example of a polyhedron defined in this way is shown in Figure 3.1 . . The shape of the region together with the fact that the objective function is linear guarantees that an optimal solution can always be found in a vertex, of which there are a finite number, provided that the polyhedron is not unbounded. Since a solution can always be found in one of these vertices, one solution method might be to check the objective function value for each vertex. One of the most commonly used algorithms for solving LP is the simplex method. The simplex method is a refinement of the idea of checking vertices which draws conclusions from how the objective function value changes when the values of the variables change. Consider now a vector of non-negative multipliers u ∈ (Rm)>, where each multiplier is associated with one row of the constraint matrix A. Since u ≥ 0, multiplying u with both sides 9 CHAPTER 3. THEORY (a) (b) Figure 3.1: The feasible set together with the domain define a convex polyhedron. Figure 3.1a . shows the feasible set defined by a linear inequalities and the domain R2 ≥0. Figure 3.1b . shows the convex polyhedron. of (3.1.1b . ) implies a new inequality uAx ≥ uc. Assuming that uA ≤ a yields the bound uc ≤ ax on the optimal solution of the model (3.1.1 . ). Since this bound is valid for all u ≥ 0 satisfying the inequality uA ≤ a, it is also valid for the u maximising uc. This new maximisation problem is referred to as the dual problem, and the original problem as the primal problem. We write it as max u uc, (3.1.2a) s. t. uA ≤ a, (3.1.2b) u ≥ 0. (3.1.2c) From the construction of the dual problem, it is apparent that the optimal value of the primal problem is bounded from below by the optimal value of the dual problem. This property is referred to as weak duality. In fact, an even stronger statement holds for LP. This is the content of the following theorem due to von Neumann, which is presented here without proof. The crucial property that is used to prove Theorem 3.1.1 . is some version of a hyperplane separation theorem, typically Farkas’ lemma. Several proofs can be found in [13 . ]. Theorem 3.1.1 (Strong Duality). If the primal problem (3.1.1 . ) has an optimal solution, then so does the dual problem (3.1.2 . ), and their optimal objective function values coincide. The LP dual can be used to perform sensitivity analysis on the primal LP model (3.1.1 . ). Let u∗ be the optimal solution to the dual problem (3.1.2 . ). If the vector c is perturbed to c+∆c, u∗ remains dual feasible, since the feasible set has not changed for the dual problem. This gives a dual objective function value u∗(c+∆c). Since u∗ is dual feasible, the optimal dual value is bounded from below by u∗(c+∆c). From Theorem 3.1.1 . it follows that the analogue is true for the optimal value of the primal problem. This gives us a lower bound ax ≥ u∗(c+∆c) on the optimal value of the perturbed primal problem. 3.1.2 Mixed Integer Linear Programming and Branch-and-Bound A general MILP can be stated as min cx, s. t. Ax ≥ a, x ≥ 0, x ∈ Zk × Rn−k, 10 CHAPTER 3. THEORY with c ∈ (Rn)>, A ∈ Rm×n, and a ∈ Rn. This class of optimisation problems contains LP, binary linear programs, and ILP as special cases. Moreover, since general MILP are known to be NP-hard, any NP-complete problem can in fact be reformulated as an MILP. (For general complexity theory, see, e.g., [1 . ].) As a consequence, any exact method for general MILP must be able to handle a very broad class of problems. In practice, this means that it is typically undesirable to use an entirely general method, and consequently more problem specific methods are usually developed. This section is concerned with branch-and-bound—one of the most prevalent design schemes for such methods. The branch-and-bound paradigm of algorithms is the backbone of many exact methods for MILP. At the core of the method is the idea of implicit enumeration of the feasible set (also referred to as the search space in the pure integer case). Consider for a moment a pure binary program in n variables x1, . . . , xn with objective function f . A divide-and-conquer approach to this problem might be to choose one of the variables, say x1, and to split the search space into two halves corresponding to the trial assignments x1 = 0 and x1 = 1. This process is called branching. If either of these smaller search spaces represents an infeasible problem, that branch need not be considered further. However, for each half that was feasible, choose a new variable, say x2, and again make trial assignments x2 = 0 and x2 = 1. A search tree is built by iterating this procedure. Iteration continues until all assignment options have been explicitly considered or deemed infeasible. After the process terminates, all feasible points can be found as leaf nodes of the search tree, along with the leaf nodes that were deemed infeasible. The optimal variable assignment x∗ = (x∗1, x ∗ 2, . . . , x ∗ n) is then the one with the lowest objective function value, f(x∗), among all these feasible leaf nodes. In general, a branching procedure like this can be performed by partitioning the variables and forming one new branch per set in the partition. Branching alone only amounts to a method for explicit enumeration of the search space; the objective function is evaluated at every feasible point. To improve on this, the method is augmented with the help of bounds on the objective function value. Every feasible point found by branching gives an upper bound on the objective function. By relaxing the condition x ∈ Zk ×Rn−k to x ∈ Rn in the node problems, LP are obtained. Solving this LP relaxation of a node problem gives a lower bound on the objective function for that branch. If at any point such a relaxed node problem is solved and the optimal solution is found to exceed the best known upper bound, then that branch need not be considered, and can be discarded from the problem. This way, parts of the search space are not considered explicitly and instead removed by bounding, and as such branch-and-bound is said to perform implicit enumeration of the search space. The general MILP case, where x may contain both discrete and continuous parts, does not represent a problem for branch-and-bound. Branching need only be performed on the discrete part of x, and since LP relaxations are solved, the continuous part is automatically handled. 3.1.3 A Branch-and-Bound Example In this section, consider a simple branch-and-bound application to the ILP min − 5x1 − 8x2, (3.1.4a) s. t. x1 + x2 ≤ 6, (3.1.4b) 5x1 + 9x2 ≤ 45, (3.1.4c) x1, x2 ∈ Z+. (3.1.4d) We will build the solution tree in a series of iterations, solve an LP relaxation in each one, and decide whether to branch or not in each individual node using the information from the LP relaxation. 11 CHAPTER 3. THEORY Iteration 0 (L0): Let L0 be the LP relaxation of problem (3.1.4 . ). The relaxation has an optimal solution x0 = (2.25, 3.75) with objective function value z0 = f(x0) = −41.25. Both x01 and x02 are fractional, so we pick (arbitrarily) x01 to branch on. We form one new LP L1 by adding the constraint x1 ≤ 2 to L0, and another new LP L2 by instead adding the constraint x1 ≥ 3 to L0. These two nodes constitute the new branches in the branches in the branching tree. Iteration 1 (L1 = L0 ∧ (x1 ≤ 2)): The solution to the LP L1 is x1 = (2, 3.888) with the objective function value z1 = −41.111. Only x12 is fractional, so it is the branching variable. We form the problem L3 by adding the constraint x2 ≥ 3, and the problem L4 by adding the constraint x2 ≥ 4. Iteration 2 (L3 = L1 ∧ (x2 ≤ 3)): The solution to the LP L3 is x3 = (2, 3) with the objective function value z3 = −34. Since x3 is feasible, there’s no reason to branch, and z∗ ≤ z3 is our current best upper bound on the optimal value of the objective function. Iteration 3 (L4 = L1 ∧ (x2 ≥ 4)): The solution to the LP L3 is x3 = (1.8, 4) with the objective function value z3 = −41. Only x31 is fractional, so it is the branching variable. We form the problem L5 by adding the constraint x1 ≤ 1 and the problem L6 by adding the constraint x1 ≥ 2. Iteration 4 (L5 = L4 ∧ (x1 ≤ 1)): The solution to the LP L5 is x5 = (1, 4.444) with the objective function value z5 = −40.555. Only x52 is fractional, so it is the branching variable. Form the problem L7 by adding the constraint x2 ≤ 4 and the problem L8 by adding the constraint x2 ≥ 5. Iteration 5 (L7 = L5 ∧ (x2 ≤ 4)): The solution to the LP L7 is x7 = (1, 4) with the objective function value z7 = −37. Since x7 is feasible, there is no need to branch, and since z7 < z3, we obtain a new best upper bound z∗ < z7 on the optimal objective function value. Iteration 6 (L8 = L5 ∧ (x2 ≥ 5)): The solution to the LP L8 is x8 = (0, 5) with the objective function value z8 = −40. Since x8 is feasible, there’s no need to branch, and since z8 < z7, we obtain a new best upper bound z∗ < z8 on the optimal objective function value. Iteration 7 (L6 = L4 ∧ (x1 ≥ 2)): The LP is infeasible. Iteration 8 (L2 = L1 ∧ (x2 ≥ 4)): The solution to the LP L2 is x2 = (3, 3) with the objective function value z2 = −39. Since z2 > z8, there is no need to branch. Figure 3.2 . depicts the branch-and-bound tree for this example graphically. The tree proves that the optimal solution is z∗ = −40, and that it is attained at x∗ = (0, 5). Only nine points are explicitly considered by the algorithm, while a complete enumeration would have to consider at least the 25 points that make up the feasible set. Remaining points are only implicitly considered by means of the bounding process, which is why the method is often described as an implicit enumeration. 12 CHAPTER 3. THEORY x1 ≤ 2 x1 ≥ 3 x2 ≤ 3 x2 ≥ 4 x1 ≤ 1 x1 ≥ 2 x2 ≤ 4 x2 ≥ 5 x = (2.25, 3.75) z = −41.25 x = (2, 3.888) z = −41.111 x = (2, 3) z = −34 x = (1.8, 4) z = −41 x = (3, 3) z = −39 x = (1, 4.444) z = −40.555 infeasible x = (1, 4) z = −37 x = (0, 5) z∗ = −40 Figure 3.2: A branch-and-bound tree for the the solution to the example problem (3.1.4 . ). Each node contains the solution to the relaxed problem and its optimal value and a visualisation of the relaxed problem in that node. The shaded regions indicate the feasible areas, the arrows represent the objective function gradient, and the optimal solutions are marked with dots. For clarity, only the relevant branch cuts are shown in each node. 13 CHAPTER 3. THEORY 3.1.4 The Benders Decomposition Algorithm The classical Benders decomposition algorithm introduced in [5 . ] is a method for solving a certain class of optimisation problems. The problem is partitioned into a master programming problem (which can be continuously linear, discrete, or nonlinear) and an LP subproblem. In other words, it is required that the subproblem can be written as an LP. For example, an MILP can be solved using this method by decomposition into a pure ILP master problem and a pure LP subproblem. These problems are then solved iteratively to arrive at the optimal solution to the full problem. Explicitly, one first solves the master problem, fixes the integer variables, and then sends them to the linear subproblem as constants. After solving the subproblem with these variables fixed, new constraints, so called cuts, are generated to feed back into the master problem, which is then re-solved. The hope is that solving these smaller problems will be considerably simpler than solving the full problem. The Benders decomposition algorithm can also be applied to an LP problem, and the goal is then to split a very large model into smaller LP that are quicker to solve. The master variables are the complicating variables, and if after fixing them the subproblem separates further into multiple smaller problems, the result is many smaller problems which are simpler to solve. Consider the problem min x,y z := ax+ by, (3.1.5a) s. t. Ax+By ≥ c, (3.1.5b) x, y ≥ 0. (3.1.5c) Any LP of the type (3.1.1 . ) can be written in this form by partitioning the variables into two groups. Let z∗ be the optimal value of this problem. The variables x are in this case continuous, but they can for example also be in {0, 1} or Z≥0. As previously stated the variables y are required to be linear. Fixing x to x̄ yields the subproblem min y z(x̄) := by + ax̄, (3.1.6a) s. t. By ≥ c−Ax̄, (3.1.6b) y ≥ 0. (3.1.6c) The terms involving x are renamed x̄ and moved to the right hand side to indicate that they are considered constants in the problem (3.1.6 . ). If the values x̄ are feasible in the original problem, the optimal value of the full problem z must be at least as small as the optimal value of (3.1.6a . ): z∗ ≤ z∗(x̄) By associating the constraint (3.1.6b . ) with a dual variable u, the dual of the subproblem can be written as min u u(c−Ax̄) + ax̄, s. t. uB ≤ b, y ≥ 0. Let u(x̄) be the optimal solution to this dual problem. Weak duality gives the bound z(x̄) ≥ u(x̄) · (c−Ax̄) + ax̄ (3.1.8) on the value of z(x̄), and since u(x̄) is feasible for any choice of x̄, the inequality (3.1.8 . ) remains valid if we replace x̄ by x. Here, it is assumed that there exists a feasible solution u to the 14 CHAPTER 3. THEORY subproblem dual, since if the subproblem is infeasible, the original problem is also infeasible for that choice of x̄. In summary, we get the bound z ≥ u(x̄) · (c−Ax) + ax (3.1.9) for the subproblem; this is the classical Benders cut. This procedure assumes that (3.1.6 . ) has a feasible solution. If instead the subproblem is infeasible or unbounded, it is possible to obtain a cut u(x̄) · (c−Ax) ≤ 0, (3.1.10) where u(x̄) is the ray that solves max v v(c−Ax̄), s. t. vB ≤ 0, v ≥ 0. The ray defined by u(x̄) is an unbounded direction and so long as inequality (3.1.10 . ) is true it will remain an unbounded direction. Let S denote the set of points x such that the system of inequalities (3.1.6b . ) and (3.1.6c . ) are consistent. With this notation, the original problem (3.1.5 . ) can be written as min x z∗(x), s. t. x ∈ S, x ≥ 0. This problem can be linearised by introducing a new variable ξ such that ξ ≥ z. Using the cuts (3.1.9 . ) and (3.1.10 . ), the problem (3.1.5 . ) can be written as min x,ξ ξ, s. t. ξ ≥ ax+ ui(c−Ax), i ∈ U, 0 ≥ vi(c−Ax), i ∈ V, x ≥ 0, where U is the set of extreme points and V is the set of extreme rays of the subproblem polytope. Instead of enumerating the full sets U and V , the Benders decomposition algorithm constructs these iteratively by generating the cuts (3.1.9 . ) and (3.1.10 . ) in the subproblem. These are then added to the master problem, which is solved again, and new trial values x̄i are sent to the subproblem. This procedure continues until the objective value of the subproblem is no longer greater than the objective value obtained from the last solution of the master problem. The general structure of the Benders decomposition algorithm is summarised in Figure 3.3 . . There are two main reasons why the Benders decomposition algorithm might be a good choice. First by partitioning the variables in a way that exploits the program structure, the resulting subproblem might be in a class of problems which are known to be simple to solve. Second, the subproblem might decouple into several smaller subproblems. These properties are interesting for more classes of optimisation problems than can be treated by the classical Benders decomposition algorithm. However, the Benders decomposition algorithm relies on the ability to obtain a bound on the optimal value of the full problem from the LP dual. If after decomposing the problem, the subproblem is not an LP, this dual does not exist. In Section 3.2 . , a new type of dual is introduced that can be written for a general optimisation problem. This is the first step to a generalised Benders decomposition algorithm. 15 CHAPTER 3. THEORY Start Master problemInfeasible problem Fix x̄ Subproblem Optimal value z Feasibility Cut Optimality CutFinished infeasible optimal (x̄, ξ̄) infeasible z > ξ̄ feasible z ≤ ξ̄ add cuts Figure 3.3: Flowchart of the general structure of the Benders decomposition algorithm. 3.2 Inference Duality Consider the general minimisation problem min x f(x), (3.2.1a) s. t. x ∈ S, (3.2.1b) x ∈ D, (3.2.1c) with domain D, the feasible set S, and objective function f : Ω→ R for some superset Ω of D and S. It need not be true that S is a subset of D, nor D of S. For example, a general LP can be described by the objective function f(x) = cx, the feasible set S = {x ∈ Rn | Ax ≥ b}, and the domain D = Rn. Choosing instead the domain D = Zk × Rn−k, for some k ≤ n, yields a general MILP. Other classes of optimisation problems can similarly be represented by appropriate choices of f , S, and D. We will in this section outline a generalisation of LP duality which is applicable in this more general setting, following the description by Hooker and Ottosson [24 . ]. Some notation is needed before the inference dual can be stated. Let P and Q be two propositions whose truth or falsehood are functions of x. This gives the following definition: Definition 3.2.1. P implies Q with respect to D (notated P D−→ Q) if Q is true for any x ∈ D for which P is true. The inference dual is defined as a maximisation problem over a proof family, where we want to find the proof which proves for x in the domain the largest lower bound of f(x) over the feasible 16 CHAPTER 3. THEORY set. Mathematically, the inference dual is written as max β β, (3.2.2a) s. t. x ∈ S D−→ f(x) ≥ β, (3.2.2b) where β can be seen as the largest lower bound on the function f(x) for any x ∈ D that can be inferred from x ∈ S. For convenience we assume that the problem has an optimal value, where we allow β∗ = −∞ as a possible value if the problem is unbounded, and β∗ =∞ if it is infeasible. With these conventions we get a form of strong duality, i.e. the problem (3.2.1 . ) has the same value as the dual problem (3.2.2 . ), which is true by definition. A solution to the inference dual is a proof that the inequality f(x) ≥ β holds. Hooker and Ottosson [24 . ] give two steps for solving the inference dual: i. Identify inference rules that are complete for the type of constraints in the problem, i.e. they can be used to infer any valid implication of the form f(x) ≥ z. ii. Use the rules to prove optimality. The exact procedure differs for each class of problem. In Section 3.2.1 . this procedure is demonstrated for a LP, where we see that this yields the standard LP dual. In Section 3.5 . this is done for a MILP by inspecting a branch-and-bound tree. 3.2.1 The Inference Dual of a Linear Program The classical LP dual is a special case of the inference dual, where the inference rule is to take a linear combination of the constraints. The proof that this is a complete inference rule is essentially the same as the classical separation lemmas for polyhedra [20 . ]. Consider a linear minimisation problem: min x cx, s. t. Ax ≥ a, x ≥ 0, and the corresponding inference dual max z z, (3.2.4a) s. t. (Ax ≥ a, x ≥ 0) Rn −−→ cx ≥ z, (3.2.4b) where Ax ≥ a and x ≥ 0 defines the feasible set and Rn is the domain of x. This constraint set can be rewritten by the following observation. Theorem 3.2.1 (Linear implications). Ax ≥ a implies cx ≥ z if and only if Ax ≥ a is infeasible or there is a real vector u ≥ 0 for which uA ≤ c and ua ≥ z. Theorem 3.2.1 . is essentially just a reformulation of the usual separation lemmas for convex polyhedra. Using Theorem 3.2.1 . , we get that (3.2.4 . ) is equivalent to finding a non-negative linear combination of Ax ≥ a that dominates cx ≥ z, i.e., uA ≤ c and ua ≥ z, and maximises z. This is the same as the classical dual for an LP: max ua, s. t. uA ≤ c, u ≥ 0. 17 CHAPTER 3. THEORY The vector u can be seen as encoding a proof of optimality, since using it we can deduce the optimal value of z. 3.3 The Logic-Based Benders Decomposition Algorithm The basic idea behind the logic-based Benders decomposition algorithm [24 . ] is the same as in the classical Benders decomposition algorithm, but in a more abstract setting. Variables in the problem are partitioned into two vectors, x and y. A subproblem containing only the variables y is obtained by fixing the variables x to some trial values x := x̄. If the solution to the subproblem reveals that the trial assignments are either infeasible or result in a non-optimal objective value for the full problem, the reason why is identified using the constraint set. New trial values for x are then generated using this information. One continues iteratively in this way until an optimal solution is found, or until the problem is found to be infeasible. The classical Benders decomposition algorithm uses the LP dual to get this information. In the logic-based Benders decomposition algorithm an inference dual (3.2.1 . ) is instead used, which is the problem of inferring a strongest possible bound from the constraint set. The solution to the dual is a proof that the bound is valid for x = x̄. Identifying for which other values of x the same reasoning holds, we get a valid bound on the objective value as a function of x. Since an inference dual is used instead of an LP dual, the logic-based Benders decomposition algorithm can be used to solve more general optimisation problems. Chu and Xia [8 . ] define a valid Benders cut by the following: i. The cut must exclude the current solution to the master problem if it is not globally feasible. ii. The cut must not exclude any globally feasible solutions. The first property guarantees finite convergence, since in the worst case scenario it results in a complete enumeration of the master variables. If the master problem is finite, the enumeration terminates in a finite number of steps. The second property guarantees optimality since no feasible solutions are removed from the problem. This means that the optimal solution remains in the problem. Consider a general optimisation problem of the form min x,y f(x, y), s. t. (x, y) ∈ S, x ∈ Dx, y ∈ Dy. The logic-based Benders decomposition algorithm goes as follows: Fix x to some trial value x̄ ∈ Dx, resulting in the subproblem min y f(x̄, y), s. t. (x̄, y) ∈ S, y ∈ Dy. The inference dual, see (3.2.2 . ), to the above subproblem, max β β, (3.3.3a) s. t. (x̄, y) ∈ S Dy−−→ f(x̄, y) ≥ β, (3.3.3b) 18 CHAPTER 3. THEORY is finding the tightest bound β on f(x̄, y) that can be inferred from (x̄, y) ∈ S. The next step is to somehow derive a function β(x) that gives a valid bound on the value of the subproblem for any value x̄ ∈ Dx. How this bound is generated differs between classes of problems, but the general idea is the following: Let β∗ be the optimal value of (3.3.3 . ), whose solution is a proof of the fact that β∗ is a lower bound on the optimal value when x = x̄. Using this same line of argument for other values of x yields a valid lower bound as a function of x. This bounding function βx̄(·) yields a Benders cut z ≥ βx̄(x) generated from the trial assignment x = x̄. The algorithm continues in the same way as the classical Benders decomposition. In each iteration the master problem has the form min z,x z̄ := z,, s. t. z ≥ βx̄i(x), i = 1, . . . ,K, x ∈ Dx, where x̄1, . . . , x̄K are the trial values obtained from earlier iterations and βx̄i(x) are the correspond- ing bounding functions. A new trial value (z̄K+1, x̄K+1) is obtained by solving the subproblem for this trial value. The algorithm terminates when the optimal value of the subproblem β∗ is equal to z̄. 3.4 Logical Clauses and Resolution This section will give a brief introduction to propositional logic so that it can be used to define the logic-based Benders decomposition algorithm for a mixed binary linear program (MBLP). For general MILP, propositional logic is not applicable, instead see [22 . ]. Propositional logic consists of formulae containing atomic propositions xj which can be either true or false. These can then be connected to create formulae with for example ‘and’, ‘or’, and ‘not’. These are written in the following way: the conjunction ‘x1 and x2’ is written as x1 ∧ x2, the disjunction ‘x1 or x2’ is written as x1 ∨ x2, and the negation ‘not x1’ is written as ¬x1. A formula F1 is said to imply another formula F2 if for any xj that make F1 true, xj also make F2 true. A literal is an atomic proposition or its negation, and a logical clause is a disjunction of literals such that the clause is true if any of the literals it contains are true. The empty clause is the clause containing no literals; it is false by necessity. We will need the following lemma from Hooker and Dawande [22 . ]: Lemma 3.4.1. A clause C1 implies clause C2 if and only if all the literals in C1 occur in C2. C1 is then said to absorb C2. A method for deriving all implications of a given set of clauses is the resolution algorithm introduced by Quine [40 . , 41 . ]. A brief overview of the resolution algorithm is given below. Algorithm 1 A brief overview of the resolution algorithm. A resolution rule is used to generate a new clause from two other clauses which contain exactly one atomic proposition xj that occurs positively in one clause (xj) and negatively in the other clause (¬xj). This new clause is called the resolvent, and consists of the disjunction of all literals from either clause except for xj and ¬xj (this is the resolution rule). Quine showed that the resolution algorithm is complete, i.e. if it begins with a set S of clauses and terminates with S′, then any clause implied by S is absorbed by some clause in S′. Then S is false if and only if S′ contains the empty clause. 19 CHAPTER 3. THEORY If after performing a resolution step, as described in Algorithm 1 . , on a set S we receive the empty set, this proves that there is no solution to the original set. As stated above the empty set is always false, so any set of clauses that imply the empty set must also be false. 3.4.1 Resolution Example Consider three atomic propositions x1, x2, and x3. Write down a set of clauses containing these propositions. The resolution method can then be used to prove that there is no solution to this set. As an example consider the set ¬x1 ∨ ¬x2 ∨ ¬x3, (a) x3, (b) x2, (c) ¬x2 ∨ ¬x3, (d) x1 ∨ x3, (e) x1 ∨ x2. (f) A resolution proof that the set (a . )–(f . ) of clauses is false can be performed as follows. Step 1. The resolution of clause (a . ) and (b . ) is (g): ¬x1 ∨ ¬x2. Step 2. The resolution of clause (d . ) and (e . ) is (h): x1 ∨ ¬x2. Step 3. The resolution of clause (g . ) and (c . ) is (i): ¬x1. Step 4. The resolution of clause (h . ) and (f . ) is (j): x1. Step 5. The resolution of clause (i . ) and (j . ) is (k): ∅. Since we get the empty clause, this proves falsehood for the set of clauses (a . )–(f . ). This procedure is visualised in Figure 3.4 . . (k . ) ∅ (i . ) ¬x1 (g . ) ¬x1 ∨ ¬x2 (a . ) ¬x1 ∨ ¬x2 ∨ ¬x3 (b . ) x3 (c . ) x2 (j . ) x1 (h . ) x1 ∨ ¬x2 (d . ) ¬x2 ∨ ¬x3 (e . ) x1 ∨ x3 (g . ) x1 ∨ x2 Figure 3.4: Visualisation of the resolution proof performed in this section. The clauses are associated with leaf nodes, and the resolvents with non-leaf nodes. The arrows show which two clauses are used to produce the corresponding resolvent. 20 CHAPTER 3. THEORY 3.5 An LBBD Method for Mixed binary Linear Programming This section broadly follows the work in Hooker and Dawande [22 . ] and Hooker [25 . ]. We consider an MBLP with n variables xj , where the binary variables correspond to j = 1, . . . , k and the continuous variables correspond to j = k + 1, . . . , n: min z := cx (3.5.1a) s. t. Ax ≥ a (3.5.1b) − x ≥ − b (3.5.1c) x ≥ 0 (3.5.1d) xj ∈ {0, 1} j = 1, . . . , k (3.5.1e) xj ∈ R j = k + 1, . . . , n (3.5.1f) Inequality (3.5.1b . ) contains the constraints of the problem and the bj is an upper bound on xj , which is taken to be 1 for the binary variables. In each node of the branch-and-bound tree a relaxed problem of the following form is solved: min z = cx (3.5.2a) s. t. Ax ≥ a (3.5.2b) Hx ≥ h (3.5.2c) − x ≥ − b (3.5.2d) x ≥ 0 (3.5.2e) The constraints Hx ≥ h correspond to the active branch cuts, which have the form of xj ≤ 0 or xj ≥ 1. Describing the variable bounds explicitly with b and the branch cuts with H and h in this way will be convenient for the sensitivity analysis performed in this section. The sensitivity analysis consists of two parts, where the first part is to recover a proof scheme that the optimal value of z is z∗, i.e. z ≥ z∗, from the solution to the primal problem. The second part of the analysis consists of fixing this proof and and investigating under which perturbations of the problem this proof remains valid. We will present two different proof schemes, which both infer bounds from the branching tree used to solve the primal problem. The first proof scheme is more easily motivated, but results in very complicated cuts. Lower bounds from the relaxed node problems are valid for the original problem. For leaf nodes this is the best value we can infer, while for a non-leaf node we get that the smallest bound from its children is a valid bound. Consider a node i ∈ I, where I are all the non-leaf nodes. Then let Ji be the children of node i and z∗i be optimal solution in node i. For a leaf node j let LBj = z∗j , and for a non-leaf node i ∈ I let the inferred bound be LBi = max {z∗i , min {LBj | j ∈ Ji}} . By starting at the bottom and going up through the tree we can recursively get a valid bound for the root node which is valid for the original MBLP, and it can be used to prove optimality after a perturbation of the right-hand side. The bounding function is then a piecewise linear function, containing nested min and max functions. For large trees, this function can be difficult to interpret, since it contains information of all the nodes. A full description of this proof scheme is presented in [26 . ]. 21 CHAPTER 3. THEORY The second proof scheme only uses information from leaf nodes, which gives a more practical way of generating the bounding function. Violated inequalities are associated with each leaf node. In feasible nodes we start with z < z∗, and in infeasible nodes with one of the violated constraints. From these we derive surrogate inequalities—see Section 3.5.1 . —which we associate with the respective nodes. These surrogate inequalities are violated in their associated leaf nodes. Logical clauses are inferred from the surrogate inequalities, and these clauses are proven to be inconsistent using a resolution proof. This is used as a refutation proof that the inequality z∗ is the optimal value, i.e. a proof that z < z∗ is false. This proof remains valid if the violated inequalities at each leaf node continue to imply the logical clauses used in the resolution proof. More details are presented in Section 3.5.2 . . In the context of logic-based Benders decomposition algorithm, this is used to generate the bounding functions βx̄(·). Due to the fact that strong duality holds for inference duals, this is a strong dual, yielding a strong cut when used in a Benders decomposition scheme. The resulting Benders cuts from using this method are in general nonlinear. However this is not a problem if the master problem is solved by branching; Hooker and Dawande [22 . ] contains a description on how this is done. The cuts can also be linearised; see [25 . ] for how this can be done. Node A Node B Node D Node E Node F Node H Node I Node G Node C x1 ≤ 2 x1 ≥ 3 x2 ≤ 3 x2 ≥ 4 x1 ≤ 1 x1 ≥ 2 x2 ≤ 4 x2 ≥ 5 −5x1 − 8x2 < −37 −5x1 − 8x2 < −34 −5x1 − 8x2 < −40 5x1 + 9x2 ≤ 45 −5x1 − 8x2 < −39 Figure 3.5: The branching tree resulting from the solution of (3.1.4 . ). Each leaf node (C, D, G, H, I) is associated with an inequality, which violates the fixed variables for that node. Each feasible leaf node (C, D, H, I) is associated with the inequality that the value of z is lower than the objective value for that node. Each infeasible leaf node (G) is associated with one of the violated constraints. If the branching tree is viewed from the bottom it represents a proof of optimality, and a solution to a suitably defined inference dual. Figure 3.5 . shows a proof of optimality for the ILP (3.1.4 . ). The branch cuts in each leaf node violate the inequality associated with that node. The resolution algorithm described in Section 3.4 . then proves that these inequalities are inconsistent. This means that a feasible solution must violate at least one of the inequalities associated with feasible nodes. This is only possible if the optimal value from the branch-and-bound tree is a lower bound on the objective function for any feasible values of x1 and x2. Solving a MIP by branching 22 CHAPTER 3. THEORY method, as in Section 3.1.2 . , solves an inference dual simultaneously. This is an analogue of the simplex method solving the primal and dual problems simultaneously [26 . ]. 3.5.1 Surrogate Inequalities Consider the tree obtained by solving (3.5.1 . ) using a branch-and-bound, and let T be the set of leaf nodes. For a leaf node t ∈ T , define a partial assignment as disjoint subsets J t1, J t0 of J = 1, . . . , k, such that J t1 is the set containing all j for which xj are fixed to 1 and J t0, analogously, those xj which are fixed to 0. If xj = 1 then the literal xj is said to be true and analogously if xj = 0 the negation ¬xj is true. This partial assignment can be uniquely associated with the weakest clause that it falsifies, i.e. , Ct = ∨ j∈Jt 0 xj ∨ ∨ j∈Jt 1 ¬xj , that one of the values of xj differ from the fixed values in the node. The set of clauses Ct for t ∈ T implies that there is an x which is not contained in any of the leaf nodes, i.e. it differs from the fixed values for all nodes. This statement is proven to be false by the resolution algorithm described in 3.4 . . Branching on xj creates two clauses that contain exactly one literal which is positive in one clause (xj) and negative in the other clause (¬xj). By starting at the bottom of the tree we can associate with the parent node the resolvent of the clauses of its children. The resolvent will contain neither xj nor ¬xj , and by continuing in this way all xj are successively removed, and the empty clause will be associated with the root node. Since the resolution algorithm is complete, it proves that the clauses Ct for all leaf nodes are inconsistent. The missing piece is how to get inequalities for each leaf node t that imply Ct. Since the set of these clauses for the leaf nodes can be proven to be inconsistent, a system of such inequalities must also be inconsistent. These inequalities are constructed such that if they are inconsistent they prove that z∗ is the optimal value of the problem (3.5.1 . ). Consider the relaxed node problem in (3.5.2 . ). Let z̄ be the best objective value found from previous nodes, and possibly z̄ =∞ if no solution has been found. In each node there are three cases: Case 1. (3.5.2 . ) is feasible and its value ẑ is not better than the current bound z̄, i.e. ẑ ≥ z̄. Case 2. (3.5.2 . ) is feasible and its value ẑ is better than the current bound z̄, i.e. ẑ < z̄. Case 3. (3.5.2 . ) is infeasible. For each of these cases we associate the node with a surrogate inequality which is violated by the fixed variables, so that they imply Ct. These inequalities are constructed from the dual solutions in each node. Associating dual variables to the relaxation (3.5.2 . ) as λ corresponding to (3.5.2b . ), µ corresponding to (3.5.2c . ), and ν corresponding to (3.5.2d . ), results in the following dual LP: max λ,µ,v λa+ µh− νb, (3.5.3a) s. t. λA+ µH − ν ≤ c, (3.5.3b) λ, µ, ν ≥ 0. (3.5.3c) Each of the Cases 1 . –3 . are considered separately and results in different surrogate inequalities. 23 CHAPTER 3. THEORY (Case 1. . ) z̄ is the minimum value of the objective function. By adding cx < z̄ to the original constraint set, the resulting constraint set becomes infeasible. This can be shown by writing down the following system: − cx ≥ ε− z̄, (3.5.4a) Ax ≥ a, (3.5.4b) Hx ≥ h, (3.5.4c) − x ≥ − b, (3.5.4d) where ε > 0 is added to remove the strict inequality. Consider a linear positive combination of inequalities (3.5.4a . –d . ) with multipliers (1, λ, µ, ν). If the above system has a feasible solution, it implies that such a combination also has a solution. The proof of inconsistency is the content of Lemma 3.5.1 . . Lemma 3.5.1. A linear positive combination of inequalities (3.5.4a . –d . ) with specified multipliers (1, λ, µ, ν) has no solution if (λ, µ, ν) is a feasible solution to the dual problem (3.5.3 . ). Proof. Assume there is a feasible solution to the inequality 1(−cx) + λ(Ax) + µ(Hx) + ν(−x) ≥ 1(ε− z̄) + λ(a) + µ(h) + ν(−b). By reordering, we obtain (λA+ µH − ν − c)x ≥ ε+ (λa+ µh− νb− z̄), (3.5.5) and since the dual solution (λ, µ, ν) satisfies the inequalities λA+ µH − ν − c ≤ 0, λa+ µh− νb− z̄ ≥ 0, the surrogate inequality (3.5.5 . ) implies 0 ≥ (λA+ µH − ν − c)x ≥ ε+ (λa+ µh− νb− z̄) ≥ ε > 0, which is false. The initial assumption must then be false, which concludes the proof. If a positive linear combination of inequalities has no solution it is implied that the system formed by these inequalities also has no solution. Combining the results above, we obtain an infeasible system: (λA− c)x ≥ λa+ ε− z̄ (3.5.6a) Hx ≥ h (3.5.6b) − x ≥ − b (3.5.6c) In branch-and-bound we only consider values of x such that Hx ≥ h and −x ≥ −b, so the inequality (3.5.6a . ) violates the fixed variables, i.e. it implies Ct. In other words, inequality (3.5.6a . ) is our surrogate inequality. (Case 2. . ) The analysis is analogous as in Case 1 . , with z̄ replaced by ẑ. The surrogate inequality thus obtained is (λA− c)x ≥ λa+ ε− ẑ. (3.5.7) 24 CHAPTER 3. THEORY (Case 3. . ) Since the primal problem (3.5.2 . ) is infeasible, there exists a non-negative vector of dual variables (λ, µ, ν) such that these dual variables prove infeasibility, that is λA+ µB − ν ≤ 0, λa+ µb− νh > 0. Like in the previous cases λAx ≥ λa, Hx ≥ h, − x ≥ − b. is infeasible. Hence the surrogate inequality is given by λAx ≥ λa. (3.5.8) The surrogate inequalities (3.5.6a . ), (3.5.7 . ), and (3.5.8 . ) all imply Ct in their respective leaf nodes. 3.5.2 Sensitivity Analysis For the logic-based Benders decomposition algorithm, we want to investigate for which perturba- tions the inequalities (3.5.6a . ), (3.5.7 . ), and (3.5.8 . ) are still valid for the refutation proof that z∗ is the optimal value, i.e. when they still imply Ct. In this algorithm only the right hand side a is perturbed. We will only consider MBLP and not general MILP. For a complete description of the general integer case, see [22 . ]. Consider thus an MBLP which has been decomposed so that the variables are split into two groups: min z = cx+ dy, s. t. Ax+By ≥ a, − x ≥ − b, x ≥ 0, xj ∈ {0, 1} , j = 1, . . . , k, y ∈ {0, 1}m . Here, x ∈ Zk × Rn−k are the subproblem variables, and y are the master variables. After fixing the master variables to y = ȳ, they can be moved to the right hand side and treated as constants: min z = cx+ dȳ, (3.5.10a) s. t. Ax ≥ a−Bȳ, (3.5.10b) − x ≥ − b, (3.5.10c) x ≥ 0, (3.5.10d) xj ∈ {0, 1} , j = 1, . . . , k. (3.5.10e) 25 CHAPTER 3. THEORY We obtain our surrogate inequalities by applying inequalities (3.5.6a . ), (3.5.7 . ), and (3.5.8 . ) to our perturbed constraint (3.5.10b . ): (Case 1. . ) (λA− c)x ≥ λ(a−Bȳ)− z̄ + ε (3.5.11a) (Case 2. . ) (λA− c)x ≥ λ(a−Bȳ)− ẑ + ε (3.5.11b) (Case 3. . ) λAx ≥ λ(a−Bȳ) (3.5.11c) The next step is to see for which other values of y these inequalities continue to imply Ct. A necessary and sufficient condition for when an inequality implies a clause of the form Ct is the content of the following lemma1 . . Lemma 3.5.2. Consider x ∈ {0, 1}k × Rn−k, a ∈ Rn, and α ∈ R. Let J = {1, . . . , n} and J0, J1 be disjoint subsets of {1, . . . , k}. Let a+j = max {0, aj}, j ∈ J . Moreover, assume that the variables x are bounded as 0 ≤ xj ≤ hj for all j ∈ J , where hj = 1 for j ∈ {1, . . . , k}. Then, the inequality ∑ j∈J ajxj ≥ α implies the clause Ct = ∨ j∈J0 xj ∨ ∨ j∈J1 ¬xj if and only if it holds that ∑ j∈J1 aj + ∑ j /∈J0 ⋃ J1 a+j hj < α. (3.5.12) Proof. We will first rewrite the first statement to show that it is equivalent to (3.5.12 . ). By contraposition ∑ j∈J ajxj ≥ α −→ Ct ⇐⇒ ¬Ct −→ ∑ j∈J ajxj < α. This is further equivalent to the inequalities ∑ j∈J ajxj < α, for all x such that Ct is false, which is true if and only if max ∑ j∈J ajxj | ¬Ct  < α. The next step is to show that the left-hand side (3.5.13 . ) is equal to the left-hand side of (3.5.12 . ). Fix x to x̄ where x̄j = 0 for j ∈ J0 and x̄j = 1 for j ∈ J1. Then the negation of Ct, i.e. , ¬Ct = ∧ j∈J0 ¬xj ∧ ∧ j∈J1 xj is always true for x = x̄. Further, it holds that∑ j∈J aj x̄j = ∑ j∈J1 aj + ∑ j /∈J0 ⋃ J1 aj x̄j . The result then follows from the fact that the maximum of ajxj is a+j hj for all j ∈ J . 1An alternate proof of this lemma can be found in Appendix C . . 26 CHAPTER 3. THEORY If the surrogate inequalities still imply Ct for a new trial assignment of y, the proof of optimality outlined in the previous section still holds. From Lemma 3.5.2 . it follows that the surrogate inequalities imply Ct for different values of y if and only if the following inequalities are satisfied: (Case 1. . ) ∑ j∈J1 (λAj − cj) + ∑ j /∈J0 ⋃ J1 (λAj − cj)+hj ≤ λ(a−By)− z̄ (3.5.13a) (Case 2. . ) ∑ j∈J1 (λAj − cj) + ∑ j /∈J0 ⋃ J1 (λAj − cj)+hj ≤ λ(a−By)− ẑ (3.5.13b) (Case 3. . ) ∑ j∈J1 λAj + ∑ j /∈J0 ⋃ J1 (λAj) +hj < λ(a−By) (3.5.13c) After inequalities for each node in the solution tree have been generated, they are combined to obtain a Benders cut. If all of the above inequalities hold in each leaf node for y it is possible to prove that z(ȳ) is still a lower bound on the value of the subproblem. For each node t ∈ T , depending on which case arose, let It denote the corresponding surrogate inequality (3.5.13a . –c . ). Let ξ be the variable introduced to the master problem by the decomposition. The Benders cut is then given by ξ ≥ f(y) + Cmin + (z(ȳ)− Cmin)1{ ∧ t∈T It}(y), where T is the set of leaf nodes, Cmin is the smallest possible value that the subproblem can attain for any fixed x and provides a default bound. This cut is then used in the logic-based Benders decomposition algorithm as described in Section 3.3 . . 27 4 Mathematical Formulation of the Problem The problem is formulated with a three-index vehicle flow formulation, first introduced by Golden, Magnanti, and Nguyen [18 . ]. As a starting point, we used the model presented by Ruffieux [44 . ], which we extended to include variable battery capacities. Definitions for the different sets, parameters, and variables used to describe the model are presented in Table 4.1 . . The full model is presented in Section 4.1 . . This model is then rewritten in Section 4.2 . to facilitate our decomposition. Table 4.1: Definitions of sets, parameters, and variables used to model the problem. Notation Description Sets Vc The set of all nodes containing customers Vr The set of all nodes containing recharge stations {0} The depot node V = {0} ∪ Vc ∪ Vr The set of all nodes in the graph K = {1, . . . ,K} The set of all vehicles A ⊆ V × V The set of all arcs in the graph N The set of all battery sizes Parameters K Size of vehicle fleet (number of vehicles) T Latest time to return to the depot [h] U Cargo storage capacity in each vehicle [kg] Qn Battery capacity for battery type n ∈ N [kWh] cfu/cel Cost of fuel/electricity [€/l]/[€/kWh] rfu/rel Consumption rate of fuel/electricity [l/h]/[kWh/h] cw n Cost of choosing battery type n ∈ N [€] tij Travelling time over arc (i, j) ∈ A [h] si Service time at node i ∈ V [h] pi Demand of cargo in node i ∈ V [kg] dij Length of arc (i, j) ∈ A [h] ei/li Earliest/latest time at which the service in node i ∈ V can start [h] MT A large enough time (MT = 2T suffices) [h] Continued on next page 29 CHAPTER 4. MATHEMATICAL FORMULATION OF THE PROBLEM Notation Description MU A large enough cargo capacity (MU = U suffices) [kg] MQ A large enough battery capacity (MQ = 2maxn∈N Qn suffices) [kWh] Variables xkij = 1 if arc (i, j) ∈ A is used by vehicle k ∈ K; = 0 otherwise τki Arrival time at node i ∈ V \ {0} for vehicle k ∈ K [h] uki Amount of cargo in vehicle k ∈ K at arrival in node i ∈ V [kg] qki Battery level for vehicle k ∈ K upon arrival in node i ∈ V [kWh] zfu,k ij /zel,k ij Fuel/electricity used on arc (i, j) ∈ A by vehicle k ∈ K [l]/[kWh] wkn = 1 if battery type n ∈ N is used by vehicle k ∈ K; = 0 otherwise 4.1 A Mixed Binary Linear Optimisation Model The three-indexed vehicle flow formulation of our VRP is given by the following: min x,τ ,u,q,zfu,,zel,,w ∑ k∈K ∑ (i, j)∈A ( cfuzfu,k ij + celzel,k ij ) + ∑ k∈K ∑ n∈N cw nw k n, (4.1a) s. t. ∑ k∈K ∑ j∈V:(i, j)∈A xkij = 1, i ∈ Vc, (4.1b) ∑ j∈V:(i, j)∈A xkij ≤ 1, k ∈ K, i ∈ V, (4.1c) ∑ j∈V:(i, j)∈A xkij = ∑ j∈V:(j, i)∈A xkji, k ∈ K, i ∈ V, (4.1d) ∑ k∈K ∑ j∈Vr∪Vc xk0j ≤ K, (4.1e) ei ≤ τki ≤ li, k ∈ K, i ∈ V, (4.1f) τki ≤ T − (si+ti0), k ∈ K, i ∈ Vr ∪ Vc, (4.1g) τkj − τki ≥ (si+tij)−MT ( 1−xkij ) , k ∈ K, i ∈ V \ {j}, j ∈ Vr ∪ Vc, (4.1h) uk0 ≤ U, k ∈ K, (4.1i) uki − ukj ≥ pixkij −MU ( 1−xkij ) , k ∈ K, i ∈ V \ {j}, j ∈ Vr ∪ Vc, (4.1j) zfu,k ij /rfu + zel,k ij /rel = dijx k ij , k ∈ K, (i, j) ∈ A, (4.1k) qki − qkj ≥ z el,k ij −MQ ( 1−xkij ) , k ∈ K, j ∈ V \ {i}, i ∈ Vc, (4.1l)∑ n∈N Qnw k n − qkj ≥ z el,k ij −MQ ( 1−xkij ) , k ∈ K, j ∈ V \ {i}, i ∈ Vr ∪ {0} , (4.1m)∑ n∈N wkn = 1, k ∈ K, (4.1n) xkij ∈ {0, 1}, k ∈ K, (i, j) ∈ A, (4.1o) wkn ∈ {0, 1}, k ∈ K, n ∈ N , (4.1p) zfu,k ij , zel,k ij ≥ 0, k ∈ K, (i, j) ∈ A, (4.1q) τki , u k i , q k i ≥ 0, k ∈ K, i ∈ V. (4.1r) 30 CHAPTER 4. MATHEMATICAL FORMULATION OF THE PROBLEM The constraints (4.1b . ) ensure that at most one vehicle visits each customer, the constraints (4.1c . ) that at most one vehicle visits each node, the constraints (4.1d . ) state that if a vehicle enters a node then it also must leave it, the constraint (4.1e . ) ensures that at most K vehicles leave the depot (node 0). Furthermore, the constraints (4.1f . ) guarantees that the time windows are followed, while the constraints (4.1g . ) and (4.1h . ) model the travel times between consecutive nodes in the routes. The cargo constraints are modelled by the constraints (4.1i . ) and (4.1j . ) which say that a vehicle can not carry more cargo than its storage capacity, and that the cargo in a vehicle decreases by the demand of each visited node. The constraints (4.1k . ) define the distance travelled on each used arc as the sum of the distances travelled using fuel and electricity, respectively. The constraints (4.1l . ) keep track of the vehicles’ battery charge levels, where constraints (4.1m . ) resets the charge to max each time it leaves the depot or visits a recharge station. The constraints (4.1n . ) ensure that each vehicle chooses exactly one battery type. There are no specific subtour elimination constraints, i.e. that each vehicle must follow a continuous route containing the depot, but these are implied by the time and cargo constraints (4.1g . –j . ). The variables zel,k ij appear in the objective function with a positive coefficient, and since problem (4.1 . ) is a minimisation problem it would at first seem like that it would be profitable to set zel,k ij to a small as possible value. The constraints (4.1l . ) would in that case not guarantee that in (4.1l . ) the battery use zel,k ij is equal to the difference in battery levels qki − qkj . This is in fact not the case since a reasonable assumption is that celrel < cfurfu, and using more electricity would always be cheaper since this uses less fuel. This forces zel,k ij to be as large as possible to minimise the use of conventional fuel, and we get equality between battery use and the difference between battery levels in (4.1l . ). Of note is that we can divide the variables and constraints in four different groups: One group with the variables xkij , one group with the time variables τki , one group with the cargo variables uki , and one group with the variables zfu,k ij , zel,k ij , and wkn. The only shared variables are xkij , which are in fact present in all four groups. This will be used when we formulate our decomposition of the model. Before we get to the decomposition, however, we will need to slightly rewrite the model, which is the contents of the next section. 4.2 A Reformulation of the Mixed Binary Linear Model The model in the previous section does not have a explicit cost related to how the route looks. Instead there are the requirement that the fuel and electricity use covers all the routes. When decomposing the problem in the next chapter, we would like there to be some information about the objective function in the master problem. This would give the master problem a reason for choosing shorter routes, instead of unnecessary long routes. The shortest route is not necessarily the optimal one, but it seems like a good place to start. Further, it is not necessary to keep track of how much fuel is used, i.e. the variables zfu,k ij , since we can assume that fuel is used for all arcs xkij that is not driven on by electricity. It is not clear that changing the model in this way does not affect the speed at we can solve the problem, but it makes it simpler when we later decompose the problem. The model is reformulated so that it no longer contains zfu,k ij . This is achieved by using the constraints (4.1k . ) to remove the variables zfu,k ij from the model. The variables zfu,k ij is only present in the objective, and the new objective now instead contains xkij and zel,k ij . This results in a explicit cost on the variables xkij which is used when solving the problem using the logic-based Benders decomposition algorithm. Renaming zel,k ij = zkij , results in the following reformulated 31 CHAPTER 4. MATHEMATICAL FORMULATION OF THE PROBLEM model: min x,τ ,u,q,z,w ∑ k∈K ∑ (i, j)∈A cfurfudijx k ij + ∑ k∈K ∑ (i, j)∈A ( cel − rfu rel c fu ) zkij + ∑ k∈K ∑ n∈N cw nw k n (4.2a) s. t.∑ k∈K ∑ j∈V:(i, j)∈A xkij = 1, i ∈ Vc, (4.2b) ∑ j∈V:(i, j)∈A xkij ≤ 1, k ∈ K, i ∈ V, (4.2c) ∑ j∈V:(i, j)∈A xkij = ∑ j∈V:(j, i)∈A xkji, k ∈ K, i ∈ V, (4.2d) ∑ k∈K ∑ j∈Vr∪Vc xk0j ≤ K, (4.2e) ei ≤ τki ≤ li, k ∈ K, i ∈ V, (4.2f) τki ≤ T − (si + ti0), k ∈ K, i ∈ Vr ∪ Vc, (4.2g) τkj − τki ≥ (si + tij)−MT ( 1− xkij ) , k ∈ K, i ∈ V \ {j}, j ∈ Vr ∪ Vc, (4.2h) uk0 ≤ U, k ∈ K, (4.2i) uki − ukj ≥ pixkij −MU ( 1− xkij ) , k ∈ K, i ∈ V \ {j}, j ∈ Vr ∪ Vc, (4.2j) zkij ≤ reldijx k ij , k ∈ K, (i, j) ∈ A, (4.2k) qki − qkj ≥ zkij −MQ ( 1− xkij ) , k ∈ K, j ∈ V \ {i}, i ∈ Vc, (4.2l)∑ n∈N Qnw k n − qkj ≥ zkij −MQ ( 1− xkij ) , k ∈ K, j ∈ V \ {i}, i ∈ Vr ∪ {0} , (4.2m)∑ n∈N wkn = 1, k ∈ K, (4.2n) xkij ∈ {0, 1}, k ∈ K, (i, j) ∈ A, (4.2o) wkn ∈ {0, 1}, k ∈ K, n ∈ N , (4.2p) zkij ≥ 0, k ∈ K, (i, j) ∈ A, (4.2q) τki , u k i , q k i ≥ 0, k ∈ K, i ∈ V. (4.2r) The cost ( cel − rfu rel c fu ) in front of zkij represents the costs incurred by using electricity instead of fuel over a arc. A reasonable assumption is that celrel < cfurfu. Hence, this cost will be negative, and it will be beneficial to use electricity over fuel whenever possible. If instead celrel ≥ cfurfu, the optimal value is attained by zkij = 0 and the problem is the same as a normal VRP. The constraints (4.1k . ) also gives a bound for zkij which is included in the model as constraints (4.2k . ). 32 5 Method A logic-based Benders decomposition algorithm was used to solve the problem defined in Chapter 4 . . The method works by splitting the optimisation model (4.2 . ) into a master problem and a subproblem. This separation was done in such a way that the variables xkij occur in the master problem and the other variables occur in the subproblem. Constraints containing only the variables xkij are the only ones left in the master problem. Considering only these parts, the following problem is obtained min x ∑ k∈K ∑ (i, j)∈A cfurfudijx k ij , (5.0.1a) s. t. ∑ k∈K ∑ j∈V:(i, j)∈A xkij = 1, i ∈ Vc, (5.0.1b) ∑ j∈V:(i, j)∈A xkij ≤ 1, k ∈ K, i ∈ V, (5.0.1c) ∑ j∈V:(i, j)∈A xkij = ∑ j∈V:(j, i)∈A xkji, k ∈ K, i ∈ V, (5.0.1d) ∑ k∈K ∑ j∈Vr∪Vc xk0j ≤ K. (5.0.1e) The problem (5.0.1 . ) will constitute the basis for the master problem, where the full master problems contains additional constraints in the form of cuts. These cuts are generated from the subproblem, which consists of three smaller subproblems. After fixing the values of the variables xkij in (4.2 . ), the model separates into three parts: one problem containing the variables τki that verify that the time window constraints are obeyed, one problem containing the variables uki that verify that the cargo capacity constraints are obeyed, and one problem containing the variables qki , zkij , and wkn which distributes where the vehicles use battery charge instead of conventional fuel. The first two are LPs, so these cuts can be generated using the LP dual. These two problems are similar in that they posses resource constraints. They differ in that cargo amount is decreasing and time is increasing, and also in that the variables representing cargo all have the same bounds whereas the time variables all have different bounds. The third subproblem is to distribute the battery charge, and is a MBLP, so it does not posses an LP dual. Instead an inference dual is used, and the solution method used for this is a branch-and-bound algorithm. Two types of Benders cuts have been formulated for the problem: classical Benders cuts from the LP subproblems, and a special cut using bounds inferred from the constraint set of the MBLP subproblem. This inference from the constraint set is achieved by inspecting the branch-and-bound tree used to solve the subproblem. A overview of the algorithm employed can 33 CHAPTER 5. METHOD be found in Section 5.5 . . The models for the different subproblems are presented in Sections 5.1 . and 5.2 . together with corresponding bounds. These bounds are then combined in Section 5.4 . to create the cut that is added to the master problem. 5.1 Linear Programming Subproblems The time subproblem and the cargo subproblem are both LP problems, which means that the same method can be used to generate cuts from both. The cuts are normal Benders feasibility cuts which are constructed by solving the LP dual. The LP dual and the corresponding feasibility cuts are defined in the following sections. Only feasibility cuts are generated since neither subproblem contains an objective function. 5.1.1 The Time Window Subproblem The time window constraints are enforced in the model by the τki variables, which appears in the constraints (4.1f . –h . ). By isolating these constraints, a time window subproblem can be formulated as min τ 0, (5.1.1a) s. t. ei ≤ τki ≤ li, k ∈ K, i ∈ V, (5.1.1b) τki ≤ T − (si + ti0), k ∈ K, i ∈ Vr ∪ Vc, (5.1.1c) τkj − τki ≥ (si + tij)−MT ( 1− x̄kij ) , k ∈ K, i ∈ V \ {j}, j ∈ Vr ∪ Vc, (5.1.1d) τki ≥ 0, k ∈ K, i ∈ V. (5.1.1e) Since the variables τki do not appear in the objective function (4.1a . ), the subproblem (5.1.1 . ) is a feasibility problem, and since the variables xkij are fixed, the subproblem is an LP. This problem separates even further since each vehicle is independent of the others, resulting in |K| subproblems. For each vehicle k ∈ K, we write down the LP dual to the corresponding subproblem: max α,β,γ,δ ∑ i∈V (αiei − βili) + ∑ i∈Vr∪Vc γi (ti0 + si − T ) + ∑ (i, j)∈A:j /∈{0} δij ( (tij + si)−MT ( 1− x̄kij )) , (5.1.2a) s. t. α0 − β0 − ∑ j∈Vr∪Vc γj ≤ 0, (5.1.2b) αi − βi − γi + ∑ j∈Vr∪Vc: (i, j)∈A δij − ∑ j:(j, i)∈A δji ≤ 0, i ∈ Vr ∪ Vc, (5.1.2c) αi, βi, γi ≥ 0, i ∈ V, (5.1.2d) γi ≥ 0, i ∈ Vr ∪ Vc, (5.1.2e) δij ≥ 0, j ∈ Vr ∪ Vc \ {i}, i ∈ V. (5.1.2f) 34 CHAPTER 5. METHOD The dual variables αi, βi are associated with the constraints (5.1.1b . ), γi with the constraint (5.1.1c . ), δij with the constraint (5.1.1d . ). Since the primal problem is a feasibility problem, the dual is either feasible or unbounded. If it is feasible the primal problem is also feasible, and no cut is generated. Let (ᾱi, β̄i, γ̄i, δ̄ij) be the direction of an unbounded ray to the dual problem so that the objective value is increasing along it. The feasibility cut, see the inequality (3.1.10 . ), from this subproblem is then given by the inequality∑ i∈V ( ᾱiei − β̄ili ) + ∑ i∈Vr∪Vc γ̄i (ti0 + si − T ) + ∑ (i, j)∈A:j /∈{0} δ̄ij ( (tij + si)−MT ( 1− xkij )) ≤ 0. (5.1.3) The inequality (5.1.3 . ) is then added to the master problem (5.0.1 . ). 5.1.2 The Cargo Subproblem The cargo capacity subproblem is a simpler variant of the time window subproblem as described in the previous section. It separates analogously, one subproblem for each k ∈ K, which are formulated as min u 0, (5.1.4a) s. t. uk0 ≤ U, (5.1.4b) uki − ukj ≥ pix̄kij −MU ( 1− x̄kij ) , i ∈ V \ {j}, j ∈ Vr ∪ Vc, (5.1.4c) uki ≥ 0, i ∈ V. (5.1.4d) It is also a feasibility problem due to the lack of objective function. There are two main differences between (5.1.4 . ) and the time window problem (5.1.1 . ). First, all variable bounds (5.1.4b . ) are equal in the cargo subproblem, and second, the cargo variable is counting down between nodes instead of counting up as in the time window problem. These differences do, however, not change how the feasibility cuts are obtained. We start by writing down the LP dual of (5.1.4 . ) as max n,m −nU + ∑ (i, j)∈A:j /∈{0} mij ( pj x̄ k ij −MU ( 1− x̄kij )) , (5.1.5a) s. t. − n+ ∑ j∈Vr∪Vc m0j ≤ 0, (5.1.5b)∑ j∈Vr∪Vc:(i, j)∈A mij − ∑ j:(j, i)∈A mji ≤ 0, i ∈ Vr ∪ Vc, (5.1.5c) n ≥ 0, (5.1.5d) mij ≥ 0, j ∈ Vr ∪ Vc \ {i}, i ∈ V. (5.1.5e) Here n corresponds to the constraint (5.1.4b . ) and mij correspond to the constraints (5.1.4c . ). By the analogous argument as for the time subproblem (5.1.2 . ), the feasibility cut corresponding to vehicle k is given by the inequality − n̄U + ∑ (i, j)∈A:j /∈{0} m̄ij ( pjx k ij −MU ( 1− xkij )) ≤ 0. (5.1.6) The inequality (5.1.6 . ) is then added to the master problem (5.0.1 . ). 35 CHAPTER 5. METHOD 5.2 The Battery Charge Distribution Subproblem The battery charge distribution subproblem is to decide which battery capacity to use for each vehicle, and to find the best way to spend the electricity charge given a specific route, i.e. for a fixed value of x. The variables wki , qki , and zkij , together with their associated constraints, are contained in one subproblem for each vehicle k ∈ K: min q,w,z ∑ (i, j)∈A ( cel − rfu rel c fu ) zkij + ∑ n∈N cw nw k n, (5.2.1a) s. t. zkij ≤ reldij x̄ k ij , (i, j) ∈ A, (5.2.1b) qki − qkj − zkij ≥ −MQ ( 1− x̄kij ) , j ∈ V \ {i}, i ∈ Vc, (5.2.1c)∑ n∈N Qnw k n − qkj − zkij ≥ −MQ ( 1− x̄kij ) , j ∈ V \ {i}, i ∈ Vr ∪ {0} , (5.2.1d)∑ n∈N wkn = 1, (5.2.1e) wkn ∈ {0, 1} , n ∈ N , (5.2.1f) zkij ≥ 0, k ∈ K, (i, j) ∈ A, (5.2.1g) qki ≥ 0 i ∈ V. (5.2.1h) This subproblem is not continuous since wkn ∈ {0, 1}, and therefore an inference dual is used to obtain Benders cuts; see Section 3.5 . . The subproblem is solved by a branch-and-bound algorithm. For a node in the branch-and-bound tree let N0 and N1 be a partial assignment corresponding to the fixed values of the variables wkn. In the node the following relaxed problem is solved: min q,w,z ∑ (i, j)∈A ( cel − rfu rel c fu ) zkij + ∑ n∈N cw nw k n, (5.2.2a) s. t. constraints (5.2.1b . –d . ),∑ n∈N wkn ≥ 1, (5.2.2b) − ∑ n∈N wkn ≥ − 1, (5.2.2c) − wkn ≥ − 1, n ∈ N , (5.2.2d) − wkn ≥ 0, n ∈ N0 (5.2.2e) wkn ≥ 1, n ∈ N1 (5.2.2f) wkn ≥ 0, n ∈ N , (5.2.2g) zkij ≥ 0, (i, j) ∈ A, (5.2.2h) qki ≥ 0, i ∈ V. (5.2.2i) The relaxed subproblem (5.2.2 . ) is feasible as long there exists a solution to the constraints (5.2.2b . –c . ), since it is always possible to choose qki = 0, zkij = 0 in the solution. Branching is done 36 CHAPTER 5. METHOD Root w1 ≥ 1 w1 ≤ 0 w2 ≥ 1 w2 ≤ 0 w3 ≥ 1 w3 ≤ 0 w4 ≥ 1 wn−1 ≤ 0 wn ≥ 1 wn ≤ 0 · · · Figure 5.1: Branching tree for the battery subproblem, where the variables wn are branched on. The structure is due to the constraint (5.2.1e . ), which forces one wn to 1, and the others to be set to 0. Consequently, if any wn is fixed to 1, this results in a feasible solution, and the branching stops. If all but one wn are fixed to zero it is immediately guaranteed that the solution is integral, and thus the bottom two nodes are never reached. by fixing the value of wn, for some n, to either 0 or 1. The structure of the solution tree from the branch-and-bound procedure is simple; see Figure 5.1 . . The branching procedure can never get to the point where all wkn have been set to zero, since a integer solution is guaranteed in the node before by the constraints (5.2.2b . ) and (5.2.2c . ). Because of this it is assured that there exists a feasible solution in each node of the branch-and-bound tree. Cuts are obtained by following the procedure outlined in Section 3.5 . . First the primal problem model (5.2.2 . ) is solved. Then we construct the necessary and sufficient conditions for the optimal objective value of (5.2.2 . ) being a valid lower bound on (5.2.1 . ) for other trial assignments of the master variable. In each node of the branch and bound tree, there is a relaxed node problem. For each of the node problems, write down the corresponding LP dual. From the optimal solution of the primal problem in each node we obtain a bound, and from the solution of the dual problems we obtain the sensitivity analysis. Exactly how this is done is the contents of the rest of this section. We use the theory that we developed in Section 3.5 . . Let µij be dual variables associated with the constraints (5.2.1b . ), πij with the constraints (5.2.1c . ) and (5.2.1d . ), σ1 with the constraint (5.2.2b . ), and σ2 with the constraint (5.2.2c . ). It is not necessary to explicitly write down how the dual variables corresponding to the branch cut constraints (5.2.2e . –f . ) and the integer relaxation constraints (5.2.2b . –d . ). Call these dual variables ψ, and as in (3.5.3 . ), let h and H be the corresponding coefficients. The dual problem for in each 37 CHAPTER 5. METHOD node is given by the following: max µij ,πij ,σ1,σ2,ψ − ∑ (i, j)∈A µijr eldij x̄ k ij − ∑ (i, j)∈A πijMQ ( 1− x̄kij ) + σ1 − σ2 + ψh, s. t. − µij − πij − σ1 + σ2 ≤ ( cel − rfu rel c fu ) , (i, j) ∈ A,∑ j∈V:(i, j)∈A (πij − πji) ≤ 0, i ∈ Vc, − ∑ j∈V:(i, j)∈A πji ≤ 0, i ∈ Vr ∪ {0} , σ1 − σ2 +Qn ∑ (i, j)∈A:i/∈Vc πji + ψH ≤ cw n , n ∈ N , µij ≥ 0, (i, j) ∈ A, µij ≥ 0, (i, j) ∈ A, σ1, σ2 ≥ 0. In (5.2.3 . ), ζ̂ denotes the optimal value of the relaxed problem (5.2.2 . ) in the node, N0 := {n ∈ N | wn is fixed to 0} , and N1 := {n ∈ N | wn is fixed to 1} . As in 3.5.2 . , define the plus operator as a+ = max {0, a}. The inequality, for x = x̄, Ix̄,t(x) from node t is obtained from the solution to (5.2.3 . ), according to Section 3.5.1 . , as Ix̄,t(x) : ∑ (i, j)∈A ( −π̄ij − µ̄ij − ( cel − rfu rel c fu ))+ · reldij + ∑ i∈Vc  ∑ j∈V:(i, j)∈A ( π̄ij − π̄ji )+ ·Qmax + ∑ i∈V\Vc  ∑ j∈V:(i, j)∈A ( −π̄ji )+ ·Qmax + ∑ n∈N1 σ̄1 − σ̄2 +Qn ∑ (i, j)∈A:i/∈Vc π̄ji − cw n  + ∑ n∈N\{N0∪N1} σ̄1 − σ̄2 +Qn ∑ (i, j)∈A:i/∈Vc π̄ji − cw n + ≤− ∑ (i, j)∈A reldij µ̄ijx k ij + ∑ (i, j)∈A MQπ̄ij ( xkij − 1 ) + σ̄1 − σ̄2 − ζ̂, (5.2.4) If the inequality (5.2.4 . ) Ix̄,t(x) is satisfied for some other route x, the surrogate inequalities (3.5.11a . –c . ) imply that the fixed variable values are violated. It is not necessary to write these surrogate inequalities down explicitly, since only the inequalities Ix̄,t(x) are needed. A resolution proof can then be used to prove that these surrogate inequalities are false. This is then a proof that the optimal value ζ̂ is a lower bound for other routes that fulfil all inequalities Ix̄,t(x), where t is a leaf node corresponding to the tree where ζ̂ was obtained. 38 CHAPTER 5. METHOD We also use the fact that the variables qki and zkij can be bounded from above as zkij ≤ reldij , (i, j) ∈ A, qki ≤ max {Qn | n ∈ N} , i ∈ V. These bounds are independent of the choice of value for xkij and can be added to the model in (5.2.2 . ) without affecting the solution. Bounds are inferred from the constraints (5.2.1b . ) for zkij and from the constraints (5.2.1c . –d . ) for qki . Before the cuts for the battery charge distribution subproblems (5.2.1 . ) can be formulated, a globally valid lower bound is needed. If the inequalities Ix̄,t(x) are not fulfilled, the refutation proof that the optimal objective value is a lower bound is no longer valid. The next section formulates two different globally valid bounds and uses them to obtain the cut for the battery charge subproblem. 5.3 Default Lower Bounds The last step is to provide a default lower bound ζmin. If the proof can no longer be used to prove that ζ̂ is a lower bound, this globally valid bound is used instead. We present bounds from two different methods. Which method results in the tightest bound depends on the problem data. The first bound is derived from the maximum amount of charge that can be used by a single vehicle. Solving a relaxed version of the battery problem, min z,w ∑ (i, j)∈A ( cel − rfu rel c fu ) zkij + ∑ n∈N cw nw k n, s. t. zkij ≤ reldijx k ij , (i, j) ∈ A, zkij ≥ 0, (i, j) ∈ A, wkn ∈ {0, 1} , n ∈ N , provides a global bound. Due to our assumptions on the fuel costs and consumption rates the coefficient in front of zkij is negative, so choosing each zkij to equal its maximum is optimal in this problem. The cost for the different battery types are all positive, so choosing wkn all zero is optimal. Since the problem above is a relaxation of the battery subproblem the optimal value of it gives the default bound ζmin := ( cel − rfu rel c fu ) rel max  ∑ (i, j)∈A dijx k ij | xkij ∈ X  , where X are the possible set of routes defined by constraints (5.0.1b . –e . ). This value is of course larger than if we set all xkij = 1, which yields the value ζmin := ( cel − rfu rel c fu ) rel ∑ (i, j)∈A dij . (5.3.2) A better bound can be found by finding the longest possible route. Let the length of this route be dmax. The lower bound is then given by ζmin := ( cel − rfu rel c fu ) reldmax. (5.3.3) 39 CHAPTER 5. METHOD Which of the two default bounds Equations (5.3.2 . ) and (5.3.3 . ) is the better depends on how difficult it is to find the longest route, compared to using a tighter bound. Depending on the problem data, these bounds may both be very poor since they calculate the cost as if it was possible to drive on electricity the whole way. If the battery capacities are small or if there are few recharge stations, a better bound is acquired by counting the number of recharge station and depot nodes. The amount of electricity used can not be greater than ζmin = Qmax(Nr +Nd), where Nr is the number of recharge stations, and Nd is the number of depots. Our model assumes that there is only one depot node. This gives the following valid lower bound on the objective function: Qmax(Nr + 1) ( cel − rfu rel c fu ) . (5.3.4) Let (ζ∗)k be the optimal objective value off the problem (5.2.1 . ) for vehicle k obtained for the fixed value x = x̄ and let T k be the set of all the leaf nodes in the solution tree for that vehicle. The set of Ix̄,t(x) represents the inequalities associated with these leaf nodes. The cut from the battery subproblem is then given by ξ ≥ ∑ k∈K ∑ (i, j)∈A ( cfurfudijx k ij ) + ∑ k∈K ξk, (5.3.5) where, for each vehicle k ∈ K, the lower bounds ξk ≥ (ζ∗)k, if ∧ t∈Tk Ix̄,t(x), ζmin, otherwise, are those that can be proven from the branch-and-bound tree. The default lower bound ζmin can be chosen as either of the two formulated above. Which is preferable depends on the problem data and on available computation time (see further in Section 6.4 . ). 5.4 The Master Problem In this section we combine everything described in this chapter to formulate the master problem, together with all of the cuts. For each vehicle k ∈ K, we obtained three separate subproblems that each yielded valid cuts. These were combined to form the Benders cut that was added to the master problem. Let Ω be the index set of previous iterations, and (β∗)k,ω be the optimal objective value obtained in iteration ω ∈ Ω of the battery subproblem for vehicle k ∈ K. The 40 CHAPTER 5. METHOD complete master problem with all of the cuts is then written as min x,ξ,ξ1,...,ξK ξ, (5.4.1a) s. t. ξ ≥ ∑ k∈K ∑ (i, j)∈A cfurfudijx k ij + ∑ k∈K ξk, (5.4.1b) ξk ≥ (ζ∗)ω,k, if ∧ t∈Tω,k Iω,t(x), ζmin, otherwise, k ∈ K, ω ∈ Ω, (5.4.1c) 0 ≥ − n̄ω,kU + ∑ (i, j)∈A: j /∈{0} m̄ω,k ij ( pjx k ij −MU ( 1− xkij )) , k ∈ K, ω ∈ Ω, (5.4.1d) 0 ≥ ∑ i∈V ( ᾱω,ki ei − β̄ω,ki li ) + ∑ i∈Vr∪Vc γ̄ω,ki (ti0 + si − T ) + ∑ (i, j)∈A: j /∈{0} δ̄ω,kij ( (tij + si)−MT ( 1− xkij )) , k ∈ K, ω ∈ Ω, (5.4.1e) 1 = ∑ k∈K ∑ (i, j)∈A xkij , i ∈ Vc, (5.4.1f) 1 ≥ ∑ j∈V:(i, j)∈A xkij , k ∈ K, i ∈ V, (5.4.1g) ∑ j∈V:(i, j)∈A xkij = ∑ j∈V:(j, i)∈A xkji, k ∈ K, i ∈ V, (5.4.1h) K ≥ ∑ k∈K ∑ i∈Vr∪Vc xk0j , (5.4.1i) xkij ∈ {0, 1} , k ∈ K, (i, j) ∈ A. (5.4.1j) The constraints (5.4.1b . –c . ) contain the optimality cuts generated from the battery subproblem, where (5.4.1c . ) define approximate values of the subproblem objective for each vehicle. Here Tω,k indexes the leaf nodes of the branch-and-bound tree from iteration ω for vehicle k, and It are the inequalities associated with these leaf nodes. If the route xk for vehicle k fulfils the logical clause∧ t∈Tk,r Iω,t(x k) for some iteration ω, where Iω,t(x) are the inequalities (5.2.4 . ), the proof that (ζ∗)ω,k is a lower bound remains a valid proof. If no bounds can be proven, ζmin (see (5.3.2 . ) and (5.3.3 . )), is used as a default lower bound. The constraints (5.4.1d . ) correspond to the feasibility cuts generated from the cargo subproblems, and the constraints (5.4.1e . ) corresponds to the feasibility cuts from the time window subproblem. Finally, the constraints (5.4.1f . –i . ) define feasible routes. To improve the performance of the logic-based Benders decomposition algorithm, a relaxed version, see (5.3.1 . ), of the battery problem is added as an initial cut: ξk ≥ ∑ (i, j)∈A ( cel − rfu rel c fu ) rfudijx k ij 41 CHAPTER 5. METHOD 5.5 Algorithm In this section we go through more of how each iteration of our LBBD algorithm functions. The full logic-based Benders decomposition algorithm is summarised in Algorithm 2 . . Algorithm 2 Logic-based Benders decomposition algorithm for the problem in Chapter 4 . 1: choose an initial route and set x̄kij to the corresponding values 2: let x̄ be the vector containing the values of x̄kij 3: set ζ∗ :=∞, ξ := −∞, where ζ∗/ξ is the current upper/lower bound on the objective value of the complete problem (4.2 . ) 4: while ζ∗ > ξ do 5: for each vehicle k ∈ K do 6: solve the battery subproblem (5.2.1 . ) for x̄ 7: for each leaf node in the solution of the battery subproblem do 8: let It be as in (5.2.4 . ) 9: end for 10: solve the time window (5.1.2 . ) and the cargo subproblem (5.1.5 . ) for x̄ 11: add Benders cuts (5.3.5 . ), (5.1.3 . ), and (5.1.6 . ) to the master problem 12: end for 13: let ζ∗ be the sum of subproblem objective values. 14: solve the master problem (5.4.1 . ) 15: if master problem infeasible then 16: original problem infeasible 17: else 18: set x̄ := optimal route obtained from the master problem 19: set ξ := objective value of route x̄. 20: end if 21: end while Algorithm 2 . starts by choosing an initial route x̄. Let ζ∗ be the current best upper bound on the objective value, which is updated each ti