Citation
Optimization methods in liquids hauling and inventory

Material Information

Title:
Optimization methods in liquids hauling and inventory
Creator:
Dyk, Wesley O. ( author )
Language:
English
Physical Description:
1 electronic file (108 pages) : ;

Subjects

Subjects / Keywords:
Petroleum industry and trade ( lcsh )
Hazardous substances -- Transportation -- United States ( lcsh )
Hazardous substances -- Transportation ( fast )
Petroleum industry and trade ( fast )
United States ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Review:
Crude oil and waste water are two liquid products of the petroleum exploration and production sector which are distributed using trucked transportation. Managing the resources to handle inventory, transportation and disposition of the products is a highly complex task. Seeking operational efficiency is aided by the use of optimization models. Initial approaches to solving these models are reported and have seen some positive results. Model enhancement and updated solution methods are required for their use in operations environments. This thesis includes discussion of the development and application of these models and methods using mathematical programming and heuristics.
Bibliography:
Includes bibliographical references,
System Details:
System requirements: Adobe Reader.
General Note:
Department of Mathematical and Statistical Sciences
Statement of Responsibility:
by Wesley O. Dyk

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
921891965 ( OCLC )
ocn921891965
Classification:
LD1193.L622 2015d D95 ( lcc )

Downloads

This item has the following downloads:


Full Text
OPTIMIZATION METHODS IN LIQUIDS HAULING AND INVENTORY
by
WESLEY O. DYK
B.S., Applied Mathematics, University of Colorado Denver, 2003
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Master of Science
Applied Mathematics
2015


This thesis for the Master of Science degree by
Wesley 0. Dyk
has been approved for the
Department of Mathematical and Statistical Sciences
by
Alexander Engau, Advisor
Stephen Billups, Chair
Burt Simon
Chrystanthos Gounaris
July 20, 2015
n


Dyk, Wesley O. (M.S., Applied Mathematics)
Optimization Methods in Liquids Hauling and Inventory
Thesis directed by Associate Professor Alexander Engau
ABSTRACT
Crude oil and waste water are two liquid products of the petroleum exploration
and production sector which are distributed using trucked transportation. Managing
the resources to handle inventory, transportation and disposition of the products is
a highly complex task. Seeking operational efficiency is aided by the use of opti-
mization models. Initial approaches to solving these models are reported and have
seen some positive results. Model enhancement and updated solution methods are
required for their use in operations environments. This thesis includes discussion of
the development and application of these models and methods using mathematical
programming and heuristics.
The form and content of this abstract are approved. I recommend its publication.
Approved: Alexander Engau
m


DEDICATION
This thesis is dedicated to my family; Jacinda, who has stood by me and supported my
efforts all through school, and our children, Dru Shaley and Tate Wesley. I wouldnt
be the man I am without them.
IV


ACKNOWLEDGMENT
This thesis would not have been possible without the guidance of Dr. Alexander
Engau. He has shown me the way through research and publication and been an
excellent teacher. I would also like to acknowledge Michael C. Maguire for supporting
my pursuit and helping me balance work and school, which have been essential to
completing this thesis.
v


TABLE OF CONTENTS
Tables.............................................................. viii
Figures ............................................................ ix
Chapter
I. INTRODUCTION......................................................... 1
1.1 General Background.............................................. 1
1.2 Liquid Inventory Management..................................... 2
1.3 Decision Process................................................ 4
1.4 Major Contributions and Findings................................ 6
1.5 Thesis Organization............................................. 7
II. LITERATURE REVIEW.................................................... 8
2.1 Modern Approaches to Inventory and Transportation Optimization . 8
2.2 Relationship Between Linear Programming and Network Flows ... 10
2.3 Transportation Problem......................................... 11
2.4 Generalized Network Flows...................................... 11
2.5 Generalized Assignment Problem................................. 13
2.6 Dispatch Scheduling Model...................................... 13
III. STOCHASTIC INVENTORY MANAGEMENT.................................... 19
3.1 Stochastic Model Description................................... 19
3.2 Solution Methods............................................... 27
3.3 Simulation..................................................... 29
3.4 Computational Results.......................................... 30
IV. BATCHED TRANSPORTATION MANAGEMENT........................... 38
4.1 Model Formulation ............................................. 38
4.2 Underlying Network Model....................................... 39
4.3 Complexity of BTM.............................................. 42
4.4 Approximization Method Using Network Flow Relaxation........... 42
vi


4.5 Goal Programming Model.......................................... 49
4.6 Solution Methods................................................ 49
4.7 Computational Results........................................... 53
V. CONCLUSION .......................................................... 60
References............................................................. 62
Appendix
A. Data Generation...................................................... 65
B. Method Execution..................................................... 68
C. General Methods ..................................................... 71
C.l Instance Generator ............................................. 71
C.2 Model Data ...................................................... 72
C.3 Method Execution................................................. 73
C. 4 Data Visualization............................................. 79
D. STIM Solutions....................................................... 82
D. l Solver......................................................... 82
D. 2 Simulation..................................................... 86
E. BTM Solvers.......................................................... 92
E. l Direct Solver.................................................. 92
E.2 Goal Programming Solver.......................................... 94
E.3 Network Flow Relaxation Solver................................... 98
vii


TABLES
Table
2.1 DSM: Index, Parameter and Variable Notation ............................. 14
3.1 STIM: Index, Parameter and Variable Notation............................. 20
4.1 BTM: Index, Parameter and Variable Notation.............................. 40
4.2 BTM: Greedy Heuristic Example 1.......................................... 51
4.3 BTM: Greedy Heuristic Example 2.......................................... 52
A.l Data Generation: Base Data for Generation of Batteries................. 66
A.2 Data Generation: Production Rate Interval Probabilities................ 66
A. 3 Data Generation: Base Data for Generation of Haulers................... 67
B. l Data Generation: Structure for the Data Generated in Every Instance . 70
viii


FIGURES
Figure
3.1 STIM: MIP Gap vs. Optimal LPR Objective ............................... 31
3.2 STIM: Wall Time Histogram ............................................. 32
3.3 MPC: Wall Time Histogram............................................... 32
3.4 MPC: Ending Inventory Simulation Sample................................ 34
3.5 MPC: Ending Inventory Simulation Scatter Plot.......................... 35
3.6 MPC: Haul Volume Simulation Scatter Plot............................... 35
3.7 MPC: SI Volume Simulation Sample....................................... 36
3.8 MPC: SI Volume Simulation Scatter Plot ................................ 36
4.1 BTM: Feasible Region Diagram........................................... 45
4.2 BTM: MIP Gap vs. Optimal LPR Objective................................. 55
4.3 BTM: Wall Time Histogram............................................... 56
4.4 BTM: Wall Time vs. Instance Size....................................... 56
4.5 GBTM: MIP Gap vs. Optimal LPR Objective............................. 57
4.6 GBTM: Wall Time Histogram........................................... 57
4.7 GBTM: Wall Time vs. Instance Size................................... 58
4.8 NFR: MIP Gap vs. Optimal LPR Objective................................. 58
4.9 NFR: Wall Time Histogram............................................... 59
4.10 NFR: Wall Time vs. Instance Size ...................................... 59
IX


CHAPTER I
INTRODUCTION
In petroleum exploration and production, production operations are the activities
that convert the capital expenditure of exploring and drilling for hydrocarbons into
a revenue stream for business. The common reference for a company that explores,
drills for and produces hydrocarbons is E&P company or also operator. This type of
company is in the upstream sector of the industry. Many individual activities fall into
the category of production operations of an E&P company. The activity of delivering
products to sales or other disposition falls under production operations. The liquid
products in many producing fields are delivered via trucked transportation.
1.1 General Background
The term field is used in multiple ways in the petroleum sector. The most specific
usage refers to the surface of the earth covering a resource-rich geologic formation,
typically many square miles in area. A geologic formation is accessed and fluids
extracted by drilling wells from the surface down to the formation. Oil wells and
gas wells both produce liquids, typically referring to the liquid hydrocarbons as con-
densate [6] when produced from a natural gas well. Wells typically produce into a
common tank battery, where products are separated by phase and saleable liquids.
These facilities local to the wells contain equipment that helps to prepare the products
for sale and disposition. The tanks in a battery provide one method of measurement
for custody transfer, which is the trigger for converting the product into revenue.
Production tanks in batteries are the main inventory vessels for oil and water in
upstream operations.
The logistics involved with selling crude oil in the absence of pipeline infrastruc-
ture involve the local inventory at tank batteries and trucking capacity. The use of
these resources can contribute significantly to the bottom line of operations. This
1


is a scenario not seen by all companies; in established liquid-rich helds, there is in-
frastructure for product transportation, namely pipelines. In developing helds, the
infrastructure does not exist at the outset and is created only if it is economically fea-
sible. Without permanent infrastructure, resource management becomes an involved
process that, if not done right, can generate waste and reduce proht.
There are many aspects of transportation in general that apply to the transporta-
tion of liquids produced in production operations. Many times, a delivery happens
between a single tank battery and market or other common delivery point. The single
truck delivery is referred to as a load. Each load hauled costs the operator an amount
determined by contract. These fees can vary by load size, distance and product. We
refer to markets and other delivery points as destinations in general. Destinations
such as disposal wells typically exist at different locations around a held, which are
used for safe disposal of waste water from production operations [6]. In specihc sce-
narios, product is picked up from two or more tank batteries before delivering the
liquids. This is referred to as a split load. Hauling from multiple locations in one
truck load is common with water, but it is most common to haul from a single lo-
cation when transporting oil and condensate. A common situation is having only a
single market for oil and condensate, dependent upon existing contracts or location.
In more complex situations, multiple destinations can receive the product from a par-
ticular tank battery. This is also true with water. In these situations where products
are dedicated to specihc destinations, the number of decisions is reduced. This is a
simplihed view of relationships between crude oil customers and their delivery points,
but will serve our purposes.
1.2 Liquid Inventory Management
Liquid inventory is a challenge for multiple reasons. As mentioned before, water
inventory diminishes the effective capacity for oil. Crude oil in inventory is cash
2


that has not been realized. The environmental challenges of oil sitting in tanks are
becoming more prevalent. Dead oil releases hydrocarbons in gaseous form that must
be injected into a gas line or burned to reduce pollution. Modern equipment is built
to handle these emissions and maximize safety of personnel performing transfers. The
combination of technology and proper procedures ensures that the transfers are safe
and sound. Within the overarching drive of integrity, maximizing hauling efficiency is
paramount. Inventory should be minimized overall, but too many tanks with partial
loads can prevent continued efficient hauling. Another property that can reduce the
options for hauling is the oil quality. Certain markets handle only a limited range of
gravity and have requirements for non-oil content below a specific level.
In addition to the procedural complications, the volume of liquid hauled in a
single load is limited by the capacity of the truck, but the capacity is additionally
restricted by the weight of the load. The units of volume used in the U.S. is barrels
(BBL), which is equal to 42 US gallons. The specific gravity of the liquid hauled,
denoted Gs below, plays a significant role in the volume that can be hauled in a
single load. The water produced with crude oil and natural gas contains minerals like
salt and is heavier than fresh water, which has a specific gravity of 1. Crude oil and
condensate are almost always lighter than water. In the U.S., the common practice is
to refer to the American Petroleum Institute (API) gravity of the substance, denoted
by Ga here, which has a relationship to specific gravity:
Ga = (141.5/Gs) 131.5,
which gives fresh water an API gravity of 10 [1]. Also note, that API gravity and
specific gravity have an inverse relationship, so that heavier oil has lower API gravity.
If condensate is as light as Ga = 120, then a single load could carry as much as 1.78
times as much condensate as water. A truck that has a volume capacity of 200 BBL
(8400 gallons) and a 50,000 lb cargo weight limit can carry approximately 200 BBL
3


of condensate lighter than
141.5/(50000/8400/8.328) 131.5 104.18 degrees API.
Note that fresh water weighs approximately 8.328 LB per US gallon. So the
maximum specihc gravity of condensate in 200 BBL is found by 50000/8400/8.328.
Then the conversion to the standard API gravity is done. The limit on carrying water
for the same trucking conditions is then
50000/8.328/42 142.95 BBL.
This shows that a truck that can ht 200 BBL of liquid has an actual capacity
that varies with the specihc gravity of that liquid and ranges from about 140 BBL to
200 BBL.
Contracts exist between operators and haulers that dehne a minimum number
of loads or barrels that will be hauled per period. Crude oil haulers typically carry
single-source loads and charge extra fees for split loads. The extra overhead involved
with measuring, verifying product quality and more complicated routing drives the
fees.
These computations are required for implementations regarding specihc products.
The data generated for use in computational results as described in Section A can be
regarded as already in common volume (net) values.
1.3 Decision Process
The decisions targeted for using data-driven methods may run into problems in im-
plementation due to complications from the order of procedures. For instance, crude
oil that must be hauled away via truck sits in tanks that is composed of water as well.
The water will settle to the bottom, but it must be drained to a specihc level before
the crude is regarded as saleable. The quality of the product has two indicators,
4


BS&W being the water and sediment, and American Petroleum Institute (API) grav-
ity being related to the specific gravity of the product. All volumes of crude in the
USA are converted to API net barrels using the measured volume, temperature, API
gravity and BS&W. This can change the effect of optimization in terms of contract
agreements if gross volumes are the only volumes considered. It can make operations
less efficient if this consideration is not regarded and crude oil from producing loca-
tions with low API gravity is taken to markets that do not accept heavy oil. The
product is subsequently rejected and the load must be taken to a different destina-
tion. It is recommended to perform calculations regarding revenue and some capacity
metrics using API net barrels. This will account for the weight and dollar value of
the product, but cant be the only consideration in handling a mix of products that
includes extremely light condensate.
These decisions are made in a variety of different ways. In small operations,
the lease operator directly contacts the hauler to request hauls. A central planning
or dispatching department is used for managing larger operations. These planners
use data at hand to decide which requests can be delayed if necessary and which
to distribute to haulers immediately. The primary motivation of the planner is to
keep inventory from building across the held as a whole, while preventing individual
wells from having to be shut in due to reaching tank capacity. The decisions of
the planner for managing inventory through hauling has two focal points where we
will concentrate. The first is in planning. Better planning of resources required for
hauling liquids will result in more efficient expenditure to arrange the resources. This
allows the prevention of wasted capacity and resource shortages. The second focus is
on how to effectively use transportation resources for minimal cost. There are some
truck routing decisions that come into consideration, even though the operator does
not usually manage truck fleets. One example is when a single location has multiple
loads to be delivered during one period. In some fleet operations, a single driver will
5


be assigned to haul all loads from a specific location. This can lead to less opportunity
for improved cost reduction, but knowledge of this situation allows the operator to
more accurately estimate travel required for loads.
1.4 Major Contributions and Findings
Two models are explored in Chapters III and IV. These two models are Stochastic
Inventory Management (STIM) and Batched Transportation Management (BTM).
Each model has different goals. The methods used for each model are tailored for
those goals, but they all use similar base data. The results desired are focused on the
usefulness of the solutions found by the distinct solution methods and the relative
performance between the methods.
The goal for STIM is around the ending inventory, total production and total
haul. The ending inventory metric is meaningless without knowing how much to-
tal production results from a plan. The relationship between the three metrics is
straightforward. BTM has an objective value that relates directly to a metric that
can show performance of a particular solution. Performance in terms of that objective
will be discussed in detail. The other two experiments will show solution performance
in relation to other solutions.
The outcomes of STIM and BTM are not directly comparable, but as they were
both born from the same master model, there are relationships that will be exposed
in data. STIM will be demonstrated primarily as a long-term planning tool with
smaller models that have tactical application. The demonstration of BTM will be as
a tactical tool for making decisions in real-time.
6


1.5 Thesis Organization
In the three chapters that follow, mathematical approaches to the decisions described
in this chapter are discussed. Chapter II covers relevant work on transportation and
inventory using mathematical optimization methods. The next two chapters cover
new work on these topics. The new work seeks to incorporate stochastic aspects
into an existing model and to formulate a precise mathematical programming model
for a specific transportation decision. Chapter III develops and tests a stochastic
model for liquid inventory management. This model is called STIM. Chapter IV
demonstrates a model for making transportation decisions and manipulates the model
to find satisficing solutions under time limits. These models are BTM, Goal BTM
(GBTM), Fixed BTM (FBTM) and Network Flow Relaxation (NFR). The two models
STIM and BTM address specific decisions around inventory management. STIM
targets the first focal point in the decision process. The goal of this model is to bring
the decision on when to haul from a specific battery to a level considering more than
just what is locally good. BTM takes on the second focal point. Each planned load
needs a destination. The particular product being transported drives the decision.
Crude oil needs to be delivered so as to satisfy market contracts and water disposal
is a limited resource. BTM can be applied in either case, but the focus is minimizing
cost, whether that cost is dollars, distance or time.
7


CHAPTER II
LITERATURE REVIEW
There are many publications relating to the mathematical optimization of trans-
portation and logistics problems. In these works, optimization can refer to cost min-
imization or throughput maximization, among others. These approaches could be
useful to find optimal solutions to the problem we are facing. Each could be benefi-
cial in the realm of tactical decision making or short and long term planning decisions.
We will see models in this thesis mostly in the space of NP, specifically NP-
complete problems. Literature under review mainly seeks to define computationally
efficient algorithms to approximate these problems. The Transportation Problem
in Section 2.3 is known to be in P, while the Generalized Assignment Problem in
Section 2.5 is iVP-complete. For a quick reminder, these classes of problems, P and
NP, are related. The problem space P is a subset of NP. Problems in NP have the
property that a proposed optimal solution can be verified as optimal in polynomial
time. Problems in P also have the property that an algorithm exists which can fold
an optimal solution in polynomial time [29].
2.1 Modern Approaches to Inventory and Transportation
Optimization
Christopherson and Dave explore the current transportation options for crude oil in
the Great Lakes Basin [20]. In their review, trucked transportation is noted as having
the highest incident risk per unit of product transported. They also observe the need
to continue this mode of transportation. With that consideration, reducing the total
number of miles required to deliver product via truck reduces environmental and
safety risks. The driving distance for transportation also has economic impact, which
benefits operators and haulers directly.
8


The scheme of optimizing inventory in regards to customer demands has been
discussed and is called the inventory routing problem (IRP) [16]. A study specific
to crude oil transportation from Shen et al. [32] seeks to meet demands by planning
transportation options and solves the MIP problem instances using Lagrangian re-
laxation. This problem is focused on meeting the demands of customers. Since our
topic is operations focused, it doesnt seem to have direct application to our topic, but
regarding batteries as customers and destinations as supplies of capacity for deliv-
ery could change that perspective, making work regarding stochastic demand inputs
relevant.
Research on transportation and processing of liquid petroleum products exists for
producers, haulers, marketers, processors and distributors. Research around refineries
with pipeline and large-batch (ship or train) transportation modes exist, as well as
some work around trucked transportation, with trucking focused on refined products.
Neiro and Pinto define a petroleum supply chain model with particular infrastructure
sub models and mass balance constraints to assist in the planning, strategic and tacti-
cal decisions around delivery of petroleum products from the perspective of chemical
engineering [28]. Cheng and Duran describe an approach to crude oil logistics that
uses simulation and optimal control [18]. They use a rolling horizon strategy to seek
a suboptimal control strategy with lower computational requirements. Abduljabbar
and Tahar build on this work with a model [34] and case study [7] to show the benefits
of those tools.
Some works around inventory optimization [9, 22, 30] use genetic algorithms to
drive decisions on inventory management in supply chains without being product
specific.
9


2.2 Relationship Between Linear Programming and Network
Flows
Network flows have the useful property that efficient methods exist for finding solu-
tions that are guaranteed to be integral as long as inputs are integral [23]. Efficient
methods for solving general LPs also exist. These methods can also guarantee inte-
gral solutions, but only when input constraint matrices are totally unimodular [8].
These two separate classes of problems are connected in that all network flows can
be represented as linear programs, but not all linear programs can be represented as
network flows. The conversion of linear programs to network flows is of interest.
Bixby and Cunningham have described an algorithm for converting LPs to net-
work flows [12] which makes the determination if a matrix is graphic or non-graphic.
This term graphic is used in reference to a matrix that is projectively equivalent to
a node-arc incidence matrix (NAIM). While Bixby and Cunningham discuss their
algorithm in matroid theory terms, Baston, Ramouni and Williams describe an algo-
rithm [10] in matrix terms for those that are less familiar with matroid theory. They
discuss tests that can be executed to prove a matrix is graphic, but explain that any
algorithm based on this approach has worse than exponential complexity. Despite
this complexity, they describe another test to demonstrate when a given matrix is
not convertible to a NAIM. This is not useful in this context because as will be seen
in Chapter IV, the conversion to a network representation is required, even if the
original problem is not maintained. For a problem with r rows and c columns in a
constraint matrix A, the algorithm proposed by Bixby and Cunningham has com-
plexity 0(rn) with n being the number of nonzero entries of A. This is efficient and
could be used to expose the underlying network structure of an LP, if not perform
the full conversion.
10


2.3 Transportation Problem
The transportation problem involves a commodity or commodities, depots and mar-
kets. Each depot has a certain supply and each market has a specified demand. An
answer to a transportation problem is one that satisfies demands from given depots
for a minimal cost. Let x E Rraxm where Xij is a quantity meeting demand at cus-
tomer i by supplier j. The formulation of the problem also involves the quantity of
demand at each consumer di and quantity of goods at each supplier Sj. A cost of
supplying goods from supplier j to customer k is Cjk- The formulation of the LP is
minimize > x / i,j (2.1a)
subject to < Sj i Vj (2.1b)
^ ^ a'xy di j V i (2.1c)
> 0 Vq j- (2. Id)
A special case of the transportation problem is the assignment problem. This
problem has a balanced number of supply and demand nodes where each node has
demand di = 1.
2.4 Generalized Network Flows
The transportation problem has demonstrated utility in industry, but is incapable
of providing solutions to related problems where the network flows have loss or gain
from one node to the next. Consider a network where flow leaving a node across
an arc does not equal the flow entering the receiving node via that arc. The flow
in and out of a node still balances as required in traditional solution methods, but
flow changing across arcs contradicts the basis of classical theory in network flows.
11


Alternative solution methods are required. A minimum cost flow problem over nodes
i E Af and arcs (i,j) E A with term Kij referred to as the gain of arc (i,j) formulated
in [11] as
minimize > CijXij X / (2.2a)
subject to KjiXji 1 M o Vz e Af (2.2b)

lij ij 0/ V(z, j) e A. (2.2c)
Equation (2.2b) is a flow conservation constraint alternative to the classical flow con-
straint. The factor indicates the flow along arc (i,j) is changed, so arriving flow at
j is different than leaving flow from i. Despite this difference from so-called ordinary
networks, duality still provides a guarantee of optimality through complementary
slackness conditions. The algorithms proposed in [11] use Lagrangian relaxation to
define a dual ascent method that iteratively updates prices from the dual problem
while maintaining flows that satisfy complementary slackness seeking primal feasi-
bility. Differences between these algorithms for ordinary and gain networks involve
the additional steps required in flow augmentation along arcs with non-unit gains
and price updates for nodes connecting these arcs. An additional step for flow aug-
mentation along a cycle could be omitted in graphs with special structure such as
bipartite graphs seen in assignment problems. These algorithms are demonstrated in
[11] as having performance improvement over classical methods for multiple classes
of problems such as assignment, transportation and transshipment problems.
Other methods exist for specific problems in the space of generalized network
flows, like one facility one commodity (OFOC). OFOC is a network flow problem
with single commodity flow between a single source and sink where flows must be
made in increments of C or 1 with costs applied in to the batch. Chopra et al.
[19] describe a polynomial algorithm for solving instances where the number of full
12


batch flows k of size C is bounded. This approach demonstrates the advantages of
concentrated research on the qualities of special instances of difficult problems.
2.5 Generalized Assignment Problem
In the generalized assignment problem (GAP), agents with specified capacity are
assigned tasks requiring general levels of service in a way to minimize the cost of
assignments. This problem is NP-complete [27], unlike the assignment problem from
network flow optimization. Starting with the same notation from the transportation
problem, GAP has the further restrictions Xij E {0,1}, the flow constraint RHS is
1 for all i and the assignment by each has a weight uy to be counted against
each agents capacity Sj. The formulation is nearly as simple as the transportation
problem;
minimize > (kjXij i,j (2.3a)
subject to WiXij < Sj i Vj (2.3b)
Xii ~ ^ j V i (2.3c)
E {0,1} Vq j- (2.3d)
Solving, or approximating, instances of GAP can include rounding methods [33]
and GPU implementations of related problems [17, 15, 24] suggest that parallelization
methods of class single instruction, multiple data (SIMD) can be used to speed up
the search for solutions.
2.6 Dispatch Scheduling Model
Engau, Moffat and Dyk present a multiperiod, multicriteria model for planning load
dispatch schedules in [21]. This model is referred to here as the dispatch scheduling
13


Table 2.1: DSM: Index, Parameter and Variable Notation
Sets Description Indices
A Set of sets of batteries which can be combined into a single load A
B Set of inventory locations or tank batteries i
C Set of haulers 3
D Set of destinations k
T Time period index set t
Parameters Units
C-ij k Revenue per barrel hauled from i G B to k G D by j G C dollars / bbl
C-k Cost per barrel for capacity constraint violations by k G D dollars / bbl
dij k Split load fee per load hauled from i G B to k G D by j G C dollars / load
hj k Full load size by hauler j and destination k bbl
h Minimum contractual haul from battery i G B (typically zero) bbl
h Minimum contractual haul using hauler j G C bbl
lk Minimum contractual supply to destination k G D bbl
Qi(t) Production rate for battery i G B during time t G T bbl / period
VjI Maximum storage capacity of battery i G B bbl
uj(t) Maximum haul capacity by j G C during t G T bbl
Uk {t) Maximum capacity at k G D during t G T bbl
Wi Parameter for battery i G B for weighting shutin avoidance unitless
V? Initial liquid inventory at battery i G B bbl
6 Shut in volume weight for multicriteria objective unitless
7 Hauled location weight for multicriteria objective unitless
Variables Domain
Vi(t) Inventory in barrels for battery i G B at the end of period t G T R+
ijk (t) Barrels of fluid hauled from ieBtofcePbyjeC during i e T R+
Djk{t) Loads hauled by hauler j G C to destination k G D during i e T
%ij k{p) Loads hauled from i G B to k G D by j G C during i e T
n{t) Time to shut in by battery i G B for t G T R+
ilk Capacity constraint overflow variable for k G D R+
14


model (DSM). It covers the decisions on when to haul loads from batteries Zijk(t),
when to split loads between batteries at cost dijk, which hauler to assign to each
load and the destination for each load. The model decisions are driven by multiple
objectives including avoidance of shutins due to exhausting tank capacity, meeting
contractual obligations, profit maximization and haul volume maximization.
To understand DSM, we begin with the concept of haul volume. The model rep-
resents the volume of liquid each battery i E B sends to destination k E D via hauler
j E C during period t E T. This volume is denoted Xijkit). The volume is trans-
ported in batches sized hjk BBL, which corresponds to the hauler and destination.
The model assumes that all loads arriving at destination k are at full capacity. The
full capacity load is enforced with constraint (2.4f), which uses the number of loads
Ujkit) sent to destination k via hauler j during period t. DSM uses the variable Zijkit)
to indicate the number of loads transported from battery to destination by a particu-
lar hauler. These loads are not necessarily full, but constraint (2.4g) ensures that the
volume does not exceed the capacity corresponding to the decision of z^kit). Along
with these variables and constraints that ensure hauling is represented consistently,
the mechanics of liquid flowing into and out of the battery must be consistent. The
inventory on site at the beginning of the planning period is a required parameter
in BBL incorporated with constraint (2.4e). The inventory at the end of each pe-
riod Vi(t) is controlled by equation (2.4d). This equation demonstrates consistency
with the liquid flowing into the battery q^t) and fluid hauled away from the battery
^2jec ^2keD xijk(t)- Each hauler j also has limits on the liquid hauled during a time
period. The constraint (2.4c) ensures that the total liquid hauled by a hauler meets
two requirements, the lower limit lj and upper limit Uj. The lower limit corresponds
to a promised amount or load count that exists in contracts or agreements, while
the upper limit corresponds to a capacity limit or other agreement. The destinations
have similar limits enforced by constraint (2.4b) with Ik and Uk being parameters
15


for lower and upper limits on fluid amount by destination. The constraints (2.4h)
ensure that the variables for load count are non-negative. The model objective (2.4a)
represents the revenue corresponding to a given hauling schedule. The parameter djk
is the revenue received for liquids delivered from battery i to destination k via hauler
j. The revenue received is reduced by the cost of the number of split loads required.
Each split load costs dijk dollars. This objective is only part of the multicriteria ob-
jectives explored. See the initial model formulation below, with notation summarized
in Table 2.1.
Mze EEEEf (c^
x,y,z,v
ieB jec keD teT
Hjk
h
'jk
Xijkit) dijkZijk(t)
subject to Ij < EE Xijk (t) < Uj (t) V j e C and teT;
ieB keD
*sEE Xijk(t) < Uk(t) Vk E D and teT;
ieB jec
'sEE Xijk{t) = qi(t) Vi(t) + Vi(t 1) Vi e B, t eT;
jec keD
(2.4a)
(2.4b)
(2.4c)
(2.4d)
0 < Vi{t) < Ui and ?y(0) = Vi e B and teT; (2.4e)
^2xijk(t) = hjkUjkit) V Ae A, i,j,k,t; (2.4f)
ieA
0 < Xijkit) < hjkZijkit) Vi, j, k, t (2.4g)
Ujkit) > 0 and Zijkit) > 0 integer. (2.4h)
Note the revenue by barrel djk term is used in the proht calculation, which can be
interchanged with cost and negated to represent water hauls, which do not generate
revenue. The system balance equation (2.4d) used in this model is common to all
calculations for change in inventory. It ensures that the ending inventory from one
period at a battery is the beginning inventory of the next period for the battery.
Model features such as split load decisions and shut in avoidance using time to sell
out are key to the design. Split loads can be controlled using the criteria (2.4f) where
all loads are assumed to be full and A is a set of batteries that can be used together
16


to build split loads. This constraint guarantees a total number of loads is honored
and the maximal number of splits can be controlled by setting a relationship between
J2ieB zijk and Ujk where Zijk is the number of loads hauled from battery i to destination
k by hauler j. Normally, since crude oil loads are split between at most two batteries,
this would be J2ieB z^k < 2but can be restricted even further to limit splits to a
fraction of total hauls to a destination by a hauler. This is not necessary, but could
be used to indicate early when desired constraints are infeasible because the model
objective guarantees 2^ takes on the smallest value where Xijk < hjkZijk is feasible
[21]. The other feature unique to this model is the method used to avoid shutins.
Using a calculation for the time to shut in by battery
Tiit)
Uj Vjit)
Qi(t)
(2.5)
with a term in the objective
(2-6)
i£B t£T
which uses a weight factor ity by battery on the time until shut in -ry (t), the model
can be made to avoid shutins. To incorporate this, (2.4a) is replaced with
Maximize s > > > >
x, y, z, v I J J J J
\ ieB jec keD teT
C-ijk
Hjk
hi
jk
Kijk(t) dijkZijk{t)
(2.7a)
££££ Xijk(t), and ££wi(i) (2.7b)
i£B j£C k£D teT i£B teT )
to seek the additional objectives beyond profit maximization. Maximizing the sum
of these variables r will result in batteries with low rates of production being favored
for hauling since rq is designed to meet at least one haul and approaches 0 as wells age. A more realistic approach to the use of this variable r in
planning is to seek a maxi-min of Tiit). This ensures that high producers maintain
a level of inventory capacity that is comfortable to operations. One drawback on
the use of this model is the deterministic approach to production operations. This
17


formulation also does not allow for any shutins. Situations where production volume
exceeds hauling capacity by enough to force the selection of wells to shut in will fail
as infeasible in this formulation. This is recognized in [21], but not addressed in a
model reformulation.
18


CHAPTER III
STOCHASTIC INVENTORY MANAGEMENT
In this chapter we describe a model that builds on the previous work from DSM
[21], which is a model for planning dispatch schedules to manage inventory for pre-
venting shut-ins, keeping sales consistent and maximizing revenue. That model does
not incorporate any stochastic components, and that work is continued in this thesis
as stochastic inventory management (STIM). Uncertainty in production and hauling
calls for a model that takes a stochastic approach to its inputs. Using the technique
of stochastic programming explored in [31], we will develop an alternative model to
handle uncertainty.
3.1 Stochastic Model Description
We will consider a model that incorporates random processes. The entropy in trans-
portation of goods arises naturally from sources such as traffic, equipment and facility
exceptions. The unknown of hydrocarbon extraction rates is one that is specific to this
industry. This model seeks to maintain efficiency by eliminating production down-
time and maximizing throughput while honoring capacity constraints. The differences
between STIM and the dispatch scheduling model presented in [21] are intended to
make computation of large instances more accessible, while maintaining the goal of
use for decisions in the near term with regard for long term planning. The notation
for STIM here builds on Table 2.1 and adds some stochastic terms given in Table 3.1.
3.1.1 Stochastic Aspects
The process of producing and selling crude oil and hauling water has stochastic com-
ponents, such as the rate of production, the time from plan to when a haul is actually
performed and the hauling capacity. Hauling capacity as a stochastic feature refers to
the number of loads a particular hauler can perform in a period from the perspective
19


Table 3.1: STIM: Index, Parameter and Variable Notation
Set Description Indices
B Set of inventory locations or tank batteries i
C Set of haulers j
D Set of destinations k
T Time period index set t
M Scenarios m
Parameters Units
rp\n') Probability of scenario m for battery i unitless
q(fn') (t) Production volume at battery i for scenario m during time t bbl
Ui Maximum storage capacity of battery i G B bbl
Uj (t) Maximum haul capacity by j G C during t G T bbl
Uk(t) Maximum capacity at k G D during t G T bbl
S Shut in volume weight for multicriteria objective unitless
7 Hauled location weight for multicriteria objective unitless
Variables Domain
g\m\t) Shut in production by location i during time t, scenario m R+
v\m\t) Inventory at battery i following time t in scenario m R+
xtjk W Haul volume from battery i during time t in scenario m R+
Zi(t) Indicator for hauling from location i during time t {ON}
20


of the operator. The focus of this model is on the random nature of the production
rate for each battery. Production decline curves from wells are modeled using deter-
ministic decline curves, but operations are subject to exceptions which are stochastic
in nature, interrupting the natural decline curve and causing unexpected drops in
production volume.
This choice of feature will represent other components implicitly. The impact
of a haul occurring later than anticipated has a similar effect to higher production
volume than anticipated. With a discrete time model, the delay would need to be
implemented by having a haul occur during a period following the planned haul.
While this is not uncommon in practice, it is limiting in that without shortening
periods and increasing the number of decision variables, only large deviations can
be represented. Having multiple possible production scenarios can help to represent
delays and early arrivals.
3.1.2 Model Development
The development of a stochastic model requires the definition of scenarios representing
different outcomes. The scenarios for outcomes are denoted m E M. Each battery
i and scenario m has an associated probability value p[m\ It is required that the
vector pi is stochastic, meaning = 1, or the expected value calculations
are incorrect. With the model choice for incorporating random events being the
production rate the inventory at the end of each period must then
be a dependent stochastic variable. Due to the inventory being stochastic, the haul
quantity x^{t) should be varied between the scenarios m E M. This is done by
adding constraints with the different scenarios of production rates, ending inventory
and haul quantity. Note the hauling decision variable Zi(t) is still deterministic. If this
hauling decision were stochastic, the model would be separable on scenario, giving
optimal outputs for each scenario, instead of a decision that seeks optimal expected
21


value considering all scenarios. Look to Table 3.1 for details on this notation. It is
important to note that even though production rates for any two distinct batteries
in a particular scenario m may have the same probability value, there is no other
relationship implied between the two in the same scenario. The stochastic variable
Qiit) for battery i is independent of the variable qj(t) for battery j. Each battery i in
each scenario has probability pf1'1, which allows for flexibility in instance specification.
Hauling capacities Uj are deterministic, possibly unique to the hauler. Inventory
capacities tq are fixed for short term planning purposes, but could be made variable
for long term decision making, from the operators point of view.
Planning hauls using a stochastic model is challenging because the dynamic na-
ture of the inventory means that a full load might be available in a particular period
in one scenario, but not in a different scenario, i.e. v[m\t) < hjk for some m G M. It
is similar in attempting to plan multiple loads. The number of loads available could
have a large range within just a few different scenarios. Consider a large producer
with a mean of 2000 BBL to be hauled per period. Under typical conditions, this
volume may only vary 10% in either direction from period to period, but this means
that the number of loads that needs to be planned to haul the liquid produced in
a period is 9, 10 or 11 (in simple 200 BBL loads). It becomes even more complex
when considering planning periods that may see new wells turn on. The time that a
new well is turned on may be known with high probability a short time in advance,
but the new production rate is more difficult to predict. Also, issues may delay the
turn on that would greatly disrupt a plan that has allocated many new resources to
hauling from a new location. The uncertainty in available inventory forces a change
to one assumption that is used in [21], that is all loads to be hauled are full. This
assumption makes possible much of the direction of the model design. The handling
of split loads follows from this assumption. Doing away with this assumption elim-
inates some structure from the model, but does not take from its ability to return
22


useful results.
One of the situations not addressed with the dispatch scheduling model is infea-
sibility of the model for unavoidable shutins. To address this we add a new term
g^it). This represents production volume that exceeds capacity of local tanks so
must be shut in. It is particular to each scenario, battery and period. Positive values
indicate that wells at a particular battery i will need to be shut in to reduce total pro-
duction volume in time t to the fraction of total production (<^"^(i) ~ 9im\t))/q[m\t)
in that scenario. This means that E^2ieBgi(t) production volume can be expected
to be lost during period t for forced shutins.
The terms discussed here are summarized in the Table 3.1. The controlling equa-
tion that drives the dynamic nature of the problem introduced in (2.4e) is changed to
accommodate stochastic elements and prevent infeasibility for cases requiring shutin
production to
TW = vT\t 1) + - y>'b(() j,k
Two more constraint sets are required to keep the decision and haul volume
variables consistent. One is required to ensure that no location has positive haul
volume if that location is not designated for hauling. The constraint
EE t) Mziit) <0 \/i e B, t eT, m e M (3.2)
jec keD
ensures this situation using a big-M method in M.Zi{t) where A4 is a value sufficiently
large enough to be greater than any other value calculated in the problem. It is only
satisfied with Zi(t) = 0 when x^(t) is also zero. So positive haul volumes require
that the particular battery decision is positive. The other requirement for maintaining
consistency is for volumes by hauler and destination. The constraint for hauler volume
by period is
xiS (t) < Uj(t) Vj G C,t eT.me M. (3.3)
ieB k£D
23


The effect of this constraint is to keep total hauled volume below the volume capacity
for a particular hauler. As well as the hauler capacity constraints, the destination
constraints are similar:
EE (t) < uk{t) \/k E D,t E T,m E M. (3.4)
i£B j£C
Each inventory variable is constrained to be non-negative and less than the upper
bound of capacity by battery tq.
The objective function for the stochastic model is to maximize the expected
value of the sum of hauls from a minimum number of locations while maintaining the
production volume. The expected haul volume is
E EEEE E4ha
i£B m£M t£T jeC keD
Production volume is maintained by subtracting ^2ieB ^9i'* W frc>m the
sum of hauled volumes, which penalizes the objective using a factor 5 for all positive
shut in volumes The number of hauled locations ^2ieB^2teT Zi(t) is implicit
to the number of split loads required to haul the total volume. This value penalizes
the objective by a factor 7. And finally, low ending inventory is a goal of operations,
so the objective is penalized with the expected value of inventory at the end of the
planning horizon. The full objective formula is
maximize ^
i£B \t£T \ m£M j£C k£D m£M )
-EErM (36)
m£M /
where A is the final period ending inventory by battery i.
In the dispatch scheduling model, the objective naturally drove the load count
to a minimal number. The objective in STIM requires a term explicitly added for
directing load count to lower values. This is the same as analyzing the cardinality of
24


haul volume vector x(t), due to (3.2). The term 7Zi(t) is meant to keep the number
of locations hauled to a minimum, considering the maximization of haul volume. The
parameter 7 can be adjusted to change the number of locations chosen for hauling.
The reasoning behind the term is to trade off the minimum load size with the choice
of hauling. In the dispatch scheduling model, the term Zijkit) is in the domain of
integers. The restriction of z^kit) to a binary variable in STIM is meant to be a small
simplification step for computational complexity.
The full model is presented here.
maximize
v,g,x,z
7 Zi(t)+ Y YY.P'
m£M jeC keD
M
(m)
Xijk
(t) Y6p'
(m) (m)
in
(*)
m£M
- £ P?)v'r)
m£M /
subject to ££ Yjk(P) uk(t) Vfc G D,t eT,m E M (3.7a)
i£B j£C
Y xi7k (*) ^ u7t) Vj e C,t eT,me M (3.7b)
i£B k£D
££ xtjk(f) ~ Mziit) < 0 V* G B, t G T, m e M (3.7c)
j£C k£D
Y\t) = Y\t 1) + qt\t) Yxm(t) gt\t) v* g 5, me M
j,k
(3.7d)
0 < v\m\t) < UiV i G B, m G M (3.7e)
vY = Ui(max{t G T}) \/i G B, m G M (3.7f)
> 0, Zi{t) G {0,1},
g(m) (^) > 0 yieBjj eC,keD,teT,meM
Also note the term g^it). This decreases the objective when decisions force
production shut in. The parameter 8 is used to indicate the importance of preventing
shut in of production volumes. This formulation does not solve the same problem as
DSM, which considers the trade off between cost of hauling split loads and inventory
25


management. This stochastic model seeks to maximize hauling and will result in a
solution that needs to be interpreted to find split loads and possibly formulated in
another model to decide how to haul the splits. This is an operation that is required in
the dispatch scheduling model as well. These subsequent problems are not necessarily
feasible, given upper bounds on load counts. Suppose that either of these models gives
a solution that includes full loads for all but hauls from four batteries, with the sum
of barrels to be hauled equal to the total number of barrels in the maximum number
of loads by hauler (which is implied in the models, but not explicitly parameterized).
Then no additional loads would be available. The haul volumes from these four
batteries are 70,80,90,160 BBL, respectively. The total of these loads is 400 BBL,
which is two full loads for hjk = 200. The four sources cant be combined in whole
to make two full loads. This is a scenario that is a feasible result from the dispatch
scheduling model and STIM. The approach proposed in [21] to solve a bin-packing
problem following a solution in this first stage would deliver an additional solution
that has a minimum number of loads required, if optimal. Although this number is
not necessarily less than or equal to the implied number of loads implied by Uj/hjk,
the uncertainty in hauling capacity mentioned earlier in this section means exceeding
the load count could be inconsequential.
3.1.3 Model Options
Features from the dispatch scheduling model [21] can be incorporated into this model
as needed. The decision on which batteries to haul can be replaced by a variable
denoting the exact number of loads to haul between locations and destinations. The
maximization of time to shut-in or maximization of time to sell-out could be brought
into the objective to prioritize different goals in inventory management. Using the
basic approach of modeling the effect of particular decisions under different scenarios
is not modified by adding multiple criteria or changing to a more complex decision.
26


3.2 Solution Methods
For large cases, STIM is difficult to solve. For 30 periods and 500 batteries, this results
in 15,000 binary variables. Applying a clever MILP solver may trim the search tree,
but the search space is so large that it will take drastic cuts to make the search
successful. Reducing the number of variables is one way to reduce the search space
directly. Reducing the number of decision variables in an instance is done here using
the method of model predictive control (MPC), discussed in Section 3.2.3.
3.2.1 MILP Solver
The use of a MILP solver on this model will still give feedback, even in large instances.
The solver (in general) will give at least a MIP gap and understanding of the size of
the tree remaining to search after a specified run time. There is the possibility of not
finding a solution in that run time, which would be an infinite MIP gap, but it will at
least show how much of the search tree is examined. It is not likely that the search
tree will have been trimmed by a significant proportion for general large instances,
but the information will still be available from a branch and bound or branch and
cut global optimization method.
3.2.2 Solution Heuristic
In DSM, the exact number of loads per battery is determined by a combination of
integer decision variables. In STIM, the decision variables are in the batteries to
haul from and quantity of product to haul per battery, hauler and destination. This
change makes STIM a cardinality problem. The minimization of J2ieBJ2teT''fZi(t)
in the objective in turn reduces the cardinality of x, due to the constraint (3.7c)
forcing Xijk(t) to 0 when Zi(t) = 0. An heuristic exists for cardinality problems
that involves repeated solutions of a relaxed version of the original problem. This
27


relaxation allows the binary variable to take any value in [0,1], subject to other
existing constraints. This is essentially a rounding heuristic on the decision to haul,
restricting the cardinality of x. It is seen in [13]. To begin this heuristic, a set
of cardinality costs must be selected. These costs are different values to apply in
the model penalizing the cardinality of x. Here is the step-by-step execution of the
heuristic.
1. Set initial objective to oo.
2. Choose a cutoff a E [0,1).
3. Build the model instance allowing fractional solutions.
4. Solve this original instance with cardinality cost 7, FAIL if infeasible.
5. Add constraints to the instance Zi(t) = 0 for each variable Zi(t) < a and Zi(t) = 1
for each Zi(t) > a.
6. Solve the new instance for the remaining variables v, g and x.
7. If objective is higher than the current objective value, then save objective and
solution.
8. Repeat Step 3-7 with next 7 from a set of cardinality costs.
This method works by fixing the cardinality at each iteration after the relaxed
solution is found. Each iterations begins with the original problem, but a new cost
for the cardinality of the solution. This means that the method does not necessarily
increase or decrease in objective value and the lowest objective should be recorded
along the way and returned as the solution after the final iteration is complete. This
algorithm is not a global search technique. It is different in at least the solution set is
not trimmed from one iteration to another; at the start of each iteration, the solution
28


set is the same. A binary integer programming technique would include branching to
search the solution space.
3.2.3 Model Predictive Control
To solve this model over a sufficiently long planning horizon, we use the method
of model predictive control (MPC) by applying a series of shifted time periods. The
MPC technique can be applied in conjunction with the solution heuristic from Section
3.2.2. This technique sets a sliding time horizon at each iteration, moving forward
one period and using the results of the execution of the previous decisions. This
technique can be applied in a planning task using simulation to determine the results
of the decisions or in a live environment, getting feedback from the system.
In a finite time horizon formulation, one choice is to set a constraint forcing the
state of the system to a specific value [14]. This has the effect of preventing the
solution in the current horizon from putting the model in a predicament in the next
horizon. Applying a cost in the objective instead of a hard constraint is an alternative
method for seeking the same goal. This is applied in the model formulation by
decreasing the objective by the expected value of the total ending inventory. Even
without MPC, this is a desired goal of the model.
3.3 Simulation
Simulation is used with STIM to drive MPC and to validate results. With MPC,
the sliding time horizon calls for inputs resulting from the previous run, which are
the outcomes of executing decisions. One method used to generate these values is
simulation. This type of simulation executes decisions and comes up with mean values
for parameters needed. Simulation can also help validate results by applying decisions
generated by the model and comparing the resulting values, such as haul volume or
29


shut in volume for STIM, to the results prescribed by the optimal solution.
As in the design of the model, the desired stochastic components for simula-
tion are the production rates for the batteries. A simulation would be more easily
designed to have the individual load pickup times be stochastic, a random variable
of time following the planned haul to when a load is picked up. Or it could be a
random variable for the number of loads picked up in a period for discrete time sim-
ulations. Alternatively, the process of simulating the inventory dynamics does not
include any representation of the trucks, tanks or even individual loads. The entire
battery capacity is treated as a single tank.
3.4 Computational Results
The data for STIM experimentation is not modified much from the generic data
discussed in Appendix A. As the model is designed to generate haul decisions, the
inventory, production rate and constraint data only needs to be translated into the
format required by the model. The model instances will feature additional variables
used to track objective, haul volume and lost production aggregates. These values
were be used to gauge and demonstrate model performance.
3.4.1 Solutions to STIM
Two solution methods have been applied to STIM.
The first is the pure MILP solver. This method was executed primarily to show
running times. The MIP gap plot in Figure 3.1 reflects the maximization goal. All
points are negative, showing a wide range of MIP gaps. Some of the instances below
5% gap may be decent performance, considering the planning application. As can
be seen in Figure 3.2, no instances ran under the cutoff time. Only a fraction of the
instances completed due to an explosion in time seen with overhead taken for building
30


2000
0
-2000
. -4000
IS
(3
Q.
s -6000
-8000
-10000
-12000
-10000 0 10000 20000 30000 40000 50000
Optimal LPR Objective
Figure 3.1: STIM: MIP Gap vs. Optimal LPR Objective
the data structures required for solving these instances.
The second method uses the heuristic algorithm and MPC. In this method, the
plot of the MIP gap is not given. The execution time of this method is quicker
than the pure implementation, shown in Figure 3.3. Since all instances of STIM
were terminated before reaching optimality conditions, MPC is used exclusively for
comparison with simulated results.
3.4.2 Simulation
Simulation applied to STIM checks the validity of the solution found with the discrete
scenarios in a model with continuous random variables. The simulation process is as
follows. The process is run for each instance and each instance generates 1000 samples.
1. Capture the output of STIM for an instance.
2. Generate base conditions, including beginning inventory ?y(0) and mean pro-
31

i 4 >
i

.




50
40
i 30h
o
U
CD
u
c
to
£ 20
10
0L
185
190
195 200
Wall Clock Time
205
210
Figure 3.2: STIM: Wall Time Histogram
13 14 15 16 17 18 19 20
Wall Clock Time
Figure 3.3: MPC: Wall Time Histogram
32


duction rate qi for each battery i.
3. For each period t, do the following
Get a new rate of production for each battery using q.\ generated from
the normal distribution for each battery i as explained in further detail in
Appendix A.
Vi(t 1) + Qi(t) are the available volume to haul.
Apply the hauling decision to the available volume and haul ay(f) up to
the least of the planned amount and available volume.
If Vi(t) = Vi(t 1) + qiit) Xi(t) is greater than the inventory cap, the
difference is the shut in volume.
Repeat for the prescribed number of periods.
This process generates a number of samples each of haul volume per period, shut
in volume per period, dry loads per period and the ending inventory for the entire
planning horizon. These data have been compared with planned output.
3.4.3 Results
The execution of the MPC method are shown here. A sample of the ending inventory
compared between the simulated and planned results is shown in Figure 3.4. This
plot is only 10% of the cases, but it shows how the simulated results closely matched
the planned results. The scatter plot in Figure 3.5 shows the simulated vs planned
ending inventory which shows nearly identical results. In the haul volume scatter
plot (Figure 3.6), each periods simulated haul volume is plotted against the planned
haul volume. The r-squared value of this data is 0.99997. Also in Figures 3.7 and 3.8,
notice a different pattern. The plot in Figure 3.8 includes the minimum, maximum
and mean of the simulated shut-in volumes by case are plotted against the planned
33


Figure 3.4: MPC: Ending Inventory Simulation Sample
shut-in volumes. Although the haul volume and ending inventory results track closely
between plan and simulation, the simulated shut in volumes are less consistent to plan
and are also less than the plan, demonstrated in the two plots.
3.4.4 Analysis of Results
The results have shown while there is still work to be done to analyze scenarios with
higher variance in the production rate, the plans and simulations have agreed so
well that it can validate the correctness of the implementations. The accuracy of
the simulated results to planned results is not to be taken as a forecast of how the
model will do in actual business use. It needs to be shown that the model will still
generate good plans in systems with more entropy. The results found by simulation of
model decisions are close to optimization results for haul volume and ending inventory
model, but less so for shut in volume. This shows that the model can be expected
34


Figure 3.5: MPC: Ending Inventory Simulation Scatter Plot
Figure 3.6: MPC: Haul Volume Simulation Scatter Plot
35


Figure 3.7: MPC: SI Volume Simulation Sample
Figure 3.8: MPC: SI Volume Simulation Scatter Plot
36


to perform well when the variance applied in planning is the same variance used in
simulation or accurate to the real-world scenario, but still requires some research on
shut in volumes.
One simulation decision that could be rethought is the decision to limit haul
volume by the planned volume. The result is simulated (expected) haul volume
is slightly lower than planned haul volume because simulated haul volume may be
reduced by experiencing a reduction in available volume to haul when production
volume is lower. A simple alternative to this approach is to allow haul volume in
the simulation to meet available volume. This would not be acceptable without
dealing with how loads will be picked up from each battery in the simulation. When
loads are not full, then there is the possibility of a partial load from one location
being topped off with the volume from another location. In a planning scenario,
this inconsistency is certainly small enough to still provide good information. In an
operational scenario, the models are perhaps good enough for business use. More
testing is required. This testing could be performed in a real world situation where
STIM presents a solution and prediction. Then legacy decision processes would be
used to carry out operations. Using simulated results from both decision sets, analysis
would be performed to compare how results would be different if the plan had been
used. This is one way to to build confidence in the model.
37


CHAPTER IV
BATCHED TRANSPORTATION MANAGEMENT
The model Batched Transportation Management (BTM) is a model that repre-
sents decisions relating to the transportation of liquids from batteries to destinations.
This decision is just one part of DSM, assuming that other decisions which lead up
to the decision on where to haul liquids have been made. Those other decisions can
be made by applying STIM or some other business process. This chapter presents
BTM, which is a model that seeks optimal decisions for delivering liquid volumes
to destinations considering cost (dollars, time, distance, etc). There are several key
elements to the model driven by the needs of decision makers. The objective of the
model is to honor destination capacities and minimize cost of hauling, which in a ref-
erence implementation is the travel distance between source inventory location and
destination. This translates to time spent on the activity and minimization of this
cost leads to efficiency gains. In many cases, the incentive to increase efficiency and
reduce costs is secondary to the ability for a method to provide a satisficing answer
[26]. Goal programming [25] is a tool from multicriteria optimization that seeks to
fold satisficing solutions to difficult problems. This modeling method assigns priori-
ties to criteria in a problem using cost coefficients and minimizes the cost of violating
these criteria. This pertains to the model in this chapter through the objective of
minimizing costs and the requirements found in the capacity constraints for haul
destinations that prove to be a challenge to meet.
4.1 Model Formulation
The mathematical model representing the decision on how much to haul to each des-
tination is given below. See Table 4.1 for the definitions of elements in the model. We
will formulate this model without time period variables. This is because the trans-
portation decision is not impacted by battery inventory, which is the complicating
38


aspect of DSM and STIM leading to multiple period analysis. This model seeks the
minimum total cost (4.1a) where Cijk is a cost for making a certain hauling assign-
ment Zijk- Each assignment Zijk results in a quantity Xijk (by constraint (4.1c)) that
represents the use of capacity whose upper bound is Uk in constraint (4.1b). The
model (4. la-4. Id) is called batched transportation management (BTM).
minimize cijkzijk (4.1a)
ieB jec keD
subject to x^k < Uk Vk e D (4.1b)
ieB jec
^ ^ y ^ %ijk V i G B (4.1c)
jec keD
%ijk$i %ijk Vi e B,j e C,k e D (4.Id)
zijk integer Vi e B,j e C,k e D (4. le)
Dropping the integrality constraint and only requiring positive values for produces
the LP relaxation of BTM, which we will refer to as LPR.
4.2 Underlying Network Model
BTM is a linear programming formulation with integer flow balance constraints that
cannot be represented using network flows, except in special cases. The necessity of
two decision variables for different aspects of a haul comes from having supply nodes
that do not fill capacity of batched deliveries alone, i.e. partial loads. BTM does
not attempt to fill individual deliveries to capacities, which is left to decision makers
(specifically the hauler) following destination assignment.
BTM can be represented as a generalized assignment problem. This problem is
known to be NP-complete. It is NP-complete to even find a feasible solution [27].
Since GAP is a generalized network flow model, this points to the underlying nature
of BTM as a network flow problem. The assignment problem has tasks of the same
39


Table 4.1: BTM: Index, Parameter and Variable Notation
Sets Description Indices
B Set of inventory locations or tank batteries %
C Set of haulers 3
D Set of destinations k
T Time period index set t
Parameters Units
Cijk Cost per load hauled from i G B to k G D by j G C dollars / load
C-k Cost per barrel for capacity constraint violations by k G D dollars / bbl
lk Minimum contractual supply to destination k G D bbl
n Total hauls required from i G B truck load
Si Volume of fluid per haul from i G B bbl
Uk Maximum capacity at k G D during i e T bbl
Variables Domain
ijk Barrels of fluid hauled from i G B to k G D by j G C R+
%ijk Loads hauled from i G B to k G D by j G C
m Capacity constraint overflow variable for k G D R+
40


weight and different costs, per assignment decision. Although the assignment prob-
lem can be solved using classical network flow methods, the generalized assignment
problem cannot. The manipulations required to formulate BTM as a GAP instance
reveal the network structure.
It can be seen in the formulation of BTM that the variable Xijk is tied directly to
Zijk- The variable Xijk must be a multiple of s*. Replacing (4.1b) with
EE SiZijk ieB jC
and (4.1c) with
Y Y zvk = A V* G B (4.3)
jec keD
allows the elimination of (4.Id). Now this structure is closer to a graphic repre-
sentation liking a transportation problem, but (4.2) still has non-graphic aspects.
Arbitrary values of s* mean BTM is a generalized network flow and preclude the use
of classic min cost flow methods to solve the problem. However, s* is bounded from
above. Let this upper bound be s. Suppose that Si = s for all batteries i E B, then
(4.2) can be rewritten
EE Zijk ieB jC
without changing the problem. The model (4.1a, 4.3, 4.4) is equivalent to a trans-
portation problem and can be solved using min cost flow methods. Call this model
FBTM. Consider FBTM with one hauler. Then (2.1a) is equivalent to (4.5a), (4.5b)
is equivalent to (2.1b) and (4.5c) is equivalent to (2.1c). The network changes nec-
essary to include multiple haulers are to add nodes for each hauler and destination
combination with arcs of cost Cijk between the battery and new node. The original
destination nodes sum the flow across arcs coming from the new nodes. This network
flow model is slightly more complex than a transportation problem (a transhipment
problem), but still has a polynomial algorithm. The full formulation of this new
41


model
minimize ^ ^ ^ cijkzijk ieB jec ken (4.5a)
^ ^ ^ ^ Zijle ^fc/S ieB jec VkeD (4.5b)
z^k = ri jec keD \/i G B (4.5c)
zijk integer Vi e B,j e C,k e D (4.5d)
only solves the original BTM problem if the load sizes all equal s.
4.3 Complexity of BTM
BTM is an iVP-complete problem. This can be seen by demonstrating how GAP
is a special case of BTM. Consider BTM with all batteries having a single load to
deliver, i.e. rj = 1 for all batteries i E B, and only a single hauler. Now it can be
seen that (4.3) is equivalent to (2.3c). The constraint (4.2) needs no modification
to be equivalent to (2.3b). With the assumption of load count and single hauler,
the objective (4.1a) is equivalent to (2.3a). Under these assumptions, the objective
and constraints between the two problems are equivalent, so GAP is a special case of
BTM. GAP is iVP-complete, therefore BTM is also iVP-complete.
4.4 Approximization Method Using Network Flow Relax-
ation
In this section we will explore an approximation method that has increasing accuracy
as load sizes s* become closer to s. Recall that when Si = s for all batteries i E B,
BTM is equivalent to FBTM and can be solved with network flow methods. The
objective of this section is to create a model that can be solved using the same methods
as FBTM, but does not eliminate solutions from the feasible region of BTM.
42


4.4.1 Model Instance Manipulation
The demonstration that specific instances of BTM are solvable using min cost flow
methods suggests that manipulations of other BTM instances could be solved in the
same way. These solutions are approximations and manipulations that are too simple
could create infeasible approximations where the original problem is feasible. For
instance, take the simple case with 5 loads, one at 150 BBL and four at 100 BBL
with total delivery capacity of 600 BBL. Assuming that every load is 150 BBL would
allow the use of min cost flow methods, but the total capacity would be less than the
total delivery requirement of 750 BBL. The manipulation in (4.6) and the resulting
new model in (4.7) creates a new problem that expands the feasibility region in BTM
and allows for the implementation of a solver with network flow methods.
Since it is assumed that not all loads to be hauled are the same size, it makes
sense to define the difference between the largest load size s and load size by battery
Si for all i E B. Let 5i = s Si. Using this value, lets build a new constraint
representing capacity limits by destination. Begin with (4.2)
iB jec
and take the inequality
which is to represent the balance between loads delivered and loads required. Then
sum these inequalities together resulting in
Since s = Si + Si\/iEB
43


Finally, dividing both sides by s produces a new constraint for our approximation
model
£ £ zijk < Uk + S6ii S,r< Vi D. (4.6)
i£B j£C
Let us call this new model network flow relaxation (NFR) after the solution
approach that is being targeted. We now have the entire formulation of NFR:
minimize ^ ^ ^ cijkzijk i£B j£C k£D (4.7a)
ZiF = ri j£C k£D Vi G B (4.7b)
V"' V"1 ^ Uk + Si£B ^ri / v / v Z%F £ i£B j£C VkeD (4.7c)
zijk integer Vi e B,j e C,k e D (4.7d)
NFR is a transportation model like FBTM, so NFR is solvable in polynomial time.
Since the capacity constraint for NFR is relaxed in a specific way, NFR has the
property that if a solution is feasible in BTM, it is also feasible in NFR. See Figure
4.1 for the set relationship for feasible solutions to NFR, BTM and the LP relaxation
of BTM (LPR). NFR does not allow fractional solutions and neither does BTM. LPR
has fractional and integer solutions. The intersection of the sets of feasible solutions
for LPR and NFR is the set of feasible solutions for BTM. This is due to LPR having
no integer solutions that are not also solutions to BTM.
In addition to the relationship between LPR and NFR, there are results for NFR
that structure its use in searching for solutions to BTM. Refer to (4.1) and (4.7) for
details regarding these proofs.
Theorem 4.4.1 Given a problem instance BTM and any feasible solution z for that
instance, z is a feasible solution for NFR.
44


Figure 4.1: BTM: Feasible Region Diagram
Proof: Let z be a feasible solution for BTM. Take the capacity constraint (4.2)
^ ^ &izijk G bT G D
ieB
and the manipulations producing (4.7c) show that z is still feasible regarding the
capacity constraints for each k G D. Therefore z is feasible for NFR.
Since s is an upper bound for Si V* G B, any assignment of loads to destinations
that honors capacity constraints uu for each k in BTM will also honor all capacity
constraints in NFR. All other constraints remain unchanged, so z is feasible for NFR.
Theorem 4.4.1 allows an exploration of the entire solution space of BTM using network
flow methods. Although, it is not true that feasible solutions (even optimal ones) to
NFR are feasible in BTM. An easy way to know when certain instances are not feasible
is using LPR. The following simple theorem demonstrates this easy check.
Theorem 4.4.2 Given a problem instance BTM, LPR and a related instance NFR,
if there exists a solution z* for NFR and J2ieB ^2keD cijkztjk ^ess than the
optimal objective for LPR, then z* is infeasible in BTM.
Proof: By contradiction. Let z* be an optimal solution to NFR for some BTM
45


instance. Let z\ be an optimal solution for LPR. Suppose that cTz* < cTz\. Suppose
also that z* is feasible for BTM. Since all z feasible for BTM are feasible for LPR, z\
must not be optimal, which is a contradiction. Thus, either cTz* > cTz\ or z* is not
feasible for BTM. Therefore, for cTz* < cTzi, z* is infeasible for BTM.
A solution to NFR that has an objective value greater than the objective value for
LPR is not necessarily feasible for BTM, so more calculations are required to check
solutions not covered by Theorem 4.4.2. But, an optimal solution with objective
value greater than the optimal objective value for LPR creates a new lower bound for
the objective value of BTM. The next theorem demonstrates the traditional use of
relaxation methods to create lower bounds for an original problem is valid for NFR.
Theorem 4.4.3 Given a feasible problem instance BTM, the related problem NFR
and an optimal solution z* for NFR, the objective value cTz* is a lower bound for the
objective cTz of BTM.
Proof: By contradiction. Suppose that there is a solution z\ to BTM with
an objective value cTzi < cTz*. By Theorem 4.4.1, z\ is feasible for NFR. Thus
cTzi is a feasible objective value for NFR, but cTz* is the optimal value of NFR.
This means cTz* < cTz for any solution z to NFR. In particular cTz* < cTz\, so
cTz* < cTz\ < cTz*, which is a contradiction. Therefore the objective value cTz*,
optimal for NFR, is a lower bound for the objective value of BTM.
FBTM is not necessarily feasible, even if BTM is feasible. Perhaps it is feasible,
with solution z*. Then z* is feasible for BTM. In terms of feasible regions, FBTM
lies entirely inside the feasible region of BTM, empty or not. The following theorem
states this explicitly.
Theorem 4.4.4 Given a problem instance BTM, the related problem FBTM and an
46


optimal point z* for FBTM. Then z* is feasible for BTM.
Proof: Let z* be an optimal solution to FBTM. Then, from the derivation
producing (4.6) we have,
sz*jk i£B j£C i£B j£C
which is the capacity constraint for BTM. Therefore, an optimal solution for FBTM
is a feasible point for BTM.
This is not only true for an optimal solution to FBTM, but any solution to FBTM.
An optimal solution to FBTM is nearly as easy to find as a feasible one, thanks to the
network flow algorithms that can solve it, so it is more interesting than any feasible
solution. An optimal solution to FBTM is not necessarily optimal in BTM, but is
useful. It may be that a solution to FBTM is good enough for a given instance or
may be a starting point for a decision maker to choose to violate some constraints to
come up with a better objective value.
4.4.2 Better Bounds
Theorem 4.4.1 shows that there is a simple manipulation of the BTM model to create
a network flow model that contains all feasible solutions for BTM. The manipulation
of the constraints is not very involved and very much a liberal manipulation. There
is a more conservative approach to manipulating the capacity constraints that still
allows for all BTM solutions to be feasible in NFR.
This method is more involved and requires a single sort, followed by m binary
search operations of size n to calculate the relaxation of capacity at each destina-
tion. The complexity of the manipulation is 0(\D\log |£>|) since binary search has
47


complexity 0(log |£>|). The method follows five steps
1. Sort the list of load sizes.
2. Choose the next destination k.
3. Use binary search on the loads to find the largest number 7 of loads that can
feasibly be assigned to the destination (call this set Gk).
4. Add ^2ieG Si to uk for the new upper limit for destination k, assigning to uk.
5. If destinations remain, repeat steps 2-5.
This algorithm is guaranteed to create a new problem that has all the feasible solutions
to the original problem. This is due to the condition that any solution z for BTM
will have a set of loads for some k E D, call it H, where for each i e H, for some
j e C zijk = 1, J2ieH Sjec Zijk < \Gk\ Vfc E D. Since the \Gk\ smallest load sizes for
i E B will have J2ieH Si < J2ieGk Si and can then feasibly fit at least as many loads
assigned to a particular destination.
Consider a simple example to demonstrate the algorithm, consider a set of loads of
sizes s = [70, 90,110], <5 = [40, 20, 0] and a destination k of size uk = 200. The binary
search technique will pick load sizes [70,90] and S = [40, 20]. Then the relaxed volume
is 200 + 60 = 260 BBL. Note any all three cannot be assigned to this destination since
70 + 90 + 110 = 270 > 200, two loads combined is less than the original capacity of
200 BBL and 2s = 220, so all assignments feasible before are feasible following the
capacity relaxation and fixing all load sizes to s.
48


4.5 Goal Programming Model
BTM can be easily extended to a goal programming model. The objective (4.1a) al-
ready seeks the minimization of cost. The constraints (4.1c) and (4.1b) can be relaxed
and brought into the objective for a goal programming approach. The flow conser-
vation constraint (4.1c) could be changed from equality to inequality and penalized
for not delivering the amount of flow r^, but it is assumed that this requirement is
already firm since it likely comes from solving a model such as STIM, which considers
production shutins. Relaxing the constraint (4.1b) by adding an overflow term rjk
to the right hand side allows for solutions to exceed the capacity constraints Uk at a
price of ck per BBL in decreased objective value, both terms available for reference
in Table 4.1. The resulting model is called GBTM and is formulated
minimize ^ ^ ^ cijkzijk + ^ ckqk i£B j£C k£D k£D (4.9a)
subject to ^ ^ xijk ^ ^ ^ ^ Xijk GG j£C k£D V i e B (4.9c)
Zijk$i Xijk Vi e B,j e C,k e D (4.9d)
zijk integer Vi e B,j e C,k e D (4.9e)
Vk> o Vk e D. (4.9f)
4.6 Solution Methods
While a problem does not have to be NP-hard to pose difficulty with finding solutions,
the nature of BTM is such that a solution may not be found without an exhaustive
search, if one exists. That being said, the solution methods for linear programs,
49


the simplex method and interior point methods can be applied and could result in
feasible, optimal solutions in the integer space. When they do not, other approaches
must be used. The methods explored are listed in this section.
4.6.1 LP Formulation
When solving general MILPs, an initial step frequently is to solve the LP relaxation.
This is no different for BTM, where the decision variables Zijk are fractional instead
of integer. This formulation may result in a solution that is in the integer domain,
is quick and easy to execute and provides an absolute lower bound on the objective
function using the dual problem. Even if the solution found is not entirely integer,
employing a rounding method may return a solution that is close to feasible. The level
of violation allowed should be left to the decision maker, which is just referring to the
violation of capacity constraints. The violation of flow should not be tolerated. This is
where a simple rounding method may not be usable. For instance, if a solution found
suggests sending less than half a particular load to three different destinations, simple
rounding would result in that load being assigned to no destination. It is possible
that each positive value in this case is the same, which calls or a more sophisticated
method to decide which of the positive values should be 1.
4.6.2 Simple Greedy Heuristic
The simple greedy heuristic uses the following algorithm. The token NEXT is an
unassigned load chosen by some mechanism that iterates over all loads. The token
QUIT is for when the algorithm successfully terminates, having found a solution.
The token FAIL is for when the algorithm terminates having failed to find a feasible
50


Table 4.2: BTM: Greedy Heuristic Example 1
Load Assigned Destination
100 200
100 200
120 220
140 250
150 None
solution.
1. If list of choices for NEXT is empty, FAIL.
2. Rank each choice for NEXT.
3. Assign the lowest cost choice to NEXT.
4. If no loads are unassigned, then QUIT.
5. Otherwise, move NEXT to an unassigned load.
This method provides no guarantee of optimality or even Ending a feasible solution if
one exists, so it is heuristic in nature. Consider this example to gain insight into the
heuristic. Let s = [100,100,120,140,150] and u = [200,220,250] for the destination
capacities. Suppose that the choice of ordering gives [100,100,120,140,150] and that
each destination costs the same between the loads, with the order of the destinations
being the preference for each load. So, the assignments are given in Table 4.2 and
the heuristic fails. But if the ordering of the loads is [140,100,120,100,150], the
assignments are given in Table 4.3 which is a successful execution of the heuristic.
51


Table 4.3: BTM: Greedy Heuristic Example 2
Load Assigned Destination
140 200
100 220
120 220
100 250
150 250
The generality of the mechanism determining NEXT allows for different rules to
be applied to find different solutions. Some examples are given.
Producer priority: the unassigned location with the highest volume to deliver
is given priority to be chosen at each iteration.
Max distance priority: the unassigned location with the furthest possible dis-
tance to travel is given priority at each iteration.
Min distance priority: the unassigned location with the nearest possible distance
to travel is given priority at each iteration.
Different combinations of these ordering mechanisms can be used to change the order
in which loads are assigned using the greedy heuristic. The producer priority ordering
seeks to give the highest producers the best choice possible. The max distance priority
ordering seeks to eliminate the worst possible trips. The min distance priority seeks
to make the decisions that seem the best for each location.
52


4.6.3 Network Flow Solution with Full Loads
Theorem 4.4.4 calls for the modification of the original problem to change all load
sizes to be full loads. This can be a wild approximation when the loads in a particular
instance are dominated by partial loads, with some full loads to be hauled. Even so,
when a feasible solution exists for the FBTM problem, it may be satisfactory for daily
use. The efficient solution methods that exist for network flow problems make it an
attractive option.
4.7 Computational Results
4.7.1 Instance Creation
To create data for use in BTM, the data from each instance described in Appendix A
needs to be translated into vectors for loads Xijk, costs and destination capacities
ttfc. The process used does not seek to create feasible instances. The load selection
method simply looks for all loads ready in inventory. The goal of experimentation
with BTM is to see how different methods behave with all types of problems. The
infeasible problems are of as much interest as the feasible ones.
The load generation is of the most interest here. Each problem instance has initial
inventory ?y(0). The smallest load size lj from the set of haulers is used to decide how
many loads were available to haul from each battery by dividing the initial inventory
by the load size and taking the floor of the result. If the locations with zero loads
following that procedure has enough inventory for 1/4 of a load, then they are selected
as well. This is to represent multiple pickups and ensure different load sizes.
53


As in other experiments, the cost for each decision was created using the G norm
of the difference of the two-vectors indicating location on the plane. The capacities
for the destinations were taken directly from the generated data.
4.7.2 Solutions to BTM
Data for each instance of BTM was generated using a single script and appropri-
ate solver. Each instance of BTM was solved using the pure MIP model, GBTM
and NFR implementation. These results of these implementations were compared to
demonstrate advantages of each method.
The plots in Figures 4.2, 4.3 and 4.4 show the performance of the MIP solver
on pure BTM. The results show that most instances completed execution in short
amounts of time. Figure 4.2 shows the MIP gap vs LPR objective. The line is the
tolerance supplied to the solver for the MIP gap to determine optimality. The points
above the line are instances that did not meet the tolerance and timed out. In Figure
4.4, the plus markers show instances deemed infeasible. Notice how some medium
sized instances reached optimality in times between the low levels and the maximum
of 180 s. There are no large instances that were executed between less than a few
seconds and the time limit. These results reinforce the difficulty of the problem.
The next set of plots, Figures 4.5, 4.6 and 4.7 show patterns that are similar,
but reveal that the goal programming approach has more of an issue with growing
complexity and execution time. The comparison is not directly comparable because
a major fraction of the large instances solved with BTM were infeasible. In GBTM,
no problem is infeasible. Note that the MIP gap used is the gap between the MIP
objective and LPR objective only in terms of the decision variables. The portion of
54


160
140
120
100
g- 80
to
0
o.
2 60
40
20
0
~20 0 2000 4000 6000 8000
Optimal LPR Objective
Figure 4.2: BTM: MIP Gap vs. Optimal LPR Objective
the objective that drove the constraint violation cost is left out for comparison to the
pure BTM results.
Finally, the NFR plots in Figures 4.8, 4.9 and 4.10 demonstrate how the new
approach can help in finding solutions in polynomial time. Notice that the infeasible
solutions still show up, but no instances took longer than 30 seconds to arrive at
optimality or infeasibility.







**<*>* -t t~
mm-1---*- ------ ---
55


300.
2501
2001
o
1501
1001
50l
50
100
Wall Clock Time
150
200
Figure 4.3: BTM: Wall Time Histogram
Figure 4.4: BTM: Wall Time vs. Instance Size
56


250
200
150
100
50
-3
000

*
t
A # :
i V i
i.
2000 4000 6000
Optimal LPR Objective
8000
10000
Figure 4.5: GBTM: MIP Gap vs. Optimal LPR Objective
Wall Clock Time
Figure 4.6: GBTM: Wall Time Plistogram
57


200
150-
100-
50
v
_5Ql-------1
-100 0
100 200 300 400 500 600
Integer Variables
Figure 4.7: GBTM: Wall Time vs. Instance Size
Figure 4.8: NFR: MIP Gap vs. Optimal LPR Objective
58


300.
2501
2001
o
1501
1001
50l
5 10 15 20
Wall Clock Time
25
30
Figure 4.9: NFR: Wall Time Histogram


t
*
t
a 41 1 ifefA*

-5i---------1---------1---------1---------1--------1---------1---------1
-100 0 100 200 300 400 500 600
Integer Variables
Figure 4.10: NFR: Wall Time vs. Instance Size
59


CHAPTER V
CONCLUSION
This thesis has demonstrated two different mathematical tools and their appli-
cation in management of liquids inventory and hauling. Stochastic inventory man-
agement (STIM) is a planning tool that helps to how resources should be arranged
and what production profile to expect. Batched transportation management (BTM)
is useful for tactical decisions for transportation and resource usage. Both problems
are NP-hard. In large implementations, heuristic and approximation methods are
required to deliver solutions that seek optimality.
STIM has been created to find hauling, production and inventory patterns seeking
optimal conditions for minimizing shut in volumes, maximizing hauling and minimiz-
ing ending inventory. It is useful for planning purposes, but also can be implemented
to find tactical decisions for hauling. The options available for this model allow it to
be used to seek different goals. The simplifications used in this model have allowed
the use of an efficient algorithm to arrive at solutions. The heuristic method applied
to STIM has been shown via simulation to provide guidance on how to manage in-
ventory through hauling. Even though the method supplied fractional solutions, the
rounding method or polishing heuristic resulted in integral decisions that seek to min-
imize the locations visited while maximizing haul volume. The patterns that result
are as expected with higher splits to haul more volume. The distribution used for
generating experimentation data has given a highly accurate planning to simulation
prediction for haul and ending inventory volumes as demonstrated by our computa-
tional results. These are key indicators for decision makers on the performance of oil
and water hauling. This demonstrates that the heuristic method can be useful for
business application of STIM.
The running times of pure BTM show a tendency to exceed time limits required of
decision makers for problem sizes equivalent to medium to large operations. GBTM
60


shows that introducing feasibility on otherwise infeasible problems only provides an
answer, while optimality is still unanswered for many instances in the time allotted.
NFR offers a solution method bounded by a reasonable polynomial. The solutions
returned by this method are not necessarily useful for decision makers outright due
to allowing capacity constraint violations, but may lead to further solution methods.
The implementation of BTM and related models GBTM and NFR has shown that
each has its application for specific cases. The pure IP implementation of BTM can be
efficient for finding results of feasible problems and frequently delivers those solutions
quickly with a method combined with branch and cut and heuristics. GBTM is not
as efficient, but does provide an answer, even if not directly feasible. The infeasibility
can be tuned to reduce capacity violations by desired destination. Feasible solutions
are delivered quickly and can be found by a direct calculation in polynomial time.
Optimal solutions have not been demonstrated for large instances due to the solver
reaching time limit. NFR offers a unique mix of the two using a network flow method
on a modified instance. It can provide approximate solutions in polynomial time.
The use of these two models in conjunction with each other will provide decision
support for tactical and planning situations. STIM is useful for planning in short and
long term and when combined with BTM for tactical decisions on resource usage, the
result is having knowledge of how to get the most out of them in the current situation.
Follow up work in this held may include better efficiencies gained through decisions
made closer to real time and an additional layer of planning in truck routing.
61


REFERENCES
[1] The engineering toolbox: Api gravity. The Engineering Toolbox, http://www.
engineeringtoolbox. com/api-gravity-d_1212. html. Accessed: 21-MAR-
2015.
[2] Gnu linear programming kit. http://www.gnu.org/software/glpk/glpk.
html.
[3] Networkx: High-productivity software for complex networks. https: //
networkx.github.io/.
[4] Pulp. http://www.coin-or.org/projects/PuLP.xml.
[5] Python language reference, version 2.7. Python Software Foundation, https:
//docs.python.org/2.7/.
[6] Oilfield Glossary. Slumberger Limited, http://glossary.oilfield.slb.com/,
2015. [Online; accessed 16-March-2015].
[7] W. K. Abduljabbar and R. M. Tahar. A case study of petroleum transportation
logistics: a decision support system based on simulation and stochastic optimal
control. African Journal of Business Management, 6(11):43504361, 2012.
[8] R. K. Ahuja, T. L. Magnanti, and J. B. Orlin. Network Flows: Theory, Algo-
rithms, and Applications. Prentice Hall, 1993.
[9] F. Altiparmak, M. Gen, L. Lin, and T. Paksoy. A genetic algorithm approach for
multi-objective optimization of supply chain networks. Computers & Industrial
Engineering, 51(1): 196215, 2006.
[10] V. J. D. Baston, M. K. Rahmouni, and H. P. Williams. The practical conversion
of linear programmes to network flow models. European Journal of Operational
Research, 50(3):325-334, 1991.
[11] D. P. Bertsekas and P. Tseng. Relaxation methods for minimum cost ordinary
and generalized network flow problems. Operations Research, 36(1):93114, 1988.
[12] R. E. Bixby and W. H. Cunningham. Converting linear programs to network
problems. Mathematics of Operations Research, 5(3):pp. 321-357, 1980.
[13] S. Boyd. LI norm methods for convex cardinality problems, http://web.
Stanford. com/class/ee364b/lectures/ll_slides .pdf. Accessed: 05-MAY-
2015.
62


[14] S. Boyd. Model predictive control, http://web.stanford.com/class/ee364b/
lectures/mpc_slides.pdf. Accessed: 05-MAY-2015.
[15] L. Bukata and P. Sucha. A gpu algorithm design for resource constrained project
scheduling problem. In Proceedings of the 2013 21st Euromicro International
Conference on Parallel, Distributed, and Network-Based Processing, PDP 13,
pages 367-374, Washington, DC, USA, 2013. IEEE Computer Society.
[16] A. Campbell, L. Clarke, A. Kleywegt, and M. Savelsbergh. The inventory routing
problem. In Fleet Management and Logistics, pages 95-113. Springer, 1998.
[17] A. Chaparala, C. Novoa, and A. Qasem. A simd solution for the quadratic assign-
ment problem with gpu acceleration. In Proceedings of the 201 f Annual Confer-
ence on Extreme Science and Engineering Discouery Enuironment, XSEDE T4,
pages 1:11:8, New York, NY, USA, 2014. ACM.
[18] L. Cheng and M. A. Duran. Logistics for world-wide crude oil transportation
using discrete event simulation and optimal control. Computers & Chemical
Engineering, 28(6-7):897-911, 2004.
[19] S. Chopra, I. Gilboa, and S. T. Sastry. Source sink flows with capacity installa-
tion in batches. Discrete Applied Mathematics, 85(3): 165192, 1998.
[20] S. Christopherson and K. Dave. Cardi reports. Technical Report 15, CaRDI
Reports, Cornell University, November 2014.
[21] A. Engau, C. Moffatt, and W. Dyk. Multicriteria modeling and tradeoff analysis
for oil load dispatch and hauling operations at Noble Energy. Optimization and
Engineering, 16(1):73 101, 2015.
[22] R. Z. Farahani and M. Elahipanah. A genetic algorithm to optimize the total
cost and service level for just-in-time distribution in a supply chain. International
Journal of Production Economics, 111 (2):229243, 2008.
[23] L. R. Ford and D. R. Fulkerson. Flows in Networks. RAND Corporation, Santa
Monica, Calif., 1962.
[24] M. E. Lalami and D. El-Baz. Gpu implementation of the branch and bound
method for knapsack problems. In Proceedings of the 2012 IEEE 26th Interna-
tional Parallel and Distributed Processing Symposium Workshops & PhD Forum,
IPDPSW T2, pages 1769-1777, Washington, DC, USA, 2012. IEEE Computer
Society.
[25] S. M. Lee et al. Goal Programming for Decision Analysis. Auerbach Philadel-
phia, 1972.
[26] W. A. Lodwick and J. Kacprzyk. Fuzzy Optimization: Recent Aduances and
Applications. Studies in Fuzziness and Soft Computing. Springer, 2010.
63


[27] S. Martello and P. Toth. Knapsack Problems: Algorithms and Computer Imple-
mentations. John Wiley & Sons, Inc., New York, NY, USA, 1990.
[28] S. M. S. Neiro and J. M. Pinto. A general modeling framework for the opera-
tional planning of petroleum supply chains. Computers & Chemical Engineering,
28(6):871-896, 2004.
[29] C. H. Papadimitriou and K. Steiglitz. Combinatorial Optimization: Algorithms
and Complexity. Dover Books on Computer Science. Dover Publications, 2013.
[30] P. Radhakrishnan, V. M. Prasad, and M. R. Gopalan. Genetic algorithm based
inventory optimization analysis in supply chain management. In Aduance Com-
puting Conference, 2009. IACC 2009. IEEE International, pages 418-422. IEEE,
2009.
[31] A. Shapiro, D. Dentcheva, et al. Lectures on Stochastic Programming: Modeling
and Theory, volume 16. SIAM, 2014.
[32] Q. Shen, F. Chu, and H. Chen. A Lagrangian relaxation approach for a multi-
mode inventory routing problem with transshipment in crude oil transportation.
Computers & Chemical Engineering, 35(10):21132123, 2011.
[33] D. B. Shmoys and E. Tardos. An approximation algorithm for the generalized
assignment problem. Mathematical Programming, 62(1-3):461474, 1993.
[34] R. M. Tahar and W. K. Abduljabbar. A novel transporting system model for oil
refinery. American Journal of Engineering and Applied Sciences, 2(1): 138143,
2010.
64


APPENDIX A. Data Generation
The bulk of the data used in experimentation has been generated. The goal of
much of the work is to show the efficiency of solvers in execution time and deter-
mining feasibility. Each instance represents a single product. It can be regarded
as representing consistent value between locations. For water, revenue is zero. For
oil, revenue can change between the destination of delivery, which is a consideration
driving the development of the multicriteria model in [21]. Otherwise, costs are rep-
resentative of either product, net oil volumes are reasonably consistent for weight and
water volumes are as well. Thus this experimentation is agnostic toward the product.
The computational methods in this thesis use a set of problem instance generated
specifically for this work. These include 100 cases each of
20 batteries, 1 hauler, 2 destinations
40 batteries, 2 haulers, 3 destinations
100 batteries, 5 haulers, 6 destinations
The batteries in each case are built using random data. A battery was initially chosen
to be of one of the types in Table A.l. After randomly choosing a type, the particular
base rate of production is chosen using a normal distribution around the types rate.
The standard deviation (a) used was 5% of the chosen type rate. Using a normal
distribution around the random base rate (ja) and the same standard deviation, the
rate values corresponding to /a 2a, fi .75a, n, n + .75a and n + 2a were assigned.
These values are referred to as qlO, q35, q50, q65 and q90 in the code. The probability
assigned to each of these cases is the same across batteries. It is calculated by using
a standard normal table. The interval for each case has a width of 1/2 standard
deviation and is centered on the distance from the mean. Using the law of total
65


Table A.l: Data Generation: Base Data for Generation of Batteries
Class Production Rate Inventory Capacity Selection Frequency
Low Producer 1 BBL 300 BBL 2/7
Low Mid Producer 1 BBL 300 BBL 2/7
High Mid Producer 100 BBL 600 BBL 2/7
High Producer 1000 BBL 2400 BBL 1/7
Table A.2: Data Generation: Production Rate Interval Probabilities
Case Interval Probability Marginal Probability
qlO (j>{n 1.75a) (j>{n 2.25a) P(ql0)/P(A) = 0.05047033
q35 4>(/i 0.5a) l.Oer) P(q35)/P(A) = 0.27098408
q50 (j>{n + 0.25a) (j>{n 0.25a) P(q50)/P(A) = 0.35709117
q65 (j){fJL + l.Oer) + 0.5a) P(q65)/P(A) = 0.27098408
q90 (j>{n + 2.25a) (j>{n + 1.75(j) P(q90)/P(A) = 0.05047033
probability, the case probabilities are calculated in Table A.2. Event A is the event
that the rate falls in one of the intervals, with probability
P(A) = P(qlO) + P(q35) + P(q50) + P(q65) + P(q90).
The inventory capacity is unmodified. The initial inventory is chosen to be a
uniform random value between 0 and the inventory capacity. Each battery was also
given coordinates, (x,y), to represent the location on the square (50, 50), (50, 50).
While it is representative of some level of operation, it has not been built around the
operations of any particular E&P company.
Each hauler was generated with a minimum and maximum load count per period,
full load size, dry load fee and split fee. The dry load fee was not used during
optimization, but during simulation to show an effect of overzealous dispatch. See
66


Table A.3: Data Generation: Base Data for Generation of Haulers
Variable Value
Lower Limit on Loads Hauled Per Period lj 0 or [/[0,40] with equal probability
Upper Limit on Loads Hauler Per Period Uj lj + U[ 0,60]
Load Size Uniform choice of {180,190, 200, 250}
Dry Load Fee Uniform choice of {75,100,125,150}
Split Load Fee Uniform choice of {75,100,125,150}
Table A.3 for the detail.
The destinations were also given coordinates in the same square. Along with co-
ordinates, a minimum and maximum delivery volume was assigned using the following
relationship. Let A,B ~ U[0,5000]. Then
h = A (A.l)
Uk = lkJr 3B, (A.2)
so that each destination has a minimum delivery (lk) up to 5000 BBL and maximum
delivery(fc) between lk and lk + 15000 BBL.
This is the whole of the base data used for generating problem instances. See
Appendix C.l for the code used.
67


APPENDIX B. Method Execution
All solver-related code created for this work has been written in Python [5].
The project PuLP [4] associated with the COIN-OR initiative allows the direct ma-
nipulation of variables, constraints and objectives to translate problem instances to
structures needed by solvers. Each model solved using LP/MILP was interfaced with
GLPK [2] through PuLP. For network flow optimization, the project NetworkX [3] for
Python was used. This project has a network simplex method that solves min cost
flow problems, satisfying all demands. A single script (see C.3 in the appendix) was
used to execute each solver for a particular instance. The script requires the instance
ID, the instance repository and output destination database. Database access was
done using anydbm in the standard Python library [5] to provide random access for
reading and writing the data. Execution of the script begins with determination of
the mode of operation, which model type or solver to run. The instance data is read
from the text hies, translated into the required structure for the requested method
and fed to the solver. If the solver is of type MILP, then the search time limit is set
to 180 seconds with a stopping criterion of 0.1% MIP gap.
Following completion of the method, either by successful optimization, certificate
of infeasibility or time limit exceeded, the resulting data is written to the destination
database. The data keys included are listed in Table B.l. Several plots and tables
were generated using this data for each method to compare performance. If additional
data points are required for a particular method, then those points are discussed in
the relevant section for that method.
For each method, three plots are generated. The first is a plot of MIP Gap vs
Optimal LP Objective. This plot is intended to demonstrate the stopping criteria of
the method in terms of the absolute lower bound. Each instance that has a feasible
68


solution is plotted as a dot and the stopping criteria of 0.1% is plotted as a solid
line. Each dot below the line represent an instance which reached tolerance or proven
optimality and the dots above the line represent those instances that reached a feasible
solution, but ran above the time limit.
The second plot is a histogram of the instance count in bins of solver duration.
This plot is straightforward and demonstrates a bimodal result where large numbers
of instances fall into categories of low processing time and high processing time.
The third plot is relatively simple as well. It plots a scatter of wall clock time
versus integer variables per instance. One additional visual representation in this plot
is the distinction of feasible and infeasible instances with feasible instances plotted
with a dot and infeasible instances plotted with crosshairs. Feasible and infeasible
instances should be mixed around the low processing time range and there should
be a higher number of feasible instances in the higher processing time range. This is
due to the number of instances that are infeasible merely because the haul volume
exceeds the delivery capacity.
69


Table B.l: Data Generation: Structure for the Data Generated in Every Instance
Key Description
clocktime processor time required to find solution (Python code only)
walltime duration of time to find a solution
objective smallest feasible objective value found
lpobjective proven minimum value for LP relaxation
haulvol total haul volume (by period, either scalar or array)
volcap total delivery capacity
loadcount total loads to be delivered
destcount count of destinations
solution solution found by solver or None, if infeasible
70


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
APPENDIX C. General Methods
C.l Instance Generator
# * coding : utf8
import numpy as np
import csv
def c r e a t e i n s t a n c e ( B C, D):
^^^^Create^a^Hauling^problem^instance^with B^- ( number o f^batter ies ) ,
( number -of^haulers ) ^and ( number of^destinations)^in^a^l00-x^l00
MMMsquareMlayout.
bbox = np arr ay ( [ ( 5 0 , 50), (50, 50)])
producers = np array([1 , 1, 10, 10, 100, 100, 1000])
prodchoice = range ( le n ( p r o d u c er s ) )
inv = np. array ( [300 , 300, 300, 300, 600, 600, 2400])
fees = np.array([75, 100, 125, 150])
loadsize = np arr ay ( [ 1 8 0 , 190, 200, 250])
chooser = np. array ([0 , 1])
def get.producer () :
coords = bbox[0] + (bbox[l] bbox [0] ) * np.random.random ((2 ,))
i = np .random, choice (prodchoice)
base = producers [ i ]
cap = inv[i]
# randomize base
sigma = base / 10.
base = np random normal(base sigma)
if base < 0: base = 0
blO = base 2 sigma
b35 = base .75 * sigma
b50 = base
b65 = base + .75 * sigma
b90 = base + 2 sigma
vO = np random random ( ) * cap
return [coords[0], coords[1], blO, b35, b50, b65, b90, cap, vO]
def get.hauler () :
d r y f ee = np random choice ( fees )
splitfee = np .random, choice ( fees )
full = np .random, choice(loadsize)
1 = i n t ( np .random, choice(chooser) * 40 * np random random ( ) )
u = 1 + int(60 np random random () )
return [1, u, full, dryfee, splitfee]
def get_destination () :
coords = bbox[0] + (bbox[l] bbox[0])
np random random (( 2 ,))
71


46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
1
1 = i n t ( np .random, choice(chooser) 5000 np random random ( ) )
u = 1 + 3 int(5000 np random random () )
return [1, u, coords [0], coords [1]]
batteries = [ get.producer () for b in range(B)]
haulers = [ get_hauler () for h in range (C) ]
destinations = [ g e t d e s t i n a t i o n ( ) for d in range(D)]
return batteries haulers destinations
def write.dat a ( header rows, fn):
with f i 1 e ( f n , wb ) as f 1 :
if __name__^= __main__ :
import getopt
import time
import sys
usage = {0} -b < battery-count >-c < hauler-count > d < mar ket count >
usage = usage, format ( sys argv [ 0] )
opts, args = getopt. getopt (sys argv [1: ] , b:c:d:)
b = None
c = None
d = None
try :
for i in range ( len ( opts )) :
if opts[i][0] = = b :
b = int(opts[i][l])
elif opts [ i] [0] = = c :
c = int(opts[i][1])
elif opts [ i ] [0] = = d :
d = int(opts[i][l])
except ValueError e:
sys e x i t ( 1)
if b is None or c is None or d is None:
sys e x i t ( 1)
batteries haulers destinations = c r e a t e i n s t a n c e ( b c, d)
bheader = [x , y , ql0 q35 , q50 , q65 q90 invcap v0]
hheader = [1, u, loadsize, d r y f e e , splitfee]
d he ad er = [ 1 , u , x , y ]
instnum = int ( 1 0 ( time t ime ( ) 14 2 9 3 8 4 6 6 9.305739))
write.dat a ( bheader batteries, modelbatteries{0}.csv. format (instnum))
write_data(hheader haulers , model h aulers {0}.csv format (instnum))
w r i t e _d at a ( d head er destinations ,
model destinations { 0 } c s v format (instnum))
time, sleep (.25)
w = csv writer ( fl )
w. writerow(header)
w. writerows (rows)
Model Data
# * coding : utf8 *
72


2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import pandas as pd
import os path
import numpy as np
def r e a d _i n s t a n c e ( i nt r e p o instnum):
bfn = os.path.join(intrepo, model batteries {0}. csv format (instnum))
hfn = os.path.join(intrepo, modelhaulers{0}.csv. format (instnum))
dfn = os.path.join(intrepo, modeldestinations{0}.csv. format (instnum))
bdf = pd.read_csv(bfn)
hdf = pd. read_csv( hfn )
ddf = pd.read_csv(dfn)
return bdf, hdf, ddf
def b uild_b t m ( bdf ddf, 1, u):
loads = bdf [bdf [v0]> 1 ]
loadcounts = np floor ( loads [ vO ] / u)
addtlloads = 1 o a d c o u n t s==0
loadsizes = 0 loads [ vO ]
loadsizes [loadcounts= = 0] = loads [ loadcounts = = 0] [ vO ]
loadsizes [ lo ad c o u nt s > 0] = u
# build costs and destinations
destinations = d d f [ u ]
costs = [ ]
for idx in ddf. index:
x = np ab s(ddf.loc[idx, x ] loads [ x ] )
y = np ab s(ddf.loc[idx, y ] loads [ y ] )
d = x + y
costs append ( np array (d) )
costs = np array(costs).transpose(). tolist ()
loads loc [:, load sizes ] = loadsizes
loads. 1 o c [ : loadcounts] = loadcounts
loads, loc [loadcounts= = 0, loadcounts] = 1
return loads costs destinations
C.3 Method Execution
# * coding : utf8
runinstancein J3TIVL solver and STINV saveoutputto JDEOVL file
from newcastle opt import (btm, nfr exc)
from newc ast le opt import c ar d i nve nt or y
from newc ast le opt import goal
from newcastle utils import modeldata
from newcastle utils simulateinv import run.sim
import numpy as np
import time
def run_btm(loads, costs, destinations):
loadsizes = loads [loadsizes]
loadcounts = loads [loadcounts]
nestloads = [[ loadsizes [i] ,]*loadcounts[i] for i in loadsizes. index]
73


19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
loads = [ ]
newcosts = []
for i, ii in zip ( range ( len ( nest loads )) loadsizes index):
loads extend( nestloads [ i ])
fo r j in range(int(loadcounts [ ii ] )) :
newcosts.append(costs[i])
b = btm BTMProblem ( lo a d s newcosts destinations r e 1 ax e d = F a 1 se )
try :
start = time. clock ( )
cstart = time time ()
solution = b.so lve ()
estop = time time ()
stop = time, clock ()
feasible = True
obj = b.objective
lpobj = b.lpobjective
except exc FailedSolution :
estop = time time ()
stop = time, clock ()
feasible = False
obj = np.inf
if b.lpobjective is None:
lpobj = np.inf
else:
lpobj = b.lpobjective
ptime = stop start
cltime = estop cstart
return sum (loads) sum (destinations), obj, lpobj, ptime, cltime, solution
def run.nfr ( loads costs, destinations):
loadsizes = loads [loadsizes]
loadcounts = loads [loadcounts]
nestloads = [[ loadsizes [i] ,]*loadcounts[i] for i in loadsizes. index]
loads = [ ]
newcosts = []
for i, ii in zip ( range ( len ( nest loads )) loadsizes index):
loads extend( nestloads [ i ])
fo r j in range(int(loadcounts [ ii ] )) :
newcosts.append(costs[i])
obj = np.inf
lpobj = np.inf
orig = btm BTMProblem ( lo ad s newcosts destinations r e 1 a x e d = Tr ue )
b = n f r BTMProblem (loads newcosts destinations)
try :
start = time. clock ( )
cstart = time time ()
orig. so lve ()
lpobj = orig.objective
solution = b.so lve ()
estop = time time ()
stop = time, clock ()
feasible = True
obj = b.objective
except exc FailedSolution :
solution = None
74


74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
estop = time time ()
stop = time, clock ()
feasible = False
ptime = stop start
cltime = estop cstart
return sum (loads) sum (destinations), obj, lpobj, ptime, cltime, solution
def run.goal ( loads costs destinations):
loadsizes = loads [ loadsizes ]
loadcounts = loads [ loadcounts ]
nestloads = [[ loadsizes [i] ,]*loadcounts[i] for i in loadsizes. index]
loads = [ ]
newcosts = []
maxcost = 0
for i, ii in zip ( range ( len ( nest loads )) loadsizes index):
loads extend( nestloads [ i ])
fo r j in range(int(loadcounts [ ii ] )) :
newcosts.append(costs[i])
if m ax co st < max (costs [ i ] ) :
max cost = max( costs [ i ] )
factors = [m axco st factor * 1000 for factor in range ( len ( dest inat ions ))]
b = goal. BTMProblem (loads newcosts destinations)
try :
start = time. clock ( )
cstart = time time ()
solution = b.solve(factors)
estop = time time ()
stop = time, clock ()
feasible = True
obj = b.objective
lpobj = b.lpobjective
except exc FailedSolution :
estop = time time ()
stop = time, clock ()
feasible = False
obj = np.inf
if b.lpobjective is None:
lpobj = np.inf
else:
lpobj = b.lpobjective
ptime = stop start
cltime = estop cstart
return sum (loads) sum (destinations), obj, lpobj, ptime, cltime, solution
def run.stim (bdf hdf, ddf, periods):
caseprob = np.array([ 0.05047033, 0.27098408, 0.35709117,
0.27098408, 0.05047033])
p = c ar d i n ve nt o r y STIMProblem ( b d f hdf, ddf, periods caseprob)
try :
start = time. clock ( )
cstart = time time ()
solution = p.so lve ()
estop = time time ()
stop = time, clock ()
75


129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
feasible
True
obj = p.objective
lpobj = p.lpobjective
siprod = p s h ut i n vo 1 ume
einv = p e n d i ng i n v e n t o r y
haulvol = p.haulvolume
hauls = p. hauls
batteryhauls = p b a 11 er y h a u 1 s
cardz = batteryhauls sum ( )
b at t ery h au 1 v o 1 u m e = p b at t ery h a u 1 v o 1 u m e
except exc FailedSolution :
estop = time time ()
stop = time, clock ()
feasible = False
obj = np.inf
if p.lpobjective is None:
lpobj = np.inf
else:
lpobj = p.lpobjective
siprod = None
einv = None
haulvol = None
hauls = None
batteryhauls = None
b at t er y h au 1 v o 1 u m e = None
print SI siprod
print HAUL haulvol
print POSSIBLEHAUL , (hdf[loadsize ] hdf [u ] ) .sum 0
print POSSIBLE JPROD , (bdf [ q90 ] .sum ())* len(periods)
print BEGINNING ^INV , bdf [ vO ] sum ( )
print HAUL^CHOICE , cardz
ptime = stop start
cltime = estop cstart
return einv siprod , obj lpobj ptime , cltime hauls , \
haulvol batteryhauls b at t er y h a u 1 vo 1 u me
len(periods)
def run_stim_mpc(bdf, hdf, ddf, periods):
^Run^a- shift ing ^time^horizon JVfPC-met hod on^STIM-
mvbdf = bdf.copyQ
window = 3
solvetimes = []
simtimes = []
shutinvol = []
objval = 0
lpobjval = 0
cardz = 0
finalhaulvol = np zeros (( bdf shape [0] hdf.shape[0] len (periods )) )
for tidx in range(len( p er iods )):
perwin = periods[tidx:tid x+window ]
einv siprod obj lpobj ptime cltime hauls , \
haulvol batteryhauls b at t er y h au 1 v o 1 u me = \
run.st im (mvbdf hdf, ddf, perwin)
solvetimes append (cltime)
onedayhauls = batteryhaulvolume [ : : 0]
76


184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
cardz += ( o ne d ay hau 1 s > 0 ) a st y p e ( np i n t ) sum ( )
shutinvo1.append(siprod [0])
finalhaulvol [: , tidx] = batteryhaulvolume [ : , 0]
startsim = time time ()
einvs haulvs sis drys prods batinv = \
run.sim (mvbdf onedayhauls)
endsim = time time ()
s i mt i me s append ( end s im startsim)
mvbdf [ vO ] = batinv
print SI shutinvol
batteryhauls = (finalhaulvol > 0) asty pe ( np int )
haulvol = finalhaulvol. sum ( axis = (0 , 1))
print HAUL haulvol
print POSSIBLE^HAUL , ( hdf [ loadsize" ] * hdf[u]). sum () * len(periods)
print POSSIblE_PROD , (bdf[q90j. sum ())* len(periods)
print BEGINNINGHNV b df [ vO ] sum ( )
print HAUL^CHOICE cardz
haullocations = (haulvol >0).astype( np .int)
minhaul = 4 hdf[ loadsize]. min () / 9 .
objval = h a u 1 v o 1 sum ( ) 2 sum ( s h u t i n v o 1 ) minhaul cardz
return einv shutinvol objval None sum (solvetimes) sum (simtimes) ,\
haullocations haulvol batteryhauls finalhaulvol
if __name__^= ..main.. :
import getopt
import time
import sys
import os path
import anydbm
import pickle
usage = {0} p [ goal | stim |btm| nfr ] b_< i n s t Mnumber>H- i _< input _p at h>_ o _ + \

usage = usage, format ( sys argv [ 0] )
opts, args = getopt.getopt(sys.argv[1:], p:b:i:o:)
b = None
0 = None
1 = None
rungoal = False
runstoch = False
runmpc = False
runbtm = True
runnfr = False
try :
for i in range ( len ( opts )) :
if opts[i][0] = = b :
b = int(opts[i][l])
elif opts [ i ] [0] = = o :
o = opt s [ i ] [1]
elif opts [ i ] [ 0 ] = = i :
ifile = opts[i][1]
elif opts [ i ] [0 ] == p :
if opts [ i ] [ 1 ] == goal :
77


239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
rungoal = True
elif opts [ i ] [1] == s t i m :
runstoch = True
elif opts [ i ] [1] == rup c :
runmpc = True
elif opts [ i ] [1 ] == n f r :
runnfr = True
except ValueError e:
print usage
sys.exit(3)
if b is None or o is None or ifile is None:
print usage
sys.exit(4)
i f not os .path. exists (o) :
print usage
sys.exit(5)
elif not os.path.exists(ifile ):
print usage
sys.exit(6)
if rungoal or runstoch or runnfr or runmpc:
runbtm = False
bdf, hdf, ddf = modeldata.read_instance(ifile b)
if rungoal or runbtm or runnfr :
u = max (hdf [ loadsize ])
1 = u / 3.
loads costs destinations = mo del d at a b u i Id _b t m ( b d f ddf, 1 u)
i f runbtm :
haulvol, volcap , o b j v a 1 , lpobjective , clocktime walltime sol = \
run.btm (loads , costs , destinations )
elif runnfr :
haulvol, volcap , o b j v a 1 , lpobjective , clocktime walltime sol = \
r u n n f r ( lo ad s , costs , destinations )
else:
haulvol, volcap , o b j v a 1 , lpobjective , clocktime walltime sol = \
run.goal ( loads , costs , destinations )
print RESULTS:_1pobj -{0},~o b j { 1} format (lpobjective objval)
db = anydbm open ( o w )
db [ instance {0} format (b ) ] = pickle, dump is({ haulvol : haulvol,
volcap: volcap , object ivevalue : obj v al ,
lpobjective : lpo bjective clocktime : clocktime ,
loadcount : len ( 1 o a d s ) , destcount : len(destinations) ,
walltime: walltime, solution: sol})
db. c1o s e ( )
elif runstoch or runmpc:
periods = range(1, 6)
if runstoch:
einv , siprod , obj , lpobj , clocktime , walltime , haulchoice , \
haulvol batteryhauls b at t er y h a u 1 vo 1 u me = \
run.stim (bdf hdf, ddf, periods)
else:
einv , siprod , obj , lpobj , clocktime , walltime , haulchoice , \
haulvol batteryhauls b at t er y h a u 1 vo 1 u me = \
run.stim.mpc (bdf hdf, ddf, periods)
print RESULTS :lpobj{0},obj{l} format (lpobj obj)
db = anydbm open ( o , w)
78


294
295
296
297
298
299
300
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
pickle dumps ( { h a u 1 v o 1 : haulvol,
db [ instance {0} format (b ) ] =
objectivevalue: obj, einv: einv, siprod: siprod,
lpobjective: lpobj, clocktime: clocktime,
walltime: walltime, haulchoice: haulchoice,
batteryhauls : batteryhauls ,
b at t ery h au 1 v o 1 u me : b at t ery h au 1 vo 1 u me } )
db. c1o s e ( )
C.4 Data Visualization
# * coding : utf8
import matplotlib
matplotlib use( Agg )
import matplotlib pyplot as pit
import numpy as np
tolerance = 0.001
def plot_mipgap_vs_lp (lpobj, obj, fn, s ho wop t co nd=Tr ue ) :
fig = pit figure ()
ax = fig.add.subplot(lll)
ax grid ()
fig set.facecolor ( w )
ax. set.ylabel ( MIP-Gap )
ax set.xlabel ( Optimal _LPR_ Objective )
ax scatter(lpobj objlpobj color=k marker = D )
yO = 0.
xO = 0.
x 1 = np .ceil (max( lpobj [ ~ np isinf (lpobj)]) / 10.) * 10
y1 = yO + tolerance (xl xO)
if showoptcond:
ax. plot ([xO xl], [yO, y1] , k )
y 1 = np .ceil (max( obj) / 10.) * 10
fig.savefig(fn)
return fig
def plot_clock(ct, fn):
fig = pit figure ()
fig set.facecolor ( w )
ax = fig.add.subplot(lll)
ax grid ()
ax. set.x label ( Wall Clock Time )
ax. set.ylabel ( InstanceCount )
ax h i s t ( ct bins=20, color=k)
fig.savefig(fn)
return fig
def p 1 o t c o m p le x i t y ( 1 c ct obj fn):
fig = pit figure ()
fig set.facecolor (w )
ax = fig.add.subplot(lll)
ax grid ()
ax set.xlabel ( Integer^Variables )
ax. set.ylabel ( Clock-Time )
79


45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
plot feasible instances in dot and infeasible in -j-
1c = np.array(lc)
c t = np. array(ct)
obj = np. array(obj )
ax sc at t er ( lc [~ np isinf(obj)], ct[~ np isinf (obj )] c= k mar ker= o)
ax sc at t er ( lc [ np isinf (obj)], c t [ np isinf (obj)], c=k, marker = + )
fig.savefig(fn)
return fig
if __name__^= ..main.. :
import getopt
import time
import sys
import os path
import anydbm
import pickle
import numpy as np
usage = { 0 } [ n ] i _< database_file >_ o _< destination_path_for_plots >
usage = usage, format ( sys argv [ 0] )
opts, args = getopt. getopt (sys argv [1:] , n i :o: )
opath = None
i n f i 1 e = None
showcondline = True
try :
for i in range ( len ( opts )) :
if opts[i][0] = = o :
opath = op t s [ i ] [1]
elif opts [ i ] [ 0 ] = = i :
infile = op t s [ i ] [ 1 ]
elif opts [ i ] [0] = = n :
showcondline = False
except ValueError e:
print usage
sys.exit(3)
i f not os .path, exists ( infile ):
print usage
sys.exit(6)
grab the data and plot each of three plots
db = anydbm open ( infile)
data = {}
for k in db.keysQ:
try :
data[k] = pickle.loads(db[k])
except:
continue
next plot the clock time histogram
c t = np array ( [ data [k ] [ walltime] for k in data])
fn = opath. format ( time_hist )
plot_clock(ct, fn)
start with the mip gap vs Ipr objective
try :
1p o b j = np.array([data[k][lpobjective] for k in data])
80


100
101
102
103
104
105
106
107
108
109
110
except KeyError :
print data[k].keys()
raise
0 b j = np array([data[k][objectivevalue] for k in data])
fn = opath. format ( mip_vs_obj )
plot_mipgap_vs_lp(lpobj, obj, fn, showcondline)
# lastly plot the time vs integer variables
1 c = np. array ([data [k][ loadcount ] * data[k] [destcount ] for k
wt = np array ( [ data [k ] [ walltime] for k in data])
fn = opath format ( time.vs.intcount)
p 1 o t c o m p le x i t y ( i c wt obj fn)
in data])
81


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
APPENDIX D. STIM Solutions
D.l Solver
# * coding : utf8
import pulp
import numpy as np
class STIMProblem ( object ) :
StochasticInventory Management
def __init__(self batteries haulers destinations periods ,
caseprob):
initialinv = [vO for vO in batteries [ vO ] ]
self, -data = {batteries: batteries,
haulers : haulers ,
destinations : destinations ,
periods : periods ,
hr ng : range (len( haulers )) ,
br ng : range ( len( batteries ) ) ,
dr ng : range ( len ( destinations)) ,
caseprob: caseprob ,
inventories : initialinv}
def clearModel ( self ) :
self, -model = {}
def buildModel(self gamma):
^Create^st ochast ic ^modeL
self clearModel ()
prob = pulp.LpProblem( STIM sen s e=p u lp Lp Maximize )
obj = self.buildVars(prob, gamma)
cnum = 0
for c in self, .model [ constraints]:
prob += c, CONS{ 0} . format ( cnum)
cnum += 1
prob += obj , stochastic_inventory management
self, .model [objective] = obj
self, .data [problem] = prob
return prob
def buildVars(self prob, gamma) :
Createthevariables .andconst raint s for theproblem
batteries = self -data [ batteries ]
haulers = self, .data [haulers]
periods = self, .data [periods]
caseprob = self, .data [caseprob]
82


46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
initialinv = self, -data [inventories]
self, .model [ constraints] = constraints = []
self, .model [ variables] = variables = []
self .model [ shutinprod ] = shutinprod = {}
self, .model [endinginv] = endinginv = []
self, .model [haulchoice] = haulchoice = {}
self._model[haulervols] = haulervols = [{} for idx in haulers, index]
for t in periods :
shutinprod [t] = []
haulchoice [t] = []
obj = 0
s e 1 f .model [ obj v al ] = objfcn = p u lp L p V ar i ab le ( OBJVAL )
q c a s e s = [ q 1 0 q 3 5 q 5 0 q 6 5 q 9 0 ]
scenarios = range ( le n ( q c a s e s ) )
xnum = 0
brng = self._data[brng]
fo r j in range ( len ( haulers )) :
for t in periods :
haulervols [j][t] = [[ None for i in brng] for m in scenarios]
# create the interperiod constraints for each battery
for i, ii in zip ( brng batteries, index):
vO = pulp. LpVariable(V{0}_0. format (i))
openinventory = {periods [0]1: [vO for m in range (5 ) ] }
variables.append(vO)
constraints .append(v0== initialinv [ i ])
endinginv.append(0)
for t in periods :
z = pulp.LpVariable(Z{0}_{l}. format ( i t ) ,
cat=pulp LpBinary ,
lowBound = 0, upBound = l)
openinventory [t] = []
variables.append(z)
haulchoice[t].append(z)
obj = gamma z
for m in scenarios:
v = pulp.LpVariable(V{0}_{l}_{2}. format ( i t m) ,
lowBound=0)
variables.append(v)
openinventory [t ] append(v)
x = n
fo r j in range ( len ( haulers )) :
xl = pulp.LpVariable(X{0}. format ( xnum ) ,
lowBound = 0)
xnum += 1
x.append(xl)
variables.append(xl)
haulervols [j ] [t ] [m] [ i] = xl
q= batteries. 1 o c [ i i qcases [m] ]
# track the production that exceeds tank inventory and haul
# this production is shut in
o = pulp.LpVariable(O{0}_{l}_{2}. format ( i t m) ,
lowBound = 0, upBound=q)
variables.append(o)
shutinprod [t ] append(caseprob [m] o)
# mod ify the objective: subtract shut in prod, add haul
83


caseprob [m] (sum(x) 20000*o)
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119 #
120
121
122
123
124
125
126
127
128
129
130
131
132
133 #
134
135
136 #
137 #
138
139
140
141
142
143
144
145
146
147
148
149
150 #
151 #
152
153
154
155
obj +=
constr = v == openinventory [t 1] [m] + q o sum (x)
constraints.append(constr)
constr = v <= batteries. 1 o c [ i i , invcap ]
constraints.append(constr)
maxhaul = batteries. 1 o c [ i i invcap] +\
batteries. 1 o c [ i i qcases [m] ]
constr = sum(x) <= z maxhaul
constraints.append(constr)
if t == periods [ 1 ] :
endinginv [ 1] += caseprob [m] * v
for j, idx in zip ( range (len( haulers)), haulers, index):
for t in periods :
for m in scenarios:
constr = sum (haulervols [j ] [t] [ m] ) <= haulers, loc [idx, u:
haulers, loc [idx, loadsize ]
constraints.append(constr)
obj = endinginv [ 1]
obj = l*sum ( [ sum ( shut inp r o d [ t ] ) for t in shut inp r o d ] )
constraints .append(o bj==o b j f c n)
return obj
def solve( seIf ) :
alpha = 1 e 3
haulers = self, .data [ haulers]
minhaul = haulers [loadsize]. min ( )
maxhaul = haulers [ loadsize ] max()
bestgamma = None
bestobj = 1*np .inf
bestlpobj = l*np. inf
def so 1 vep a ir (gamma) :
prob = self.buildModel( gamma)
opts = [--mip g ap , 0.001, -tmlim , 180]
solve the LP relaxation explicitly first
opts = [--nomip ]
solver = pulp GLPK-CMD ( msg=F alse p at h = /usr / b in / g Ip s o l ,
o p ti o ns = o pt s )
solver = pulp .PYGLPK(msg=False mip=False)
status = prob.solve(solver)
if pulp. LpStatusf status] != Optimal:
raise F ai 1 e d S o 1 u t io n ( s t a t u s )
lpobj = self .model [ objval ] value ()
add the cardinality fix constraints
znum = 0
for t in self .model [ haulchoice ] :
for z in self, .model [ haulchoice ] [ t ] :
if z value () < alpha:
prob += z <= 0 Z {0} format ( znum )
znum += 1
solver = pulp GLPK_CMD ( msg=F alse p at h = /usr / b in / g Ip s o l ,
o p ti o ns = o pt s )
solver = pulp PYGLPK( msg=F al se mip=False)
status = prob.solve(solver)
if pulp.LpStatusfstatus] != Optimal:
raise F ai 1 e d S o 1 u t io n ( s t a t u s )
\


156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
objective = self -model [ objval ] value ()
return lpobj objective
fo r gamma i n np linspace (minhaul / 4., 2* maxhaul , 10):
self.lpsolution = []
lpobj objective = s o 1 v e p a i r (gamma)
if lpobj > bestlpobj:
bestlpobj = lpobj
if objective > bestobj:
bestgamma = gamma
bestobj = objective
if gamma != bestgamma:
# need to run the model again
lpobj objective = s o 1 v e p a i r ( bestgamma )
self.objective = bestobj
self.lpobjective = bestlpobj
v a 1 s = [ vv value () for vv in self, .model [ variables]]
self, solution = vals
return vals
(property
def shutinvolume (self):
v o 1 s = [ ]
shutinvars = self .model [ shutinprod ]
for t in self .data [ periods ] :
v = 0
for svar in shutinvars [t ] :
v += svar value ()
vols.append(v)
return vols
property
def endinginventory( self ) :
return sum ( [v. value () for v in self, .model [ endinginv ] ] )
property
def haulvolume (self):
periods = self, -data [periods]
volumes = np zeros ( len (periods ) np float32 )
haulervols = self, .model [ haulervols ]
caseprob = self.-data[caseprob]
haulers = self._data[hrng]
batteries = self, .data [brng]
for j in haulers :
for idx, t in z ip ( range ( len ( periods ) ) periods):
for m in range(len(caseprob )) :
for i in batteries:
vol = haulervols [j ] [t] [m] [ i]
volumes [ idx ] += caseprob [m] vol value ()
return volumes
property
def hauls(self):
periods = self, -data [periods]
haulcount = np zeros (len(periods ) np float32 )
haulchoice = self .model [ haulchoice ]
85


211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
1
2
3
4
5
6
7
8
9
10
11
12
13
for i in self, -data [brng]:
for idx, t in z ip ( range ( len ( periods ) ) periods):
haulcount [idx] += np ceil ( haulchoice [t ] [ i ] value ())
return haulcount
(property
def batteryhauls (self):
Report-haul-indicator by ^ battery^byperiod
brng = self._data[brng]
periods = self, -data [periods]
haulchoice = self .model [ haulchoice ]
haulind = np zeros (( len(brng) len (periods )) )
for i in self..data[brng]:
for idx, t in z ip ( range ( len ( periods ) ) periods):
haulind [i, idx] = np ceil ( haulchoice [t ] [ i ] value ())
return haulind
property
def batteryhaulvolume(self):
Report-haul-indicator by ^ battery^by-period-
brng = self..data[brng]
hrng = self._data[hrng]
periods = self, -data [periods]
haulervols = self, .model [ haulervols ]
caseprob = self.-data[caseprob]
haulers = self..data[hrng]
batteries = self, .data [brng]
shp = (len(brng), len(hrng), len (periods ))
haulvol = np zeros (shp np float32 )
for j in haulers :
for idx, t in zip( range ( len ( periods ) ) periods):
for m in range(len(caseprob )) :
for i in batteries:
vol = haulervols [j][t][m][i]. value ()
haulvol [i, j, idx] += caseprob [m] * vol
return haulvol
D.2 Simulation
# * coding : utf8
import anydbm
import pickle
import matplotlib
matplotlib use( Agg )
import matplotlib pyplot as pit
import numpy as np
from newcastle utils import modeldata
import os path
def get_rate(mu, sigma):
Get^rates^distributed^normally ,-dont allow-zero
86


14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
rate = np .random. normal( mu, sigma, siz e=mu shape)
rate [ rate <0] = 0
return rate
def get.shutin.volume (rate vO invcap, haulvol):
Determine_theexcess_volumefromtherate and availablesapce_
sivolume = vO + rate haulvol invcap
sivolume [ sivolume <0] = 0
return sivolume
def get.haul.volume (rate vO hauldecision maxhaul minhaul):
^Use^the^haul^f lag ^and^ availab le-volume-to^get^total^haul^
a v a i 1 ab 1 e v o 1 u m e = hauldecision (vO + rate)
shortvolume = a v a i 1 ab le v o 1 u me maxhaul
shortvolume [ shortvolume >0] = 0
haulvol = maxhaul + shortvolume
dryload = ( h au 1 v o 1 0)
haulvol [ hau 1 vo 1 return maxhaul + shortvolume dryload astype (np int)
def g et _e n d i n g i n v e n t o r y ( r at e vO haulvolume):
Usethestandardrelationshiptodetermineendinginventory
return vO + rate haulvolume
def simulate.hauling (mu, sigma, invcap, haulvolume, vO):
q = get_rate(mu, sigma)
haulind = ( h au lvo lu me > 0 ) a st y p e ( np i n t )
hv dryload = get_haul_volume(q, v0, haulind, haulvolume, 45)
sivol = get.shutin.volume (q vO invcap hv)
actprod = q sivol
einv = get.ending.inventory(q, vO hv)
einv [einv>invcap ] = invcap [ einv>invcap ]
return einv hv, sivol dryload actprod
def test.sim () :
# test expected ending inventory
cap = np. array ([300. , 300.])
haulvol = np. array ([180. , 185.])
mu = np array ( [30 , 60.])
sigma = 0.05 * mu
vO = np array ( [170. , 200.])
einv = [ ]
m = 1000
for i in range(m):
inv hv sivol dryload actprod = \
simulate.hauling (mu, sigma, cap haulvol vO)
einv.append(inv)
meaninv = np array ( einv ) mean ( axis=0)
print Planned :{0}Simulated : { 1} format ( vO + mu haulvol meaninv )
def run.sim (bdf haulvol s d f a c t = 0.0 5 ) :
87


69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
^Execute^simulation^:
try :
numper = haulvol shape [2]
except IndexError:
numper = 1
haulvol = haulvol. reshape ((haulvol. shape [0], haulvol shape [1] numper ) )
periods = range ( numper )
cap = bdf[invcap]
mu = bdf [ q50 ]
stddev = sdfact mu
mrng = range(1000)
sisample = np zeros ((1000, numper ) )
haulsample = np . z e r o s ( ( 1 0 0 0 numper))
d r y co u nt s a mp le = np z e r o s ( ( 1 0 0 0 numper))
prodsample = np . z e r o s ( ( 1 0 0 0 numper))
einvsample = np . zeros ((1000 ,))
einvbat = np zeros ((bdf. shape [0] ,))
for m in mrng :
inv = bdf [ vO ]
for t in periods :
inv hv sivol dryload actprod = \
simulate-hauling (mu, stddev c ap haulvol [: : t ] sum (axis=l), inv)
sisample [m, t] = sum (sivol)
haulsample [m, t] = sum(hv)
drycountsample [m, t] = sum ( dryload )
prodsample [m, t] = sum ( actprod )
einvsample [m] = sum (inv)
einvbat [:] += inv / 1000
return einvsample haulsample sisample drycountsample prodsample ,\
einvbat
def plot_einv(sim, plan, size=20):
^Plot^sim-and-plan^einv^
s i n v = np array ( [sim [k] [ einv ] mean (axis=0) for k in
p in v = np array ([plan[k] [einv] for k in sim keys ( ) ] )
figl = pit. figure ()
figl set.facecolor (w )
ax = figl. add_subplot (111)
ax b ar ( np arange (1 size +1) sinv [ : size ] ,
#hatch=/,
1abe1= Si m ul ated color=k)
ax b ar(0.5 + np arange (1 size+1), pinv [ : size] ,
#hatch =\\ ,
label=Plan color=w)
ax.set_xlabel(Instance)
ax set.ylabel ( Barrels )
ax set_xticks ([])
ax grid ()
ax. legend ( )
fig2 = pit. figure ()
fig2 set.facecolor (w )
ax = fig2.add_subplot(lll)
ax set_xlabel(Plan^Barrels)
sim keys ( ) ] )
88


124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
ax set.ylabel ( Simulated-Barrels )
ax.scatter(pinv, sinv, facecolors=none, edgecolors=k)
ax grid ()
return figl fig2
def p 1 o t h a u 1 ( sim plan):
Plotsim-vsplan-haulvol-
fig = pit figure ()
fig .set.facecolor(w)
# flatten the data series
keys = sim keys ()
periodcount = pc = sim [keys [ 0] ] [ haul ] shape [ 1]
casecount = len(keys)
n = casecount periodcount
simhaul = np zeros ((n,) np float32)
planhaul = np zeros ((n,) np float32 )
for i in range(casecount):
k = keys[i]
simhaul [ i*pc : i*p c+p c] = sim [k][ haul], mean ( axis=0)
planhaul [ i *pc : i*p c+p c] = plan [k ] [ batteryhaulvolume ] sum (axis=(0, 1))
ax = fig. add.subplot (111)
ax set.xlabel ( Plan Barrels)
ax set.ylabel ( Simulated-Barrels )
ax grid ()
ax.scatter(planhaul, simhaul, facecolors=none, edgecolors=k)
return fig
def plot_si(sim, plan, size=40):
_Plot sim-and-plan-si volumes-
skeys = sim .keys()
N = len(skeys)
M = sim. values () [0] [ si ] shape [0]
figl = pit. figure ()
figl set.facecolor (w )
ax = figl. add.subplot (111)
s i = np. zeros ((M, N) np. f1o at 3 2 )
p 1 a n s i = np array ( [sum (plan[k][siprod]) for k in skeys])
for n in range(N):
k = skeys[n]
s i [ : n] = s i m [ k ] [ s i ] sum ( axis=l)
ax .boxplot(si [: , : size])
ax plot ( range (1, size+1), plansi [ : size ] , k )
ax.set.xlabel(Instance)
ax set.ylabel ( Barrels )
ax set.xticks ([])
ax grid ()
fig2 = pit. figure ()
fig2 set.facecolor (w )
ax = fig2. add.subplot (111)
ax set.xlabel (Plan Barrels)
ax set.ylabel ( Simulated-Barrels )
89


179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
ax scatter (plansi
marker:
ax scatter (plansi
marker:
ax scatter (plansi
marker:
ax grid ()
ax. legend ( )
return figl fig2
s i max ( axis = 0) facecolors=none
o label = max )
s i mean ( axis = 0) facecolors=none
+ label- mean )
s i min ( axis = 0) facecolors=none
D , 1 a b e 1= min )
e d g e c o 1 o r s= k ,
e d g e c o 1 o r s= k ,
e d g e c o 1 o r s= k ,
def p 1 o t _d r y ( sim ) :
Plotsimdryloads
fig = pit figure ()
fig .set.facecolor(w)
def make.plots (repo dest, plandbname, simdbname) :
^Organize^results -for-plots and^ generate ^them^
output = anydbm open ( simdb name )
plan = anydbm open ( plandbname )
indata = {}
simdata = {}
plandata = {}
for k in output keys ( ) :
junk, instnum = k. split ()
simdata [k] = pickle loads (output [ k ] )
plandata [k] = pickle loads (plan [ k ] )
bdf, hdf, ddf = modeldata.read.instance (repo instnum)
indata [k] = {bdf: bdf, ddf: ddf, hd f : hdf}
print Plotting^hauls .
fig = plot-haul ( simdata plandata)
fn = os.path.join(dest, stim_sim_haul.eps)
fig.savefig(fn)
print Plotting^ending^inventory .
figl fig2 = plot.einv (simdata plandata, size=30)
fn = os.path.join(dest, stim_sim_einv_bar.eps)
figl.savefig(fn)
fn = os.path.join(dest, stim_sim_einv_scatter.eps)
fig2.savefig(fn)
print Plottng_si_volumes ...
figl fig2 = plot.si (simdata plandata, 30)
fn = os.path.join(dest, stim_sim_si_boxplot.eps)
figl.savefig(fn)
fn = os.path.join(dest, stim_sim_si_scatter.eps)
fig2.savefig(fn)
if __name__^= __main__ :
import getopt
import time
import sys
usage = {0} i _< database_file >_ o _< destination_path_for_plots >_ r _< instance_repo >
usage = usage, format ( sys argv [ 0] )
90


234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
opts, args = getopt.getopt(sys.argv[l:], i:o:r:)
opath = None
i n f i 1 e = None
repo = None
try :
for i in range ( len ( opts )) :
if opts[i][0] = = o :
opath = opts[i][1]
elif opts [ i ] [ 0 ] = = i :
infile = opts [ i ] [ 1 ]
elif opts [ i ] [ 0 ] = = r :
repo = opts[i][1]
except ValueError e:
print usage
sys.exit(3)
i f not os .path, exists (repo):
print usage
sys.exit(4)
i f not os .path, exists ( infile ):
print usage
sys.exit(6)
db = anydbm open ( infile)
if os.path.exists(os.path.join(opath, output.db)):
fig = w
else:
fig = c
output = anydbm open (os.path.join(opath, output, db), fig)
data = {}
for k in db:
d at a [k] = pickle. loads(db[k])
junk, inst = k. split ( )
bdf, hdf, ddf = modeldata.read.instance (repo inst)
einvsample haulsample sisample d ry c o unt s amp le prodsample ,\
batinv = run.sim (bdf data [k ] [ batteryhaulvolume ] )
output [k] = pickle, dumps ( { e i n v : einvsample, haul: haulsample,
si: sisample, dry: drycountsample,
prod : prodsample batinv : batinv})
print . ,
output.sync()
output. close()
91


Full Text

PAGE 1

OPTIMIZATIONMETHODSINLIQUIDSHAULINGANDINVENTORY by WESLEYO.DYK B.S.,AppliedMathematics,UniversityofColoradoDenver,2003 Athesissubmittedtothe FacultyoftheGraduateSchoolofthe UniversityofColoradoinpartialfulllment oftherequirementsforthedegreeof MasterofScience AppliedMathematics 2015

PAGE 2

ThisthesisfortheMasterofSciencedegreeby WesleyO.Dyk hasbeenapprovedforthe DepartmentofMathematicalandStatisticalSciences by AlexanderEngau,Advisor StephenBillups,Chair BurtSimon ChrystanthosGounaris July20,2015 ii

PAGE 3

Dyk,WesleyO.M.S.,AppliedMathematics OptimizationMethodsinLiquidsHaulingandInventory ThesisdirectedbyAssociateProfessorAlexanderEngau ABSTRACT Crudeoilandwastewateraretwoliquidproductsofthepetroleumexploration andproductionsectorwhicharedistributedusingtruckedtransportation.Managing theresourcestohandleinventory,transportationanddispositionoftheproductsis ahighlycomplextask.Seekingoperationaleciencyisaidedbytheuseofoptimizationmodels.Initialapproachestosolvingthesemodelsarereportedandhave seensomepositiveresults.Modelenhancementandupdatedsolutionmethodsare requiredfortheiruseinoperationsenvironments.Thisthesisincludesdiscussionof thedevelopmentandapplicationofthesemodelsandmethodsusingmathematical programmingandheuristics. Theformandcontentofthisabstractareapproved.Irecommenditspublication. Approved:AlexanderEngau iii

PAGE 4

DEDICATION Thisthesisisdedicatedtomyfamily;Jacinda,whohasstoodbymeandsupportedmy eortsallthroughschool,andourchildren,DruShaleyandTateWesley.Iwouldn't bethemanIamwithoutthem. iv

PAGE 5

ACKNOWLEDGMENT ThisthesiswouldnothavebeenpossiblewithouttheguidanceofDr.Alexander Engau.Hehasshownmethewaythroughresearchandpublicationandbeenan excellentteacher.IwouldalsoliketoacknowledgeMichaelC.Maguireforsupporting mypursuitandhelpingmebalanceworkandschool,whichhavebeenessentialto completingthisthesis. v

PAGE 6

TABLEOFCONTENTS Tables........................................viii Figures.......................................ix Chapter I.INTRODUCTION...............................1 1.1GeneralBackground...........................1 1.2LiquidInventoryManagement.....................2 1.3DecisionProcess.............................4 1.4MajorContributionsandFindings...................6 1.5ThesisOrganization...........................7 II.LITERATUREREVIEW...........................8 2.1ModernApproachestoInventoryandTransportationOptimization.8 2.2RelationshipBetweenLinearProgrammingandNetworkFlows...10 2.3TransportationProblem.........................11 2.4GeneralizedNetworkFlows.......................11 2.5GeneralizedAssignmentProblem....................13 2.6DispatchSchedulingModel.......................13 III.STOCHASTICINVENTORYMANAGEMENT..............19 3.1StochasticModelDescription......................19 3.2SolutionMethods............................27 3.3Simulation................................29 3.4ComputationalResults.........................30 IV.BATCHEDTRANSPORTATIONMANAGEMENT............38 4.1ModelFormulation...........................38 4.2UnderlyingNetworkModel.......................39 4.3ComplexityofBTM...........................42 4.4ApproximizationMethodUsingNetworkFlowRelaxation......42 vi

PAGE 7

4.5GoalProgrammingModel........................49 4.6SolutionMethods............................49 4.7ComputationalResults.........................53 V.CONCLUSION................................60 References ......................................62 Appendix A.DataGeneration................................65 B.MethodExecution...............................68 C.GeneralMethods...............................71 C.1InstanceGenerator...........................71 C.2ModelData...............................72 C.3MethodExecution............................73 C.4DataVisualization............................79 D.STIMSolutions................................82 D.1Solver...................................82 D.2Simulation................................86 E.BTMSolvers..................................92 E.1DirectSolver...............................92 E.2GoalProgrammingSolver........................94 E.3NetworkFlowRelaxationSolver....................98 vii

PAGE 8

TABLES Table 2.1DSM:Index,ParameterandVariableNotation..............14 3.1STIM:Index,ParameterandVariableNotation..............20 4.1BTM:Index,ParameterandVariableNotation..............40 4.2BTM:GreedyHeuristicExample1.....................51 4.3BTM:GreedyHeuristicExample2.....................52 A.1DataGeneration:BaseDataforGenerationofBatteries.........66 A.2DataGeneration:ProductionRateIntervalProbabilities.........66 A.3DataGeneration:BaseDataforGenerationofHaulers..........67 B.1DataGeneration:StructurefortheDataGeneratedinEveryInstance.70 viii

PAGE 9

FIGURES Figure 3.1STIM:MIPGapvs.OptimalLPRObjective...............31 3.2STIM:WallTimeHistogram........................32 3.3MPC:WallTimeHistogram.........................32 3.4MPC:EndingInventorySimulationSample................34 3.5MPC:EndingInventorySimulationScatterPlot..............35 3.6MPC:HaulVolumeSimulationScatterPlot................35 3.7MPC:SIVolumeSimulationSample....................36 3.8MPC:SIVolumeSimulationScatterPlot.................36 4.1BTM:FeasibleRegionDiagram.......................45 4.2BTM:MIPGapvs.OptimalLPRObjective................55 4.3BTM:WallTimeHistogram.........................56 4.4BTM:WallTimevs.InstanceSize.....................56 4.5GBTM:MIPGapvs.OptimalLPRObjective...............57 4.6GBTM:WallTimeHistogram........................57 4.7GBTM:WallTimevs.InstanceSize....................58 4.8NFR:MIPGapvs.OptimalLPRObjective................58 4.9NFR:WallTimeHistogram.........................59 4.10NFR:WallTimevs.InstanceSize.....................59 ix

PAGE 10

CHAPTERI INTRODUCTION Inpetroleumexplorationandproduction,productionoperationsaretheactivities thatconvertthecapitalexpenditureofexploringanddrillingforhydrocarbonsinto arevenuestreamforbusiness.Thecommonreferenceforacompanythatexplores, drillsforandproduceshydrocarbonsis E&Pcompany oralso operator .Thistypeof companyisinthe upstream sectoroftheindustry.Manyindividualactivitiesfallinto thecategoryofproductionoperationsofanE&Pcompany.Theactivityofdelivering productstosalesorotherdispositionfallsunderproductionoperations.Theliquid productsinmanyproducingeldsaredeliveredviatruckedtransportation. 1.1GeneralBackground Theterm eld isusedinmultiplewaysinthepetroleumsector.Themostspecic usagereferstothesurfaceoftheearthcoveringaresource-richgeologicformation, typicallymanysquaremilesinarea.Ageologic formation isaccessedanduids extractedbydrillingwellsfromthesurfacedowntotheformation.Oilwellsand gaswellsbothproduceliquids,typicallyreferringtotheliquidhydrocarbonsascondensate[6]whenproducedfromanaturalgaswell.Wellstypicallyproduceintoa common tankbattery ,whereproductsareseparatedbyphaseandsaleableliquids. Thesefacilitieslocaltothewellscontainequipmentthathelpstopreparetheproducts forsaleanddisposition.Thetanksinabatteryprovideonemethodofmeasurement forcustodytransfer,whichisthetriggerforconvertingtheproductintorevenue. Productiontanksinbatteriesarethemaininventoryvesselsforoilandwaterin upstreamoperations. Thelogisticsinvolvedwithsellingcrudeoilintheabsenceofpipelineinfrastructureinvolvethelocalinventoryattankbatteriesandtruckingcapacity.Theuseof theseresourcescancontributesignicantlytothebottomlineofoperations.This 1

PAGE 11

isascenarionotseenbyallcompanies;inestablishedliquid-richelds,thereisinfrastructureforproducttransportation,namelypipelines.Indevelopingelds,the infrastructuredoesnotexistattheoutsetandiscreatedonlyifitiseconomicallyfeasible.Withoutpermanentinfrastructure,resourcemanagementbecomesaninvolved processthat,ifnotdoneright,cangeneratewasteandreduceprot. Therearemanyaspectsoftransportationingeneralthatapplytothetransportationofliquidsproducedinproductionoperations.Manytimes,adeliveryhappens betweenasingletankbatteryandmarketorothercommondeliverypoint.Thesingle truckdeliveryisreferredtoasa load .Eachloadhauledcoststheoperatoranamount determinedbycontract.Thesefeescanvarybyloadsize,distanceandproduct.We refertomarketsandotherdeliverypointsas destinations ingeneral.Destinations suchas disposalwells typicallyexistatdierentlocationsaroundaeld,whichare usedforsafedisposalofwastewaterfromproductionoperations[6].Inspecicscenarios,productispickedupfromtwoormoretankbatteriesbeforedeliveringthe liquids.Thisisreferredtoasa splitload .Haulingfrommultiplelocationsinone truckloadiscommonwithwater,butitismostcommontohaulfromasinglelocationwhentransportingoilandcondensate.Acommonsituationishavingonlya singlemarketforoilandcondensate,dependentuponexistingcontractsorlocation. Inmorecomplexsituations,multipledestinationscanreceivetheproductfromaparticulartankbattery.Thisisalsotruewithwater.Inthesesituationswhereproducts are dedicated tospecicdestinations,thenumberofdecisionsisreduced.Thisisa simpliedviewofrelationshipsbetweencrudeoilcustomersandtheirdeliverypoints, butwillserveourpurposes. 1.2LiquidInventoryManagement Liquidinventoryisachallengeformultiplereasons.Asmentionedbefore,water inventorydiminishestheeectivecapacityforoil.Crudeoilininventoryiscash 2

PAGE 12

thathasnotbeenrealized.Theenvironmentalchallengesofoilsittingintanksare becomingmoreprevalent.Deadoilreleaseshydrocarbonsingaseousformthatmust beinjectedintoagaslineorburnedtoreducepollution.Modernequipmentisbuilt tohandletheseemissionsandmaximizesafetyofpersonnelperformingtransfers.The combinationoftechnologyandproperproceduresensuresthatthetransfersaresafe andsound.Withintheoverarchingdriveofintegrity,maximizinghaulingeciencyis paramount.Inventoryshouldbeminimizedoverall,buttoomanytankswithpartial loadscanpreventcontinuedecienthauling.Anotherpropertythatcanreducethe optionsforhaulingistheoilquality.Certainmarketshandleonlyalimitedrangeof gravityandhaverequirementsfornon-oilcontentbelowaspeciclevel. Inadditiontotheproceduralcomplications,thevolumeofliquidhauledina singleloadislimitedbythecapacityofthetruck,butthecapacityisadditionally restrictedbytheweightoftheload.TheunitsofvolumeusedintheU.S.isbarrels BBL,whichisequalto42USgallons.Thespecicgravityoftheliquidhauled, denoted G S below,playsasignicantroleinthevolumethatcanbehauledina singleload.Thewaterproducedwithcrudeoilandnaturalgascontainsmineralslike saltandisheavierthanfreshwater,whichhasaspecicgravityof1.Crudeoiland condensatearealmostalwayslighterthanwater.IntheU.S.,thecommonpracticeis torefertotheAmericanPetroleumInstituteAPIgravityofthesubstance,denoted by G A here,whichhasarelationshiptospecicgravity: G A = : 5 =G S )]TJ/F15 11.9552 Tf 11.955 0 Td [(131 : 5 ; whichgivesfreshwateranAPIgravityof10[1].Alsonote,thatAPIgravityand specicgravityhaveaninverserelationship,sothatheavieroilhaslowerAPIgravity. Ifcondensateisaslightas G A =120,thenasingleloadcouldcarryasmuchas1.78 timesasmuchcondensateaswater.Atruckthathasavolumecapacityof200BBL gallonsanda50,000lbcargoweightlimitcancarryapproximately200BBL 3

PAGE 13

ofcondensatelighterthan 141 : 5 = = 8400 = 8 : 328 )]TJ/F15 11.9552 Tf 11.955 0 Td [(131 : 5 104 : 18degreesAPI : Notethatfreshwaterweighsapproximately8.328LBperUSgallon.Sothe maximumspecicgravityofcondensatein200BBLisfoundby50000 = 8400 = 8 : 328. ThentheconversiontothestandardAPIgravityisdone.Thelimitoncarryingwater forthesametruckingconditionsisthen 50000 = 8 : 328 = 42 142 : 95BBL : Thisshowsthatatruckthatcant200BBLofliquidhasanactualcapacity thatvarieswiththespecicgravityofthatliquidandrangesfromabout140BBLto 200BBL. Contractsexistbetweenoperatorsandhaulersthatdeneaminimumnumber ofloadsorbarrelsthatwillbehauledperperiod.Crudeoilhaulerstypicallycarry single-sourceloadsandchargeextrafeesforsplitloads.Theextraoverheadinvolved withmeasuring,verifyingproductqualityandmorecomplicatedroutingdrivesthe fees. Thesecomputationsarerequiredforimplementationsregardingspecicproducts. ThedatageneratedforuseincomputationalresultsasdescribedinSectionAcanbe regardedasalreadyincommonvolumenetvalues. 1.3DecisionProcess Thedecisionstargetedforusingdata-drivenmethodsmayrunintoproblemsinimplementationduetocomplicationsfromtheorderofprocedures.Forinstance,crude oilthatmustbehauledawayviatrucksitsintanksthatiscomposedofwateraswell. Thewaterwillsettletothebottom,butitmustbedrainedtoaspeciclevelbefore thecrudeisregardedassaleable.Thequalityoftheproducthastwoindicators, 4

PAGE 14

BS&Wbeingthewaterandsediment,andAmericanPetroleumInstituteAPIgravitybeingrelatedtothespecicgravityoftheproduct.Allvolumesofcrudeinthe USAareconvertedtoAPInetbarrelsusingthemeasuredvolume,temperature,API gravityandBS&W.Thiscanchangetheeectofoptimizationintermsofcontract agreementsifgrossvolumesaretheonlyvolumesconsidered.Itcanmakeoperations lessecientifthisconsiderationisnotregardedandcrudeoilfromproducinglocationswithlowAPIgravityistakentomarketsthatdonotacceptheavyoil.The productissubsequentlyrejectedandtheloadmustbetakentoadierentdestination.Itisrecommendedtoperformcalculationsregardingrevenueandsomecapacity metricsusingAPInetbarrels.Thiswillaccountfortheweightanddollarvalueof theproduct,butcan'tbetheonlyconsiderationinhandlingamixofproductsthat includesextremelylightcondensate. Thesedecisionsaremadeinavarietyofdierentways.Insmalloperations, theleaseoperatordirectlycontactsthehaulertorequesthauls.Acentralplanning ordispatchingdepartmentisusedformanaginglargeroperations.Theseplanners usedataathandtodecidewhichrequestscanbedelayedifnecessaryandwhich todistributetohaulersimmediately.Theprimarymotivationoftheplanneristo keepinventoryfrombuildingacrosstheeldasawhole,whilepreventingindividual wellsfromhavingtobeshutinduetoreachingtankcapacity.Thedecisionsof theplannerformanaginginventorythroughhaulinghastwofocalpointswherewe willconcentrate.Therstisinplanning.Betterplanningofresourcesrequiredfor haulingliquidswillresultinmoreecientexpendituretoarrangetheresources.This allowsthepreventionofwastedcapacityandresourceshortages.Thesecondfocusis onhowtoeectivelyusetransportationresourcesforminimalcost.Therearesome truckroutingdecisionsthatcomeintoconsideration,eventhoughtheoperatordoes notusuallymanagetruckeets.Oneexampleiswhenasinglelocationhasmultiple loadstobedeliveredduringoneperiod.Insomeeetoperations,asingledriverwill 5

PAGE 15

beassignedtohaulallloadsfromaspeciclocation.Thiscanleadtolessopportunity forimprovedcostreduction,butknowledgeofthissituationallowstheoperatorto moreaccuratelyestimatetravelrequiredforloads. 1.4MajorContributionsandFindings TwomodelsareexploredinChaptersIIIandIV.ThesetwomodelsareStochastic InventoryManagementSTIMandBatchedTransportationManagementBTM. Eachmodelhasdierentgoals.Themethodsusedforeachmodelaretailoredfor thosegoals,buttheyallusesimilarbasedata.Theresultsdesiredarefocusedonthe usefulnessofthesolutionsfoundbythedistinctsolutionmethodsandtherelative performancebetweenthemethods. ThegoalforSTIMisaroundtheendinginventory,totalproductionandtotal haul.Theendinginventorymetricismeaninglesswithoutknowinghowmuchtotalproductionresultsfromaplan.Therelationshipbetweenthethreemetricsis straightforward.BTMhasanobjectivevaluethatrelatesdirectlytoametricthat canshowperformanceofaparticularsolution.Performanceintermsofthatobjective willbediscussedindetail.Theothertwoexperimentswillshowsolutionperformance inrelationtoothersolutions. TheoutcomesofSTIMandBTMarenotdirectlycomparable,butastheywere bothbornfromthesamemastermodel,therearerelationshipsthatwillbeexposed indata.STIMwillbedemonstratedprimarilyasalong-termplanningtoolwith smallermodelsthathavetacticalapplication.ThedemonstrationofBTMwillbeas atacticaltoolformakingdecisionsinreal-time. 6

PAGE 16

1.5ThesisOrganization Inthethreechaptersthatfollow,mathematicalapproachestothedecisionsdescribed inthischapterarediscussed.ChapterIIcoversrelevantworkontransportationand inventoryusingmathematicaloptimizationmethods.Thenexttwochapterscover newworkonthesetopics.Thenewworkseekstoincorporatestochasticaspects intoanexistingmodelandtoformulateaprecisemathematicalprogrammingmodel foraspecictransportationdecision.ChapterIIIdevelopsandtestsastochastic modelforliquidinventorymanagement.ThismodeliscalledSTIM.ChapterIV demonstratesamodelformakingtransportationdecisionsandmanipulatesthemodel tondsatiscingsolutionsundertimelimits.ThesemodelsareBTM,GoalBTM GBTM,FixedBTMFBTMandNetworkFlowRelaxationNFR.Thetwomodels STIMandBTMaddressspecicdecisionsaroundinventorymanagement.STIM targetstherstfocalpointinthedecisionprocess.Thegoalofthismodelistobring thedecisiononwhentohaulfromaspecicbatterytoalevelconsideringmorethan justwhatislocallygood.BTMtakesonthesecondfocalpoint.Eachplannedload needsadestination.Theparticularproductbeingtransporteddrivesthedecision. Crudeoilneedstobedeliveredsoastosatisfymarketcontractsandwaterdisposal isalimitedresource.BTMcanbeappliedineithercase,butthefocusisminimizing cost,whetherthatcostisdollars,distanceortime. 7

PAGE 17

CHAPTERII LITERATUREREVIEW Therearemanypublicationsrelatingtothemathematicaloptimizationoftransportationandlogisticsproblems.Intheseworks,optimizationcanrefertocostminimizationorthroughputmaximization,amongothers.Theseapproachescouldbe usefultondoptimalsolutionstotheproblemwearefacing.Eachcouldbebenecialintherealmoftacticaldecisionmakingorshortandlongtermplanningdecisions. Wewillseemodelsinthisthesismostlyinthespaceof NP ,specically NP completeproblems.Literatureunderreviewmainlyseekstodenecomputationally ecientalgorithmstoapproximatetheseproblems.TheTransportationProblem inSection2.3isknowntobein P ,whiletheGeneralizedAssignmentProblemin Section2.5is NP -complete.Foraquickreminder,theseclassesofproblems, P and NP ,arerelated.Theproblemspace P isasubsetof NP .Problemsin NP havethe propertythataproposedoptimalsolutioncanbeveriedasoptimalinpolynomial time.Problemsin P alsohavethepropertythatanalgorithmexistswhichcannd anoptimalsolutioninpolynomialtime[29]. 2.1ModernApproachestoInventoryandTransportation Optimization ChristophersonandDaveexplorethecurrenttransportationoptionsforcrudeoilin theGreatLakesBasin[20].Intheirreview,truckedtransportationisnotedashaving thehighestincidentriskperunitofproducttransported.Theyalsoobservetheneed tocontinuethismodeoftransportation.Withthatconsideration,reducingthetotal numberofmilesrequiredtodeliverproductviatruckreducesenvironmentaland safetyrisks.Thedrivingdistancefortransportationalsohaseconomicimpact,which benetsoperatorsandhaulersdirectly. 8

PAGE 18

Theschemeofoptimizinginventoryinregardstocustomerdemandshasbeen discussedandiscalledtheinventoryroutingproblemIRP[16].Astudyspecic tocrudeoiltransportationfromShenetal.[32]seekstomeetdemandsbyplanning transportationoptionsandsolvestheMIPprobleminstancesusingLagrangianrelaxation.Thisproblemisfocusedonmeetingthedemandsofcustomers.Sinceour topicisoperationsfocused,itdoesn'tseemtohavedirectapplicationtoourtopic,but regardingbatteriesascustomersanddestinationsassuppliesofcapacityfordeliverycouldchangethatperspective,makingworkregardingstochasticdemandinputs relevant. Researchontransportationandprocessingofliquidpetroleumproductsexistsfor producers,haulers,marketers,processorsanddistributors.Researcharoundreneries withpipelineandlarge-batchshiportraintransportationmodesexist,aswellas someworkaroundtruckedtransportation,withtruckingfocusedonrenedproducts. NeiroandPintodeneapetroleumsupplychainmodelwithparticularinfrastructure submodelsandmassbalanceconstraintstoassistintheplanning,strategicandtacticaldecisionsarounddeliveryofpetroleumproductsfromtheperspectiveofchemical engineering[28].ChengandDurandescribeanapproachtocrudeoillogisticsthat usessimulationandoptimalcontrol[18].Theyusearollinghorizonstrategytoseek asuboptimalcontrolstrategywithlowercomputationalrequirements.Abduljabbar andTaharbuildonthisworkwithamodel[34]andcasestudy[7]toshowthebenets ofthosetools. Someworksaroundinventoryoptimization[9,22,30]usegeneticalgorithmsto drivedecisionsoninventorymanagementinsupplychainswithoutbeingproduct specic. 9

PAGE 19

2.2RelationshipBetweenLinearProgrammingandNetwork Flows Networkowshavetheusefulpropertythatecientmethodsexistforndingsolutionsthatareguaranteedtobeintegralaslongasinputsareintegral[23].Ecient methodsforsolvinggeneralLPsalsoexist.Thesemethodscanalsoguaranteeintegralsolutions,butonlywheninputconstraintmatricesaretotallyunimodular[8]. Thesetwoseparateclassesofproblemsareconnectedinthatallnetworkowscan berepresentedaslinearprograms,butnotalllinearprogramscanberepresentedas networkows.Theconversionoflinearprogramstonetworkowsisofinterest. BixbyandCunninghamhavedescribedanalgorithmforconvertingLPstonetworkows[12]whichmakesthedeterminationifamatrixis graphic ornon-graphic. Thistermgraphicisusedinreferencetoamatrixthatisprojectivelyequivalentto anode-arcincidencematrixNAIM.WhileBixbyandCunninghamdiscusstheir algorithminmatroidtheoryterms,Baston,RamouniandWilliamsdescribeanalgorithm[10]inmatrixtermsforthosethatarelessfamiliarwithmatroidtheory.They discussteststhatcanbeexecutedtoproveamatrixisgraphic,butexplainthatany algorithmbasedonthisapproachhas"worsethanexponential"complexity.Despite thiscomplexity,theydescribeanothertesttodemonstratewhenagivenmatrixis not convertibletoaNAIM.Thisisnotusefulinthiscontextbecauseaswillbeseen inChapterIV,theconversiontoanetworkrepresentationisrequired,evenifthe originalproblemisnotmaintained.Foraproblemwith r rowsand c columnsina constraintmatrix A ,thealgorithmproposedbyBixbyandCunninghamhascomplexity O rn with n beingthenumberofnonzeroentriesof A .Thisisecientand couldbeusedtoexposetheunderlyingnetworkstructureofanLP,ifnotperform thefullconversion. 10

PAGE 20

2.3TransportationProblem Thetransportationprobleminvolvesacommodityorcommodities,depotsandmarkets.Eachdepothasacertainsupplyandeachmarkethasaspecieddemand.An answertoatransportationproblemisonethatsatisesdemandsfromgivendepots foraminimalcost.Let x 2 R n m where x ij isaquantitymeetingdemandatcustomer i bysupplier j .Theformulationoftheproblemalsoinvolvesthequantityof demandateachconsumer d i andquantityofgoodsateachsupplier s j .Acostof supplyinggoodsfromsupplier j tocustomer k is c jk .TheformulationoftheLPis minimize x X i;j c ij x ij .1a subjectto X i x ij s j 8 j .1b X j x ij = d i 8 i .1c x ij 0 8 i;j: .1d Aspecialcaseofthetransportationproblemistheassignmentproblem.This problemhasabalancednumberofsupplyanddemandnodeswhereeachnodehas demand d i =1. 2.4GeneralizedNetworkFlows Thetransportationproblemhasdemonstratedutilityinindustry,butisincapable ofprovidingsolutionstorelatedproblemswherethenetworkowshavelossorgain fromonenodetothenext.Consideranetworkwhereowleavinganodeacross anarcdoesnotequaltheowenteringthereceivingnodeviathatarc.Theow inandoutofanodestillbalancesasrequiredintraditionalsolutionmethods,but owchangingacrossarcscontradictsthebasisofclassicaltheoryinnetworkows. 11

PAGE 21

Alternativesolutionmethodsarerequired.Aminimumcostowproblemovernodes i 2N andarcs i;j 2A withterm K ij referredtoasthegainofarc i;j formulated in[11]as minimize x X i;j c ij x ij .2a subjectto X j : j;i 2A K ji x ji )]TJ/F25 11.9552 Tf 20.5 11.358 Td [(X j : i;j 2A x ij =0 8 i 2N .2b l ij x ij u ij 8 i;j 2A : .2c Equation.2bisaowconservationconstraintalternativetotheclassicalowconstraint.Thefactor K ij indicatestheowalongarc i;j ischanged,soarrivingowat j isdierentthanleavingowfrom i .Despitethisdierencefromso-calledordinary networks,dualitystillprovidesaguaranteeofoptimalitythroughcomplementary slacknessconditions.Thealgorithmsproposedin[11]useLagrangianrelaxationto deneadualascentmethodthatiterativelyupdatespricesfromthedualproblem whilemaintainingowsthatsatisfycomplementaryslacknessseekingprimalfeasibility.Dierencesbetweenthesealgorithmsforordinaryandgainnetworksinvolve theadditionalstepsrequiredinowaugmentationalongarcswithnon-unitgains andpriceupdatesfornodesconnectingthesearcs.Anadditionalstepforowaugmentationalongacyclecouldbeomittedingraphswithspecialstructuresuchas bipartitegraphsseeninassignmentproblems.Thesealgorithmsaredemonstratedin [11]ashavingperformanceimprovementoverclassicalmethodsformultipleclasses ofproblemssuchasassignment,transportationandtransshipmentproblems. Othermethodsexistforspecicproblemsinthespaceofgeneralizednetwork ows,likeonefacilityonecommodityOFOC.OFOCisanetworkowproblem withsinglecommodityowbetweenasinglesourceandsinkwhereowsmustbe madeinincrementsof C or1withcostsappliedintothebatch.Chopraetal. [19]describeapolynomialalgorithmforsolvinginstanceswherethenumberoffull 12

PAGE 22

batchows k ofsize C isbounded.Thisapproachdemonstratestheadvantagesof concentratedresearchonthequalitiesofspecialinstancesofdicultproblems. 2.5GeneralizedAssignmentProblem InthegeneralizedassignmentproblemGAP,agentswithspeciedcapacityare assignedtasksrequiringgenerallevelsofserviceinawaytominimizethecostof assignments.ThisproblemisNP-complete[27],unliketheassignmentproblemfrom networkowoptimization.Startingwiththesamenotationfromthetransportation problem,GAPhasthefurtherrestrictions x ij 2f 0 ; 1 g ,theowconstraintRHSis 1forall i andtheassignmentbyeach x ij hasaweight w i tobecountedagainst eachagent'scapacity s j .Theformulationisnearlyassimpleasthetransportation problem; minimize x X i;j c ij x ij .3a subjectto X i w i x ij s j 8 j .3b X j x ij =1 8 i .3c x ij 2f 0 ; 1 g8 i;j: .3d Solving,orapproximating,instancesofGAPcanincluderoundingmethods[33] andGPUimplementationsofrelatedproblems[17,15,24]suggestthatparallelization methodsofclass"singleinstruction,multipledata"SIMDcanbeusedtospeedup thesearchforsolutions. 2.6DispatchSchedulingModel Engau,MoatandDykpresentamultiperiod,multicriteriamodelforplanningload dispatchschedulesin[21].Thismodelisreferredtohereasthedispatchscheduling 13

PAGE 23

Table2.1:DSM:Index,ParameterandVariableNotation SetsDescriptionIndices A Setofsetsofbatterieswhichcanbecombinedintoasingleload A B Setofinventorylocationsortankbatteries i C Setofhaulers j D Setofdestinations k T Timeperiodindexset t Parameters Units c ijk Revenueperbarrelhauledfrom i 2 B to k 2 D by j 2 C dollars/bbl c k Costperbarrelforcapacityconstraintviolationsby k 2 D dollars/bbl d ijk Splitloadfeeperloadhauledfrom i 2 B to k 2 D by j 2 C dollars/load h jk Fullloadsizebyhauler j anddestination k bbl l i Minimumcontractualhaulfrombattery i 2 B typicallyzerobbl l j Minimumcontractualhaulusinghauler j 2 C bbl l k Minimumcontractualsupplytodestination k 2 D bbl q i t Productionrateforbattery i 2 B duringtime t 2 T bbl/period u i Maximumstoragecapacityofbattery i 2 B bbl u j t Maximumhaulcapacityby j 2 C during t 2 T bbl u k t Maximumcapacityat k 2 D during t 2 T bbl w i Parameterforbattery i 2 B forweightingshutinavoidanceunitless v 0 i Initialliquidinventoryatbattery i 2 B bbl Shutinvolumeweightformulticriteriaobjectiveunitless Hauledlocationweightformulticriteriaobjectiveunitless Variables Domain v i t Inventoryinbarrelsforbattery i 2 B attheendofperiod t 2 T R + x ijk t Barrelsofuidhauledfrom i 2 B to k 2 D by j 2 C during t 2 T R + y jk t Loadshauledbyhauler j 2 C todestination k 2 D during t 2 T Z + z ijk t Loadshauledfrom i 2 B to k 2 D by j 2 C during t 2 T Z + i t Timetoshutinbybattery i 2 B for t 2 T R + k Capacityconstraintoverowvariablefor k 2 D R + 14

PAGE 24

modelDSM.Itcoversthedecisionsonwhentohaulloadsfrombatteries z ijk t whentosplitloadsbetweenbatteriesatcost d ijk ,whichhaulertoassigntoeach loadandthedestinationforeachload.Themodeldecisionsaredrivenbymultiple objectivesincludingavoidanceofshutinsduetoexhaustingtankcapacity,meeting contractualobligations,protmaximizationandhaulvolumemaximization. TounderstandDSM,webeginwiththeconceptofhaulvolume.Themodelrepresentsthevolumeofliquideachbattery i 2 B sendstodestination k 2 D viahauler j 2 C duringperiod t 2 T .Thisvolumeisdenoted x ijk t .Thevolumeistransportedinbatchessized h jk BBL,whichcorrespondstothehauleranddestination. Themodelassumesthatallloadsarrivingatdestination k areatfullcapacity.The fullcapacityloadisenforcedwithconstraint.4f,whichusesthenumberofloads y jk t senttodestination k viahauler j duringperiod t .DSMusesthevariable z ijk t toindicatethenumberofloadstransportedfrombatterytodestinationbyaparticularhauler.Theseloadsarenotnecessarilyfull,butconstraint.4gensuresthatthe volumedoesnotexceedthecapacitycorrespondingtothedecisionof z ijk t .Along withthesevariablesandconstraintsthatensurehaulingisrepresentedconsistently, themechanicsofliquidowingintoandoutofthebatterymustbeconsistent.The inventoryonsiteatthebeginningoftheplanningperiod v 0 i isarequiredparameter inBBLincorporatedwithconstraint.4e.Theinventoryattheendofeachperiod v i t iscontrolledbyequation.4d.Thisequationdemonstratesconsistency withtheliquidowingintothebattery q i t anduidhauledawayfromthebattery P j 2 C P k 2 D x ijk t .Eachhauler j alsohaslimitsontheliquidhauledduringatime period.Theconstraint.4censuresthatthetotalliquidhauledbyahaulermeets tworequirements,thelowerlimit l j andupperlimit u j .Thelowerlimitcorresponds toapromisedamountorloadcountthatexistsincontractsoragreements,while theupperlimitcorrespondstoacapacitylimitorotheragreement.Thedestinations havesimilarlimitsenforcedbyconstraint.4bwith l k and u k beingparameters 15

PAGE 25

forlowerandupperlimitsonuidamountbydestination.Theconstraints.4h ensurethatthevariablesforloadcountarenon-negative.Themodelobjective.4a representstherevenuecorrespondingtoagivenhaulingschedule.Theparameter c ijk istherevenuereceivedforliquidsdeliveredfrombattery i todestination k viahauler j .Therevenuereceivedisreducedbythecostofthenumberofsplitloadsrequired. Eachsplitloadcosts d ijk dollars.Thisobjectiveisonlypartofthemulticriteriaobjectivesexplored.Seetheinitialmodelformulationbelow,withnotationsummarized inTable2.1. Maximize x;y;z;v X i 2 B X j 2 C X k 2 D X t 2 T c ijk + d ijk h jk x ijk t )]TJ/F19 11.9552 Tf 11.956 0 Td [(d ijk z ijk t .4a subjectto l j X i 2 B X k 2 D x ijk t u j t 8 j 2 C and t 2 T ;.4b l k X i 2 B X j 2 C x ijk t u k t 8 k 2 D and t 2 T ;.4c l i X j 2 C X k 2 D x ijk t = q i t )]TJ/F19 11.9552 Tf 11.956 0 Td [(v i t + v i t )]TJ/F15 11.9552 Tf 11.955 0 Td [(1 8 i 2 B;t 2 T ;.4d 0 v i t u i and v i = v 0 i 8 i 2 B and t 2 T ;.4e X i 2 A x ijk t = h jk y jk t 8 A 2A ;i;j;k;t ;.4f 0 x ijk t h jk z ijk t 8 i;j;k;t .4g y jk t 0and z ijk t 0integer : .4h Notetherevenuebybarrel c ijk termisusedintheprotcalculation,whichcanbe interchangedwithcostandnegatedtorepresentwaterhauls,whichdonotgenerate revenue.Thesystembalanceequation.4dusedinthismodeliscommontoall calculationsforchangeininventory.Itensuresthattheendinginventoryfromone periodatabatteryisthebeginninginventoryofthenextperiodforthebattery. Modelfeaturessuchassplitloaddecisionsandshutinavoidanceusingtimetosell outarekeytothedesign.Splitloadscanbecontrolledusingthecriteria.4fwhere allloadsareassumedtobefulland A isasetofbatteriesthatcanbeusedtogether 16

PAGE 26

tobuildsplitloads.Thisconstraintguaranteesatotalnumberofloadsishonored andthemaximalnumberofsplitscanbecontrolledbysettingarelationshipbetween P i 2 B z ijk and y jk where z ijk isthenumberofloadshauledfrombattery i todestination k byhauler j .Normally,sincecrudeoilloadsaresplitbetweenatmosttwobatteries, thiswouldbe P i 2 B z ijk 2 y jk ,butcanberestrictedevenfurthertolimitsplitstoa fractionoftotalhaulstoadestinationbyahauler.Thisisnotnecessary,butcould beusedtoindicateearlywhendesiredconstraintsareinfeasiblebecausethemodel objectiveguarantees z ijk takesonthesmallestvaluewhere x ijk h jk z ijk isfeasible [21].Theotherfeatureuniquetothismodelisthemethodusedtoavoidshutins. Usingacalculationforthetimetoshutinbybattery i t = u i )]TJ/F19 11.9552 Tf 11.955 0 Td [(v i t q i t .5 withatermintheobjective X i 2 B X t 2 T w i i t .6 whichusesaweightfactor w i bybatteryonthetimeuntilshutin i t ,themodel canbemadetoavoidshutins.Toincorporatethis,.4aisreplacedwith Maximize x;y;z;v X i 2 B X j 2 C X k 2 D X t 2 T c ijk + d ijk h jk x ijk t )]TJ/F19 11.9552 Tf 11.955 0 Td [(d ijk z ijk t ; .7a X i 2 B X j 2 C X k 2 D X t 2 T x ijk t ; and X i 2 B X t 2 T w i i t .7b toseektheadditionalobjectivesbeyondprotmaximization.Maximizingthesum ofthesevariables willresultinbatterieswithlowratesofproductionbeingfavored forhaulingsince u i isdesignedtomeetatleastonehauland q i t asymptotically approaches0aswellsage.Amorerealisticapproachtotheuseofthisvariable in planningistoseekamaxi-minof i t .Thisensuresthathighproducersmaintain alevelofinventorycapacitythatiscomfortabletooperations.Onedrawbackon theuseofthismodelisthedeterministicapproachtoproductionoperations.This 17

PAGE 27

formulationalsodoesnotallowforanyshutins.Situationswhereproductionvolume exceedshaulingcapacitybyenoughtoforcetheselectionofwellstoshutinwillfail asinfeasibleinthisformulation.Thisisrecognizedin[21],butnotaddressedina modelreformulation. 18

PAGE 28

CHAPTERIII STOCHASTICINVENTORYMANAGEMENT InthischapterwedescribeamodelthatbuildsonthepreviousworkfromDSM [21],whichisamodelforplanningdispatchschedulestomanageinventoryforpreventingshut-ins,keepingsalesconsistentandmaximizingrevenue.Thatmodeldoes notincorporateanystochasticcomponents,andthatworkiscontinuedinthisthesis asstochasticinventorymanagementSTIM.Uncertaintyinproductionandhauling callsforamodelthattakesastochasticapproachtoitsinputs.Usingthetechnique ofstochasticprogrammingexploredin[31],wewilldevelopanalternativemodelto handleuncertainty. 3.1StochasticModelDescription Wewillconsideramodelthatincorporatesrandomprocesses.Theentropyintransportationofgoodsarisesnaturallyfromsourcessuchastrac,equipmentandfacility exceptions.Theunknownofhydrocarbonextractionratesisonethatisspecictothis industry.Thismodelseekstomaintaineciencybyeliminatingproductiondowntimeandmaximizingthroughputwhilehonoringcapacityconstraints.Thedierences betweenSTIMandthedispatchschedulingmodelpresentedin[21]areintendedto makecomputationoflargeinstancesmoreaccessible,whilemaintainingthegoalof usefordecisionsintheneartermwithregardforlongtermplanning.Thenotation forSTIMherebuildsonTable2.1andaddssomestochastictermsgiveninTable3.1. 3.1.1StochasticAspects Theprocessofproducingandsellingcrudeoilandhaulingwaterhasstochasticcomponents,suchastherateofproduction,thetimefromplantowhenahaulisactually performedandthehaulingcapacity.Haulingcapacityasastochasticfeaturerefersto thenumberofloadsaparticularhaulercanperforminaperiodfromtheperspective 19

PAGE 29

Table3.1:STIM:Index,ParameterandVariableNotation SetDescriptionIndices B Setofinventorylocationsortankbatteries i C Setofhaulers j D Setofdestinations k T Timeperiodindexset t M Scenarios m ParametersUnits p m i Probabilityofscenario m forbattery i unitless q m i t Productionvolumeatbattery i forscenario m duringtime t bbl u i Maximumstoragecapacityofbattery i 2 B bbl u j t Maximumhaulcapacityby j 2 C during t 2 T bbl u k t Maximumcapacityat k 2 D during t 2 T bbl Shutinvolumeweightformulticriteriaobjectiveunitless Hauledlocationweightformulticriteriaobjectiveunitless VariablesDomain g m i t Shutinproductionbylocation i duringtime t ,scenario m R + v m i t Inventoryatbattery i followingtime t inscenario m R + x m ijk t Haulvolumefrombattery i duringtime t inscenario m R + z i t Indicatorforhaulingfromlocation i duringtime t f 0 ; 1 g 20

PAGE 30

oftheoperator.Thefocusofthismodelisontherandomnatureoftheproduction rateforeachbattery.Productiondeclinecurvesfromwellsaremodeledusingdeterministicdeclinecurves,butoperationsaresubjecttoexceptionswhicharestochastic innature,interruptingthenaturaldeclinecurveandcausingunexpecteddropsin productionvolume. Thischoiceoffeaturewillrepresentothercomponentsimplicitly.Theimpact ofahauloccurringlaterthananticipatedhasasimilareecttohigherproduction volumethananticipated.Withadiscretetimemodel,thedelaywouldneedtobe implementedbyhavingahauloccurduringaperiodfollowingtheplannedhaul. Whilethisisnotuncommoninpractice,itislimitinginthatwithoutshortening periodsandincreasingthenumberofdecisionvariables,onlylargedeviationscan berepresented.Havingmultiplepossibleproductionscenarioscanhelptorepresent delaysandearlyarrivals. 3.1.2ModelDevelopment Thedevelopmentofastochasticmodelrequiresthedenitionofscenariosrepresenting dierentoutcomes.Thescenariosforoutcomesaredenoted m 2 M .Eachbattery i andscenario m hasanassociatedprobabilityvalue p m i .Itisrequiredthatthe vector p i isstochastic,meaning P m 2 M p m i =1,ortheexpectedvaluecalculations areincorrect.Withthemodelchoiceforincorporatingrandomeventsbeingthe productionrate q m i t ,theinventoryattheendofeachperiod v m i t mustthen beadependentstochasticvariable.Duetotheinventorybeingstochastic,thehaul quantity x m ijk t shouldbevariedbetweenthescenarios m 2 M .Thisisdoneby addingconstraintswiththedierentscenariosofproductionrates,endinginventory andhaulquantity.Notethehaulingdecisionvariable z i t isstilldeterministic.Ifthis haulingdecisionwerestochastic,themodelwouldbeseparableonscenario,giving optimaloutputsforeachscenario,insteadofadecisionthatseeksoptimalexpected 21

PAGE 31

valueconsideringallscenarios.LooktoTable3.1fordetailsonthisnotation.Itis importanttonotethateventhoughproductionratesforanytwodistinctbatteries inaparticularscenario m mayhavethesameprobabilityvalue,thereisnoother relationshipimpliedbetweenthetwointhesamescenario.Thestochasticvariable q i t forbattery i isindependentofthevariable q j t forbattery j .Eachbattery i in eachscenariohasprobability p m i ,whichallowsforexibilityininstancespecication. Haulingcapacities u j aredeterministic,possiblyuniquetothehauler.Inventory capacities u i arexedforshorttermplanningpurposes,butcouldbemadevariable forlongtermdecisionmaking,fromtheoperator'spointofview. Planninghaulsusingastochasticmodelischallengingbecausethedynamicnatureoftheinventorymeansthatafullloadmightbeavailableinaparticularperiod inonescenario,butnotinadierentscenario,i.e. v m i t
PAGE 32

usefulresults. Oneofthesituationsnotaddressedwiththedispatchschedulingmodelisinfeasibilityofthemodelforunavoidableshutins.Toaddressthisweaddanewterm g m i t .Thisrepresentsproductionvolumethatexceedscapacityoflocaltanksso mustbeshutin.Itisparticulartoeachscenario,batteryandperiod.Positivevalues indicatethatwellsataparticularbattery i willneedtobeshutintoreducetotalproductionvolumeintime t tothefractionoftotalproduction q m i t )]TJ/F19 11.9552 Tf 10.855 0 Td [(g m i t =q m i t inthatscenario.Thismeansthat E P i 2 B g i t productionvolumecanbeexpected tobelostduringperiod t forforcedshutins. ThetermsdiscussedherearesummarizedintheTable3.1.Thecontrollingequationthatdrivesthedynamicnatureoftheproblemintroducedin.4eischangedto accommodatestochasticelementsandpreventinfeasibilityforcasesrequiringshutin productionto v m i t = v m i t )]TJ/F15 11.9552 Tf 11.955 0 Td [(1+ q m i t )]TJ/F25 11.9552 Tf 11.956 11.357 Td [(X j;k x m ijk t )]TJ/F19 11.9552 Tf 11.956 0 Td [(g m i t 8 m 2 M;i 2 B: .1 Twomoreconstraintsetsarerequiredtokeepthedecisionandhaulvolume variablesconsistent.Oneisrequiredtoensurethatnolocationhaspositivehaul volume x m ijk t ,ifthatlocationisnotdesignatedforhauling.Theconstraint X j 2 C X k 2 D x m ijk t )-222(M z i t 0 8 i 2 B;t 2 T;m 2 M .2 ensuresthissituationusingabig-Mmethodin M z i t where M isavaluesuciently largeenoughtobegreaterthananyothervaluecalculatedintheproblem.Itisonly satisedwith z i t =0when x m ijk t isalsozero.Sopositivehaulvolumesrequire thattheparticularbatterydecisionispositive.Theotherrequirementformaintaining consistencyisforvolumesbyhauleranddestination.Theconstraintforhaulervolume byperiodis X i 2 B X k 2 D x m ijk t u j t 8 j 2 C;t 2 T;m 2 M: .3 23

PAGE 33

Theeectofthisconstraintistokeeptotalhauledvolumebelowthevolumecapacity foraparticularhauler.Aswellasthehaulercapacityconstraints,thedestination constraintsaresimilar: X i 2 B X j 2 C x m ijk t u k t 8 k 2 D;t 2 T;m 2 M: .4 Eachinventoryvariableisconstrainedtobenon-negativeandlessthantheupper boundofcapacitybybattery u i Theobjectivefunctionforthestochasticmodelistomaximizetheexpected valueofthesumofhaulsfromaminimumnumberoflocationswhilemaintainingthe productionvolume.Theexpectedhaulvolumeis X i 2 B X m 2 M p m i X t 2 T X j 2 C X k 2 D x m ijk t : .5 Productionvolumeismaintainedbysubtracting P i 2 B P t 2 T P m 2 M g m i t fromthe sumofhauledvolumes,whichpenalizestheobjectiveusingafactor forallpositive shutinvolumes g m i t .Thenumberofhauledlocations P i 2 B P t 2 T z i t isimplicit tothenumberofsplitloadsrequiredtohaulthetotalvolume.Thisvaluepenalizes theobjectivebyafactor .Andnally,lowendinginventoryisagoalofoperations, sotheobjectiveispenalizedwiththeexpectedvalueofinventoryattheendofthe planninghorizon.Thefullobjectiveformulais maximize X i 2 B X t 2 T )]TJ/F19 11.9552 Tf 9.299 0 Td [(z i t + X m 2 M X j 2 C X k 2 D p m i x m ijk t )]TJ/F25 11.9552 Tf 14.379 11.357 Td [(X m 2 M p m i g m i t )]TJ/F25 11.9552 Tf 13.715 11.357 Td [(X m 2 M p m i ^ v m i .6 where^ v i isthenalperiodendinginventorybybattery i Inthedispatchschedulingmodel,theobjectivenaturallydrovetheloadcount toaminimalnumber.TheobjectiveinSTIMrequiresatermexplicitlyaddedfor directingloadcounttolowervalues.Thisisthesameasanalyzingthecardinalityof 24

PAGE 34

haulvolumevector x t ,dueto.2.Theterm z i t ismeanttokeepthenumber oflocationshauledtoaminimum,consideringthemaximizationofhaulvolume.The parameter canbeadjustedtochangethenumberoflocationschosenforhauling. Thereasoningbehindthetermistotradeotheminimumloadsizewiththechoice ofhauling.Inthedispatchschedulingmodel,theterm z ijk t isinthedomainof integers.Therestrictionof z ijk t toabinaryvariableinSTIMismeanttobeasmall simplicationstepforcomputationalcomplexity. Thefullmodelispresentedhere. maximize v;g;x;z X i 2 B X t 2 T )]TJ/F19 11.9552 Tf 9.299 0 Td [(z i t + X m 2 M X j 2 C X k 2 D p m i x m ijk t )]TJ/F25 11.9552 Tf 14.379 11.358 Td [(X m 2 M p m i g m i t )]TJ/F25 11.9552 Tf 13.716 11.357 Td [(X m 2 M p m i ^ v m i subjectto X i 2 B X j 2 C x m ijk t u k t 8 k 2 D;t 2 T;m 2 M .7a X i 2 B X k 2 D x m ijk t u j t 8 j 2 C;t 2 T;m 2 M .7b X j 2 C X k 2 D x m ijk t )-222(M z i t 0 8 i 2 B;t 2 T;m 2 M .7c v m i t = v m i t )]TJ/F15 11.9552 Tf 11.956 0 Td [(1+ q m i t )]TJ/F25 11.9552 Tf 11.955 11.358 Td [(X j;k x m ijk t )]TJ/F19 11.9552 Tf 11.955 0 Td [(g m i t 8 i 2 B;m 2 M .7d 0 v m i t u i 8 i 2 B;m 2 M .7e ^ v m i = v i max f t 2 T g 8 i 2 B;m 2 M .7f x m ijk t 0 ;z i t 2f 0 ; 1 g ; g m i t 0 8 i 2 B;j 2 C;k 2 D;t 2 T;m 2 M Alsonotetheterm g m i t .Thisdecreasestheobjectivewhendecisionsforce productionshutin.Theparameter isusedtoindicatetheimportanceofpreventing shutinofproductionvolumes.Thisformulationdoesnotsolvethesameproblemas DSM,whichconsidersthetradeobetweencostofhaulingsplitloadsandinventory 25

PAGE 35

management.Thisstochasticmodelseekstomaximizehaulingandwillresultina solutionthatneedstobeinterpretedtondsplitloadsandpossiblyformulatedin anothermodeltodecidehowtohaulthesplits.Thisisanoperationthatisrequiredin thedispatchschedulingmodelaswell.Thesesubsequentproblemsarenotnecessarily feasible,givenupperboundsonloadcounts.Supposethateitherofthesemodelsgives asolutionthatincludesfullloadsforallbuthaulsfromfourbatteries,withthesum ofbarrelstobehauledequaltothetotalnumberofbarrelsinthemaximumnumber ofloadsbyhaulerwhichisimpliedinthemodels,butnotexplicitlyparameterized. Thennoadditionalloadswouldbeavailable.Thehaulvolumesfromthesefour batteriesare70 ; 80 ; 90 ; 160BBL,respectively.Thetotaloftheseloadsis400BBL, whichistwofullloadsfor h jk =200.Thefoursourcescan'tbecombinedinwhole tomaketwofullloads.Thisisascenariothatisafeasibleresultfromthedispatch schedulingmodelandSTIM.Theapproachproposedin[21]tosolveabin-packing problemfollowingasolutioninthisrststagewoulddeliveranadditionalsolution thathasaminimumnumberofloadsrequired,ifoptimal.Althoughthisnumberis notnecessarilylessthanorequaltotheimpliednumberofloadsimpliedby u j =h jk theuncertaintyinhaulingcapacitymentionedearlierinthissectionmeansexceeding theloadcountcouldbeinconsequential. 3.1.3ModelOptions Featuresfromthedispatchschedulingmodel[21]canbeincorporatedintothismodel asneeded.Thedecisiononwhichbatteriestohaulcanbereplacedbyavariable denotingtheexactnumberofloadstohaulbetweenlocationsanddestinations.The maximizationoftimetoshut-inormaximizationoftimetosell-outcouldbebrought intotheobjectivetoprioritizedierentgoalsininventorymanagement.Usingthe basicapproachofmodelingtheeectofparticulardecisionsunderdierentscenarios isnotmodiedbyaddingmultiplecriteriaorchangingtoamorecomplexdecision. 26

PAGE 36

3.2SolutionMethods Forlargecases,STIMisdiculttosolve.For30periodsand500batteries,thisresults in15,000binaryvariables.ApplyingacleverMILPsolvermaytrimthesearchtree, butthesearchspaceissolargethatitwilltakedrasticcutstomakethesearch successful.Reducingthenumberofvariablesisonewaytoreducethesearchspace directly.Reducingthenumberofdecisionvariablesinaninstanceisdonehereusing themethodofmodelpredictivecontrolMPC,discussedinSection3.2.3. 3.2.1MILPSolver TheuseofaMILPsolveronthismodelwillstillgivefeedback,eveninlargeinstances. ThesolveringeneralwillgiveatleastaMIPgapandunderstandingofthesizeof thetreeremainingtosearchafteraspeciedruntime.Thereisthepossibilityofnot ndingasolutioninthatruntime,whichwouldbeaninniteMIPgap,butitwillat leastshowhowmuchofthesearchtreeisexamined.Itisnotlikelythatthesearch treewillhavebeentrimmedbyasignicantproportionforgenerallargeinstances, buttheinformationwillstillbeavailablefromabranchandboundorbranchand cutglobaloptimizationmethod. 3.2.2SolutionHeuristic InDSM,theexactnumberofloadsperbatteryisdeterminedbyacombinationof integerdecisionvariables.InSTIM,thedecisionvariablesareinthebatteriesto haulfromandquantityofproducttohaulperbattery,hauleranddestination.This changemakesSTIMacardinalityproblem.Theminimizationof P i 2 B P t 2 T z i t intheobjectiveinturnreducesthecardinalityof x ,duetotheconstraint.7c forcing x ijk t to0when z i t =0.Anheuristicexistsforcardinalityproblems thatinvolvesrepeatedsolutionsofarelaxedversionoftheoriginalproblem.This 27

PAGE 37

relaxationallowsthebinaryvariabletotakeanyvaluein[0 ; 1],subjecttoother existingconstraints.Thisisessentiallyaroundingheuristiconthedecisiontohaul, restrictingthecardinalityof x .Itisseenin[13].Tobeginthisheuristic,aset ofcardinalitycostsmustbeselected.Thesecostsaredierentvaluestoapplyin themodelpenalizingthecardinalityof x .Hereisthestep-by-stepexecutionofthe heuristic. 1.Setinitialobjectiveto 2.Chooseacuto 2 [0 ; 1. 3.Buildthemodelinstanceallowingfractionalsolutions. 4.Solvethisoriginalinstancewithcardinalitycost ,FAILifinfeasible. 5.Addconstraintstotheinstance z i t =0foreachvariable z i t < and z i t =1 foreach z i t 6.Solvethenewinstancefortheremainingvariables v g and x 7.Ifobjectiveishigherthanthecurrentobjectivevalue,thensaveobjectiveand solution. 8.RepeatStep3-7withnext fromasetofcardinalitycosts. Thismethodworksbyxingthecardinalityateachiterationaftertherelaxed solutionisfound.Eachiterationsbeginswiththeoriginalproblem,butanewcost forthecardinalityofthesolution.Thismeansthatthemethoddoesnotnecessarily increaseordecreaseinobjectivevalueandthelowestobjectiveshouldberecorded alongthewayandreturnedasthesolutionafterthenaliterationiscomplete.This algorithmisnotaglobalsearchtechnique.Itisdierentinatleastthesolutionsetis nottrimmedfromoneiterationtoanother;atthestartofeachiteration,thesolution 28

PAGE 38

setisthesame.Abinaryintegerprogrammingtechniquewouldincludebranchingto searchthesolutionspace. 3.2.3ModelPredictiveControl Tosolvethismodeloverasucientlylongplanninghorizon,weusethemethod ofmodelpredictivecontrolMPCbyapplyingaseriesofshiftedtimeperiods.The MPCtechniquecanbeappliedinconjunctionwiththesolutionheuristicfromSection 3.2.2.Thistechniquesetsaslidingtimehorizonateachiteration,movingforward oneperiodandusingtheresultsoftheexecutionofthepreviousdecisions.This techniquecanbeappliedinaplanningtaskusingsimulationtodeterminetheresults ofthedecisionsorinaliveenvironment,gettingfeedbackfromthesystem. Inanitetimehorizonformulation,onechoiceistosetaconstraintforcingthe stateofthesystemtoaspecicvalue[14].Thishastheeectofpreventingthe solutioninthecurrenthorizonfromputtingthemodelinapredicamentinthenext horizon.Applyingacostintheobjectiveinsteadofahardconstraintisanalternative methodforseekingthesamegoal.Thisisappliedinthemodelformulationby decreasingtheobjectivebytheexpectedvalueofthetotalendinginventory.Even withoutMPC,thisisadesiredgoalofthemodel. 3.3Simulation SimulationisusedwithSTIMtodriveMPCandtovalidateresults.WithMPC, theslidingtimehorizoncallsforinputsresultingfromthepreviousrun,whichare theoutcomesofexecutingdecisions.Onemethodusedtogeneratethesevaluesis simulation.Thistypeofsimulationexecutesdecisionsandcomesupwithmeanvalues forparametersneeded.Simulationcanalsohelpvalidateresultsbyapplyingdecisions generatedbythemodelandcomparingtheresultingvalues,suchashaulvolumeor 29

PAGE 39

shutinvolumeforSTIM,totheresultsprescribedbytheoptimalsolution. Asinthedesignofthemodel,thedesiredstochasticcomponentsforsimulationaretheproductionratesforthebatteries.Asimulationwouldbemoreeasily designedtohavetheindividualloadpickuptimesbestochastic,arandomvariable oftimefollowingtheplannedhaultowhenaloadispickedup.Oritcouldbea randomvariableforthenumberofloadspickedupinaperiodfordiscretetimesimulations.Alternatively,theprocessofsimulatingtheinventorydynamicsdoesnot includeanyrepresentationofthetrucks,tanksorevenindividualloads.Theentire batterycapacityistreatedasasingletank. 3.4ComputationalResults ThedataforSTIMexperimentationisnotmodiedmuchfromthegenericdata discussedinAppendixA.Asthemodelisdesignedtogeneratehauldecisions,the inventory,productionrateandconstraintdataonlyneedstobetranslatedintothe formatrequiredbythemodel.Themodelinstanceswillfeatureadditionalvariables usedtotrackobjective,haulvolumeandlostproductionaggregates.Thesevalues werebeusedtogaugeanddemonstratemodelperformance. 3.4.1SolutionstoSTIM TwosolutionmethodshavebeenappliedtoSTIM. TherstisthepureMILPsolver.Thismethodwasexecutedprimarilytoshow runningtimes.TheMIPgapplotinFigure3.1reectsthemaximizationgoal.All pointsarenegative,showingawiderangeofMIPgaps.Someoftheinstancesbelow 5%gapmaybedecentperformance,consideringtheplanningapplication.Ascan beseeninFigure3.2,noinstancesranunderthecutotime.Onlyafractionofthe instancescompletedduetoanexplosionintimeseenwithoverheadtakenforbuilding 30

PAGE 40

Figure3.1:STIM:MIPGapvs.OptimalLPRObjective thedatastructuresrequiredforsolvingtheseinstances. ThesecondmethodusestheheuristicalgorithmandMPC.Inthismethod,the plotoftheMIPgapisnotgiven.Theexecutiontimeofthismethodisquicker thanthepureimplementation,showninFigure3.3.SinceallinstancesofSTIM wereterminatedbeforereachingoptimalityconditions,MPCisusedexclusivelyfor comparisonwithsimulatedresults. 3.4.2Simulation SimulationappliedtoSTIMchecksthevalidityofthesolutionfoundwiththediscrete scenariosinamodelwithcontinuousrandomvariables.Thesimulationprocessisas follows.Theprocessisrunforeachinstanceandeachinstancegenerates1000samples. 1.CapturetheoutputofSTIMforaninstance. 2.Generatebaseconditions,includingbeginninginventory v i andmeanpro31

PAGE 41

Figure3.2:STIM:WallTimeHistogram Figure3.3:MPC:WallTimeHistogram 32

PAGE 42

ductionrate q i foreachbattery i 3.Foreachperiod t ,dothefollowing Getanewrateofproductionforeachbatteryusing q i generatedfrom thenormaldistributionforeachbattery i asexplainedinfurtherdetailin AppendixA. v i t )]TJ/F15 11.9552 Tf 11.955 0 Td [(1+ q i t aretheavailablevolumetohaul. Applythehaulingdecisiontotheavailablevolumeandhaul x i t upto theleastoftheplannedamountandavailablevolume. If v i t = v i t )]TJ/F15 11.9552 Tf 12.584 0 Td [(1+ q i t )]TJ/F19 11.9552 Tf 12.584 0 Td [(x i t isgreaterthantheinventorycap,the dierenceistheshutinvolume. Repeatfortheprescribednumberofperiods. Thisprocessgeneratesanumberofsampleseachofhaulvolumeperperiod,shut involumeperperiod,dryloadsperperiodandtheendinginventoryfortheentire planninghorizon.Thesedatahavebeencomparedwithplannedoutput. 3.4.3Results TheexecutionoftheMPCmethodareshownhere.Asampleoftheendinginventory comparedbetweenthesimulatedandplannedresultsisshowninFigure3.4.This plotisonly10%ofthecases,butitshowshowthesimulatedresultscloselymatched theplannedresults.ThescatterplotinFigure3.5showsthesimulatedvsplanned endinginventorywhichshowsnearlyidenticalresults.Inthehaulvolumescatter plotFigure3.6,eachperiod'ssimulatedhaulvolumeisplottedagainsttheplanned haulvolume.Ther-squaredvalueofthisdatais0.99997.AlsoinFigures3.7and3.8, noticeadierentpattern.TheplotinFigure3.8includestheminimum,maximum andmeanofthesimulatedshut-involumesbycaseareplottedagainsttheplanned 33

PAGE 43

Figure3.4:MPC:EndingInventorySimulationSample shut-involumes.Althoughthehaulvolumeandendinginventoryresultstrackclosely betweenplanandsimulation,thesimulatedshutinvolumesarelessconsistenttoplan andarealsolessthantheplan,demonstratedinthetwoplots. 3.4.4AnalysisofResults Theresultshaveshownwhilethereisstillworktobedonetoanalyzescenarioswith highervarianceintheproductionrate,theplansandsimulationshaveagreedso wellthatitcanvalidatethecorrectnessoftheimplementations.Theaccuracyof thesimulatedresultstoplannedresultsisnottobetakenasaforecastofhowthe modelwilldoinactualbusinessuse.Itneedstobeshownthatthemodelwillstill generategoodplansinsystemswithmoreentropy.Theresultsfoundbysimulationof modeldecisionsareclosetooptimizationresultsforhaulvolumeandendinginventory model,butlesssoforshutinvolume.Thisshowsthatthemodelcanbeexpected 34

PAGE 44

Figure3.5:MPC:EndingInventorySimulationScatterPlot Figure3.6:MPC:HaulVolumeSimulationScatterPlot 35

PAGE 45

Figure3.7:MPC:SIVolumeSimulationSample Figure3.8:MPC:SIVolumeSimulationScatterPlot 36

PAGE 46

toperformwellwhenthevarianceappliedinplanningisthesamevarianceusedin simulationoraccuratetothereal-worldscenario,butstillrequiressomeresearchon shutinvolumes. Onesimulationdecisionthatcouldberethoughtisthedecisiontolimithaul volumebytheplannedvolume.Theresultissimulatedexpectedhaulvolume isslightlylowerthanplannedhaulvolumebecausesimulatedhaulvolumemaybe reducedbyexperiencingareductioninavailablevolumetohaulwhenproduction volumeislower.Asimplealternativetothisapproachistoallowhaulvolumein thesimulationtomeetavailablevolume.Thiswouldnotbeacceptablewithout dealingwithhowloadswillbepickedupfromeachbatteryinthesimulation.When loadsarenotfull,thenthereisthepossibilityofapartialloadfromonelocation beingtoppedowiththevolumefromanotherlocation.Inaplanningscenario, thisinconsistencyiscertainlysmallenoughtostillprovidegoodinformation.Inan operationalscenario,themodelsareperhapsgoodenoughforbusinessuse.More testingisrequired.Thistestingcouldbeperformedinarealworldsituationwhere STIMpresentsasolutionandprediction.Thenlegacydecisionprocesseswouldbe usedtocarryoutoperations.Usingsimulatedresultsfrombothdecisionsets,analysis wouldbeperformedtocomparehowresultswouldbedierentiftheplanhadbeen used.Thisisonewaytotobuildcondenceinthemodel. 37

PAGE 47

CHAPTERIV BATCHEDTRANSPORTATIONMANAGEMENT ThemodelBatchedTransportationManagementBTMisamodelthatrepresentsdecisionsrelatingtothetransportationofliquidsfrombatteriestodestinations. ThisdecisionisjustonepartofDSM,assumingthatotherdecisionswhichleadup tothedecisiononwheretohaulliquidshavebeenmade.Thoseotherdecisionscan bemadebyapplyingSTIMorsomeotherbusinessprocess.Thischapterpresents BTM,whichisamodelthatseeksoptimaldecisionsfordeliveringliquidvolumes todestinationsconsideringcostdollars,time,distance,etc.Thereareseveralkey elementstothemodeldrivenbytheneedsofdecisionmakers.Theobjectiveofthe modelistohonordestinationcapacitiesandminimizecostofhauling,whichinareferenceimplementationisthetraveldistancebetweensourceinventorylocationand destination.Thistranslatestotimespentontheactivityandminimizationofthis costleadstoeciencygains.Inmanycases,theincentivetoincreaseeciencyand reducecostsissecondarytotheabilityforamethodtoprovideasatiscinganswer [26].Goalprogramming[25]isatoolfrommulticriteriaoptimizationthatseeksto ndsatiscingsolutionstodicultproblems.Thismodelingmethodassignsprioritiestocriteriainaproblemusingcostcoecientsandminimizesthecostofviolating thesecriteria.Thispertainstothemodelinthischapterthroughtheobjectiveof minimizingcostsandtherequirementsfoundinthecapacityconstraintsforhaul destinationsthatprovetobeachallengetomeet. 4.1ModelFormulation Themathematicalmodelrepresentingthedecisiononhowmuchtohaultoeachdestinationisgivenbelow.SeeTable4.1forthedenitionsofelementsinthemodel.We willformulatethismodelwithouttimeperiodvariables.Thisisbecausethetransportationdecisionisnotimpactedbybatteryinventory,whichisthecomplicating 38

PAGE 48

aspectofDSMandSTIMleadingtomultipleperiodanalysis.Thismodelseeksthe minimumtotalcost.1awhere c ijk isacostformakingacertainhaulingassignment z ijk .Eachassignment z ijk resultsinaquantity x ijk byconstraint.1cthat representstheuseofcapacitywhoseupperboundis u k inconstraint.1b.The model.1a-4.1discalledbatchedtransportationmanagementBTM. minimize x;z X i 2 B X j 2 C X k 2 D c ijk z ijk .1a subjectto X i 2 B X j 2 C x ijk u k 8 k 2 D .1b X j 2 C X k 2 D x ijk = r i s i 8 i 2 B .1c z ijk s i = x ijk 8 i 2 B;j 2 C;k 2 D .1d z ijk integer 8 i 2 B;j 2 C;k 2 D .1e Droppingtheintegralityconstraintandonlyrequiringpositivevaluesfor z ijk produces theLPrelaxationofBTM,whichwewillrefertoasLPR. 4.2UnderlyingNetworkModel BTMisalinearprogrammingformulationwithintegerowbalanceconstraintsthat cannotberepresentedusingnetworkows,exceptinspecialcases.Thenecessityof twodecisionvariablesfordierentaspectsofahaulcomesfromhavingsupplynodes thatdonotllcapacityofbatcheddeliveriesalone,i.e.partialloads.BTMdoes notattempttollindividualdeliveriestocapacities,whichislefttodecisionmakers specicallythehaulerfollowingdestinationassignment. BTMcanberepresentedasageneralizedassignmentproblem.Thisproblemis knowntobeNP-complete.ItisNP-completetoevenndafeasiblesolution[27]. SinceGAPisageneralizednetworkowmodel,thispointstotheunderlyingnature ofBTMasanetworkowproblem.Theassignmentproblemhastasksofthesame 39

PAGE 49

Table4.1:BTM:Index,ParameterandVariableNotation SetsDescriptionIndices B Setofinventorylocationsortankbatteries i C Setofhaulers j D Setofdestinations k T Timeperiodindexset t ParametersUnits c ijk Costperloadhauledfrom i 2 B to k 2 D by j 2 C dollars/load c k Costperbarrelforcapacityconstraintviolationsby k 2 D dollars/bbl l k Minimumcontractualsupplytodestination k 2 D bbl r i Totalhaulsrequiredfrom i 2 B truckload s i Volumeofuidperhaulfrom i 2 B bbl u k Maximumcapacityat k 2 D during t 2 T bbl VariablesDomain x ijk Barrelsofuidhauledfrom i 2 B to k 2 D by j 2 C R + z ijk Loadshauledfrom i 2 B to k 2 D by j 2 C Z + k Capacityconstraintoverowvariablefor k 2 D R + 40

PAGE 50

weightanddierentcosts,perassignmentdecision.Althoughtheassignmentproblemcanbesolvedusingclassicalnetworkowmethods,thegeneralizedassignment problemcannot.ThemanipulationsrequiredtoformulateBTMasaGAPinstance revealthenetworkstructure. ItcanbeseenintheformulationofBTMthatthevariable x ijk istieddirectlyto z ijk .Thevariable x ijk mustbeamultipleof s i .Replacing.1bwith X i 2 B X j 2 C s i z ijk u k 8 k 2 D .2 and.1cwith X j 2 C X k 2 D z ijk = r i 8 i 2 B .3 allowstheeliminationof.1d.Nowthisstructureisclosertoagraphicrepresentationlikingatransportationproblem,but.2stillhasnon-graphicaspects. Arbitraryvaluesof s i meanBTMisageneralizednetworkowandprecludetheuse ofclassicmincostowmethodstosolvetheproblem.However, s i isboundedfrom above.Letthisupperboundbe^ s .Supposethat s i =^ s forallbatteries i 2 B ,then .2canberewritten X i 2 B X j 2 C z ijk u k = ^ s 8 k 2 D .4 withoutchangingtheproblem.Themodel.1a,4.3,4.4isequivalenttoatransportationproblemandcanbesolvedusingmincostowmethods.Callthismodel FBTM.ConsiderFBTMwithonehauler.Then.1aisequivalentto.5a,.5b isequivalentto.1band.5cisequivalentto.1c.Thenetworkchangesnecessarytoincludemultiplehaulersaretoaddnodesforeachhauleranddestination combinationwitharcsofcost c ijk betweenthebatteryandnewnode.Theoriginal destinationnodessumtheowacrossarcscomingfromthenewnodes.Thisnetwork owmodelisslightlymorecomplexthanatransportationproblematranshipment problem,butstillhasapolynomialalgorithm.Thefullformulationofthisnew 41

PAGE 51

model, minimize z X i 2 B X j 2 C X k 2 D c ijk z ijk .5a X i 2 B X j 2 C z ijk u k = ^ s 8 k 2 D .5b X j 2 C X k 2 D z ijk = r i 8 i 2 B .5c z ijk integer 8 i 2 B;j 2 C;k 2 D .5d onlysolvestheoriginalBTMproblemiftheloadsizesallequal^ s 4.3ComplexityofBTM BTMisan NP -completeproblem.ThiscanbeseenbydemonstratinghowGAP isaspecialcaseofBTM.ConsiderBTMwithallbatterieshavingasingleloadto deliver,i.e. r i =1forallbatteries i 2 B ,andonlyasinglehauler.Nowitcanbe seenthat.3isequivalentto.3c.Theconstraint.2needsnomodication tobeequivalentto.3b.Withtheassumptionofloadcountandsinglehauler, theobjective.1aisequivalentto.3a.Undertheseassumptions,theobjective andconstraintsbetweenthetwoproblemsareequivalent,soGAPisaspecialcaseof BTM.GAPis NP -complete,thereforeBTMisalso NP -complete. 4.4ApproximizationMethodUsingNetworkFlowRelaxation Inthissectionwewillexploreanapproximationmethodthathasincreasingaccuracy asloadsizes s i becomecloserto^ s .Recallthatwhen s i =^ s forallbatteries i 2 B BTMisequivalenttoFBTMandcanbesolvedwithnetworkowmethods.The objectiveofthissectionistocreateamodelthatcanbesolvedusingthesamemethods asFBTM,butdoesnoteliminatesolutionsfromthefeasibleregionofBTM. 42

PAGE 52

4.4.1ModelInstanceManipulation ThedemonstrationthatspecicinstancesofBTMaresolvableusingmincostow methodssuggeststhatmanipulationsofotherBTMinstancescouldbesolvedinthe sameway.Thesesolutionsareapproximationsandmanipulationsthataretoosimple couldcreateinfeasibleapproximationswheretheoriginalproblemisfeasible.For instance,takethesimplecasewith5loads,oneat150BBLandfourat100BBL withtotaldeliverycapacityof600BBL.Assumingthateveryloadis150BBLwould allowtheuseofmincostowmethods,butthetotalcapacitywouldbelessthanthe totaldeliveryrequirementof750BBL.Themanipulationin.6andtheresulting newmodelin.7createsanewproblemthatexpandsthefeasibilityregioninBTM andallowsfortheimplementationofasolverwithnetworkowmethods. Sinceitisassumedthatnotallloadstobehauledarethesamesize,itmakes sensetodenethedierencebetweenthelargestloadsize^ s andloadsizebybattery s i forall i 2 B .Let i =^ s )]TJ/F19 11.9552 Tf 12.944 0 Td [(s i .Usingthisvalue,let'sbuildanewconstraint representingcapacitylimitsbydestination.Beginwith.2 X i 2 B X j 2 C s i z ijk u k 8 k 2 D andtaketheinequality X i 2 B X j 2 C i z ijk X i 2 B i r i 8 k 2 D whichistorepresentthebalancebetweenloadsdeliveredandloadsrequired.Then sumtheseinequalitiestogetherresultingin X i 2 B X j 2 C s i + i z ijk u k + X i 2 B i r i 8 k 2 D: Since^ s = s i + i 8 i 2 B X i 2 B X j 2 C ^ sz ijk u k + X i 2 B i r i 8 k 2 D: 43

PAGE 53

Finally,dividingbothsidesby^ s producesanewconstraintforourapproximation model X i 2 B X j 2 C z ijk u k + P i 2 B i r i ^ s 8 k 2 D: .6 LetuscallthisnewmodelnetworkowrelaxationNFRafterthesolution approachthatisbeingtargeted.WenowhavetheentireformulationofNFR: minimize z X i 2 B X j 2 C X k 2 D c ijk z ijk .7a X j 2 C X k 2 D z ijk = r i 8 i 2 B .7b X i 2 B X j 2 C z ijk u k + P i 2 B i r i ^ s 8 k 2 D .7c z ijk integer 8 i 2 B;j 2 C;k 2 D .7d NFRisatransportationmodellikeFBTM,soNFRissolvableinpolynomialtime. SincethecapacityconstraintforNFRisrelaxedinaspecicway,NFRhasthe propertythatifasolutionisfeasibleinBTM,itisalsofeasibleinNFR.SeeFigure 4.1forthesetrelationshipforfeasiblesolutionstoNFR,BTMandtheLPrelaxation ofBTMLPR.NFRdoesnotallowfractionalsolutionsandneitherdoesBTM.LPR hasfractionalandintegersolutions.Theintersectionofthesetsoffeasiblesolutions forLPRandNFRisthesetoffeasiblesolutionsforBTM.ThisisduetoLPRhaving nointegersolutionsthatarenotalsosolutionstoBTM. InadditiontotherelationshipbetweenLPRandNFR,thereareresultsforNFR thatstructureitsuseinsearchingforsolutionstoBTM.Referto.1and.7for detailsregardingtheseproofs. Theorem4.4.1 GivenaprobleminstanceBTMandanyfeasiblesolution z forthat instance, z isafeasiblesolutionforNFR. 44

PAGE 54

NFR LPR BTM Figure4.1:BTM:FeasibleRegionDiagram Proof: Let z beafeasiblesolutionforBTM.Takethecapacityconstraint4.2 X i 2 B s i z ijk u k 8 k 2 D andthemanipulationsproducing.7cshowthat z isstillfeasibleregardingthe capacityconstraintsforeach k 2 D .Therefore z isfeasibleforNFR. Since^ s isanupperboundfor s i 8 i 2 B ,anyassignmentofloadstodestinations thathonorscapacityconstraints u k foreach k inBTMwillalsohonorallcapacity constraintsinNFR.Allotherconstraintsremainunchanged,so z isfeasibleforNFR. Theorem4.4.1allowsanexplorationoftheentiresolutionspaceofBTMusingnetwork owmethods.Although,itisnottruethatfeasiblesolutionsevenoptimalonesto NFRarefeasibleinBTM.Aneasywaytoknowwhencertaininstancesarenotfeasible isusingLPR.Thefollowingsimpletheoremdemonstratesthiseasycheck. Theorem4.4.2 GivenaprobleminstanceBTM,LPRandarelatedinstanceNFR, ifthereexistsasolution z forNFRand P i 2 B P j 2 C P k 2 D c ijk z ijk islessthanthe optimalobjectiveforLPR,then z isinfeasibleinBTM. Proof: Bycontradiction.Let z beanoptimalsolutiontoNFRforsomeBTM 45

PAGE 55

instance.Let z 1 beanoptimalsolutionforLPR.Supposethat c T z
PAGE 56

optimalpoint z forFBTM.Then z isfeasibleforBTM. Proof: Let z beanoptimalsolutiontoFBTM.Then,fromthederivation producing.6wehave, X i 2 B X j 2 C s i z ijk X i 2 B X j 2 C ^ sz ijk u k 8 k 2 D; .8 whichisthecapacityconstraintforBTM.Therefore,anoptimalsolutionforFBTM isafeasiblepointforBTM. ThisisnotonlytrueforanoptimalsolutiontoFBTM,butanysolutiontoFBTM. AnoptimalsolutiontoFBTMisnearlyaseasytondasafeasibleone,thankstothe networkowalgorithmsthatcansolveit,soitismoreinterestingthananyfeasible solution.AnoptimalsolutiontoFBTMisnotnecessarilyoptimalinBTM,butis useful.ItmaybethatasolutiontoFBTMisgoodenoughforagiveninstanceor maybeastartingpointforadecisionmakertochoosetoviolatesomeconstraintsto comeupwithabetterobjectivevalue. 4.4.2BetterBounds Theorem4.4.1showsthatthereisasimplemanipulationoftheBTMmodeltocreate anetworkowmodelthatcontainsallfeasiblesolutionsforBTM.Themanipulation oftheconstraintsisnotveryinvolvedandverymuchaliberalmanipulation.There isamoreconservativeapproachtomanipulatingthecapacityconstraintsthatstill allowsforallBTMsolutionstobefeasibleinNFR. Thismethodismoreinvolvedandrequiresasinglesort,followedby m binary searchoperationsofsize n tocalculatetherelaxationofcapacityateachdestination.Thecomplexityofthemanipulationis O j D j log j B j sincebinarysearchhas 47

PAGE 57

complexity O log j B j .Themethodfollowsvesteps 1.Sortthelistofloadsizes. 2.Choosethenextdestination k 3.Usebinarysearchontheloadstondthelargestnumber ofloadsthatcan feasiblybeassignedtothedestinationcallthisset G k 4.Add P i 2 G k i to u k forthenewupperlimitfordestination k ,assigningto u k 5.Ifdestinationsremain,repeatsteps2-5. Thisalgorithmisguaranteedtocreateanewproblemthathasallthefeasiblesolutions totheoriginalproblem.Thisisduetotheconditionthatanysolution z forBTM willhaveasetofloadsforsome k 2 D ,callit H ,whereforeach i 2 H ,forsome j 2 Cz ijk =1, P i 2 H P j 2 C z ijk j G k j8 k 2 D .Sincethe j G k j smallestloadsizesfor i 2 B willhave P i 2 H i P i 2 G k i andcanthenfeasiblytatleastasmanyloads assignedtoaparticulardestination. Considerasimpleexampletodemonstratethealgorithm,considerasetofloadsof sizes s =[70 ; 90 ; 110], =[40 ; 20 ; 0]andadestination k ofsize u k =200.Thebinary searchtechniquewillpickloadsizes[70 ; 90]and =[40 ; 20].Thentherelaxedvolume is200+60=260BBL.Noteanyallthreecannotbeassignedtothisdestinationsince 70+90+110=270 > 200,twoloadscombinedislessthantheoriginalcapacityof 200BBLand2^ s =220,soallassignmentsfeasiblebeforearefeasiblefollowingthe capacityrelaxationandxingallloadsizesto^ s 48

PAGE 58

4.5GoalProgrammingModel BTMcanbeeasilyextendedtoagoalprogrammingmodel.Theobjective.1aalreadyseekstheminimizationofcost.Theconstraints.1cand.1bcanberelaxed andbroughtintotheobjectiveforagoalprogrammingapproach.Theowconservationconstraint.1ccouldbechangedfromequalitytoinequalityandpenalized fornotdeliveringtheamountofow r i s i ,butitisassumedthatthisrequirementis alreadyrmsinceitlikelycomesfromsolvingamodelsuchasSTIM,whichconsiders productionshutins.Relaxingtheconstraint.1bbyaddinganoverowterm k totherighthandsideallowsforsolutionstoexceedthecapacityconstraints u k ata priceof c k perBBLindecreasedobjectivevalue,bothtermsavailableforreference inTable4.1.TheresultingmodeliscalledGBTMandisformulated minimize x;z X i 2 B X j 2 C X k 2 D c ijk z ijk + X k 2 D c k k .9a subjectto X i 2 B X j 2 C x ijk u k + k 8 k 2 D .9b X j 2 C X k 2 D x ijk = r i s i 8 i 2 B .9c z ijk s i = x ijk 8 i 2 B;j 2 C;k 2 D .9d z ijk integer 8 i 2 B;j 2 C;k 2 D .9e k 0 8 k 2 D: .9f 4.6SolutionMethods WhileaproblemdoesnothavetobeNP-hardtoposedicultywithndingsolutions, thenatureofBTMissuchthatasolutionmaynotbefoundwithoutanexhaustive search,ifoneexists.Thatbeingsaid,thesolutionmethodsforlinearprograms, 49

PAGE 59

thesimplexmethodandinteriorpointmethodscanbeappliedandcouldresultin feasible,optimalsolutionsintheintegerspace.Whentheydonot,otherapproaches mustbeused.Themethodsexploredarelistedinthissection. 4.6.1LPFormulation WhensolvinggeneralMILPs,aninitialstepfrequentlyistosolvetheLPrelaxation. ThisisnodierentforBTM,wherethedecisionvariables z ijk arefractionalinstead ofinteger.Thisformulationmayresultinasolutionthatisintheintegerdomain, isquickandeasytoexecuteandprovidesanabsolutelowerboundontheobjective functionusingthedualproblem.Evenifthesolutionfoundisnotentirelyinteger, employingaroundingmethodmayreturnasolutionthatisclosetofeasible.Thelevel ofviolationallowedshouldbelefttothedecisionmaker,whichisjustreferringtothe violationofcapacityconstraints.Theviolationofowshouldnotbetolerated.Thisis whereasimpleroundingmethodmaynotbeusable.Forinstance,ifasolutionfound suggestssendinglessthanhalfaparticularloadtothreedierentdestinations,simple roundingwouldresultinthatloadbeingassignedtonodestination.Itispossible thateachpositivevalueinthiscaseisthesame,whichcallsoramoresophisticated methodtodecidewhichofthepositivevaluesshouldbe1. 4.6.2SimpleGreedyHeuristic Thesimplegreedyheuristicusesthefollowingalgorithm.ThetokenNEXTisan unassignedloadchosenbysomemechanismthatiteratesoverallloads.Thetoken QUITisforwhenthealgorithmsuccessfullyterminates,havingfoundasolution. ThetokenFAILisforwhenthealgorithmterminateshavingfailedtondafeasible 50

PAGE 60

Table4.2:BTM:GreedyHeuristicExample1 Load AssignedDestination 100 200 100 200 120 220 140 250 150 None solution. 1.IflistofchoicesforNEXTisempty,FAIL. 2.RankeachchoiceforNEXT. 3.AssignthelowestcostchoicetoNEXT. 4.Ifnoloadsareunassigned,thenQUIT. 5.Otherwise,moveNEXTtoanunassignedload. Thismethodprovidesnoguaranteeofoptimalityorevenndingafeasiblesolutionif oneexists,soitisheuristicinnature.Considerthisexampletogaininsightintothe heuristic.Let s =[100 ; 100 ; 120 ; 140 ; 150]and u =[200 ; 220 ; 250]forthedestination capacities.Supposethatthechoiceoforderinggives[100 ; 100 ; 120 ; 140 ; 150]andthat eachdestinationcoststhesamebetweentheloads,withtheorderofthedestinations beingthepreferenceforeachload.So,theassignmentsaregiveninTable4.2and theheuristicfails.Butiftheorderingoftheloadsis[140 ; 100 ; 120 ; 100 ; 150],the assignmentsaregiveninTable4.3whichisasuccessfulexecutionoftheheuristic. 51

PAGE 61

Table4.3:BTM:GreedyHeuristicExample2 Load AssignedDestination 140 200 100 220 120 220 100 250 150 250 ThegeneralityofthemechanismdeterminingNEXTallowsfordierentrulesto beappliedtonddierentsolutions.Someexamplesaregiven. Producerpriority:theunassignedlocationwiththehighestvolumetodeliver isgivenprioritytobechosenateachiteration. Maxdistancepriority:theunassignedlocationwiththefurthestpossibledistancetotravelisgivenpriorityateachiteration. Mindistancepriority:theunassignedlocationwiththenearestpossibledistance totravelisgivenpriorityateachiteration. Dierentcombinationsoftheseorderingmechanismscanbeusedtochangetheorder inwhichloadsareassignedusingthegreedyheuristic.Theproducerpriorityordering seekstogivethehighestproducersthebestchoicepossible.Themaxdistancepriority orderingseekstoeliminatetheworstpossibletrips.Themindistancepriorityseeks tomakethedecisionsthatseemthebestforeachlocation. 52

PAGE 62

4.6.3NetworkFlowSolutionwithFullLoads Theorem4.4.4callsforthemodicationoftheoriginalproblemtochangeallload sizestobefullloads.Thiscanbeawildapproximationwhentheloadsinaparticular instancearedominatedbypartialloads,withsomefullloadstobehauled.Evenso, whenafeasiblesolutionexistsfortheFBTMproblem,itmaybesatisfactoryfordaily use.Theecientsolutionmethodsthatexistfornetworkowproblemsmakeitan attractiveoption. 4.7ComputationalResults 4.7.1InstanceCreation TocreatedataforuseinBTM,thedatafromeachinstancedescribedinAppendixA needstobetranslatedintovectorsforloads x ijk ,costs c ijk anddestinationcapacities u k .Theprocessuseddoesnotseektocreatefeasibleinstances.Theloadselection methodsimplylooksforallloadsreadyininventory.Thegoalofexperimentation withBTMistoseehowdierentmethodsbehavewithalltypesofproblems.The infeasibleproblemsareofasmuchinterestasthefeasibleones. Theloadgenerationisofthemostinteresthere.Eachprobleminstancehasinitial inventory v i .Thesmallestloadsize l j fromthesetofhaulersisusedtodecidehow manyloadswereavailabletohaulfromeachbatterybydividingtheinitialinventory bytheloadsizeandtakingtheooroftheresult.Ifthelocationswithzeroloads followingthatprocedurehasenoughinventoryfor1/4ofaload,thentheyareselected aswell.Thisistorepresentmultiplepickupsandensuredierentloadsizes. 53

PAGE 63

Asinotherexperiments,thecostforeachdecisionwascreatedusingthe ` 1 norm ofthedierenceofthetwo-vectorsindicatinglocationontheplane.Thecapacities forthedestinationsweretakendirectlyfromthegenerateddata. 4.7.2SolutionstoBTM DataforeachinstanceofBTMwasgeneratedusingasinglescriptandappropriatesolver.EachinstanceofBTMwassolvedusingthepureMIPmodel,GBTM andNFRimplementation.Theseresultsoftheseimplementationswerecomparedto demonstrateadvantagesofeachmethod. TheplotsinFigures4.2,4.3and4.4showtheperformanceoftheMIPsolver onpureBTM.Theresultsshowthatmostinstancescompletedexecutioninshort amountsoftime.Figure4.2showstheMIPgapvsLPRobjective.Thelineisthe tolerancesuppliedtothesolverfortheMIPgaptodetermineoptimality.Thepoints abovethelineareinstancesthatdidnotmeetthetoleranceandtimedout.InFigure 4.4,theplusmarkersshowinstancesdeemedinfeasible.Noticehowsomemedium sizedinstancesreachedoptimalityintimesbetweenthelowlevelsandthemaximum of180s.Therearenolargeinstancesthatwereexecutedbetweenlessthanafew secondsandthetimelimit.Theseresultsreinforcethedicultyoftheproblem. Thenextsetofplots,Figures4.5,4.6and4.7showpatternsthataresimilar, butrevealthatthegoalprogrammingapproachhasmoreofanissuewithgrowing complexityandexecutiontime.Thecomparisonisnotdirectlycomparablebecause amajorfractionofthelargeinstancessolvedwithBTMwereinfeasible.InGBTM, noproblemisinfeasible.NotethattheMIPgapusedisthegapbetweentheMIP objectiveandLPRobjectiveonlyintermsofthedecisionvariables.Theportionof 54

PAGE 64

Figure4.2:BTM:MIPGapvs.OptimalLPRObjective theobjectivethatdrovetheconstraintviolationcostisleftoutforcomparisontothe pureBTMresults. Finally,theNFRplotsinFigures4.8,4.9and4.10demonstratehowthenew approachcanhelpinndingsolutionsinpolynomialtime.Noticethattheinfeasible solutionsstillshowup,butnoinstancestooklongerthan30secondstoarriveat optimalityorinfeasibility. 55

PAGE 65

Figure4.3:BTM:WallTimeHistogram Figure4.4:BTM:WallTimevs.InstanceSize 56

PAGE 66

Figure4.5:GBTM:MIPGapvs.OptimalLPRObjective Figure4.6:GBTM:WallTimeHistogram 57

PAGE 67

Figure4.7:GBTM:WallTimevs.InstanceSize Figure4.8:NFR:MIPGapvs.OptimalLPRObjective 58

PAGE 68

Figure4.9:NFR:WallTimeHistogram Figure4.10:NFR:WallTimevs.InstanceSize 59

PAGE 69

CHAPTERV CONCLUSION Thisthesishasdemonstratedtwodierentmathematicaltoolsandtheirapplicationinmanagementofliquidsinventoryandhauling.StochasticinventorymanagementSTIMisaplanningtoolthathelpstohowresourcesshouldbearranged andwhatproductionproletoexpect.BatchedtransportationmanagementBTM isusefulfortacticaldecisionsfortransportationandresourceusage.Bothproblems areNP-hard.Inlargeimplementations,heuristicandapproximationmethodsare requiredtodeliversolutionsthatseekoptimality. STIMhasbeencreatedtondhauling,productionandinventorypatternsseeking optimalconditionsforminimizingshutinvolumes,maximizinghaulingandminimizingendinginventory.Itisusefulforplanningpurposes,butalsocanbeimplemented tondtacticaldecisionsforhauling.Theoptionsavailableforthismodelallowitto beusedtoseekdierentgoals.Thesimplicationsusedinthismodelhaveallowed theuseofanecientalgorithmtoarriveatsolutions.Theheuristicmethodapplied toSTIMhasbeenshownviasimulationtoprovideguidanceonhowtomanageinventorythroughhauling.Eventhoughthemethodsuppliedfractionalsolutions,the roundingmethodorpolishingheuristicresultedinintegraldecisionsthatseektominimizethelocationsvisitedwhilemaximizinghaulvolume.Thepatternsthatresult areasexpectedwithhighersplitstohaulmorevolume.Thedistributionusedfor generatingexperimentationdatahasgivenahighlyaccurateplanningtosimulation predictionforhaulandendinginventoryvolumesasdemonstratedbyourcomputationalresults.Thesearekeyindicatorsfordecisionmakersontheperformanceofoil andwaterhauling.Thisdemonstratesthattheheuristicmethodcanbeusefulfor businessapplicationofSTIM. TherunningtimesofpureBTMshowatendencytoexceedtimelimitsrequiredof decisionmakersforproblemsizesequivalenttomediumtolargeoperations.GBTM 60

PAGE 70

showsthatintroducingfeasibilityonotherwiseinfeasibleproblemsonlyprovidesan answer,whileoptimalityisstillunansweredformanyinstancesinthetimeallotted. NFRoersasolutionmethodboundedbyareasonablepolynomial.Thesolutions returnedbythismethodarenotnecessarilyusefulfordecisionmakersoutrightdue toallowingcapacityconstraintviolations,butmayleadtofurthersolutionmethods. TheimplementationofBTMandrelatedmodelsGBTMandNFRhasshownthat eachhasitsapplicationforspeciccases.ThepureIPimplementationofBTMcanbe ecientforndingresultsoffeasibleproblemsandfrequentlydeliversthosesolutions quicklywithamethodcombinedwithbranchandcutandheuristics.GBTMisnot asecient,butdoesprovideananswer,evenifnotdirectlyfeasible.Theinfeasibility canbetunedtoreducecapacityviolationsbydesireddestination.Feasiblesolutions aredeliveredquicklyandcanbefoundbyadirectcalculationinpolynomialtime. Optimalsolutionshavenotbeendemonstratedforlargeinstancesduetothesolver reachingtimelimit.NFRoersauniquemixofthetwousinganetworkowmethod onamodiedinstance.Itcanprovideapproximatesolutionsinpolynomialtime. Theuseofthesetwomodelsinconjunctionwitheachotherwillprovidedecision supportfortacticalandplanningsituations.STIMisusefulforplanninginshortand longtermandwhencombinedwithBTMfortacticaldecisionsonresourceusage,the resultishavingknowledgeofhowtogetthemostoutoftheminthecurrentsituation. Followupworkinthiseldmayincludebetterecienciesgainedthroughdecisions madeclosertorealtimeandanadditionallayerofplanningintruckrouting. 61

PAGE 71

REFERENCES [1]Theengineeringtoolbox:Apigravity.TheEngineeringToolbox, http://www. engineeringtoolbox.com/api-gravity-d_1212.html .Accessed:21-MAR2015. [2]Gnulinearprogrammingkit. http://www.gnu.org/software/glpk/glpk. html [3]Networkx:High-productivitysoftwareforcomplexnetworks. https:// networkx.github.io/ [4]Pulp. http://www.coin-or.org/projects/PuLP.xml [5]Pythonlanguagereference,version2.7.PythonSoftwareFoundation, https: //docs.python.org/2.7/ [6]OileldGlossary.SlumbergerLimited, http://glossary.oilfield.slb.com/ 2015.[Online;accessed16-March-2015]. [7]W.K.AbduljabbarandR.M.Tahar.Acasestudyofpetroleumtransportation logistics:adecisionsupportsystembasedonsimulationandstochasticoptimal control. AfricanJournalofBusinessManagement ,6:4350{4361,2012. [8]R.K.Ahuja,T.L.Magnanti,andJ.B.Orlin. NetworkFlows:Theory,Algorithms,andApplications .PrenticeHall,1993. [9]F.Altiparmak,M.Gen,L.Lin,andT.Paksoy.Ageneticalgorithmapproachfor multi-objectiveoptimizationofsupplychainnetworks. Computers&Industrial Engineering ,51:196{215,2006. [10]V.J.D.Baston,M.K.Rahmouni,andH.P.Williams.Thepracticalconversion oflinearprogrammestonetworkowmodels. EuropeanJournalofOperational Research ,50:325{334,1991. [11]D.P.BertsekasandP.Tseng.Relaxationmethodsforminimumcostordinary andgeneralizednetworkowproblems. OperationsResearch ,36:93{114,1988. [12]R.E.BixbyandW.H.Cunningham.Convertinglinearprogramstonetwork problems. MathematicsofOperationsResearch ,5:pp.321{357,1980. [13]S.Boyd.L1normmethodsforconvexcardinalityproblems. http://web. stanford.com/class/ee364b/lectures/l1_slides.pdf .Accessed:05-MAY2015. 62

PAGE 72

[14]S.Boyd.Modelpredictivecontrol. http://web.stanford.com/class/ee364b/ lectures/mpc_slides.pdf .Accessed:05-MAY-2015. [15]L.BukataandP. Sucha.Agpualgorithmdesignforresourceconstrainedproject schedulingproblem.In Proceedingsofthe201321stEuromicroInternational ConferenceonParallel,Distributed,andNetwork-BasedProcessing ,PDP'13, pages367{374,Washington,DC,USA,2013.IEEEComputerSociety. [16]A.Campbell,L.Clarke,A.Kleywegt,andM.Savelsbergh.Theinventoryrouting problem.In FleetManagementandLogistics ,pages95{113.Springer,1998. [17]A.Chaparala,C.Novoa,andA.Qasem.Asimdsolutionforthequadraticassignmentproblemwithgpuacceleration.In Proceedingsofthe2014AnnualConferenceonExtremeScienceandEngineeringDiscoveryEnvironment ,XSEDE'14, pages1:1{1:8,NewYork,NY,USA,2014.ACM. [18]L.ChengandM.A.Duran.Logisticsforworld-widecrudeoiltransportation usingdiscreteeventsimulationandoptimalcontrol. Computers&Chemical Engineering ,28-7:897{911,2004. [19]S.Chopra,I.Gilboa,andS.T.Sastry.Sourcesinkowswithcapacityinstallationinbatches. DiscreteAppliedMathematics ,85:165{192,1998. [20]S.ChristophersonandK.Dave.Cardireports.TechnicalReport15,CaRDI Reports,CornellUniversity,November2014. [21]A.Engau,C.Moatt,andW.Dyk.Multicriteriamodelingandtradeoanalysis foroilloaddispatchandhaulingoperationsatNobleEnergy. Optimizationand Engineering ,16:73{101,2015. [22]R.Z.FarahaniandM.Elahipanah.Ageneticalgorithmtooptimizethetotal costandservicelevelforjust-in-timedistributioninasupplychain. International JournalofProductionEconomics ,111:229{243,2008. [23]L.R.FordandD.R.Fulkerson. FlowsinNetworks .RANDCorporation,Santa Monica,Calif.,1962. [24]M.E.LalamiandD.El-Baz.Gpuimplementationofthebranchandbound methodforknapsackproblems.In Proceedingsofthe2012IEEE26thInternationalParallelandDistributedProcessingSymposiumWorkshops&PhDForum IPDPSW'12,pages1769{1777,Washington,DC,USA,2012.IEEEComputer Society. [25]S.M.Leeetal. GoalProgrammingforDecisionAnalysis .AuerbachPhiladelphia,1972. [26]W.A.LodwickandJ.Kacprzyk. FuzzyOptimization:RecentAdvancesand Applications .StudiesinFuzzinessandSoftComputing.Springer,2010. 63

PAGE 73

[27]S.MartelloandP.Toth. KnapsackProblems:AlgorithmsandComputerImplementations .JohnWiley&Sons,Inc.,NewYork,NY,USA,1990. [28]S.M.S.NeiroandJ.M.Pinto.Ageneralmodelingframeworkfortheoperationalplanningofpetroleumsupplychains. Computers&ChemicalEngineering 28:871{896,2004. [29]C.H.PapadimitriouandK.Steiglitz. CombinatorialOptimization:Algorithms andComplexity .DoverBooksonComputerScience.DoverPublications,2013. [30]P.Radhakrishnan,V.M.Prasad,andM.R.Gopalan.Geneticalgorithmbased inventoryoptimizationanalysisinsupplychainmanagement.In AdvanceComputingConference,2009.IACC2009.IEEEInternational ,pages418{422.IEEE, 2009. [31]A.Shapiro,D.Dentcheva,etal. LecturesonStochasticProgramming:Modeling andTheory ,volume16.SIAM,2014. [32]Q.Shen,F.Chu,andH.Chen.ALagrangianrelaxationapproachforamultimodeinventoryroutingproblemwithtransshipmentincrudeoiltransportation. Computers&ChemicalEngineering ,35:2113{2123,2011. [33]D.B.ShmoysandE.Tardos.Anapproximationalgorithmforthegeneralized assignmentproblem. MathematicalProgramming ,62-3:461{474,1993. [34]R.M.TaharandW.K.Abduljabbar.Anoveltransportingsystemmodelforoil renery. AmericanJournalofEngineeringandAppliedSciences ,2:138{143, 2010. 64

PAGE 74

APPENDIXA.DataGeneration Thebulkofthedatausedinexperimentationhasbeengenerated.Thegoalof muchoftheworkistoshowtheeciencyofsolversinexecutiontimeanddeterminingfeasibility.Eachinstancerepresentsasingleproduct.Itcanberegarded asrepresentingconsistentvaluebetweenlocations.Forwater,revenueiszero.For oil,revenuecanchangebetweenthedestinationofdelivery,whichisaconsideration drivingthedevelopmentofthemulticriteriamodelin[21].Otherwise,costsarerepresentativeofeitherproduct,netoilvolumesarereasonablyconsistentforweightand watervolumesareaswell.Thusthisexperimentationisagnostictowardtheproduct. Thecomputationalmethodsinthisthesisuseasetofprobleminstancegenerated specicallyforthiswork.Theseinclude100caseseachof 20batteries,1hauler,2destinations 40batteries,2haulers,3destinations 100batteries,5haulers,6destinations Thebatteriesineachcasearebuiltusingrandomdata.Abatterywasinitiallychosen tobeofoneofthetypesinTableA.1.Afterrandomlychoosingatype,theparticular baserateofproductionischosenusinganormaldistributionaroundthetype'srate. Thestandarddeviation usedwas5%ofthechosentyperate.Usinganormal distributionaroundtherandombaserate andthesamestandarddeviation,the ratevaluescorrespondingto )]TJ/F15 11.9552 Tf 11.627 0 Td [(2 )]TJ/F19 11.9552 Tf 11.627 0 Td [(: 75 , + : 75 and +2 wereassigned. Thesevaluesarereferredtoasq10,q35,q50,q65andq90inthecode.Theprobability assignedtoeachofthesecasesisthesameacrossbatteries.Itiscalculatedbyusing astandardnormaltable.Theintervalforeachcasehasawidthof1/2standard deviationandiscenteredonthedistancefromthemean.Usingthelawoftotal 65

PAGE 75

TableA.1:DataGeneration:BaseDataforGenerationofBatteries Class ProductionRateInventoryCapacitySelectionFrequency LowProducer 1BBL300BBL2 = 7 LowMidProducer 1BBL300BBL2 = 7 HighMidProducer 100BBL600BBL2 = 7 HighProducer 1000BBL2400BBL1 = 7 TableA.2:DataGeneration:ProductionRateIntervalProbabilities Case IntervalProbabilityMarginalProbability q10 )]TJ/F15 11.9552 Tf 11.955 0 Td [(1 : 75 )]TJ/F19 11.9552 Tf 11.956 0 Td [( )]TJ/F15 11.9552 Tf 11.955 0 Td [(2 : 25 P q 10 =P A =0 : 05047033 q35 )]TJ/F15 11.9552 Tf 11.955 0 Td [(0 : 5 )]TJ/F19 11.9552 Tf 11.956 0 Td [( )]TJ/F15 11.9552 Tf 11.955 0 Td [(1 : 0 P q 35 =P A =0 : 27098408 q50 +0 : 25 )]TJ/F19 11.9552 Tf 11.955 0 Td [( )]TJ/F15 11.9552 Tf 11.955 0 Td [(0 : 25 P q 50 =P A =0 : 35709117 q65 +1 : 0 )]TJ/F19 11.9552 Tf 11.955 0 Td [( +0 : 5 P q 65 =P A =0 : 27098408 q90 +2 : 25 )]TJ/F19 11.9552 Tf 11.955 0 Td [( +1 : 75 P q 90 =P A =0 : 05047033 probability,thecaseprobabilitiesarecalculatedinTableA.2.Event A istheevent thattheratefallsinoneoftheintervals,withprobability P A = P q 10+ P q 35+ P q 50+ P q 65+ P q 90 : Theinventorycapacityisunmodied.Theinitialinventoryischosentobea uniformrandomvaluebetween0andtheinventorycapacity.Eachbatterywasalso givencoordinates, x;y ,torepresentthelocationonthesquare )]TJ/F15 11.9552 Tf 9.298 0 Td [(50 ; )]TJ/F15 11.9552 Tf 9.298 0 Td [(50 ; ; 50. Whileitisrepresentativeofsomelevelofoperation,ithasnotbeenbuiltaroundthe operationsofanyparticularE&Pcompany. Eachhaulerwasgeneratedwithaminimumandmaximumloadcountperperiod, fullloadsize,dryloadfeeandsplitfee.Thedryloadfeewasnotusedduring optimization,butduringsimulationtoshowaneectofoverzealousdispatch.See 66

PAGE 76

TableA.3:DataGeneration:BaseDataforGenerationofHaulers Variable Value LowerLimitonLoadsHauledPerPeriod l j 0or U [0 ; 40]withequalprobability UpperLimitonLoadsHaulerPerPeriod u j l j + U [0 ; 60] LoadSize Uniformchoiceof f 180 ; 190 ; 200 ; 250 g DryLoadFee Uniformchoiceof f 75 ; 100 ; 125 ; 150 g SplitLoadFee Uniformchoiceof f 75 ; 100 ; 125 ; 150 g TableA.3forthedetail. Thedestinationswerealsogivencoordinatesinthesamesquare.Alongwithcoordinates,aminimumandmaximumdeliveryvolumewasassignedusingthefollowing relationship.Let A;B U [0 ; 5000].Then l k = A A.1 u k = l k +3 B; A.2 sothateachdestinationhasaminimumdelivery l k upto5000BBLandmaximum delivery u k between l k and l k +15000BBL. Thisisthewholeofthebasedatausedforgeneratingprobleminstances.See AppendixC.1forthecodeused. 67

PAGE 77

APPENDIXB.MethodExecution Allsolver-relatedcodecreatedforthisworkhasbeenwritteninPython[5]. TheprojectPuLP[4]associatedwiththeCOIN-ORinitiativeallowsthedirectmanipulationofvariables,constraintsandobjectivestotranslateprobleminstancesto structuresneededbysolvers.EachmodelsolvedusingLP/MILPwasinterfacedwith GLPK[2]throughPuLP.Fornetworkowoptimization,theprojectNetworkX[3]for Pythonwasused.Thisprojecthasanetworksimplexmethodthatsolvesmincost owproblems,satisfyingalldemands.AsinglescriptseeC.3intheappendixwas usedtoexecuteeachsolverforaparticularinstance.Thescriptrequirestheinstance ID,theinstancerepositoryandoutputdestinationdatabase.Databaseaccesswas doneusing anydbm inthestandardPythonlibrary[5]toproviderandomaccessfor readingandwritingthedata.Executionofthescriptbeginswithdeterminationof themodeofoperation,whichmodeltypeorsolvertorun.Theinstancedataisread fromthetextles,translatedintotherequiredstructurefortherequestedmethod andfedtothesolver.IfthesolverisoftypeMILP,thenthesearchtimelimitisset to180secondswithastoppingcriterionof0.1%MIPgap. Followingcompletionofthemethod,eitherbysuccessfuloptimization,certicate ofinfeasibilityortimelimitexceeded,theresultingdataiswrittentothedestination database.ThedatakeysincludedarelistedinTableB.1.Severalplotsandtables weregeneratedusingthisdataforeachmethodtocompareperformance.Ifadditional datapointsarerequiredforaparticularmethod,thenthosepointsarediscussedin therelevantsectionforthatmethod. Foreachmethod,threeplotsaregenerated.TherstisaplotofMIPGapvs OptimalLPObjective.Thisplotisintendedtodemonstratethestoppingcriteriaof themethodintermsoftheabsolutelowerbound.Eachinstancethathasafeasible 68

PAGE 78

solutionisplottedasadotandthestoppingcriteriaof0.1%isplottedasasolid line.Eachdotbelowthelinerepresentaninstancewhichreachedtoleranceorproven optimalityandthedotsabovethelinerepresentthoseinstancesthatreachedafeasible solution,butranabovethetimelimit. Thesecondplotisahistogramoftheinstancecountinbinsofsolverduration. Thisplotisstraightforwardanddemonstratesabimodalresultwherelargenumbers ofinstancesfallintocategoriesoflowprocessingtimeandhighprocessingtime. Thethirdplotisrelativelysimpleaswell.Itplotsascatterofwallclocktime versusintegervariablesperinstance.Oneadditionalvisualrepresentationinthisplot isthedistinctionoffeasibleandinfeasibleinstanceswithfeasibleinstancesplotted withadotandinfeasibleinstancesplottedwithcrosshairs.Feasibleandinfeasible instancesshouldbemixedaroundthelowprocessingtimerangeandthereshould beahighernumberoffeasibleinstancesinthehigherprocessingtimerange.Thisis duetothenumberofinstancesthatareinfeasiblemerelybecausethehaulvolume exceedsthedeliverycapacity. 69

PAGE 79

TableB.1:DataGeneration:StructurefortheDataGeneratedinEveryInstance Key Description clocktime processortimerequiredtondsolutionPythoncodeonly walltime durationoftimetondasolution objective smallestfeasibleobjectivevaluefound lpobjective provenminimumvalueforLPrelaxation haulvol totalhaulvolumebyperiod,eitherscalarorarray volcap totaldeliverycapacity loadcount totalloadstobedelivered destcount countofdestinations solution solutionfoundbysolverorNone,ifinfeasible 70

PAGE 80

APPENDIXC.GeneralMethods C.1InstanceGenerator 1 # )61()]TJ/F47 5.9776 Tf 19.179 0 Td [(coding:utf )]TJ/F47 5.9776 Tf 5.558 0 Td [(8 )61()]TJ/F18 5.9776 Tf -100.56 -11.577 Td [(2 import numpyasnp 3 import csv 4 5 def create instanceB,C,D: 6""" 7 Create a Hauling problem instance with B number of batteries, 8 C number of haulers and D number of destinations in a 100 x 100 9 square layout. 10 11 """ 12bbox=np.array[ )]TJ/F18 5.9776 Tf 6.522 0 Td [(50, )]TJ/F18 5.9776 Tf 6.24 0 Td [(50,50,50] 13producers=np.array[1,1,10,10,100,100,1000] 14prodchoice= range len producers 15inv=np.array[300,300,300,300,600,600,2400] 16fees=np.array[75,100,125,150] 17loadsize=np.array[180,190,200,250] 18chooser=np.array[0,1] 19 def get producer: 20coords=bbox[0]+bbox[1] )]TJ/F18 5.9776 Tf 10.073 0 Td [(bbox[0] np.random.random2, 21i=np.random.choiceprodchoice 22base=producers[i] 23cap=inv[i] 24 #randomizebase 25sigma=base/10. 26base=np.random.normalbase,sigma 27 if base < 0:base=0 28b10=base )]TJ/F18 5.9776 Tf 10.047 0 Td [(2 sigma 29b35=base )]TJ/F18 5.9776 Tf 10.617 0 Td [(.75 sigma 30b50=base 31b65=base+.75 sigma 32b90=base+2 sigma 33v0=np.random.random cap 34 return [coords[0],coords[1],b10,b35,b50,b65,b90,cap,v0] 35 36 def get hauler: 37dryfee=np.random.choicefees 38splitfee=np.random.choicefees 39full=np.random.choiceloadsize 40l= int np.random.choicechooser 40 np.random.random 41u=l+ int 60 np.random.random 42 return [l,u,full,dryfee,splitfee] 43 44 def get destination: 45coords=bbox[0]+bbox[1] )]TJ/F18 5.9776 Tf 10.073 0 Td [(bbox[0] np.random.random2, 71

PAGE 81

46l= int np.random.choicechooser 5000 np.random.random 47u=l+3 int 5000 np.random.random 48 return [l,u,coords[0],coords[1]] 49 50batteries=[get producer for b inrange B] 51haulers=[get hauler for h inrange C] 52destinations=[get destination for d inrange D] 53 return batteries,haulers,destinations 54 55 def write dataheader,rows,fn: 56with file fn,"wb"asf1: 57w=csv.writerf1 58w.writerowheader 59w.writerowsrows 60 61 if name ==" main ": 62 import getopt 63 import time 64 import sys 65 66usage=" f 0 g )]TJ/F18 5.9776 Tf 5.471 0 Td [(b < battery count > )]TJ/F18 5.9776 Tf 5.858 0 Td [(c < hauler count > )]TJ/F18 5.9776 Tf 5.47 0 Td [(d < market count > 67usage=usage. format sys.argv[0] 68 69opts,args=getopt.getoptsys.argv[1:],"b:c:d:" 70b=None 71c=None 72d=None 73 try : 74 for i inrange len opts: 75 if opts[i][0]==" )]TJ/F18 5.9776 Tf 5.47 0 Td [(b": 76b= int opts[i][1] 77 elif opts[i][0]==" )]TJ/F18 5.9776 Tf 5.857 0 Td [(c": 78c= int opts[i][1] 79 elif opts[i][0]==" )]TJ/F18 5.9776 Tf 5.47 0 Td [(d": 80d= int opts[i][1] 81 except ValueError,e: 82sys.exit )]TJ/F18 5.9776 Tf 6.151 0 Td [(1 83 if b is None or c is None or d is None: 84sys.exit )]TJ/F18 5.9776 Tf 6.151 0 Td [(1 85 86batteries,haulers,destinations=create instanceb,c,d 87bheader=["x","y","q10","q35","q50","q65","q90","invcap","v0"] 88hheader=["l","u","loadsize","dryfee","splitfee"] 89dheader=["l","u","x","y"] 90instnum= int 10 time.time )]TJ/F18 5.9776 Tf 10.531 0 Td [(1429384669.305739 91write databheader,batteries,"model )]TJ/F18 5.9776 Tf 6.122 0 Td [(batteries )-77(f 0 g .csv". format instnum 92write datahheader,haulers,"model )]TJ/F18 5.9776 Tf 5.936 0 Td [(haulers )-77(f 0 g .csv". format instnum 93write datadheader,destinations, 94"model )]TJ/F18 5.9776 Tf 6.078 0 Td [(destinations )-77(f 0 g .csv". format instnum 95time.sleep.25 C.2ModelData 1 # )61()]TJ/F47 5.9776 Tf 19.179 0 Td [(coding:utf )]TJ/F47 5.9776 Tf 5.558 0 Td [(8 )61()]TJ0 g 0 G/F15 11.9552 Tf 333.35 -18 Td [(72

PAGE 82

2 import pandasaspd 3 import os.path 4 import numpyasnp 5 6 def read instanceintrepo,instnum: 7bfn=os.path.joinintrepo,"model )]TJ/F18 5.9776 Tf 6.122 0 Td [(batteries )-77(f 0 g .csv". format instnum 8hfn=os.path.joinintrepo,"model )]TJ/F18 5.9776 Tf 5.935 0 Td [(haulers )-77(f 0 g .csv". format instnum 9dfn=os.path.joinintrepo,"model )]TJ/F18 5.9776 Tf 6.078 0 Td [(destinations )-77(f 0 g .csv". format instnum 10bdf=pd.read csvbfn 11hdf=pd.read csvhfn 12ddf=pd.read csvdfn 13 return bdf,hdf,ddf 14 15 def build btmbdf,ddf,l,u: 16loads=bdf[bdf["v0"] > l] 17loadcounts=np.floorloads["v0"]/u 18addtlloads=loadcounts==0 19loadsizes=0 loads["v0"] 20loadsizes[loadcounts==0]=loads[loadcounts==0]["v0"] 21loadsizes[loadcounts > 0]=u 22 #buildcostsanddestinations 23destinations=ddf["u"] 24costs=[] 25 for idx in ddf.index: 26x=np. abs ddf.loc[idx,"x"] )]TJ/F18 5.9776 Tf 10.607 0 Td [(loads["x"] 27y=np. abs ddf.loc[idx,"y"] )]TJ/F18 5.9776 Tf 10.607 0 Td [(loads["y"] 28d=x+y 29costs.appendnp.arrayd 30costs=np.arraycosts.transpose.tolist 31loads.loc[:,"loadsizes"]=loadsizes 32loads.loc[:,"loadcounts"]=loadcounts 33loads.loc[loadcounts==0,"loadcounts"]=1 34 return loads,costs,destinations C.3MethodExecution 1 # )61()]TJ/F47 5.9776 Tf 19.179 0 Td [(coding:utf )]TJ/F47 5.9776 Tf 5.558 0 Td [(8 )61()]TJ/F18 5.9776 Tf -100.56 -11.576 Td [(2""" 3run instance in BTM solver and STINV, save output to DBM file 4 5""" 6 7 from newcastle.opt import btm,nfr,exc 8 from newcastle.opt import cardinventory 9 from newcastle.opt import goal 10 from newcastle.utils import modeldata 11 from newcastle.utils.simulateinv import run sim 12 import numpyasnp 13 import time 14 15 def run btmloads,costs,destinations: 16loadsizes=loads["loadsizes"] 17loadcounts=loads["loadcounts"] 18nestloads=[[loadsizes[i],] loadcounts[i] for i in loadsizes.index] 73

PAGE 83

19loads=[] 20newcosts=[] 21 for i,ii inzip range len nestloads,loadsizes.index: 22loads.extendnestloads[i] 23 for j inrange int loadcounts[ii]: 24newcosts.appendcosts[i] 25b=btm.BTMProblemloads,newcosts,destinations,relaxed=False 26 try : 27start=time.clock 28cstart=time.time 29solution=b.solve 30cstop=time.time 31stop=time.clock 32feasible=True 33obj=b.objective 34lpobj=b.lpobjective 35 except exc.FailedSolution: 36cstop=time.time 37stop=time.clock 38feasible=False 39obj=np.inf 40 if b.lpobjective is None: 41lpobj=np.inf 42 else : 43lpobj=b.lpobjective 44ptime=stop )]TJ/F18 5.9776 Tf 10.802 0 Td [(start 45cltime=cstop )]TJ/F18 5.9776 Tf 10.801 0 Td [(cstart 46 returnsum loads, sum destinations,obj,lpobj,ptime,cltime,solution 47 48 def run nfrloads,costs,destinations: 49loadsizes=loads["loadsizes"] 50loadcounts=loads["loadcounts"] 51nestloads=[[loadsizes[i],] loadcounts[i] for i in loadsizes.index] 52loads=[] 53newcosts=[] 54 for i,ii inzip range len nestloads,loadsizes.index: 55loads.extendnestloads[i] 56 for j inrange int loadcounts[ii]: 57newcosts.appendcosts[i] 58obj=np.inf 59lpobj=np.inf 60orig=btm.BTMProblemloads,newcosts,destinations,relaxed=True 61b=nfr.BTMProblemloads,newcosts,destinations 62 try : 63start=time.clock 64cstart=time.time 65orig.solve 66lpobj=orig.objective 67solution=b.solve 68cstop=time.time 69stop=time.clock 70feasible=True 71obj=b.objective 72 except exc.FailedSolution: 73solution=None 74

PAGE 84

74cstop=time.time 75stop=time.clock 76feasible=False 77ptime=stop )]TJ/F18 5.9776 Tf 10.802 0 Td [(start 78cltime=cstop )]TJ/F18 5.9776 Tf 10.801 0 Td [(cstart 79 returnsum loads, sum destinations,obj,lpobj,ptime,cltime,solution 80 81 82 def run goalloads,costs,destinations: 83loadsizes=loads["loadsizes"] 84loadcounts=loads["loadcounts"] 85nestloads=[[loadsizes[i],] loadcounts[i] for i in loadsizes.index] 86loads=[] 87newcosts=[] 88maxcost=0 89 for i,ii inzip range len nestloads,loadsizes.index: 90loads.extendnestloads[i] 91 for j inrange int loadcounts[ii]: 92newcosts.appendcosts[i] 93 if maxcost < max costs[i]: 94maxcost= max costs[i] 95factors=[maxcost factor 1000 for factor inrange len destinations] 96b=goal.BTMProblemloads,newcosts,destinations 97 try : 98start=time.clock 99cstart=time.time 100solution=b.solvefactors 101cstop=time.time 102stop=time.clock 103feasible=True 104obj=b.objective 105lpobj=b.lpobjective 106 except exc.FailedSolution: 107cstop=time.time 108stop=time.clock 109feasible=False 110obj=np.inf 111 if b.lpobjective is None: 112lpobj=np.inf 113 else : 114lpobj=b.lpobjective 115ptime=stop )]TJ/F18 5.9776 Tf 10.801 0 Td [(start 116cltime=cstop )]TJ/F18 5.9776 Tf 10.801 0 Td [(cstart 117 returnsum loads, sum destinations,obj,lpobj,ptime,cltime,solution 118 119 def run stimbdf,hdf,ddf,periods: 120caseprob=np.array[0.05047033,0.27098408,0.35709117, 1210.27098408,0.05047033] 122p=cardinventory.STIMProblembdf,hdf,ddf,periods,caseprob 123 try : 124start=time.clock 125cstart=time.time 126solution=p.solve 127cstop=time.time 128stop=time.clock 75

PAGE 85

129feasible=True 130obj=p.objective 131lpobj=p.lpobjective 132siprod=p.shutinvolume 133einv=p.endinginventory 134haulvol=p.haulvolume 135hauls=p.hauls 136batteryhauls=p.batteryhauls 137cardz=batteryhauls. sum 138batteryhaulvolume=p.batteryhaulvolume 139 except exc.FailedSolution: 140cstop=time.time 141stop=time.clock 142feasible=False 143obj=np.inf 144 if p.lpobjective is None: 145lpobj=np.inf 146 else : 147lpobj=p.lpobjective 148siprod=None 149einv=None 150haulvol=None 151hauls=None 152batteryhauls=None 153batteryhaulvolume=None 154 print "SI",siprod 155 print "HAUL",haulvol 156 print "POSSIBLE HAUL",hdf["loadsize"] hdf["u"]. sum len periods 157 print "POSSIBLE PROD",bdf["q90"]. sum len periods 158 print "BEGINNING INV",bdf["v0"]. sum 159 print "HAUL CHOICE",cardz 160ptime=stop )]TJ/F18 5.9776 Tf 10.802 0 Td [(start 161cltime=cstop )]TJ/F18 5.9776 Tf 10.801 0 Td [(cstart 162 return einv,siprod,obj,lpobj,ptime,cltime,hauls, n 163haulvol,batteryhauls,batteryhaulvolume 164 165 def run stim mpcbdf,hdf,ddf,periods: 166""" Run a shifting time horizon MPC method on STIM """ 167 168mvbdf=bdf.copy 169window=3 170solvetimes=[] 171simtimes=[] 172shutinvol=[] 173objval=0 174lpobjval=0 175cardz=0 176finalhaulvol=np.zerosbdf.shape[0],hdf.shape[0], len periods 177 for tidx inrange len periods: 178perwin=periods[tidx:tidx+window] 179einv,siprod,obj,lpobj,ptime,cltime,hauls, n 180haulvol,batteryhauls,batteryhaulvolume= n 181run stimmvbdf,hdf,ddf,perwin 182solvetimes.appendcltime 183onedayhauls=batteryhaulvolume[:,:,0] 76

PAGE 86

184cardz+=onedayhauls > 0.astypenp. int sum 185shutinvol.appendsiprod[0] 186finalhaulvol[:,:,tidx]=batteryhaulvolume[:,:,0] 187startsim=time.time 188einvs,haulvs,sis,drys,prods,batinv= n 189run simmvbdf,onedayhauls 190endsim=time.time 191simtimes.appendendsim )]TJ/F18 5.9776 Tf 10.668 0 Td [(startsim 192mvbdf["v0"]=batinv 193 print "SI",shutinvol 194batteryhauls=finalhaulvol > 0.astypenp. int 195haulvol=finalhaulvol. sum axis=0,1 196 print "HAUL",haulvol 197 print "POSSIBLE HAUL",hdf["loadsize"] hdf["u"]. sum len periods 198 print "POSSIblE PROD",bdf["q90"]. sum len periods 199 print "BEGINNING INV",bdf["v0"]. sum 200 print "HAUL CHOICE",cardz 201haullocations=haulvol > 0.astypenp. int 202minhaul=4 hdf["loadsize"]. min /9. 203objval=haulvol. sum )]TJ/F18 5.9776 Tf 10.047 0 Td [(2 sum shutinvol )]TJ/F18 5.9776 Tf 10.273 0 Td [(minhaul cardz 204 return einv,shutinvol,objval,None, sum solvetimes, sum simtimes, n 205haullocations,haulvol,batteryhauls,finalhaulvol 206 207 if name ==" main ": 208 import getopt 209 import time 210 import sys 211 import os.path 212 import anydbm 213 import pickle 214 215usage=" f 0 g )]TJ/F18 5.9776 Tf 5.471 0 Td [(p [goal j stim j btm j nfr] )]TJ/F18 5.9776 Tf 5.47 0 Td [(b < inst number > )]TJ/F18 5.9776 Tf 6.438 0 Td [(i < input path > )]TJ/F18 5.9776 Tf 5.664 0 Td [(o "+ n 216" < output file > 217usage=usage. format sys.argv[0] 218 219opts,args=getopt.getoptsys.argv[1:],"p:b:i:o:" 220 221b=None 222o=None 223i=None 224rungoal=False 225runstoch=False 226runmpc=False 227runbtm=True 228runnfr=False 229 try : 230 for i inrange len opts: 231 if opts[i][0]==" )]TJ/F18 5.9776 Tf 5.47 0 Td [(b": 232b= int opts[i][1] 233 elif opts[i][0]==" )]TJ/F18 5.9776 Tf 5.664 0 Td [(o": 234o=opts[i][1] 235 elif opts[i][0]==" )]TJ/F18 5.9776 Tf 6.438 0 Td [(i": 236ifile=opts[i][1] 237 elif opts[i][0]==" )]TJ/F18 5.9776 Tf 5.47 0 Td [(p": 238 if opts[i][1]=="goal": 77

PAGE 87

239rungoal=True 240 elif opts[i][1]=="stim": 241runstoch=True 242 elif opts[i][1]=="mpc": 243runmpc=True 244 elif opts[i][1]=="nfr": 245runnfr=True 246 except ValueError,e: 247 print usage 248sys.exit3 249 if b is None or o is None or ifile is None: 250 print usage 251sys.exit4 252 ifnot os.path.existso: 253 print usage 254sys.exit5 255 elifnot os.path.existsifile: 256 print usage 257sys.exit6 258 if rungoal or runstoch or runnfr or runmpc: 259runbtm=False 260bdf,hdf,ddf=modeldata.read instanceifile,b 261 if rungoal or runbtm or runnfr: 262u= max hdf["loadsize"] 263l=u/3. 264loads,costs,destinations=modeldata.build btmbdf,ddf,l,u 265 if runbtm: 266haulvol,volcap,objval,lpobjective,clocktime,walltime,sol= n 267run btmloads,costs,destinations 268 elif runnfr: 269haulvol,volcap,objval,lpobjective,clocktime,walltime,sol= n 270run nfrloads,costs,destinations 271 else : 272haulvol,volcap,objval,lpobjective,clocktime,walltime,sol= n 273run goalloads,costs,destinations 274 print "RESULTS: lpobj f 0 g obj f 1 g ". format lpobjective,objval 275db=anydbm. open o,'w' 276db["instance )-16(f 0 g ". format b]=pickle.dumps f "haulvol":haulvol, 277"volcap":volcap,"objectivevalue":objval, 278"lpobjective":lpobjective,"clocktime":clocktime, 279"loadcount": len loads,"destcount": len destinations, 280"walltime":walltime,"solution":sol g 281db.close 282 elif runstoch or runmpc: 283periods= range 1,6 284 if runstoch: 285einv,siprod,obj,lpobj,clocktime,walltime,haulchoice, n 286haulvol,batteryhauls,batteryhaulvolume= n 287run stimbdf,hdf,ddf,periods 288 else : 289einv,siprod,obj,lpobj,clocktime,walltime,haulchoice, n 290haulvol,batteryhauls,batteryhaulvolume= n 291run stim mpcbdf,hdf,ddf,periods 292 print "RESULTS: lpobj f 0 g obj f 1 g ". format lpobj,obj 293db=anydbm. open o,'w' 78

PAGE 88

294db["instance )-16(f 0 g ". format b]=pickle.dumps f "haulvol":haulvol, 295"objectivevalue":obj,"einv":einv,"siprod":siprod, 296"lpobjective":lpobj,"clocktime":clocktime, 297"walltime":walltime,"haulchoice":haulchoice, 298"batteryhauls":batteryhauls, 299"batteryhaulvolume":batteryhaulvolume g 300db.close C.4DataVisualization 1 # )61()]TJ/F47 5.9776 Tf 19.179 0 Td [(coding:utf )]TJ/F47 5.9776 Tf 5.558 0 Td [(8 )61()]TJ/F18 5.9776 Tf -100.56 -11.576 Td [(2 import matplotlib 3matplotlib.use'Agg' 4 import matplotlib.pyplotasplt 5 import numpyasnp 6 7tolerance=0.001 8 9 def plot mipgap vs lplpobj,obj,fn,showoptcond=True: 10fig=plt.figure 11ax=fig.add subplot111 12ax.grid 13fig.set facecolor'w' 14ax.set ylabel"MIP Gap" 15ax.set xlabel"Optimal LPR Objective" 16ax.scatterlpobj,obj )]TJ/F18 5.9776 Tf 5.896 0 Td [(lpobj,color='k',marker='D' 17y0=0. 18x0=0. 19x1=np.ceil max lpobj[~np.isinflpobj]/10. 10 20y1=y0+tolerance x1 )]TJ/F18 5.9776 Tf 10.104 0 Td [(x0 21 if showoptcond: 22ax.plot[x0,x1],[y0,y1],'k )]TJ/F18 5.9776 Tf 6.59 0 Td [(' 23y1=np.ceil max obj/10. 10 24fig.savefigfn 25 return fig 26 27 def plot clockct,fn: 28fig=plt.figure 29fig.set facecolor'w' 30ax=fig.add subplot111 31ax.grid 32ax.set xlabel'Wall Clock Time' 33ax.set ylabel'Instance Count' 34ax.histct,bins=20,color='k' 35fig.savefigfn 36 return fig 37 38 def plot complexitylc,ct,obj,fn: 39fig=plt.figure 40fig.set facecolor"w" 41ax=fig.add subplot111 42ax.grid 43ax.set xlabel"Integer Variables" 44ax.set ylabel"Clock Time" 79

PAGE 89

45 #plotfeasibleinstancesindotandinfeasiblein+ 46lc=np.arraylc 47ct=np.arrayct 48obj=np.arrayobj 49ax.scatterlc[~np.isinfobj],ct[~np.isinfobj],c='k',marker='o' 50ax.scatterlc[np.isinfobj],ct[np.isinfobj],c='k',marker='+' 51fig.savefigfn 52 return fig 53 54 if name ==" main ": 55 import getopt 56 import time 57 import sys 58 import os.path 59 import anydbm 60 import pickle 61 import numpyasnp 62 63usage=" f 0 g [ )]TJ/F18 5.9776 Tf 6.155 0 Td [(n] )]TJ/F18 5.9776 Tf 6.439 0 Td [(i < database file > )]TJ/F18 5.9776 Tf 5.664 0 Td [(o < destination path for plots > 64usage=usage. format sys.argv[0] 65 66opts,args=getopt.getoptsys.argv[1:],"ni:o:" 67 68opath=None 69infile=None 70showcondline=True 71 try : 72 for i inrange len opts: 73 if opts[i][0]==" )]TJ/F18 5.9776 Tf 5.663 0 Td [(o": 74opath=opts[i][1] 75 elif opts[i][0]==" )]TJ/F18 5.9776 Tf 6.438 0 Td [(i": 76infile=opts[i][1] 77 elif opts[i][0]==" )]TJ/F18 5.9776 Tf 5.47 0 Td [(n": 78showcondline=False 79 except ValueError,e: 80 print usage 81sys.exit3 82 ifnot os.path.existsinfile: 83 print usage 84sys.exit6 85 #grabthedataandploteachofthreeplots 86db=anydbm. open infile 87data= fg 88 for k in db.keys: 89 try : 90data[k]=pickle.loadsdb[k] 91 except : 92 continue 93 #nextplottheclocktimehistogram 94ct=np.array[data[k]["walltime"] for k in data] 95fn=opath. format "time hist" 96plot clockct,fn 97 #startwiththemipgapvslprobjective 98 try : 99lpobj=np.array[data[k]["lpobjective"] for k in data] 80

PAGE 90

100 except KeyError: 101 print data[k].keys 102 raise 103obj=np.array[data[k]["objectivevalue"] for k in data] 104fn=opath. format "mip vs obj" 105plot mipgap vs lplpobj,obj,fn,showcondline 106 #lastly,plotthetimevsintegervariables 107ic=np.array[data[k]["loadcount"] data[k]["destcount"] for k in data] 108wt=np.array[data[k]["walltime"] for k in data] 109fn=opath. format "time vs intcount" 110plot complexityic,wt,obj,fn 81

PAGE 91

APPENDIXD.STIMSolutions D.1Solver 1 # )61()]TJ/F47 5.9776 Tf 19.179 0 Td [(coding:utf )]TJ/F47 5.9776 Tf 5.558 0 Td [(8 )61()]TJ/F18 5.9776 Tf -100.56 -11.577 Td [(2 import pulp 3 import numpyasnp 4 5 class STIMProblem object : 6""" Stochastic Inventory Management """ 7 8 def init self,batteries,haulers,destinations,periods, 9caseprob: 10initialinv=[v0 for v0 in batteries["v0"]] 11self. data= f "batteries":batteries, 12"haulers":haulers, 13"destinations":destinations, 14"periods":periods, 15"hrng": range len haulers, 16"brng": range len batteries, 17"drng": range len destinations, 18"caseprob":caseprob, 19"inventories":initialinv g 20 21 def clearModelself: 22self. model= fg 23 24 def buildModelself,gamma: 25""" Create stochastic model """ 26 27self.clearModel 28prob=pulp.LpProblem"STIM",sense=pulp.LpMaximize 29obj=self.buildVarsprob,gamma 30cnum=0 31 for c in self. model["constraints"]: 32prob+=c,"CONS f 0 g ". format cnum 33cnum+=1 34prob+=obj,"stochastic inventory management" 35self. model["objective"]=obj 36self. data["problem"]=prob 37 return prob 38 39 def buildVarsself,prob,gamma: 40""" Create the variables and constraints for the problem """ 41 42batteries=self. data["batteries"] 43haulers=self. data["haulers"] 44periods=self. data["periods"] 45caseprob=self. data["caseprob"] 82

PAGE 92

46initialinv=self. data["inventories"] 47self. model["constraints"]=constraints=[] 48self. model["variables"]=variables=[] 49self. model["shutinprod"]=shutinprod= fg 50self. model["endinginv"]=endinginv=[] 51self. model["haulchoice"]=haulchoice= fg 52self. model["haulervols"]=haulervols=[ fg for idx in haulers.index] 53 for t in periods: 54shutinprod[t]=[] 55haulchoice[t]=[] 56obj=0 57self. model["objval"]=objfcn=pulp.LpVariable"OBJVAL" 58qcases=["q10","q35","q50","q65","q90"] 59scenarios= range len qcases 60xnum=0 61brng=self. data["brng"] 62 for j inrange len haulers: 63 for t in periods: 64haulervols[j][t]=[[None for i in brng] for m in scenarios] 65 #createtheinterperiodconstraintsforeachbattery 66 for i,ii inzip brng,batteries.index: 67v0=pulp.LpVariable"V f 0 g 0". format i 68openinventory= f periods[0] )]TJ/F18 5.9776 Tf 6.623 0 Td [(1:[v0 for m inrange 5] g 69variables.appendv0 70constraints.appendv0==initialinv[i] 71endinginv.append0 72 for t in periods: 73z=pulp.LpVariable"Z f 0 g f 1 g ". format i,t, 74cat=pulp.LpBinary, 75lowBound=0,upBound=1 76openinventory[t]=[] 77variables.appendz 78haulchoice[t].appendz 79obj )]TJ/F18 5.9776 Tf 5.112 0 Td [(=gamma z 80 for m in scenarios: 81v=pulp.LpVariable"V f 0 g f 1 g f 2 g ". format i,t,m, 82lowBound=0 83variables.appendv 84openinventory[t].appendv 85x=[] 86 for j inrange len haulers: 87x1=pulp.LpVariable"X f 0 g ". format xnum, 88 lowBound=0 89xnum+=1 90x.appendx1 91variables.appendx1 92haulervols[j][t][m][i]=x1 93q=batteries.loc[ii,qcases[m]] 94 #tracktheproductionthatexceedstankinventoryandhaul 95 #thisproductionis"shutin" 96o=pulp.LpVariable"O f 0 g f 1 g f 2 g ". format i,t,m, 97lowBound=0,upBound=q 98variables.appendo 99shutinprod[t].appendcaseprob[m] o 100 #modifytheobjective:subtractshutinprod,addhaul 83

PAGE 93

101obj+=caseprob[m] sum x )]TJ/F18 5.9776 Tf 10.285 0 Td [(20000 o 102constr=v==openinventory[t )]TJ/F18 5.9776 Tf 6.54 0 Td [(1][m]+q )]TJ/F18 5.9776 Tf 10.048 0 Td [(o )]TJ/F48 5.9776 Tf 9.289 0 Td [(sum x 103constraints.appendconstr 104constr=v < =batteries.loc[ii,"invcap"] 105constraints.appendconstr 106maxhaul=batteries.loc[ii,"invcap"]+ n 107batteries.loc[ii,qcases[m]] 108constr= sum x < =z maxhaul 109constraints.appendconstr 110 if t==periods[ )]TJ/F18 5.9776 Tf 6.642 0 Td [(1]: 111endinginv[ )]TJ/F18 5.9776 Tf 6.409 0 Td [(1]+=caseprob[m] v 112 for j,idx inzip range len haulers,haulers.index: 113 for t in periods: 114 for m in scenarios: 115constr= sum haulervols[j][t][m] < =haulers.loc[idx,"u"] n 116haulers.loc[idx,"loadsize"] 117constraints.appendconstr 118obj )]TJ/F18 5.9776 Tf 5.111 0 Td [(=endinginv[ )]TJ/F18 5.9776 Tf 6.409 0 Td [(1] 119 #obj= )]TJ/F47 5.9776 Tf 5.75 0 Td [(1 sum[sumshutinprod[t]fortinshutinprod] 120constraints.appendobj==objfcn 121 return obj 122 123 def solveself: 124alpha=1e )]TJ/F18 5.9776 Tf 5.542 0 Td [(3 125haulers=self. data["haulers"] 126minhaul=haulers["loadsize"]. min 127maxhaul=haulers["loadsize"]. max 128bestgamma=None 129bestobj= )]TJ/F18 5.9776 Tf 5.741 0 Td [(1 np.inf 130bestlpobj= )]TJ/F18 5.9776 Tf 5.741 0 Td [(1 np.inf 131 def solvepairgamma: 132prob=self.buildModelgamma 133 #opts=[" )-10()]TJ/F47 5.9776 Tf 12.221 0 Td [(mipgap","0.001"," )110()]TJ/F47 5.9776 Tf 11.207 0 Td [(tmlim","180"] 134 #solvetheLPrelaxationexplicitlyfirst 135opts=[" )184()]TJ/F18 5.9776 Tf 10.215 0 Td [(nomip"] 136 #solver=pulp.GLPK CMDmsg=False,path="/usr/bin/glpsol", 137 #options=opts 138solver=pulp.PYGLPKmsg=False,mip=False 139status=prob.solvesolver 140 if pulp.LpStatus[status]!="Optimal": 141 raise FailedSolutionstatus 142lpobj=self. model["objval"].value 143 #addthecardinalityfixconstraints 144znum=0 145 for t in self. model["haulchoice"]: 146 for z in self. model["haulchoice"][t]: 147 if z.value < alpha: 148prob+=z < =0,"Z )-16(f 0 g ". format znum 149znum+=1 150 #solver=pulp.GLPK CMDmsg=False,path="/usr/bin/glpsol", 151 #options=opts 152solver=pulp.PYGLPKmsg=False,mip=False 153status=prob.solvesolver 154 if pulp.LpStatus[status]!='Optimal': 155 raise FailedSolutionstatus 84

PAGE 94

156objective=self. model["objval"].value 157 return lpobj,objective 158 for gamma in np.linspaceminhaul/4.,2 maxhaul,10: 159self.lpsolution=[] 160lpobj,objective=solvepairgamma 161 if lpobj > bestlpobj: 162bestlpobj=lpobj 163 if objective > bestobj: 164bestgamma=gamma 165bestobj=objective 166 if gamma!=bestgamma: 167 #needtorunthemodelagain 168lpobj,objective=solvepairbestgamma 169self.objective=bestobj 170self.lpobjective=bestlpobj 171vals=[vv.value for vv in self. model["variables"]] 172self.solution=vals 173 return vals 174 175@property 176 def shutinvolumeself: 177vols=[] 178shutinvars=self. model["shutinprod"] 179 for t in self. data["periods"]: 180v=0 181 for svar in shutinvars[t]: 182v+=svar.value 183vols.appendv 184 return vols 185 186@property 187 def endinginventoryself: 188 returnsum [v.value for v in self. model["endinginv"]] 189 190@property 191 def haulvolumeself: 192periods=self. data["periods"] 193volumes=np.zeros len periods,np.float32 194haulervols=self. model["haulervols"] 195caseprob=self. data["caseprob"] 196haulers=self. data["hrng"] 197batteries=self. data["brng"] 198 for j in haulers: 199 for idx,t inzip range len periods,periods: 200 for m inrange len caseprob: 201 for i in batteries: 202vol=haulervols[j][t][m][i] 203volumes[idx]+=caseprob[m] vol.value 204 return volumes 205 206@property 207 def haulsself: 208periods=self. data["periods"] 209haulcount=np.zeros len periods,np.float32 210haulchoice=self. model["haulchoice"] 85

PAGE 95

211 for i in self. data["brng"]: 212 for idx,t inzip range len periods,periods: 213haulcount[idx]+=np.ceilhaulchoice[t][i].value 214 return haulcount 215 216@property 217 def batteryhaulsself: 218""" Report haul indicator by battery by period """ 219 220brng=self. data["brng"] 221periods=self. data["periods"] 222haulchoice=self. model["haulchoice"] 223haulind=np.zeros len brng, len periods 224 for i in self. data["brng"]: 225 for idx,t inzip range len periods,periods: 226haulind[i,idx]=np.ceilhaulchoice[t][i].value 227 return haulind 228 229@property 230 def batteryhaulvolumeself: 231""" Report haul indicator by battery by period """ 232 233brng=self. data["brng"] 234hrng=self. data["hrng"] 235periods=self. data["periods"] 236haulervols=self. model["haulervols"] 237caseprob=self. data["caseprob"] 238haulers=self. data["hrng"] 239batteries=self. data["brng"] 240shp= len brng, len hrng, len periods 241haulvol=np.zerosshp,np.float32 242 for j in haulers: 243 for idx,t inzip range len periods,periods: 244 for m inrange len caseprob: 245 for i in batteries: 246vol=haulervols[j][t][m][i].value 247haulvol[i,j,idx]+=caseprob[m] vol 248 return haulvol D.2Simulation 1 # )61()]TJ/F47 5.9776 Tf 19.179 0 Td [(coding:utf )]TJ/F47 5.9776 Tf 5.558 0 Td [(8 )61()]TJ/F18 5.9776 Tf -100.559 -11.576 Td [(2 3 import anydbm 4 import pickle 5 import matplotlib 6matplotlib.use'Agg' 7 import matplotlib.pyplotasplt 8 import numpyasnp 9 from newcastle.utils import modeldata 10 import os.path 11 12 def get ratemu,sigma: 13""" Get rates distributed normally, don't allow zero """ 86

PAGE 96

14 15rate=np.random.normalmu,sigma,size=mu.shape 16rate[rate < 0]=0 17 return rate 18 19 def get shutin volumerate,v0,invcap,haulvol: 20""" Determine the excess volume from the rate and available sapce """ 21 22sivolume=v0+rate )]TJ/F18 5.9776 Tf 10.588 0 Td [(haulvol )]TJ/F18 5.9776 Tf 10.447 0 Td [(invcap 23sivolume[sivolume < 0]=0 24 return sivolume 25 26 def get haul volumerate,v0,hauldecision,maxhaul,minhaul: 27""" Use the haul flag and available volume to get total haul """ 28 29availablevolume=hauldecision v0+rate 30shortvolume=availablevolume )]TJ/F18 5.9776 Tf 10.103 0 Td [(maxhaul 31shortvolume[shortvolume > 0]=0 32haulvol=maxhaul+shortvolume 33dryload=haulvol < minhaul&hauldecision > 0 34haulvol[haulvol < minhaul]=0 35 return maxhaul+shortvolume,dryload.astypenp. int 36 37 def get ending inventoryrate,v0,haulvolume: 38""" Use the standard relationship to determine ending inventory """ 39 40 return v0+rate )]TJ/F18 5.9776 Tf 10.328 0 Td [(haulvolume 41 42 def simulate haulingmu,sigma,invcap,haulvolume,v0: 43q=get ratemu,sigma 44haulind=haulvolume > 0.astypenp. int 45hv,dryload=get haul volumeq,v0,haulind,haulvolume,45 46sivol=get shutin volumeq,v0,invcap,hv 47actprod=q )]TJ/F18 5.9776 Tf 10.898 0 Td [(sivol 48einv=get ending inventoryq,v0,hv 49einv[einv > invcap]=invcap[einv > invcap] 50 return einv,hv,sivol,dryload,actprod 51 52 def test sim: 53 #testexpectedendinginventory 54cap=np.array[300.,300.] 55haulvol=np.array[180.,185.] 56mu=np.array[30.,60.] 57sigma=0.05 mu 58v0=np.array[170.,200.] 59einv=[] 60m=1000 61 for i inrange m: 62inv,hv,sivol,dryload,actprod= n 63simulate haulingmu,sigma,cap,haulvol,v0 64einv.appendinv 65meaninv=np.arrayeinv.meanaxis=0 66 print "Planned: f 0 g Simulated: f 1 g ". format v0+mu )]TJ/F18 5.9776 Tf 10.488 0 Td [(haulvol,meaninv 67 68 def run simbdf,haulvol,sdfact=0.05: 87

PAGE 97

69""" Execute simulation """ 70 71 try : 72numper=haulvol.shape[2] 73 except IndexError: 74numper=1 75haulvol=haulvol.reshapehaulvol.shape[0],haulvol.shape[1],numper 76periods= range numper 77cap=bdf["invcap"] 78mu=bdf["q50"] 79stddev=sdfact mu 80mrng= range 1000 81sisample=np.zeros1000,numper 82haulsample=np.zeros1000,numper 83drycountsample=np.zeros1000,numper 84prodsample=np.zeros1000,numper 85einvsample=np.zeros1000, 86einvbat=np.zerosbdf.shape[0], 87 for m in mrng: 88inv=bdf["v0"] 89 for t in periods: 90inv,hv,sivol,dryload,actprod= n 91simulate haulingmu,stddev,cap,haulvol[:,:,t]. sum axis=1,inv 92sisample[m,t]= sum sivol 93haulsample[m,t]= sum hv 94drycountsample[m,t]= sum dryload 95prodsample[m,t]= sum actprod 96einvsample[m]= sum inv 97einvbat[:]+=inv/1000 98 return einvsample,haulsample,sisample,drycountsample,prodsample, n 99einvbat 100 101 def plot einvsim,plan,size=20: 102""" Plot sim and plan einv """ 103 104sinv=np.array[sim[k]["einv"].meanaxis=0 for k in sim.keys] 105pinv=np.array[plan[k]["einv"] for k in sim.keys] 106fig1=plt.figure 107fig1.set facecolor"w" 108ax=fig1.add subplot111 109ax.barnp.arange1,size+1,sinv[:size], 110 #hatch="/", 111label="Simulated",color='k' 112ax.bar0.5+np.arange1,size+1,pinv[:size], 113 #hatch=" nn ", 114label="Plan",color='w' 115ax.set xlabel"Instance" 116ax.set ylabel"Barrels" 117ax.set xticks[] 118ax.grid 119ax.legend 120fig2=plt.figure 121fig2.set facecolor"w" 122ax=fig2.add subplot111 123ax.set xlabel"Plan Barrels" 88

PAGE 98

124ax.set ylabel"Simulated Barrels" 125ax.scatterpinv,sinv,facecolors='none',edgecolors='k' 126ax.grid 127 return fig1,fig2 128 129 130 def plot haulsim,plan: 131""" Plot sim vs plan haulvol """ 132 133fig=plt.figure 134fig.set facecolor"w" 135 #flattenthedataseries 136keys=sim.keys 137periodcount=pc=sim[keys[0]]["haul"].shape[1] 138casecount= len keys 139n=casecount periodcount 140simhaul=np.zerosn,,np.float32 141planhaul=np.zerosn,,np.float32 142 for i inrange casecount: 143k=keys[i] 144simhaul[i pc:i pc+pc]=sim[k]["haul"].meanaxis=0 145planhaul[i pc:i pc+pc]=plan[k]["batteryhaulvolume"]. sum axis=0,1 146ax=fig.add subplot111 147ax.set xlabel"Plan Barrels" 148ax.set ylabel"Simulated Barrels" 149ax.grid 150ax.scatterplanhaul,simhaul,facecolors='none',edgecolors='k' 151 return fig 152 153 154 def plot sisim,plan,size=40: 155""" Plot sim and plan si volumes """ 156 157skeys=sim.keys 158N= len skeys 159M=sim.values[0]["si"].shape[0] 160fig1=plt.figure 161fig1.set facecolor"w" 162ax=fig1.add subplot111 163si=np.zerosM,N,np.float32 164plansi=np.array[ sum plan[k]["siprod"] for k in skeys] 165 for n inrange N: 166k=skeys[n] 167si[:,n]=sim[k]["si"]. sum axis=1 168ax.boxplotsi[:,:size] 169ax.plot range 1,size+1,plansi[:size],"k )]TJ/F18 5.9776 Tf 5.435 0 Td [(" 170ax.set xlabel"Instance" 171ax.set ylabel"Barrels" 172ax.set xticks[] 173ax.grid 174fig2=plt.figure 175fig2.set facecolor"w" 176ax=fig2.add subplot111 177ax.set xlabel"Plan Barrels" 178ax.set ylabel"Simulated Barrels" 89

PAGE 99

179ax.scatterplansi,si. max axis=0,facecolors='none',edgecolors='k', 180marker='o',label="max" 181ax.scatterplansi,si.meanaxis=0,facecolors='none',edgecolors='k', 182marker='+',label="mean" 183ax.scatterplansi,si. min axis=0,facecolors='none',edgecolors='k', 184marker='D',label="min" 185ax.grid 186ax.legend 187 return fig1,fig2 188 189 def plot drysim: 190""" Plot sim dry loads""" 191 192fig=plt.figure 193fig.set facecolor"w" 194 195 def make plotsrepo,dest,plandbname,simdbname: 196""" Organize results for plots and generate them """ 197 198output=anydbm. open simdbname 199plan=anydbm. open plandbname 200indata= fg 201simdata= fg 202plandata= fg 203 for k in output.keys: 204junk,instnum=k.split" )]TJ/F18 5.9776 Tf 5.664 0 Td [(" 205simdata[k]=pickle.loadsoutput[k] 206plandata[k]=pickle.loadsplan[k] 207bdf,hdf,ddf=modeldata.read instancerepo,instnum 208indata[k]= f "bdf":bdf,"ddf":ddf,"hdf":hdf g 209 print "Plotting hauls..." 210fig=plot haulsimdata,plandata 211fn=os.path.joindest,"stim sim haul.eps" 212fig.savefigfn 213 print "Plotting ending inventory..." 214fig1,fig2=plot einvsimdata,plandata,size=30 215fn=os.path.joindest,"stim sim einv bar.eps" 216fig1.savefigfn 217fn=os.path.joindest,"stim sim einv scatter.eps" 218fig2.savefigfn 219 print "Plottng si volumes..." 220fig1,fig2=plot sisimdata,plandata,30 221fn=os.path.joindest,"stim sim si boxplot.eps" 222fig1.savefigfn 223fn=os.path.joindest,"stim sim si scatter.eps" 224fig2.savefigfn 225 226 if name ==" main ": 227 import getopt 228 import time 229 import sys 230 231usage=" f 0 g )]TJ/F18 5.9776 Tf 6.439 0 Td [(i < database file > )]TJ/F18 5.9776 Tf 5.664 0 Td [(o < destination path for plots > )]TJ/F18 5.9776 Tf 6.051 0 Td [(r < instance repo > 232usage=usage. format sys.argv[0] 233 90

PAGE 100

234opts,args=getopt.getoptsys.argv[1:],"i:o:r:" 235 236opath=None 237infile=None 238repo=None 239 try : 240 for i inrange len opts: 241 if opts[i][0]==" )]TJ/F18 5.9776 Tf 5.663 0 Td [(o": 242opath=opts[i][1] 243 elif opts[i][0]==" )]TJ/F18 5.9776 Tf 6.438 0 Td [(i": 244infile=opts[i][1] 245 elif opts[i][0]==" )]TJ/F18 5.9776 Tf 6.051 0 Td [(r": 246repo=opts[i][1] 247 except ValueError,e: 248 print usage 249sys.exit3 250 ifnot os.path.existsrepo: 251 print usage 252sys.exit4 253 ifnot os.path.existsinfile: 254 print usage 255sys.exit6 256db=anydbm. open infile 257 if os.path.existsos.path.joinopath,"output.db": 258flg="w" 259 else : 260flg="c" 261output=anydbm. open os.path.joinopath,"output.db",flg 262data= fg 263 for k in db: 264data[k]=pickle.loadsdb[k] 265junk,inst=k.split" )]TJ/F18 5.9776 Tf 5.664 0 Td [(" 266bdf,hdf,ddf=modeldata.read instancerepo,inst 267einvsample,haulsample,sisample,drycountsample,prodsample, n 268batinv=run simbdf,data[k]["batteryhaulvolume"] 269output[k]=pickle.dumps f "einv":einvsample,"haul":haulsample, 270"si":sisample,"dry":drycountsample, 271"prod":prodsample,"batinv":batinv g 272 print ".", 273output.sync 274output.close 91

PAGE 101

APPENDIXE.BTMSolvers E.1DirectSolver 1 # )61()]TJ/F47 5.9776 Tf 19.179 0 Td [(coding:utf )]TJ/F47 5.9776 Tf 5.558 0 Td [(8 )61()]TJ/F18 5.9776 Tf -100.56 -11.577 Td [(2 import pulpaslp 3 import numpyasnp 4 from newcastle.opt.exc import FailedSolution,BTMLogicError 5 6 class BTMProblem object : 7"""Standard BTM Problem 8 9 Arguments: 10 loads )184()]TJETq1 0 0 1 161.519 483.169 cm[]0 d 0 J 0.772 w 0 0 m 0.398 0 l SQq1 0 0 1 161.918 482.982 cm[]0 d 0 J 0.398 w 0 0 m 2.192 0 l SQq1 0 0 1 164.11 483.169 cm[]0 d 0 J 0.772 w 0 0 m 0.398 0 l SQBT/F18 5.9776 Tf 166.195 482.782 Td [(iterable of loads load sizes 11 costs )184()]TJETq1 0 0 1 161.519 471.592 cm[]0 d 0 J 0.772 w 0 0 m 0.398 0 l SQq1 0 0 1 161.918 471.405 cm[]0 d 0 J 0.398 w 0 0 m 2.192 0 l SQq1 0 0 1 164.11 471.592 cm[]0 d 0 J 0.772 w 0 0 m 0.398 0 l SQBT/F18 5.9776 Tf 166.195 471.206 Td [(iterable of iterables assignment costs by load per destination 12 destinations )184()]TJETq1 0 0 1 192.204 460.015 cm[]0 d 0 J 0.772 w 0 0 m 0.398 0 l SQq1 0 0 1 192.603 459.829 cm[]0 d 0 J 0.398 w 0 0 m 2.192 0 l SQq1 0 0 1 194.794 460.015 cm[]0 d 0 J 0.772 w 0 0 m 0.398 0 l SQBT/F18 5.9776 Tf 196.88 459.629 Td [(iterable of capacities 13 14 Keyword arguments: 15 relaxed )184()]TJETq1 0 0 1 170.286 425.285 cm[]0 d 0 J 0.772 w 0 0 m 0.398 0 l SQq1 0 0 1 170.685 425.099 cm[]0 d 0 J 0.398 w 0 0 m 2.192 0 l SQq1 0 0 1 172.877 425.285 cm[]0 d 0 J 0.772 w 0 0 m 0.398 0 l SQBT/F18 5.9776 Tf 174.183 424.899 Td [(True for solving LP relaxation of BTM, False for solving 16 integer problem 17 18 """ 19 20 model=None 21solution=None 22lpobjective=None 23 24 def init self,loads,costs,destinations,relaxed=True: 25self.loads=loads 26self.costs=costs 27self.destinations=destinations 28self. relaxed=relaxed 29self.dcrng= range len self.destinations 30self.lcrng= range len self.loads 31 32 def clearModelself: 33self. model=None 34 35 def buildModelself: 36""" """ 37 38self.clearModel 39 #shapeofvariablesisnxm 40dcrng=self.dcrng 41lcrng=self.lcrng 42prob=lp.LpProblem"BTM",lp.LpMinimize 43constraints,variables,obj,objfcn=self.buildVarsprob 44constraints,variables=self.buildConstraintsprob,constraints, 45 variables 92

PAGE 102

46prob+=obj,"assignment cost minimization" 47self. model= f "variables":variables,"constraints":constraints, 48"objective":obj,"objval":objfcn g 49 return prob 50 51 def buildVarsself,prob: 52constraints=[] 53variables=[] 54 #createavariableforeachload )]TJ/F47 5.9776 Tf 6.251 0 Td [(destinationassignment 55 #buildtheobjectivewhiledefiningvariables 56dcrng=self.dcrng 57lcrng=self.lcrng 58costs=self.costs 59objfcn=lp.LpVariable"OBJVAL",0 60obj=0 61 if self. relaxed: 62vardom=lp.LpContinuous 63 else : 64vardom=lp.LpBinary 65 for i in lcrng: 66loadvars=[None for j in dcrng] 67constr=0 68 for j in dcrng: 69v=lp.LpVariable"V f 0 g f 1 g ". format i,j,cat=vardom, 70lowBound=0,upBound=1 71prob+=v 72loadvars[j]=v 73constr+=v 74obj+=costs[i][j] v 75prob+=constr==1,"FC f 0 g ". format i 76variables.appendloadvars 77constraints.appendconstr 78constraints.appendobjfcn==obj 79prob+=constraints[ )]TJ/F18 5.9776 Tf 6.531 0 Td [(1],"ObjVarSet" 80 return constraints,variables,obj,objfcn 81 82 def buildConstraintsself,prob,constraints,variables: 83dcrng=self.dcrng 84lcrng=self.lcrng 85loads=self.loads 86destinations=self.destinations 87 #buildthecapacityconstraintsatthedestinations 88 for j in dcrng: 89constr=0 90 for i in lcrng: 91v=variables[i][j] 92constr+=loads[i] v 93 #throwtheoverloadvariableontopofthe 94msg="CAP f 0 g ". format j 95prob+=constr < =destinations[j],msg 96constraints.appendconstr 97 return constraints,variables 98 99 def solveself: 100prob=self.buildModel 93

PAGE 103

101opts=[" )184()]TJ/F18 5.9776 Tf 10.27 0 Td [(mipgap","0.001"," )184()]TJ/F18 5.9776 Tf 10.345 0 Td [(tmlim","180"] 102 ifnot self. relaxed: 103newopts=[" )184()]TJ/F18 5.9776 Tf 10.216 0 Td [(nomip"] 104solver=lp.GLPK CMDmsg=True,path="/usr/bin/glpsol", 105options=newopts 106status=prob.solvesolver 107 if lp.LpStatus[status]!="Optimal": 108 raise FailedSolutionstatus 109lpvals=[] 110 for vv in self. model["variables"]: 111lpvals.extend[v.value for v in vv] 112self.lpsolution=lpvals 113self.lpobjective=self. model["objval"].value 114 if self. relaxed and )184()]TJ/F18 5.9776 Tf 10.215 0 Td [(nomip" notin opts: 115opts.append" )184()]TJ/F18 5.9776 Tf 10.216 0 Td [(nomip" 116solver=lp.GLPK CMDmsg=True,path="/usr/bin/glpsol",options=opts 117status=prob.solvesolver 118 if lp.LpStatus[status]!='Optimal': 119 raise FailedSolutionstatus 120vals=[] 121 for vv in self. model["variables"]: 122vals.extend[v.value for v in vv] 123self.solution=vals 124 return vals 125 126@property 127 def objectiveself: 128""" 129 Return the value of the objective without any constraint penalties 130 131 """ 132 133 ifnot self.solution: 134 raise BTMLogicError"No solution has been obtained" 135val=self. model["objval"].value 136 print "OBJECTIVE VALUE: ",val 137 return val 138 139 #@property 140 #defchoicebinsself: 141 #"""Returnthecountsofloadsassignedtotheirpreferreddest""" 142 # 143 #dc=lenself.destinations 144 #lc=lenself.loads 145 #sol=np.arrayself.solution[:lc dc] 146 #sol=sol.reshapelc,dc 147 #returnnp.sumsol,axis=0 E.2GoalProgrammingSolver 1 # )61()]TJ/F47 5.9776 Tf 19.179 0 Td [(coding:utf )]TJ/F47 5.9776 Tf 5.558 0 Td [(8 )61()]TJ/F18 5.9776 Tf -100.56 -11.577 Td [(2 import pulpaslp 3 import numpyasnp 4 from exc import BTMError,BTMLogicError,FailedSolution 94

PAGE 104

5 6 class BTMProblem object : 7 model=None 8solution=None 9 10 def init self,loads,costs,destinations: 11self.loads=loads 12self.costs=costs 13self.destinations=destinations 14 15 def clearModelself: 16self. model=None 17 18 def buildModelself,factors: 19""" """ 20 21self.clearModel 22loads=self.loads 23destinations=self.destinations 24costs=self.costs 25 #shapeofvariablesisnxm 26variables=[] 27constraints=[] 28dcrng= range len self.destinations 29lcrng= range len self.loads 30 #theoverloadvariables 31overload=[lp.LpVariable"O f 0 g ". format j,0 for j in dcrng] 32prob=lp.LpProblem"BTM",lp.LpMinimize 33 #createavariableforeachload )]TJ/F47 5.9776 Tf 6.251 0 Td [(destinationassignment 34 #buildtheobjectivewhiledefiningvariables 35objfcn=lp.LpVariable"OBJVAL",0 36obj=0 37 for i in lcrng: 38loadvars=[None for j in dcrng] 39constr=0 40 for j in dcrng: 41v=lp.LpVariable"V f 0 g f 1 g ". format i,j,cat=lp.LpBinary 42prob+=v 43loadvars[j]=v 44constr+=v 45obj+=costs[i][j] v 46prob+=constr==1,"FC f 0 g ". format i 47variables.appendloadvars 48constraints.appendconstr 49constraints.appendobjfcn==obj 50prob+=constraints[ )]TJ/F18 5.9776 Tf 6.531 0 Td [(1],"ObjVarSet" 51 #buildthecapacityconstraintsatthedestinations 52 for j in dcrng: 53constr=0 54 for i in lcrng: 55v=variables[i][j] 56constr+=loads[i] v 57 #throwtheoverloadvariableontopofthe 58msg="CAP f 0 g ". format j 59prob+=constr < =destinations[j]+overload[j],msg 95

PAGE 105

60obj+=factors[j] overload[j] 61constraints.appendconstr 62prob+=obj,"assignment cost minimization" 63self. model= f "variables":variables,"constraints":constraints, 64"objective":obj,"factors":factors, 65"overload":overload,"objval":objfcn g 66 return prob 67 68 def solveself,factors,nomip=False: 69prob=self.buildModelfactors 70 #alwaysruntheLPRelaxationfirst 71 if 1: 72newopts=[" )184()]TJ/F18 5.9776 Tf 10.216 0 Td [(nomip"] 73solver=lp.GLPK CMDmsg=True,path="/usr/bin/glpsol", 74options=newopts 75status=prob.solvesolver 76 if lp.LpStatus[status]!="Optimal": 77 raise FailedSolutionstatus 78lpvals=[] 79 for vv in self. model["variables"]: 80lpvals.extend[v.value for v in vv] 81self.lpsolution=lpvals 82self.lpobjective=self. model["objval"].value 83opts=[" )184()]TJ/F18 5.9776 Tf 10.27 0 Td [(mipgap","0.001"," )184()]TJ/F18 5.9776 Tf 10.345 0 Td [(tmlim","180"] 84 if nomip and )184()]TJ/F18 5.9776 Tf 10.216 0 Td [(nomip" notin opts: 85opts.append" )184()]TJ/F18 5.9776 Tf 10.216 0 Td [(nomip" 86solver=lp.GLPK CMDmsg=True,path="/usr/bin/glpsol",options=opts 87status=prob.solvesolver 88 if lp.LpStatus[status]!='Optimal': 89 raise FailedSolutionstatus 90vals=[] 91 for vv in self. model["variables"]: 92vals.extend[v.value for v in vv] 93 for v in self. model["overload"]: 94vals.appendv.value 95self.solution=vals 96 return vals 97 98@property 99 def objectiveself: 100""" 101 Return the value of the objective without any constraint penalties 102 103 """ 104 105 ifnot self.solution: 106 raise BTMLogicError"No solution has been obtained" 107val=self. model["objval"].value 108 print "OBJECTIVE VALUE: ",val 109 return val 110 111 class GraduatedOverflowBTMBTMProblem: 112 113 def buildModelself,factors: 114""" """ 96

PAGE 106

115 116self.clearModel 117loads=self.loads 118destinations=self.destinations 119costs=self.costs 120 #shapeofvariablesisnxm 121variables=[] 122constraints=[] 123dcrng= range len self.destinations 124lcrng= range len self.loads 125 #theoverloadvariables 126 #createavariableforeachlevelofoverload )]TJ/F47 5.9776 Tf 10.037 0 Td [(oneload,twoloads... 127factors=np.arrayfactors 128n,m=factors.shape 129nmrng= range n m 130overload=[lp.LpVariable"O f 0 g ". format i,0,150 for i in nmrng] 131prob=lp.LpProblem"BTM",lp.LpMinimize 132 #createavariableforeachload )]TJ/F47 5.9776 Tf 6.251 0 Td [(destinationassignment 133 #buildtheobjectivewhiledefiningvariables 134objfcn=lp.LpVariable"OBJVAL",0 135obj=0 136 for i in lcrng: 137loadvars=[None for j in dcrng] 138constr=0 139 for j in dcrng: 140v=lp.LpVariable"V f 0 g f 1 g ". format i,j,cat=lp.LpBinary 141prob+=v 142loadvars[j]=v 143constr+=v 144obj+=costs[i][j] v 145prob+=constr==1,"FC f 0 g ". format i 146variables.appendloadvars 147constraints.appendconstr 148 #thisconstrainttrackstheactualobjectivevalue 149constraints.appendobjfcn==obj 150prob+=constraints[ )]TJ/F18 5.9776 Tf 6.531 0 Td [(1],"ObjVarSet" 151 #buildthecapacityconstraintsatthedestinations 152 for j in dcrng: 153constr=0 154 for i in lcrng: 155v=variables[i][j] 156constr+=loads[i] v 157 for i inrange m: 158obj+=factors[j,i] overload[j m+i] 159 #throwtheoverloadvariableontopofthe 160msg="CAP f 0 g ". format j 161prob+=constr < =destinations[j]+ n 162 sum [overload[j m+i] for i inrange m],msg 163 #constraints.appendconstr 164prob+=obj,"assignment cost minimization" 165self. model= f "variables":variables,"constraints":constraints, 166"objective":obj,"factors":factors, 167"overload":overload,"objval":objfcn g 168 return prob 97

PAGE 107

E.3NetworkFlowRelaxationSolver 1 # )61()]TJ/F47 5.9776 Tf 19.179 0 Td [(coding:utf )]TJ/F47 5.9776 Tf 5.558 0 Td [(8 )61()]TJ/F18 5.9776 Tf -100.56 -11.576 Td [(2 3 import networkxasnx 4 import numpyasnp 5 from newcastle.opt.exc import FailedSolution,BTMLogicError 6 7scale=1000000 8 9 class BTMProblem object : 10""" NFR Implementation of BTM Solver """ 11 12 def init self,loads,costs,destinations: 13self. model= f "loads":loads,"costs":costs, 14"destinations":destinations, 15"flowcost":None,"flowdict":None, 16"assignment":None,"graph":None g 17self.buildModel 18 19 def clearModelself: 20self. model= f "loads":self. model["loads"], 21"costs":self. model["costs"], 22"destinations":self. model["destinations"], 23"flowcost":None,"flowdict":None, 24"assignment":None,"graph":None g 25 26 def buildModelself: 27"""Calculate the relaxation of delivery capacity and assemble 28 the nodes and edges 29 30 """ 31 32loads=self. model["loads"] 33costs=self. model["costs"] 34destinations=self. model["destinations"] 35lcrng= range len loads 36dcrng= range len destinations 37maxsize= max loads 38destxtra=[] 39destnames=[] 40g=nx.DiGraph 41g.add node"t",demand= len lcrng 42loadcopy=[load for load in loads] 43loadcopy.sort 44 for k in dcrng: 45 #findthelargestnumberofloadsthatfitthedestination 46 #capacityusethoseloadstocalculatetherelaxationfor 47 #thatfacilitysummaxsize )]TJ/F47 5.9776 Tf 10.74 0 Td [(loads[i] 48lo=0 49hi= len loadcopy 50mid=hi+lo//2 51 while lo < mid and mid < hi: 52 ifsum loads[:mid+1] > destinations[k]: 98

PAGE 108

53hi=mid 54 else : 55lo=mid 56mid=hi+lo//2 57xtra=mid+1 maxsize )]TJ/F48 5.9776 Tf 9.289 0 Td [(sum loads[:mid+1] 58destxtra.appendxtra 59 #addanodeforthedestination 60nodename="D )-16(f 0 g ". format k 61destnames.appendnodename 62g.add nodenodename 63g.add edgenodename,"t",weight=0, 64capacity= int np.ceildestinations[k]+xtra/maxsize 65 for i in lcrng: 66nodename="B )-16(f 0 g ". format i 67g.add nodenodename,demand= )]TJ/F18 5.9776 Tf 5.699 0 Td [(1 68 for k in dcrng: 69g.add edgenodename,destnames[k], 70weight= int scale costs[i][k] 71self. model["graph"]=g 72 73 def solveself: 74g=self. model["graph"] 75 if g is None: 76 raise RuntimeError"Model must be built before solving" 77lcrng= range len self. model["loads"] 78dcrng= range len self. model["destinations"] 79 try : 80flowcost,flowdict=nx.network simplexg 81 #flowdict=nx.min cost flowg 82 #flowcost=nx.cost of flowg,flowdict 83 except nx.NetworkXUnfeasible,e: 84 print e 85 raise FailedSolution"Model is infeasible" 86self. model["flowcost"]=flowcost 87self. model["flowdict"]=flowdict 88solution=[flowdict["B )-16(f 0 g ". format u]["D )-16(f 0 g ". format v] n 89 for u in lcrng for v in dcrng] 90 return solution 91 92@property 93 def objectiveself: 94 return self. model["flowcost"]/ float scale 99