MA30087/50087: Optimisation methods of operational research

Brief notes to accompany the lectures.^{1}

Prepared by A. E. Kyprianou^{2} Department of Mathematical Sciences The University of Bath Claverton Down BA2 7AY

Suggested further reading: Elementary Linear Programming With Applications, (Second Edition). B. Kolman and R. E. Beck, Academic Press 1980. ISBN 0-12-417910-X. This book goes at rather a slow pace but spells the theory and practice out very clearly.

Other books:

- Understanding and Using Linear Programming Series: Universitext J. Matoušek, B. Gärtner, Springer (2007).

- B.D. Bunday, Basic Optimization Methods, Hodder Arnold (1984) .

- Bazaran, J.J. Harvis and H.D. Shara’i, Linear Programming and Network Flows, Wiley (2005).

- V. Chvatal, Linear Programming, Freeman (1983).

- D. G. Luenberger, Linear and Nonlinear Programming, Addison-Wesley (1984).

- F. S. Hillier, G. J. Lieberman, Introduction to Operations Research. McGraw-Hill (2001), 7th Edition.

- A. M. Glicksman, An Introduction to Linear Programming and the Theory of Games. Dover.

A very interesting source of material: The internet - try googling ‘linear programming’ and look for lecture notes given at other universities!!

1Corrections are strongly encouraged to be sent to a.kyprianou@bath.ac.uk

2Also partly based on the older notes of M.E. Brigden. Thanks also to M.E. Brigden for supplying me with a comprehensive set of notes with additional examples and exercises.

Contents

1 Introduction: aims of the course 4

2 The linear programming problem 5 2.1 A Simple Example . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Another example . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 A third example . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.4 General statement of the linear programming problem . . . . . 7 2.5 Slack variables and canonical form . . . . . . . . . . . . . . . . 8

3 Geometrical considerations 11 3.1 Constraints, objective functions and hyperplanes . . . . . . . . 11 3.2 Feasible solutions and convexity . . . . . . . . . . . . . . . . . 18 3.3 Extreme points . . . . . . . . . . . . . . . . . . . . . . . . . . 20

4 The fundamental theorem of linear programming 21 4.1 Basic solutions . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.2 The fundamental theorem . . . . . . . . . . . . . . . . . . . . 23

5 The simplex method 26 5.1 Simplex algorithm in action: example 1 . . . . . . . . . . . . . 26 5.2 Simplex algorithm in action: example 2 . . . . . . . . . . . . . 30 5.3 Simplex algorithm in action: example 3 . . . . . . . . . . . . . 32 5.4 The theory behind the simplex method . . . . . . . . . . . . . 34 5.5 Degeneracy and cycling . . . . . . . . . . . . . . . . . . . . . . 36 5.6 The two-phase method . . . . . . . . . . . . . . . . . . . . . . 37 5.7 An example of the two-phase method . . . . . . . . . . . . . . 39

6 Duality 42 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 6.2 Symmetric dual . . . . . . . . . . . . . . . . . . . . . . . . . . 43 6.3 The asymmetric dual . . . . . . . . . . . . . . . . . . . . . . . 44 6.4 Primal to the Dual via the diet problem . . . . . . . . . . . . 45 6.5 The Duality Theorem . . . . . . . . . . . . . . . . . . . . . . . 46 6.6 Complementary slackness . . . . . . . . . . . . . . . . . . . . 54 7 The transportation problem 59

7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 7.2 The canonical form . . . . . . . . . . . . . . . . . . . . . . . . 59 7.3 Dual of the transportation problem . . . . . . . . . . . . . . . 60 7.4 Properties of the solution to the transportation problem … 63 7.5 Solving the transportation problem . . . . . . . . . . . . . . . 65 7.6 Initial feasible solutions . . . . . . . . . . . . . . . . . . . . . . 67 7.7 A worked example . . . . . . . . . . . . . . . . . . . . . . . . 68 7.8 Improvement at each iteration . . . . . . . . . . . . . . . . . . 70 7.9 Degeneracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 7.10 Pricing out . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

8 Optimisation over Networks 77 8.1 Capacitated Networks . . . . . . . . . . . . . . . . . . . . . . 77 8.2 Flows in Networks . . . . . . . . . . . . . . . . . . . . . . . . 78 8.3 Maximal flow problem and labelling algorithm . . . . . . . . . 79 8.4 A worked example . . . . . . . . . . . . . . . . . . . . . . . . 81 8.5 Maximal flow and minimal cuts . . . . . . . . . . . . . . . . . 83 8.6 Example revisited . . . . . . . . . . . . . . . . . . . . . . . . . 86 8.7 A second worked example . . . . . . . . . . . . . . . . . . . . 87

1 Introduction: aims of the course

The principle aim of this course will be to cover the basic ideas and techniques which lie behind modern-day linear programming. One may think of the latter as a field of applied mathematics which concerns itself with resource allocation making use of classical elements of Linear Algebra. Specifically the reader, who should be well versed in an understanding of basic Linear Algebra, will see that many of the results and techniques presented boil down to an application of concepts such as

linear independence and matrix inversion,

convexity in Euclidian spaces and

geometric interpretation of linear equations.

The course is divided into essentially two parts. First we consider the classical linear programming problem and show that there exists a generic method for analysing such problems with the so called the simplex method. In the second part of the course we look at a more complicated class of linear programming problems which come under the heading of the transport problem and optimisation on networks. We consider methods of solution to the latter class of problems which can be thought of as variants of the simplex method.

Some of the proofs in this course are rather lengthy but none are particularly demanding in view of the volume of mathematics that 3rd year students in their first semester will have seen thus far. Although you may be expecting this course to be more ‘applied’ and therfore ‘easier’, please remember that even applications require some degree of intellectual justification for them to be of any value. Therefore the algorithms alone will be meaningless to you unless you make the effort to understand the mathematics behind them.

2 The linear programming problem

2.1 A Simple Example

Let us begin by with an example. Betta Machine Products plc makes two different products. Each requires casting, machining and assembly time. Products are made at a fixed cost and likewise sold at a fixed cost. This information is given in the table below.

Casting Machining Assembly Cost Selling Price
Product 1  5  5  1  25  40
Product 2  8  4  3  30  50

Times are given in minutes per unit. Costs and Prices are in pounds. Each week there are 16,000 mins of casting time, 14,000 mins of machining time and 5,000 mins of assembly time available and there is no limit to numbers that can be sold. The objective of this company is to maximise the difference between total revenue and total cost.

We may now formulate this optimisation problem in a purely mathematical context. Define x_{1} units of product 1 and x_{2} units of product 2 each week. Our objective is to:

maximise z = 15x_{1} + 20x_{2} subject to: 5x_{1} + 8x_{2} ≤ 16, 000 5x_{1} + 4x_{2} ≤ 14, 000 x_{1} + 3x_{2} ≤ 5, 000 x_{1} ≥ 0, x_{2} ≥ 0.

This is a linear programming problem with two variables and three constraints. (Strictly speaking there are in fact five constraints but, as we shall shortly see, positivity of variables is a usual requirement). Many practical problems have literally thousands of variables and constraints.

2.2 Another example

We wish to feed cattle, meeting nutritional requirements at minimum cost. A farmer has two feeds A and B at her disposal. There are certain nutritional requirements which stipulate that cattle should receive minimum quantities of carbohydrate, vitamins and protein. The table below gives the number of units per kilo of the nutrients.

Food  Carbo. Vit. Pro. Cost (pence/kilo)
A  10  2  2  40
B  5  3  8  80
Daily Requirement 10  3  4

The mathematical formulation of the above problem is as follows. Use x_{1} kilos of food A and x_{2} kilos of B per day. Our objective is to:

minimise z = 40x_{1} + 80x_{2} subject to: 10x_{1} + 5x_{2} ≥ 10 2x_{1} + 3x_{2} ≥ 3 2x_{1} + 8x_{2} ≥ 4 x_{1}, x_{2} ≥ 0.

This is an example of a diet problem. Note that where as in the previous example we were asked to maximise a linear function, in this problem we are asked to minimise a linear function. None the less, as we see in the next section, both problems fit into the same framework.

2.3 A third example

The R.H.Lawn Products Co. has available 80 metric tons of nitrate and 50 metric tons of phosphate to use in producing its three types of fertiliser during the coming week. The mixture ratios and profit figures are given in the table below. Determine how the current inventory should be used to maximise profit.

Metric tons / 1000 bags  Profit
Nitrate  Phosphate (dollars / 1000 bags)
Regular lawn 4  2  300
Super lawn 4  3  500
Garden  2  2  400

The above problem is again a maximisation problem. In this case we let x_{1}, x_{2}, x_{3} denote the quantities of the number of bags of the three products in the left hand column of the given table respectively. Our objective is to:

maximise z = 300x_{1} + 500x_{2} + 400x_{3} subject to: 4x_{1} + 4x_{2} + 2x_{3} ≤ 80 2x_{1} + 3x_{2} + 2x_{3} ≤ 50 x_{1}, x_{2}, x_{3} ≥ 0.

Thus we reach a linear programming problem with three variables and two constraints.

2.4 General statement of the linear programming problem

In general we may identify a linear programming problem to take the following standard form. Suppose that A is a matrix belonging to ℝ^{m×n} where m, n ≥ 1, b is a vector belonging to ℝ^{m} and c is a vector belonging to ℝ^{n} and 0_{n} is an n-dimensional vector filled with zeros. Our objective is to find a vector x ∈ ℝ^{n} that will:

maximise z := c · x subject to: Ax ≤ b x ≥ 0_{n}.

(1)

This is a compact way of saying:

maximise z := c_{1}x_{1} + c_{2}x_{2} + ⋯ + c_{n}x_{n} subject to: A_{11}x_{1} + A_{12}x_{2} ⋯ A_{1n}x_{n} ≤ b_{1} A_{21}x_{1} + A_{22}x_{2} ⋯ A_{2n}x_{n} ≤ b_{2}

⋮

A_{m1}x_{1} + A_{m2}x_{2} ⋯ A_{mn}x_{n} ≤ b_{m}

x_{j} ≥ 0, j = 1, 2, ..., n.

Notes.

1. In general there are no restrictions on the dimension of the matrix A so that either m ≤ n or n ≤ m.

2. A ‘minimisation’ problem (i.e. a problem which requires one to minimise a linear function of the form c · x) can easily be converted to a maximisation problem by changing c to −c.

3. A problem which appears to require Ax ≥ b (i.e. the inequality goes the wrong way) can reduced to the standard form by writing it in the form −Ax ≤ −b.

4. A problem which appears to require that Ax = b (i.e. an equality instead of an inequality) can be re-expressed in terms of the requirements that Ax ≤ b and −Ax ≤ −b. However, as we shall see later when we define the canonical form of a linear programming problem, having an equality in place of an inequality is an advantage.

5. Sometimes one may find that instead of x ≥ 0, one has unconstrained values for the vector x. In that case, one may always write x = x^{+} − x^{−} where x^{±} ≥ 0 and then one may re-write the conditions on x as conditions on the concatenation of x^{+} and x^{−}.

We conclude this section with a definition of two types of solutions.

Definition 2.1 Consider the standard form of the linear programming problem (2). We say that x ∈ ℝ^{n} is a feasible solution if it respects all the constraints of the specified problem. If x is feasible and maximises the value of z then we call it an optimal solution. (We will also refer to it on occasion as an optimal feasible solution^{3}).

Note, we have no reason to believe at this point in time that any given linear programming problem has a feasible solution, let alone an optimal feasible solution. Moreover, even if it does have an optimal solution, then it may not necessarily be unique.

2.5 Slack variables and canonical form

If instead of Ax ≤ b we have Ax = b in the formulation of the linear programming problem then we say that it is in canonical form. In other words, a linear programming problem in canonical form requires one to find a vector x ∈ ℝ^{n} that will:

3Sometimes we will drop the word ‘feasible’ as clearly an optimal solution must necessarily be feasible!

maximise z := c · x subject to: Ax = b x ≥ 0_{n}.

(2)

The way to convert a problem from standard form into canonical form is to introduce slack variables. That is to say, if A is a matrix belonging to A^{m×n} then we introduce m new variables, say u_{1}, ..., u_{m} so that the i-th constraint becomes

A_{i1}x_{1} + ⋯ + A_{in}x_{n} + u_{i} = b_{i}

for i = 1, ..., m. If we now introduce the new vector belonging to ℝ^{n+m} which is the concatenation of x and u^{T} = (u_{1}, ..., u_{m}), say

x′ = ( \begin{matrix}x \\ u \end{matrix} ) .

Similarly we introduce the vector c′ which is the concatenation of c and 0_{m},

c′ = ( \begin{matrix}c \\ 0_{m} \end{matrix} )

and define a new ℝ^{m×(m+n)} matrix

A′ = (A|I_{m})

where I_{m} is the m-dimensional identity matrix, then the linear programming problem in standard form can be re-expressed as a linear programming problem in linear form by identifying it as having objective to find a vector x′ in ℝ^{n+m} that will:

maximise z′ := c′ · x′ subject to: A′x′ = b x′ ≥ 0_{m+n}.

As with the standard form of the linear programming problem, we call a feasible solution for the canonical form to be any solution which satisfies Ax = b and x ≥ 0_{n}. Further, an optimal solution for the canonical form is any feasible solution which optimises the value of objective function.

An important observation is that if one is presented with a linear programming problem in standard form, solving the associated linear programming problem in canonical form will produce the same optimal value of the objective. Further, restricting the optimal solution of the canonical problem to the original variables will also give an optimal solution to the problem in standard form. In principle this is all intuitively obvious, however for a proof, see the accompanying exercise sheet.

Let us conclude this chapter with a note on notation. In the above discussion on the introduction of slack variables, we started with n variables and m constraints which meant that we needed to introduce m slack variables. Hence we ended up with a matrix of dimension m × (m + n). However, more often than not, when introducing a linear programming problem already in canonical form for the purpose of theoretical discussion, we shall write the constraints simply as Ax = b where A is a matrix of dimension m×n. That is to say, we shall not bother indicating the number of slack variables that have been introduced (if al all) and use m and n as generic terms for the dimensions of the matrix A as it appears in the given linear programming problem.

3 Geometrical considerations

3.1 Constraints, objective functions and hyperplanes

There is a very natural geometrical interpretation of a system of linear equations which are presented in the form Ax ≤ b where A is a prespecified ℝ^{m×n} matrix, x is an ℝ^{n} vector and b is a prespecified ℝ^{m} vector. Suppose for i = 1, ..., m we write r_{i} as the vector corresponding to the i-th row of A then the system Ax ≤ b may otherwise be written

r_{i} · x ≤ b_{i} for i = 1, ..., m

where b_{i} is the i-th entry of b.

At this point it is instructive to recall that in n-dimensional Euclidian space one may describe any hyperplane in the the form

n · x = d

where n is a vector of unit length orthogonal to the hyperplane and d is the orthogonal distance of the hyperplane from the origin. Indeed the intuition behind this formula is that any vector x which belongs to such a hyperplane can be decomposed into the sum of two vectors. The first vector moves from the origin onto the hyperplane at the closest position to the origin and the second vector moves from the aforementioned point on the hyperplane to x. That is to say

x = dn + y

where obviously the first vector on the right hand side is the shortest straight line from the origin to the hyperplane and the vector y is the relative position of x to the closest point on the hyperplane to the origin and must be orthogonal to n. See Figure 1.

Hence n · x = dn · n + n · y = d + 0. Note in particular that when each of the vectors is taken in two dimensional Euclidian space then

n · x = \binom{n_{1}}{n_{2}} · \binom{x_{1}}{x_{2}} = n_{1}x_{1} + x_{2}x_{2} = d

which is the equation of a line. If we are working in three-dimensional Euclidian space then

n · x = ( \begin{matrix}n_{1} \\ n_{2} \\ n_{3} \end{matrix} ) · ( \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ) = n_{1}x_{1} + n_{2}x_{2} + n_{3}x_{3} = d

y  x n  d

Figure 1: A visial representation of the plane n · x = d.

which is the equation of a plane.

If instead we now consider the equation

n · x ≤ d

then in light of the above intuitive explanation one sees that the latter equation describes a closed half space which lies ‘to one side’ of the hyperplane n · x = d. Indeed, if a vector x satisfies x · n < d then we can imagine that it lies on another hyperplane with the same normal n which is at some other orthogonal distance, say D, from the origin; i.e. n · x = D. Then if, for example, d > 0 and D < d then x lies on a hyperplane parallel to n · x = d, which is contained in the same half-space as the origin is and which is an orthogonal distance D from the origin.

Returning to the linear programming problem in standard form (1) we now see that if the vector x satisfies Ax ≤ b then it must lie within the intersection of a number of closed half spaces whose boundaries are described by the hyperplanes

r_{i} · x ≤ b_{i} for i = 1, ..., m.

Indeed, even the additional positivity constraint x ≥ 0_{n} is also a requirement that x lies to one side of a hyperplane. It says that x lies in the positive orthant of ℝ^{n} which is a quick way of saying that

( −1 ; 0 ; ⋮ ; 0 )·( x_{1} ; x_{2} ; ⋮ ; x_{n} ) ≤ 0, ( 0 ; −1 ; ⋮ ; 0 )·( x_{1} ; x_{2} ; ⋮ ; x_{n} ) ≤ 0, ⋯ and ( 0 ; 0 ; ⋮ ; −1 )·( x_{1} ; x_{2} ; ⋮ ; x_{n} ) ≤ 0.

This interpretation also has meaning when we look at the objective function of a linear programming problem too. Indeed in the standard form we are required to maximise the function c · x. If we write

z = c · x

then we have the usual constraints Ax ≤ b and x ≥ 0_{n} then are required to: find a vector x which belongs to a hyperplane which is orthogonal to c and whose distance from the origin is as large as as possible^{4} so that x is still contained in the region which lies in the intersection of the closed half spaces specified by the constraints. We call the hyperplane given by the objective equation z = c · x the objective hyperplane or objective function.

To some extent it is easier to visualise what is being said above with a diagram. This can be done however only in dimensions 2 and 3 for obvious reasons.^{5} Below are several examples of systems of constraints in either two or three variables together with a plot of the region which lies in the intersection of the closed half spaces specified by the constraints.

Example 3.1 This is a system of constraints in two dimensions,

x_{1} + x_{2} ≤ 5
2x_{1} + x_{2} ≤ 8
x_{1} ≥ 0, x_{2} ≥ 0.

As mentioned earlier, a hyperplane in two dimensions is nothing more than a line. If we take the first constraint and change the inequality to an equality then it reads

x_{2} = 5 − x_{1}.

4Even if this is negative.

5I have not yet found a way of drawing diagrams in 4 or more dimensional space which can be put effectively on paper.

2  4  6 4 2 6 8 x_{2}  x_{1} x_{1} + x_{2} = 5 2x_{1} + x_{2} = 8

Figure 2: The region described by the equations x_{1} + x_{2} ≤ 5, 2x_{1} + x_{2} ≤ 8 and x_{1}, x_{2} ≥ 0.

That is to say it describes the unique line which passes through the points (x_{1}, x_{2}) = (0, 5) and (x_{1}, x_{2}) = (5, 0) (and therefore has gradient −1). Or, said another way, it describes the line which is an orthogonal distance 5/\sqrt{2} for the origin with orthogonal unit vector

\binom{\frac{1}{\sqrt{2}}}{\frac{1}{\sqrt{2}}.}

When we put the inequality back in so that it reads

x_{2} ≤ 5 − x_{1}

then it describes all pairs (x_{1}, x_{2}) which are ‘below’ the line x_{2} = 5 − x_{1}. A similar analysis of the second constraint requires us to consider all points which are ‘below’ the points described by the line x_{2} = 8 − 2x_{1}. Finally taking account of the fact that both variables must be positive, we should only consider points which are ‘to the right’ of the line described by the vertical axis x_{1} = 0 and ‘above’ the line described by the horizontal axis x_{2} = 0. We thus arrive at the representation in Figure 2 for points satisfying the given constraints (shaded in grey).

−3x_{1} + 2x_{2} = 6 x^{1} + x^{2} = 3 2  4 2 4  x_{1} x_{2}

Figure 3: The region described by the equations x_{1} + x_{2} ≥ 3, −3x_{1} + 2x_{2} ≤ 6 and x_{1}, x_{2} ≥ 0.

Example 3.2 A similar analysis of the constraints

x_{1} + x_{2} ≥ 3
−3x_{1} + 2x_{2} ≤ 6
x_{1} ≥ 0, x_{2} ≥ 0

yields Figure 3 for the points satisfying the constraints (shaded in grey).

Example 3.3 Now consider the following system of constraints

6x_{1} + 4x_{2} + 9x_{3} ≤ 36
2x_{1} + 5x_{2} + 4x_{3} ≤ 20
x_{1} ≥ 0, x_{2} ≥ 0, x_{3} ≥ 0.

Take for example the first constraint. Changing the inequality to an equality it may be written in the form 6x_{1} + 4x_{2} + 9x_{3} = 36 or equivalently in the form

( \frac{6}{\sqrt{133}} ; \frac{4}{\sqrt{133}} ; \frac{9}{\sqrt{133}} ) · ( \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ) = \frac{36}{\sqrt{133}}

suggesting that it is a plane whose orthogonal distance from the origin is 36/\sqrt{133} which has a unit normal vector equal to the first vector in the inner product above. It might seem difficult to imagine how one can easily draw this plane, however looking again at the equation of this plane in the form 6x_{1} + 4x_{2} + 9x_{3} = 36 one easily sees that it passes through the points (x_{1}, x_{2}, x_{3}) = (6, 0, 0), (0, 9, 0) and (0, 0, 4) which is enough information to see how this plane appears in the positive orthant x_{1}, x_{2}, x_{3} ≥ 0. This appears as one of the planes in Figure 4. Returning to the original constraint 6x_{1} + 4x_{2} + 9x_{3} ≤ 36, we are required to look for points (x_{1}, x_{2}, x_{3}) which are on the ‘same side’ of the plane 6x_{1} + 4x_{2} + 9x_{3} = 36 as the origin is.

The other plane in the diagram in Figure 4 represents the second constraint. Note that, for example, the constraint x_{1} ≥ 0 corresponds to the requirement that points (x_{1}, x_{2}, x_{3}) lie ‘on the positive side’ of the (x_{2}, x_{3}) plane. All points satisfying the constraints are again shaded in grey.

9 x^{1}  x_{2} x_{3} 5 4 6 10

Figure 4: The region described by the equations 6x_{1} + 4x_{2} + 9x_{3} ≤ 36, 2x_{1} + 5x_{2} + 4x_{2} ≤ 20 and x_{1} ≥ 0, x_{2} ≥ 0, x_{3} ≥ 0.

Example 3.4 In our final example we leave the reader to examine the following system of equations,

x_{1} + x_{2} + x_{3} ≤ 1
2x_{1} + 5x_{2} + 3x_{2} ≤ 4
4x_{1} + x_{2} + 3x_{3} ≤ 2
x_{1} ≥ 0, x_{2} ≥ 0, x_{3} ≥ 0.

We leave it as an exercise to associate them with the sketch in Figure 5.

x^{1}  x_{3}  x_{2}

Figure 5: The region described by the equations x_{1} + x_{2} + x_{3} ≤ 1, 2x_{1} + 5x_{2} + 3x_{2} ≤ 4, 4x_{1} + x_{2} + 3x_{3} ≤ 2 and x_{1} ≥ 0, x_{2} ≥ 0, x_{3} ≥ 0.

Examples 3.1 and 3.2 examples above in the context of a linear programming problem, would require an objective function of the form z = c_{1}x_{1}+c_{2}x_{2}. For a fixed objective value z we thus have an equation of a straight line which passes through the points (x_{1}, x_{2}) = (0, z/c_{2}) and (z/c_{1}, 0). Maximising the value of z constitutes increasing the value of z, thereby moving the line in the direction of the vector

\binom{c_{1}}{c_{2},}

until the value z^{∗} for which there is no intersection between the line and region described by the constraints for all z > z^{∗}.

The second two examples in the context of a linear programming problem would require an objective function taking the form z = c_{1}x_{1} + c_{2}x_{2} + c_{3}x_{3} which is the equation of a plane. Analogously to the two dimensional case, maximising the value of z requires us to move the plane in the direction of its orthogonal such that z increases until we reach a value z^{∗} beyond which there is no longer contact with the region described by the constraints.

How should one understand the canonical form of the linear programming problem in the above geometrical context? The introduction of slack variables increases the dimension of the problem and one looks for a vector x which solves a system of m linear equations in the positive orthant. Therefore any solution, written in terms of the original and slack variables, will lie at the intersection of the m specified hyperplanes. In general as there are less equations than variables, in dimension three and above, there must be an infinite number of intersection points whose intersection with the positive orthant describe the feasible region. This is difficult to see in higher dimensions, but note that two planes in three dimensional space with distinct normal vectors intersect on a line; that is to say common solutions to both planes are points on a line.

Alternatively, since an equation of the form r · x = d may be written as the two constraints r · x ≤ d and r · x ≥ d, then the feasible region of a linear programming problem in canonical form is still the intersection of a finite number of closed half spaces. There is a special definition for the latter.

Definition 3.5 The intersection of a finite number of closed half spaces is a closed polyhedron.

Note that the intersection of a finite number of closed spaces is still closed^{6} which qualifies the use of the word ‘closed’ in the definition.

3.2 Feasible solutions and convexity

For linear programming problems in either canonical or standard form, let us define F, the feasible domain, to be the set of feasible solutions. From the previous section we argued that F consists of the intersection of a finite number of closed half spaces and hence constitutes a closed polyhedron.

In this section our objective is to give a mathematical description of the set F. Indeed we shall see that F is a closed convex polyhedron. To this end, let us recall the notion of a convex set in ℝ^{n}.

Definition 3.6 A non-empty set S in ℝ^{n} is a convex set if for all x_{1}, x_{2} ∈ S and λ ∈ (0, 1) it holds that

λx_{1} + (1 − λ)x_{2} ∈ S.

6Suppose that C_{1}, ..., C_{n} are closed spaces. If {x_{n} : n ≥ 1} is a sequence of points in ⋂_{i=1}^{n} C_{i} with an accumulation point, say x. Then for each i = 1, ..., n we have {x_{n} : n ≥ 1} is a sequence of points in C_{i}. Since, for each i = 1, ..., n, we are given that C_{i} is closed then it follows that x ∈ C_{i}. Hence x ∈ ⋂_{i=1}^{n} C_{i} and thus ⋂_{i=1}^{n} C_{i} is also closed.

Roughly speaking what this definition means that S is a convex set if, when choosing any two points in that set, the line joining those points is also contained in S. Here are some relevant examples of convex sets.

Example 3.7 Suppose that for some given vector n in ℝ^{n} and constant d ≥ 0, S = {x ∈ ℝ^{n} : n · x = d} (i.e. S is a hyperplane in ℝ^{n}). Then S is a convex set. To see why note that if x_{1}, x_{2} ∈ S then since for λ ∈ (0, 1)

n · (λx_{1} + (1 − λ)x_{2}) = λn · x_{1} + (1 − λ)n · x_{2}
= λd + (1 − λ)d
= d

it follows that λx_{1} + (1 − λ)x_{2} ∈ S.

Example 3.8 Suppose now that we take S = {x ∈ ℝ^{n} : n · x ≤ d} where the quantities n and d are as in the previous example (i.e. S is the closed half space with boundary given by the hyperplane in the previous example). Then again S is convex. In this case the proof of this fact is almost the same as before. Indeed as before we may write

n · (λx_{1} + (1 − λ)x_{2}) = λn · x_{1} + (1 − λ)n · x_{2}
≤ λd + (1 − λ)d
= d

showing that λx_{1} + (1 − λ)x_{2} ∈ S whenever x_{1}, x_{2} ∈ S and λ ∈ (0, 1).

Example 3.9 Suppose that S_{1}, ..., S_{n} are convex sets in ℝ^{m} such that ⋂_{i=1}^{m} S_{i} ≠ ∅. Then ⋂_{i=1}^{m} S_{i} is also a convex set. In this case, one supposes that x_{1}, x_{2} ∈ ⋂_{i=1}^{n} S_{i} then in particular x_{1}, x_{2} ∈ S_{i} for each i = 1, 2, ..., m and hence by convexity of each of the latter sets, λx_{1} + (1 − λ)x_{2} ∈ S_{i} for any λ ∈ (0, 1) and i = 1, 2, ..., m. Consequently, λx_{1} + (1 − λ)x_{2} ∈ ⋂_{i=1}^{n} S_{i} as required.

The last three examples can be used to prove in particular the following result.

Theorem 3.10 Consider either the standard or the canonical form of the linear programming problem, (1) and (2). Then the set F of feasible solutions is a closed convex polyhedron.

Proof. By the examples above and the earlier observation that F is the intersection of a finite number closed half spaces, the statement of the theorem follows.

3.3 Extreme points

Let us conclude this chapter on the geometry of the linear programming problem by hinting at where we should expect to find optimal solutions in the feasible region F.

Definition 3.11 Given a non-empty convex set S, we say that x ∈ S is an extreme point in S if x is not an interior point of any line segment in S. That is to say, there exists no two vectors x_{1}, x_{2} ∈ S and λ ∈ (0, 1) such that

x = λx_{1} + (1 − λ)x_{2}.

A good example to visualise extreme points is to think of the vertices of any Platonic Solid or indeed any n-dimensional simplex^{7}. The following result is intuitively obvious and we exclude its proof.

Lemma 3.12 Any closed convex polyhedron generated by the intersection of a finite number of half spaces (and hence any feasible region) has at most a finite number of extreme points.

The importance of extreme points is that one should look for optimal solutions to linear programming problems by considering the extreme points of the feasible region as the following theorem confirms. The proof is omitted.

Theorem 3.13 Consider a linear programming problem either in standard or canonical form. Let F be the space of feasible solutions. Then one of the following three scenarios necessarily holds.

(i) If F ≠ ∅ and is bounded^{8} and there is an optimal solution which occurs at an extreme point.

(ii) If F ≠ ∅ and not bounded, but an optimal solution exists which occurs at an extreme point.

(iii) There exists no optimal solution in which case F is either unbounded or F = ∅.

7Look them up on Google!

8One may understand bounded in this context to mean that F ⊂ B(r) where B(r) = {x ∈ ℝ^{n} : ||x|| < r} for some 0 < r < ∞. In other words, there exists a big enough hyper-sphere of radius r, centred at the origin such that F fits inside this sphere.

4 The fundamental theorem of linear programming

4.1 Basic solutions

Recall that, by the introduction of slack variables, the linear programming problem can be expressed in canonical form. That is to say one has a given matrix A ∈ ℝ^{m×n} and given vector b ∈ ℝ^{m} and the objective is to find a vector x ∈ ℝ^{n} that will:

maximise z := c · x subject to: Ax = b x ≥ 0_{n}.

(3)

Let us assume that c_{1}, ..., c_{n} are the columns of A so that we may write the constraint Ax = b as

x_{1}c_{1} + ⋯ + x_{n}c_{n} = b

where x_{i} is the i-th entry in the vector x.

Note we may assume without loss of generality that m ≤ n and that rank(A) = m. Otherwise when m > n there are more linear equations than unknowns and by a series of linear row operations, one may use row reduction^{9} to decrease the constraints Ax = b to a system of linear equations whose rank can be no greater than the number of unknowns. Also if rank(A) < m then one may again reduce the constraints to a linear system where the rank is equal to the number of equations.

The assumption that rank(A) = m implies that there are precisely m columns of A which are linearly independent. For the sake of notational convenience it will be henceforth assumed without loss of generality that the last m columns of A are the ones which are linearly independent. Under such circumstances we may write A = (A^{∘}|B) where B ∈ ℝ^{m×m} and rank(B) = m so that in particular B is invertible.

Definition 4.1 Consider the canonical form of the linear programming problem (3). A basic solution to Ax = b is a vector x such that if I = {i : x_{i} ≠ 0} (that is to say the set of indices for which x has a non-zero entry) then the columns of A corresponding to the index set I are linearly independent. Note that necessarily |I| ≤ rank(A). If x is a basic solution and x ≥ 0_{n} then we say that x is a basic feasible solution.

9Row reduction means taking linear combinations of the rows of A to produce as many rows of zeros as possible. Recall that the latter is only possible if there is linear dependence between rows.

Another way of understanding what a basic solution is as follows. Suppose we decompose the matrix A = (A^{∘}|B) as above and consider solutions of the form

x = \binom{0_{n−m}}{x_{B}}  (4)

where x_{B} ∈ ℝ^{m}. In particular Bx_{B} = b and since B is invertible, then x_{B} is uniquely identifiable as x_{B} = B^{−1}b. Here x is a basic solution.

In the next theorem it is proved that basic solutions are necessarily extreme points. This conclusion suggests that since optimal solutions occur at extreme points, one may try to find an optimal solution in the class of basic solutions. The reader should naturally be concerned that the optimal solution may of course occur at an extreme point which is not a basic solution. However we shall dispel this concern in the next section. For now we conclude this section with the aforementioned theorem.

Theorem 4.2 Suppose that x is a basic feasible solution of the linear programming problem (3). Then x is an extreme point of the feasible region.

Proof. Without loss of generality we may arrange the given basic solution in the form (4). Let us assume for contradiction that x is not an extreme point so that we may write it as

x = λv + (1 − λ)w

where λ ∈ (0, 1) and v and w belong to the feasible region of the linear programming problem (3). Let us write in the usual way v_{j} and w_{j} for the j-th elements of the vectors v and w. Since x_{1}, .., x_{n−m} are all equal to zero, it follows that for j = 1, ..., n − m,

0 = λv_{j} + (1 − λ)w_{j}.

Since λ ∈ (0, 1) and v_{j}, w_{j} ≥ 0 by assumption, it follows that v_{j} = w_{j} = 0 for j = 1, ..., n − m.

Recalling that v is a feasible solution, we may now write it in the form

v = \binom{0_{n−m}}{v_{B}}

where Bv_{B} = b. Consequently we have

0_{m} = Bv_{B} − Bx_{B} = (v_{n−m+1} − x_{n−m+1})c_{n−m+1} + ⋯ + (v_{n} − x_{n})c_{n}

and hence by the linear independence of c_{n−m+1}, ..., c_{n} it follows that v_{j} = x_{j} for j = n − m + 1, ..., n. However this contradicts the assumption that x was not an extreme point (i.e. the assumption that λ ∈ (0, 1)).

4.2 The fundamental theorem

Thus far we have seen that a bounded optimal solution, when it exists, can always be found at an extreme point and that all basic solutions are extreme points. This makes the class of basic solutions a good class in which to look for solutions. As mentioned earlier however, one may worry that the optimal solution may be an extreme point (or indeed any other point) which cannot be expressed as a basic solution. The fundamental theorem dismays this concern.

Theorem 4.3 Consider the linear programming problem in canonical form (3). Then if there is a finite optimal feasible solution then there is an optimal basic feasible solution. (In other words, by Theorem 4.2, an extreme point corresponding to a finite optimal solution can be chosen to be a basic feasible solution).

Proof. Suppose that x is an optimal feasible solution and without loss of generality we may write it in the form

x = \binom{0_{n−k}}{y}

where n − k is the number of zeros entries in x and y > 0_{k}. Next partition A = (H|G) where the columns of H are equal to the first n − k columns of A and the columns of G are equal to the last k columns of A. Consequently we have

Ax = Gy = b.

If it so happens that the columns of G are linearly independent (in which case it must necessarily be the case that k ≤ m)^{10} and the proof is complete.

Assume then that the columns of G are not linearly independent so that

α_{1}c_{n−k+1} + ⋯ + α_{k}c_{n} = 0_{m}  (5)

implies that for at least one j = 1, ..., k, we have α_{j} ≠ 0. Without loss of generality we can assume that α_{j} > 0.^{11} Now writing the equation Gy = b in the form

y_{1}c_{n−k+1} + ⋯ y_{k}c_{n} = b,

since (5) holds, we have for any ε > 0 that

(y_{1} − ϵα_{1})c_{n−k+1} + ⋯ + (y_{k} − ϵα_{k})c_{k} = b.

Now let us define a new vector

x′_{ϵ} = x − ϵα

where the first n − k entries of α are zero and the remaining k entries are equal to α_{1}, ..., α_{k} respectively. It follows that Ax′_{ϵ} = b.

Note that for ϵ sufficiently small, it can be arranged that all the entries of x′_{ϵ} are positive. Indeed let

δ_{1} = max_{α_{j}<0} \frac{y_{j}}{α_{j}} or − ∞ if all α_{j} ≥ 0

δ_{2} = min_{α_{j}>0} \frac{y_{j}}{α_{j}}

and note that it is always true that δ_{1} < 0 < δ_{2} and that x′_{ϵ} ≥ 0 for all ϵ ∈ [δ_{1}, δ_{2}]. In other words, x′_{ϵ} ≥ 0 is a feasible solution for ϵ ∈ [δ_{1}, δ_{2}]. On account of the optimality of x we have for ϵ ∈ [δ_{1}, δ_{2}] that

c · x ≥ c · x′_{ϵ} = c · x − ϵc · α.

In other words

ϵc · α ≥ 0

10Otherwise one would have k > m vectors in ℝ^{m} which necessarily must have linear dependency.

11Note that at least one of the α_{j} ≠ 0. If that α_{j} < 0 then multiplying (5) by −1 then relabelling each α_{i} by −α_{i} we have the existence of at least one α_{j} > 0.

for all ϵ ∈ [δ_{1}, δ_{2}]. This can only happen if c · α = 0 and so for ϵ ∈ [δ_{1}, δ_{2}] we have shown that x′_{ϵ} is also an optimal feasible solution! Now choose ϵ = δ_{2} and note that by doing so, the number of zero terms in the vector x′_{δ_{2}} has increased by at least one over and above the number of zeros in x.

We may now take return to the very beginning of this proof working with x′_{δ_{2}} instead of x and apply the same reasoning to deduce that either we have found an optimal basic feasible solution, or we can find another optimal feasible solution which has at least one more zero in its entries. In the latter case, an iterative procedure will eventually bring us to an optimal basic feasible solution.

Let us conclude this chapter by reiterating the main drive of the discussion.

For a given linear programming problem, an optimal solution (if it exists) will occur at an extreme point and by the fundamental theorem of linear programming this extreme point may be chosen to correspond to a basic feasible solution. We therefore need a mechanism which looks systematically through the basic feasible solutions in order to identify one which gives an optimal value.

Thus enter the simplex method!

5 The simplex method

In 1947 George Dantzig carried out an analysis of military programming and realised that many of the activities of large organisations could be viewed as linear programming problems. During a major international economic conference in 1949, Dantzig presented a method of solution for these problems which became known as the simplex method. It also became even clearer at the conference how wide were the possible applications.

The simplex method is an iterative procedure for solving high-dimensional^{12} linear programming problems. The simplex method gives an algorithm by which one may search for an optimal basic feasible solution by proceeding from a vertex of the feasible region to an adjacent vertex, moving in a direction that improves the objective function.

The algorithm is quite straight-forward, but it is quite difficult to describe in a precise way, mainly because of the notation required to cope with different possible choices of basic feasible solutions and how to move between them.

We shall present the algorithm in this chapter, initially with three concrete examples, followed by some theoretical considerations for a general set-up and conclude with some discussion on degenerate solutions and extra care that one needs to take in certain situations.

5.1 Simplex algorithm in action: example 1

Let us consider the following explicit problem. Our objective is to:

maximise z = 5x_{1} + 6x_{2} + 4x_{3} subject to: x_{1} + x_{2} + x_{3} ≤ 10 3x_{1} + 2x_{2} + 4x_{3} ≤ 21 3x_{1} + 2x_{2} ≤ 15 x_{1} ≥ 0, x_{2} ≥ 0, x_{3} ≥ 0.

12High dimensional essentially means of dimension three or above since these problems cannot be illustrated graphically.

Step 1. Let us introduce slack variables u_{1}, u_{2}, u_{3} allowing us to replace inequalities by equalities so that the constraints become

\begin{aligned}x_{1} + x_{2} + x_{3} & +u_{1} & & & = 10 \\ 3x_{1} + 2x_{2} + 4x_{3} & & +u_{2} & & = 21 \\ 3x_{1} + 2x_{2} & & & +u_{3} & = 15 \end{aligned}

and the new variables u_{1} ≥ 0, u_{2} ≥ 0, u_{3} ≥ 0. We may now take as our first basic feasible solution

(x_{1}, x_{2}, x_{3}, u_{1}, u_{2}, u_{3}) = (0, 0, 0, 10, 21, 15).

Take a moment and look back at the previous chapter to see why this is a basic feasible solution! Note that with this choice of basic feasible solution, our objective function z = 0. The reason for this is that our basic feasible solution has not included any of the variables x_{1}, x_{2}, x_{3} and clearly including at least one of these variables into the basic feasible solution would increase the objective function from it current value.

Step 2. Recall the conclusion of the previous chapter, that it suffices to look through all the basic feasible solutions in order to find the solution to a linear programming problem when it is finite. Further from the previous step, we can already see that moving to another basic feasible solution which includes one of the variables x_{1}, x_{2}, x_{3} would increase the objective function. Keeping the constraints tight, the simplest thing to do is consider introducing the variable x_{1} into the basis and removing one of the basis variables u_{1}, u_{2} or u_{3}. From the three main constraints, our options are:

(i) x_{1} : 0 → 10 and u_{1} : 10 → 0
or (ii) x_{1} : 0 → 7 and u_{2} : 21 → 0
or (iii) x_{1} : 0 → 5 and u_{3} : 15 → 0.

Whichever of these options we choose, we must respect the feasibility of the solution. Below we see what the values of the basis variables would become in the three scenarios above:

(i) (x_{1}, x_{2}, x_{3}, u_{1}, u_{2}, u_{3}) = (10, 0, 0, 0, −9, −15)
(ii) (x_{1}, x_{2}, x_{3}, u_{1}, u_{2}, u_{3}) = (7, 0, 0, 3, 0, −6)
(iii) (x_{1}, x_{2}, x_{3}, u_{1}, u_{2}, u_{3}) = (5, 0, 0, 5, 6, 0).

Clearly only the third case respects the constraints of non-negativity and so we move to a new basis containing the variables (x_{1}, u_{1}, u_{2}).

Step 3. The next step is to reformulate the linear programming problem in terms of the new basis variables so that they appear only once in each constraint and have unit coefficients. A little linear algebra^{13} gives us the following

\frac{1}{3}x_{2} +x_{3} +u_{1}  −\frac{1}{3}u_{3} = 5
4x_{3}  +u_{2} −u_{3} = 6
x_{1} +\frac{2}{3}x_{2}  \frac{1}{3}u_{3} = 5.

Further, by writing the objective function in terms of the non-basic elements, we see that

z = 5(5 − \frac{1}{3}u_{3} − \frac{2}{3}x_{2}) + 6x_{2} + 4x_{3} = 25 − \frac{5}{3}u_{3} + \frac{8}{3}x_{2} + 4x_{3}.

This tells us immediately two things. Firstly that with the new basis, the objective function is valued at z = 25. Secondly, if we are to change to a new basic feasible solution, it would only make sense to introduce the variables x_{2} or x_{3} as otherwise introducing u_{3} would decrease the value of z.

Step 4. Let us now consider the two alternatives of introducing x_{3} or x_{2} into the basis.

If we introduce x_{3} then our options are:

(i) x_{3} : 0 → 5 and u_{1} : 5 → 0
or (ii) x_{3} : 0 → \frac{3}{2} and u_{2} : 6 → 0
or (iii) x_{3} ∈ ℝ  and x_{1} = 5.

In order to respect the positivity of all basis variables as well as maximising the contribution to the objective function, the best choice from the above is (ii). In that case z increases by 6 units.

If we introduce x_{2} then our options are:

(i) x_{2} : 0 → 15 and u_{1} : 5 → 0
or (ii) x_{2} ∈ ℝ  and u_{2} = 6
or (iii) x_{2} : 0 → \frac{15}{2} and x_{1} : 5 → 0.

In order to respect the positivity of all basis variables as well as maximising the contribution to the objective function, the best choice from the above is (iii). In that case z increases by 20.

13We did row 1 - \frac{1}{3} row 3 and row 2 - row 1.

The conclusion of these calculations is that we should move to a new basic feasible solution

(x_{1}, x_{2}, x_{3}, u_{1}, u_{2}, u_{3}) = (0, \frac{15}{2}, 0, \frac{5}{2}, 6, 0).

giving value z = 45 to the objective function.

Step 5. Again writing the constraints so that the basic variables appear only once in each constraint with coefficient 1, and in particular keeping x_{2} in the third row, we have^{14}:

−\frac{1}{2}x_{1}  +x_{3} +u_{1}  −\frac{1}{2}u_{3} = \frac{5}{2}
4x_{3}  +u_{2} −u_{3} = 6
\frac{3}{2}x_{1} +x_{2}  +\frac{1}{2}u_{3} = \frac{15}{2}

Further, writing z in terms of the non-basic variables we obtain

z = 15 − \frac{5}{3}u_{3} + \frac{8}{3}(\frac{15}{2} − \frac{8}{2}x_{1} − \frac{1}{2}u_{3}) + 4x_{3} = 45 − 4x_{1} + 4x_{3} − 3u_{3}.

The only way we can improve on the value of the objective function is to include the variable x_{3} into the basis and this becomes the final step of the algorithm.

Step 6. Introducing x_{3}, our options are:

(i) x_{3} : 0 → \frac{5}{2} and u_{1} : \frac{5}{2} → 0
or (ii) x_{3} : 0 → \frac{3}{2} and u_{2} : 6 → 0
or (iii) x_{3} ∈ ℝ  and x_{2} = \frac{15}{2}.

The first option would not work as it would force u_{2} to take a negative value. The second and third options are essentially the same when one takes account of the need for positive values. In conclusion the new basic feasible solution

(x_{1}, x_{2}, x_{3}, u_{1}, u_{2}, u_{3}) = (0, \frac{15}{2}, \frac{3}{2}, 1, 0, 0).

Increasing the value of z by 6 to 51.

Step 7. A quick calculation allows us to see that

z = 51 − 4x_{1} − 2u_{3} − u_{2}

14We did row operations: row1 → row1 - \frac{1}{2}row3 and row3 → \frac{3}{2}row3.

and there is no point to introduce any further variables into the basis as this would only decrease the value of the objective problem.

Step 8. The fundamental theorem of linear programming tells us that we must have found the optimal solution since we have found the optimal value of z amongst all the basic feasible solutions.

5.2 Simplex algorithm in action: example 2

In this second example, we present our computations in a more sophisticated way by making use of a standard tool called a tableau. Let us briefly introduce the latter. Suppose that we are given a linear programming problem with m equations in n variables in standard form which we then extend to canonical form with the introduction of slack variables which form the basis of an initial basic feasible solution. We can represent all this information in a tableau of the form given below.

x_{1} x_{2} ⋯ x_{n} u_{1} u_{2} ⋯ u_{m}
u_{1} a_{11} a_{12} ⋯ a_{1n} 1 0 ⋯ 0 b_{1}
u_{2} a_{21} a_{22} ⋯ a_{2n} 0 1 ⋯ 0 b_{2}
⋮  ⋮  ⋮  ⋮  ⋮  ⋮ ⋮ ⋮  ⋮ ⋮
u_{m} a_{m1} a_{m2} ⋯ a_{mn} 0 0 ⋯ 1 b_{m}
−c_{1} −c_{2} ⋯ −c_{n} 0 0 ⋯ 0 0

The left most column describes the basis variables and the bottom row gives the negative of the coefficients of the objective function in terms of the nonbasic variables. The bottom right entry is the value of the objective function with the current choice of basis. The rest is self-explanatory.

Let us now state an example of a linear programming problem and, following similar logic to the example given in the previous section, we shall produce a series of tableaus which describe how the optimal solution is obtained.

Consider now the linear programming problem. Our objective is to:

maximise z = 8x_{1} + 9x_{2} + 5x_{3} subject to: x_{1} + x_{2} + 2x_{3} ≤ 2 2x_{1} + 3x_{2} + 4x_{3} ≤ 3 6x_{1} + 6x_{2} + 2x_{3} ≤ 8 x_{1} ≥ 0, x_{2} ≥ 0, x_{3} ≥ 0.

Moving to the tableau for the canonical form we have

x_{1} x_{2} x_{3} u_{1} u_{2} u_{3}
u_{1} 1 1 2 1 0 0 2
u_{2} 2 3 4 0 1 0 3
u_{3} 6 6 2 0 0 1 8
−8 −9 −5 0 0 0 0

Following the logic of the previous example, we note that non-basic variables which have a negative coefficient in the final row will contribute a positive increase in the value of the objective if introduced to the basis. For each such variable x_{j} we compute θ_{j} = min_{i}{b_{i}/a_{ij} : b_{j}/a_{ij} ≥ 0}. Of these x_{j} introduce the one for which o_{j}θ_{j}, the increase in the objective, is the greatest where o_{j} is the coefficient of x_{j} when the objective function is written in terms of the non-basic variables and x_{j}. (In other words o_{j} is minus the entry in the j-th column of the last row and o_{j}).

Pivoting about the entry in bold in the above tableau so that x_{1} is introduced into the basis and u_{3} is removed from the basis (note that one will have chosen this change of basis variable on by considering feasibility) one obtains the following.

x_{1} x_{2}  x_{3} u_{1} u_{2}  u_{3}
u_{1} 0 0 5/3 1 0 −1/6 2/3
u_{2} 0 1 10/3 0 1 −1/3 1/3
x_{1} 1 1 1/3 0 0 1/6 4/3
0 −1 −7/3 0 0 4/3 32/3

Pivoting again about the entry shown in bold we obtain:

x_{1} x_{2}  x_{3} u_{1} u_{2}  u_{3}
u_{1} 0 0 5/3 1 0 −1/6 2/3
x_{2} 0 1 10/3 0 1 −1/3 1/3
x_{1} 1 0 −3 0 −1 1/2 1
0 0  1 0 1  1 11

Since all the coefficients of the non-basic variables in the bottom row are now positive, this means that an increase in their value would cause a decrease in the value of the objective function which is now identifiable as

z = 11 − x_{3} − u_{2} − u_{3}.

As an exercise, go back to the first example and see if you can summarise the key steps in a sequence of tableaus too. Alternatively consider the linear programming problem in the next section whose solution is given simply as a sequence of tableaus.

5.3 Simplex algorithm in action: example 3

Here is a third example of the simplex algorithm where we discover at the end that there is in fact no finite optimal solution. The problem is to:

maximise z = 2x_{1} + 3x_{2} + x_{3} + x_{4} subject to: x_{1} − x_{2} − x_{3} ≤ 2 −2x_{1} + 5x_{2} − 3x_{3} − 3x_{4} ≤ 10 2x_{1} − 5x_{2} + 3x_{4} ≤ 5 x_{1} ≥ 0, x_{2} ≥ 0, x_{3} ≥ 0, x_{4} ≥ 0.

Introducing slack variables u_{1}, u_{2}, u_{3} we move straight to to the first tableau and from which we produce two further tableaus. As usual we have pivoted about the entry indicated in bold. The details are left to the reader.

x_{1} x_{2} x_{3} x_{4} u_{1} u_{2} u_{3}
u_{1} 1 −1 −1 0 1 0 0 2
u_{2} −2 5 −3 −3 0 1 0 10
u_{3} 2 −5 0 3 0 0 1 5
−2 −3 −1 −1 0 0 0 0

x_{1} x_{2}  x_{3}  x_{4} u_{1} u_{2} u_{3}
u_{1}  3/5 0 −8/5 −3/5 1 1/5 0 4
x_{2} −2/5 1 −3/5 −3/5 0 1/5 0 2
u_{3} 0 0 −3 0 0 1 1 15
−16/5 0 −14/5 −14/5 0 3/5 0 6

x_{1} x_{2}  x_{3} x_{4}  u_{1} u_{2} u_{3}
x_{1} 1 0 −8/3 −1 5/3 1/3 0 20/3
x_{2} 0 1 −5/3 −1 2/3 1/3 0 14/3
u_{3} 0 0 −3 0 0 1 1 15
0 0 −34/3 −6 16/3 5/3 0 82/3

In this last tableau all the coefficients of x_{3} are negative and hence it would be easy to create a solution with the four variables, x_{1}, x_{2}, x_{3} and u_{3} by setting x_{3} = M > 0 and then taking

x_{1} = \frac{20}{3} + \frac{8}{3}M
x_{2} = \frac{14}{3} + \frac{5}{3}M
u_{3} = 15 + 3M.

Note that all three of the above variables are, like x_{3}, strictly positive and hence we have constructed a feasible (but not necessarily basic) solution; and this reasoning is valid for any M > 0. Note also from the last row of the final tableau above that since one may write the objective as

z = \frac{34}{3}x_{3} + 6x_{4} − \frac{16}{3}u_{1} + \frac{5}{3}u_{2} + \frac{82}{3}

and hence including x_{3} in the solution increases the objective by 34/3 for every unit increase in x_{3}. The preceding analysis would seem to indicate that we can set x_{3} = M > 0 for any arbitrary large M and hence z can be made arbitrarily large with a feasible solution. In other words there is no bounded optimal solution.

5.4 The theory behind the simplex method

Given a linear programming problem in canonical form

maximise z := c · x subject to: Ax = b x ≥ 0_{n}.

where A ∈ ℝ^{m×n} and we are at some point in the simplex algorithm were we are in possession of a basic feasible solution of the form

x = \binom{0_{n−m}}{x_{B}}

where x_{B} > 0_{m} represents the basic variables corresponding to the decomposition of the matrix A = (A^{∘}|B). Arranging the constraints so that the basic variables appear in only one row and with unit coefficient we obtain the tableau

B^{−1}A^{∘}  ⋮ I_{m}  x_{B}
c_{B}^{T}B^{−1}A^{∘} − (c^{∘})^{T} ⋮ 0_{m}^{T}  z = c_{B}^{T}B^{−1}b

The last row gives the coefficients in the objective function when it is written in terms of the non-basic variables. To see where this comes from, suppose that we remove one basic variable from the current solution and introduce another variable giving a new solution

y = \binom{y^{∘}}{y_{B}}

where y^{∘} are the first n−m entries and y_{B} are the last m entries corresponding to the old basis variables. Since

Ay = A^{∘}y^{∘} + By_{B} = b

it follows that

y_{B} = B^{−1}b − B^{−1}A^{∘}y^{∘}.

If we also write the objective vector in obvious notation

c = \binom{c^{∘}}{c_{B}}

then since the objective function valued with the initial basic feasible solution z(x) = c_{B} · x_{B} and Bx_{B} = b so that z(x) = c_{B} · B^{−1}b it follows that the new value of the objective function

z(y) = c^{∘} · y^{∘} + c_{B} · y_{B}
= c^{∘} · y^{∘} + c_{B} · (B^{−1}b − B^{−1}A^{∘}y^{∘})
= z(x) + ((c^{∘})^{T} − c_{B}^{T}B^{−1}A^{∘})y^{∘}.

The last equality shows us that the coefficients in the objective function when written in terms of non-basic variables which fill y^{∘} are the entries of the vector (c^{∘})^{T} − c_{B}^{T}B^{−1}A^{∘}. This explains where the final row of the tableau given above comes from.^{15}

For i = 1, ..., n − m, write o_{i} for the i-th element of the vector (c^{∘})^{T} − c_{B}^{T}B^{−1}A^{∘} so that

z(y) = z(x) + ∑_{i=1}^{n−m} o_{i}y_{i}.

If we see that o_{j} > 0 for some j = 1, ..., n − m then it would be profitable to take our new basis y so that

y^{∘} = (0, ⋯, y_{j}, 0, ⋯, 0) where y_{j} > 0

and remove one of the elements of the basis variables in x_{B}. This would ensure that we have improved on the value of the objective function since then z(y) = z(x) + o_{j}y_{j} > z(x). It is of course necessary to do so in a way that y remains in the feasible region. (This corresponds to considering how one may introduce the variable y_{j} into each of the m equations given by Ay = b whilst respecting the requirement that y ≥ 0_{n}).

The essence of the algorithm is thus updating the basis of the feasible solution in such a way that one is always introducing a variable from the non-basic variables to create a new basic feasible solution such that the value of the objective function increases. One keeps a tab on which variables to introduce by always writing the objective function in terms of non-basic variables. Eventually, all the coefficients o_{i} will be negative in which case one has obtained the optimal basic feasible solution. The fundamental theorem of linear programming then tells us that we have found an optimal solution.

15Don’t forget that one negates the coefficients before filling them in the tableau.

5.5 Degeneracy and cycling

There is always the possibility that, when introducing a new variable into the basis, one of the basic variables coincidentally becomes equal to zero. We define such basic feasible solutions as degenerate. Recall that a basic feasible solution is an extreme point of the convex polyhedron describing the feasible region and the simplex method is an algorithm which systematically moves through extreme points of the feasible region.

It may so happen that the simplex method brings you to an extreme point which lies on the the hyperplane x_{i} = 0 for some i = 1, ..., n. (Recall that there are n constraining hyperplanes to the feasible region which are given by the condition x ≥ 0_{n}). In that case, the basic feasible solution will be such that x_{i} belongs to the basis and yet x_{i} = 0.

As an exercise, one may draw a picture and consider the extreme points of a linear programming problem with the constraints

x_{1} − x_{2} ≤ 2
2x_{1} + x_{2} ≤ 4
−3x_{1} + 2x_{2} ≤ 6
x_{1}, x_{2} ≥ 0.

For this reason, one can say that although the dimension of the basis is always equal to m =rank(A) it may happen that the number of zeros in a basic feasible solution may be greater than n − m.

Note that the analysis in the previous section concerning the introduction of a non-basic variable y_{j} to improve the objective value assumed that x_{B} > 0_{m} and y_{j} > 0. The analysis is not necessarily the same if the solution x_{B} is degenerate.

Because of the phenomenon of degeneracy, it is possible to slip into a cycle with the simplex method as we now explain. Suppose one has a basic feasible solution which is degenerate. In that case, recalling that the simplex tableau may always be rearranged in the form

B^{−1}A^{∘}  ⋮ I_{m}  x_{B}
c_{B}^{T}B^{−1}A^{∘} − (c^{∘})^{T} ⋮ 0_{m}^{T}  z = c_{B}^{T}B^{−1}b

then a degenerate solution implies that one of the entries of x_{B} is equal to zero. Then removing the corresponding element of the basis and introducing a new basic variable will have the effect that the new basic variable will also be zero. In other words, after pivoting, the next basic feasible solution will be degenerate too. Further, the objective will remain at the same value after pivoting since the new basic variable, being equal to zero, contributes no value of the objective function. Consequently, the simplex algorithm will move a degenerate solution to a degenerate solution without increasing the objective value and ultimately one ends up cycling between degenerate solutions. Most linear programming problems do not exhibit this phenomena however.

5.6 The two-phase method

There is a particular point that we need to be aware of in what we have discussed thus far. In the three examples of Sections 5.1, 5.2 and 5.3 we were able to specify an initial basic feasible solution by simply setting the slack variables equal to the entries of the vector b. The above theoretical considerations also show how one moves from one basic feasible solution to another under the assumption that one has a basic feasible solution in the first place.

What would happen if it is not immediately obvious how to form an initial basic feasible solution?

Here is a very likely scenario in which one might find this to be an issue. Suppose that at least one of the entries of the vector b is strictly negative. Naively following what we have done in the two examples would tell us to set the slack variable as equal to that value. However one should of course not forget about the constraints of positivity which would then be violated! So how does one set up the initial basic feasible solution from which to proceed with the simplex method? This is where the two-phase method can prove to be useful.

Phase I. The first step of the two-phase method is to introduce even more^{16} variables which are known as artificial variables. The way to do this is via the following series of actions.

(i) Arrange the linear programming problem in standard form so that the constraints are given as Ax ≤ b (even if the entries in the vector b are negative) and x ≥ 0_{n}. We assume as usual that A ∈ ℝ^{m×n}.

16In this context ‘even more’ means over and above the slack variables

(ii) Introduce slack variables in the usual way.

(iii) For each row that has a negative entry in b, multiply by −1 producing constraints of a new linear programming problem, say A′x′ = b′ and x′ ≥ 0_{m+n}, where x′ includes the slack variables and now b′ ≥ 0_{m}. Note at this point, if i is a row which has been multiplied through by −1, the slack variables u_{i} has coefficient −1 and hence cannot be set as equal to b′_{i} = −b_{i} as part of a basis solution.

(iv) Introduce to each row a new artificial variable, say w_{i} for i = 1, ..., m.

The principle action of phase I is to solve the following auxiliary linear programming problem.

Let A′′ = (A′|I_{m}),

x′′ = \binom{x′}{w}

and write 1_{m} for the vector in ℝ^{m} which has all of its entries equal to one.

maximise z′′ = −1_{m} · w subject to: A′′x′′ = b′ x′′ ≥ 0_{n+2m}.

The importance of this linear programming problem is that since it is clear that the last condition also implies that w ≥ 0_{m}, then the maximal value of z′′ is in fact zero. Further, if we solve it using the simplex method, then the solution which will achieve this maximal value of zero will necessarily be a basic solution and it will necessarily have w = 0_{m}. In other words, solving the above auxiliary linear programming problem will provide a basic feasible solution to the system

\begin{matrix}A′x′ = b′ \\ x′ ≥ 0_{n+m}, \end{matrix}

which is exactly what we are after to proceed to phase II.

Phase II. Taking the solution of the auxiliary linear programming problem we may now return to the original linear programming problem in canonical form using it as an initial basic solution from which to proceed in the usual way with the simplex method. One should take care however, to immediately write the original objective function in terms of the the non-basic variables at the end of phase I in order to proceed to the first tableau.

Roughly speaking, the two-phase method introduces an additional set of variables for the purpose of establishing an initial basic feasible solution in the easiest way of ‘setting x_{B} = b’ (once things have been set up in the right way of course) and then ‘pushing’ this basic feasible solution (with the help of the simplex algorithm) onto another basic feasible solution in terms of the original set of variables.

5.7 An example of the two-phase method

These theoretical considerations of the the two-phase method are perhaps best illustrated with an example. Consider the problem:

maximise z = 3x_{1} + x_{2} subject to: x_{1} + x_{2} ≤ 6 4x_{1} − x_{2} ≥ 8 2x_{1} + x_{2} = 8 x_{1}, x_{2} ≥ 0

This problem has a mixture of constraint types^{17} and after a slack variable x_{3} is added to the first constraint, a surplus variable x_{4} subtracted from the second, we have the following problem

maximise z = 3x_{1} + x_{2} subject to: x_{1} + x_{2} + x_{3} = 6 4x_{1} − x_{2} − x_{4} = 8 2x_{1} + x_{2} = 8 x_{1}, x_{2}, x_{3}, x_{4} ≥ 0

and there is no obvious basic feasible solution.

In order to overcome this, non negative artificial variables, x_{5} and x_{6} are added to the left hand side of the second and third constraints. Note that this is in slight contradiction with the method described in the previous section which dictates that we should add an artificial variable to the first constraint too. However a little thought regarding how we proceed in this example reveals that in the end this will make no difference. This creates the auxiliary linear programming problem,

17Note that this problem is slightly different to the ones discussed up until now in the sense that some of the constraints are slack and some are already tight. The reader will learn with experience that one only need introduce as many slack variables as there are slack constraints and proceed in the usual way with the simplex method.

maximise z′′ = −x_{5} − x_{6}

subject to:

x_{1} + x_{2} + x_{3} = 6 4x_{1} − x_{2} − x_{4} + x_{5} = 8 2x_{1} + x_{2} + x_{6} = 8 x_{1}, ..., x_{6} ≥ 0

In order for the augmented system to correspond to the original constraint system, the artificial variables must be zero. This is achieved by using the Simplex method to minimise their sum (or maximise the negative sum) and this is what we called Phase I. If this terminates with all artificial variables at zero, a basic feasible solution to the original constraint system will have been obtained. The Simplex method is then used to maximise the original objective in Phase II of the procedure.

Of course, in this example it would be a simple matter to use one of the constraints to eliminate x_{2}, or x_{1} from the other constraints and objective, leaving a trivial problem. In general however, this would not be at all convenient with many variables and many constraints. This simple example shows how the general procedure would be used.

Phase I

x_{1} x_{2} x_{3} x_{4} x_{5} x_{6}
x_{3} 1 1 1 0 0 0 6
x_{5} 4 -1 0 -1 1 0 8
x_{6} 2 1 0 0 0 1 8
-6 0 0 1 0 0 -16

Now x_{1} enters and x_{5} leaves the set of basic variables. We shall not bother to transform the x_{5} column for the next tableau as the artificial variable will not be considered to enter the basis again.

x_{1} x_{2} x_{3} x_{4} x_{5} x_{6}
x_{3} 0 5/4 1 1/4 - 0 4
x_{1} 1 -1/4 0 -1/4 - 0 2
x_{6} 0 3/2 0 1/2 - 1 4
0 -3/2 0 -1/2 - 0 -4

Now x_{2} enters and x_{6} leaves the set of basic variables and we reach the final tableau for phase I. In this tableau we have added an extra row for the coefficients of the objective function z in terms of the basis variables which conclude phase I and which initiate phase II.

x_{1} x_{2} x_{3} x_{4} x_{5} x_{6}
x_{3} 0 0 1 -1/6 - - 2/3
x_{1} 1 0 0 -1/6 - - 8/3
x_{2} 0 1 0 1/3 - - 8/3
z′′ 0 0 0 0 - -  0
z 0 0 0 -1/6 - - 32/3

The last row indicates that we can improve the value of z by introducing x_{4} into the basis which leads to the following tableau (now for the original problem in canonical form - ie we have now entered phase II).

x_{1} x_{2} x_{3} x_{4}
x_{3} 0 1/2 1 0 2
x_{1} 1 1/2 0 0 4
x_{4} 0 3 0 1 8
0 1/2 0 0 12

From this final table we see that there are no further variables which can be introduced into the basis in order to improve the value function and we have reached an optimal solution with the basis x_{1} = 4, x_{2} = 0, x_{3} = 2, x_{4} = 8 and z = 12.

6 Duality

6.1 Introduction

Recall the Betta Machine Product problem.

Maximise z = 15x_{1} + 20x_{2} subject to: 5x_{1} + 8x_{2} ≤ 16, 000 5x_{1} + 4x_{2} ≤ 14, 000 x_{1} + 3x_{2} ≤ 5, 000 x_{1} ≥ 0, x_{2} ≥ 0.

The variables x_{1} and x_{2} represent quantities of products 1 and 2 whilst the constraints are concerned with the available numbers of minutes of time for casting, machining and assembly respectively.

Multiplying the first constraint by 3 gives:

15x_{1} + 24x_{2} ≤ 48, 000.

Comparing this with the equation for z shows that we have found an upper bound for z as we must have z ≤ 48, 000. Could we use a similar trick to find a tighter bound? By adding twice the first constraint to the second we obtain

15x_{1} + 20x_{2} ≤ 46, 000

from which we can deduce that z ≤ 46, 000.

How can we find the best possible bound on z ? Suppose we follow the above technique and look for linear combinations of the three constraints that would give us the best upper bound. Let us add multiples of the linear constraints together using the multipliers w_{1}, w_{2} and w_{3}. We would need

5w_{1} + 5w_{2} + w_{3} ≥ 15
8w_{1} + 4w_{2} + 3w_{3} ≥ 20
w_{1}, w_{2}, w_{3} ≥ 0

giving an upper bound

16, 000w_{1} + 14, 000w_{2} + 5, 000w_{3}.

Note that the positivity of the variables w_{1}, w_{2} and w_{3} is necessary to preserve the inequality in each of the original constraints.

Hence, if we wish to find the smallest bound possible, we have to solve another linear programming problem; namely

minimise 16, 000w_{1} + 14, 000w_{2} + 5, 000w_{3} subject to: 5w_{1} + 5w_{2} + w_{3} ≥ 15 8w_{1} + 4w_{2} + 3w_{3} ≥ 20 w_{1}, w_{2}, w_{3} ≥ 0

This is called the dual problem. To distinguish the two problems, the original is called the primal problem.

6.2 Symmetric dual

Consider the standard maximisation problem (P): given, A ∈ ℝ^{m×n}, b ∈ ℝ^{m} and c ∈ ℝ^{n},

maximise c · x subject to: Ax ≤ b x ≥ 0_{n}.

Its symmetric dual is the standard minimisation problem (D): given A^{T} ∈ ℝ^{n×m}, c ∈ ℝ^{n} and b ∈ ℝ^{m},

minimise b · w subject to: A^{T}w ≥ c w ≥ 0_{m}.

Remarks

1. The primal has n variables and m constraints whilst the dual has m variables and n constraints. In the last chapter we had n ≥ m and m =rank(A) because the number of variables, n, included slack and surplus variables. In the present situation slack and surplus variables have not been counted so either n or m could be the larger.

2. Any problem may be put into form (P) or (D). To do this you may need to multiply by −1 or to write an equation as two inequalities. Note it is not required here that the elements of b should be non negative.

3. It is a very important point that the dual of the dual is the primal.

4. To appreciate why the dual problem we have identified is symmetric it is necessary to consider the asymmetric dual problem given below.

6.3 The asymmetric dual

Suppose that instead of starting with the linear programming problem in standard form, we started with the linear programming problem in canonical form. That is to say: given, A ∈ ℝ^{m×n}, b ∈ ℝ^{m} and c ∈ ℝ^{n},

maximise c · x subject to: Ax = b x ≥ 0_{n}.

Following previous advice, we can write this linear programming problem in the form

maximise c · x subject to: \binom{A}{−A}x ≤ \binom{b}{−b} x ≥ 0_{n}.

Now following the formulation for the dual but partitioning the dual vector

ℝ^{2m} ∋ y = \binom{u}{v},

we have dual problem minimise \binom{b}{−b} · \binom{u}{v} subject to: (A^{T}| − A^{T})\binom{u}{v} ≥ c u, v ≥ 0.

For convenience we can write z = u − v and now the dual can be more conveniently written as

minimise b · z ; subject to: ; A^{T}z ≥ c

noting in particular that z ∈ ℝ^{m}. Hence the dual of the canonical linear programming problem has no positivity constraint and hence is asymmetric to its primal.

6.4 Primal to the Dual via the diet problem

Here is an example of how we may interpret duality by considering a classic diet problem.

Suppose that there are nutrients N_{i}, i = 1, ..., m and foods F_{j}, j = 1, ..., n with food costs c_{j} per unit of food F_{j}. Suppose further that each food F_{j} contains A_{ij} units of of nutrient N_{i}. The diet problem is to minimise the cost of the diet subject to the constraints that the diet provides the minimal nutritional requirements which are specified by the quantities b_{i}; i = 1, ..., m for each nutrient N_{i}. If (x_{1}, ..., x_{n}) represents the quantity of the respective foods F_{1}, ..., F_{n} in the diet (and clearly we must have x_{j} ≥ 0 for all j = 1, ..., n), then the cost of the diet becomes ∑_{j} c_{j}x_{j}. On the other hand, to respect the minimum requirement of nutrient N_{i} we must also have that ∑_{j} A_{ij}x_{j} ≥ b_{i}. In other words the problem becomes:

minimise c · x ; subject to ; Ax ≥ b ; x ≥ 0_{n}.

The above we may consider as our primal problem.

Now suppose that there is a company which, instead of selling the foods F_{1}, ..., F_{n}, sells the raw nutrients (for example in the form of pills) at respective costs y_{1}, ..., y_{m}. Since food F_{j} contains A_{ij} units of nutrient of N_{i}, then ∑_{i} y_{i}A_{ij} is the equivalent cost of food F_{j} for j = 1, ..., n. The manufacturer of the raw nutrients wishes to maximise the profits of selling the recommended dietary intake, ∑_{i} b_{i}y_{i}, whilst competing with the primal diet. In other words, the equivalent cost of food F_{j} should not be more expensive than the primal cost of food F_{j} so that ∑_{i} y_{i}A_{ij} ≤ c_{j}.

From the perspective of the manufacturer of the raw nutrients the problem is to:

maximise b · y subject to y^{T}A ≤ c^{T} y ≥ 0_{m}.

We see that the manufacturer of the raw nutrients is thus concerned with the dual problem.

Intuitively we would expect that any feasible solution to the above problem, which would specify a profit for a particular diet of the raw nutrients N_{1}, ..., N_{m}, should not be greater than any feasible solution to the primal problem which is the cost of a diet using the foods F_{1}, ..., F_{n}. That is to say, we would expect b · y ≤ c · x where x and y are feasible solutions to the primal and dual respectively. Moreover in that case, it would seem sensible that both have optimal solutions if there exist feasible solutions, x and y, such that b·y = c·x. This is the beginning of duality theory which consumes the remainder of this chapter.

6.5 The Duality Theorem

We shall henceforth continue our analysis for the case of the symmetric duality. That is to say we shall analyse the relationship between the primal (P) and dual (D) given in Section 6.2.

As we have seen in the preliminary example, the dual objective gives a bound for the primal objective (and vice versa). This is exactly what the Weak Duality Theorem says.

Theorem 6.1 (Weak Duality Theorem) Suppose that x_{0} is a feasible solution to the primal problem (P) and y_{0} is a feasible solution to its dual (D). Then

c · x_{0} ≤ b · y_{0}.

Proof. Since Ax_{0} ≤ b and y_{0} ≥ 0_{m} it follows that

y_{0} · Ax_{0} ≤ y_{0} · b

Similarly, since A^{T}y_{0} ≥ c and x_{0} ≥ 0_{n} it follows that x_{0} · A^{T}y_{0} ≥ x_{0} · c; or transposing both sides, that

y_{0} · Ax_{0} ≥ c · x_{0}.

Comparing the established inequalities we see immediately that

c · x_{0} ≤ b · y_{0}

as required.

Since in the proof above, x_{0} and y_{0} are abitrary, it follows that

sup_{x∈F (P)} c · x ≤ inf_{y∈F (D)} b · y  (6)

where F (P) and F (D) are the feasible regions of (P) and (D) respectively. Two consequences follow immediately from this observation which we give below.

Corollary 6.2 If (P) has a feasible solution but no bounded optimal solution then (D) has no feasible solution. Since (P) is the dual of (D) it also follows that if (D) has a feasible solution but no bounded optimal solution then (P) has no feasible solution.

Proof. When (P) has an unbounded optimal solution then (6) shows that, should it be the case that F (D) ≠ ∅, inf_{y∈F (D)} b · y = ∞ which is a contradiction. It must therefore follow that F (D) = ∅.

Corollary 6.3 Suppose there exist feasible solutions for (P) and (D), say x_{0} and y_{0} respectively, such that c · x_{0} = y_{0} · b then necessarily x_{0} and y_{0} are optimal solutions for (P) and (D) respectively.

Proof. Note that for any pair of solutions x_{0} and y_{0} we have

c · x_{0} ≤ sup_{x∈F (P)} c · x ≤ inf_{y∈F (D)} b · y ≤ b · y_{0}

and hence if c · x_{0} = y_{0} · b then all of the above inequalities are necessarily equalities an in particular we see immediately that x_{0} and y_{0} are optimal.

The main result of this section is the following.

Theorem 6.4 (Duality Theorem) The primal problem (P) has a finite optimal solution if and only if the dual problem (D) does too in which case the optimal value of their objective functions are the same.

The way we will prove the Duality Theorem is to consider first four auxiliary results. The first two of these four are classic results coming from convex analysis.

The first result makes the somewhat intuitively obvious statement that for any point which lies outside of a closed convex set, it is possible to pass a hyperplane between that point and the convex set.

Lemma 6.5 (Separating Hyperplane Lemma) Suppose that C is a closed convex set in ℝ^{n} and that b ∉ C. Then there exists y ∈ ℝ^{n} such that y · b < y · z for all z ∈ C.

C  b z_{0}  y 0

Figure 6: A diagrammatic representation of elements of the Separating Hyperplane Lemma.

Proof. Let \overline{B}(b, r) be a closed hyper-sphere whose centre is at b and radius r is sufficiently large that C ∩ \overline{B}(b, r) ≠ ∅. Now consider the function f : ℝ^{n} → ℝ defined by

f(z) = {(b − z) · (b − z)}^{1/2} = ||b − z||

for z ∈ ℝ^{n}. As f is a continuous function and C ∩\overline{B}(b, r) is a closed bounded convex domain^{18} then f is bounded and there exists a z_{0} ∈ C ∩ \overline{B}(b, r) such that ^{19}

r ≥ f(z_{0}) = inf\{f(z) : z ∈ C ∩ \overline{B}(b, r)\} =: δ > 0.

Note that the first inequality is a consequence of the fact that z ∈ C ∩ \overline{B}(b, r) ⊆ \overline{B}(b, r) can never be further than distance r from b and the last inequality is due to the fact that b ∉ C. Hence for all z ∈ C ∩ \overline{B}(b, r)

δ = ||b − z_{0}|| ≤ ||b − z||  (7)

18Recall that the intersection of two closed domains is closed and that the non-empty intersection of two convex domains is convex

19Google Weierstrauss’ Theorem!

and otherwise for z ∈ C\\overline{B}(b, r)

δ = ||b − z_{0}|| ≤ r ≤ ||b − z||.

In other words (7) holds for all z ∈ C (even when C is an unbounded domain). As C is a convex domain, it follows in particular that for each λ ∈ (0, 1) and z ∈ C

(b − z_{0}) · (b − z_{0}) ≤ (b − λz − (1 − λ)z_{0}) · (b − λz − (1 − λ)z_{0}).

This may otherwise be written as a quadratic inequality in λ,

0 ≤ λ^{2}||z_{0} − z||^{2} + 2λ(b − z_{0}) · (z_{0} − z)

and since λ ∈ (0, 1) we may divide through by λ and then let λ ↓ 0 to discover that necessarily (b − z_{0}) · (z_{0} − z) ≥ 0. Note we have used the fact that ||z_{0} − z)||^{2} ≥ 0. We may now write

(b − z_{0}) · b = (b − z_{0}) · (b − z_{0}) + (b − z_{0}) · z_{0}
= ||b − z_{0}||^{2} + (b − z_{0}) · z_{0}
> (b − z_{0}) · z_{0}  (8)
≥ (b − z_{0}) · z

for all z ∈ C. Now take y = z_{0} − b and we see that

y · b < y · z

for all z ∈ C as required. Note that the strict inequality follows from (8).

The second result is slighly less intuitively obvious and says that the range of a linear transformation, when restricted in its domain to the positive orthant, is a closed convex domain. Convexity is the easy part, being closed is difficut!

Lemma 6.6 Suppose as usual that A ∈ ℝ^{m×n}. Define

C = {z = Ax|x ≥ 0_{n}}

Then C is closed and convex.

Proof. To check that C is convex, suppose that z_{1}, z_{2} ∈ C and λ ∈ (0, 1). Then λz_{1} + (1 − λ)z_{2} ∈ C because

λz_{1} + (1 − λ)z_{2} = A(λx_{1} + (1 − λ)x_{2})

and λx_{1} + (1 − λ)x_{2} ≥ 0_{n} (the latter is a consequence of the fact that the positive orthant is a convex space).

To show that C is closed, suppose that {z_{n} : n ≥ 1} is a sequence in C which has some limit z. We need to show that z ∈ C. To do this, note that by definition there exists a sequence of vectors {x_{n} : n ≥ 1} such that z_{n} = Ax_{n}. If it so happens that x_{n} has a limit, say x, then since the positive orthant is a closed space, it must follow that x ≥ 0_{n}. Then by continuity of the mapping A : ℝ^{n} ↦ ℝ^{m} it follows that z = Ax.

If the sequence of vectors x_{n} does not converge then we can look for a convergent subsequence instead. We do this by assuming without loss of generality that x_{n} is a basic feasible solution and hence can be written in the form

x_{n} = \binom{0}{x_{B_{n}}}

where x_{B_{n}} is the basic part which induces a partition A = (A_{n}^{∘}|B_{n}) which depends on n. As there can only be a finite number of variables we can include in the basis, along the sequence x_{n} we will see the same basis being used infinitely often. For any particular basis which partitions A = (A^{∘}|B) we can therefore identify a subsequence {x_{B_{n_{k}}} : k ≥ 1} of the basic feasible solutions x_{n} which may be written in terms of that basis. Since B is inveritble we thus have x_{B_{n_{k}}} = B^{−1}z_{n_{k}} and hence the existence of the limit z and continuity of the inverse mapping implies the existence of a limit of the x_{n_{k}} to some x satisfying

x = \binom{0}{B^{−1}z}.

As before, we necessarily have x in the positive orthant as it is a closed set and clearly z = Ax. Hence C contains all its limit points and thus, by definition, is closed.

The remaining two lemmas give complementary existence results using the Separating Hyperplane Lemma. The first is used to prove the second and the second is the principle result which is instrumental in proving the Duality Theorem. We assume that the quantities A and b are given as in the primal linear programming problem.

Lemma 6.7 (Farkas’ Lemma) Either

(i) Ax = b and x ≥ 0_{n} has a solution or

(ii) y^{T}A ≥ 0_{n}^{T} and y · b < 0 has a solution

but not both.

Proof. Suppose both (i) and (ii) hold. Then, on the one hand from (i), since x ≥ 0_{n}, we have y^{T}Ax ≥ 0. On the other hand from (ii) we have y^{T}Ax = y · b < 0. Hence (i) and (ii) cannot be true at the same time. In particular when (i) is true then (ii) is false. To complete the proof it suffices to show that when (i) is false then (ii) is true.

Suppose now that (i) is false and recall the definition

C = {z = Ax|x ≥ 0_{n}}

which was shown to be closed and convex in the previous Lemma. Note that by assumption b ∉ C. Then by The Separating Hyperplane Lemma we have the existence of y such that y · b < y · z for all z ∈ C. Taking x = 0_{n} we see then that y · b < y · A0_{n} = 0. The proof is thus concluded by showing that y^{T}A ≥ 0_{n}^{T}.

To this end, suppose that (y^{T}A)_{i} =: λ_{i} < 0. Then let

x = \frac{1}{λ_{i}} (y · b)e_{i} ≥ 0_{n}

where e_{i} is the vector in R^{n} with zero entries except in the i-th row where it has a unit entry. Now note that

y^{T}Ax = y · b

but on account of the definition of y from the Separating Hyperplane Lemma and the positivity of this particular choice of x, we also have

y · b < y · Ax = y^{T}Ax.

This constitutes a contradiction and hence all elements of y^{T}A are positive as required.

Lemma 6.8 Either (i) Ax ≤ b and x ≥ 0_{n} has a solution, or

(ii) y^{T}A ≥ 0_{n}^{T}, y · b < 0 and y ≥ 0_{m} has a solution

but not both.

Proof. If both (i) and (ii) hold, then

0 ≤ y^{T}Ax = y · Ax ≤ y · b < 0

which is a contradiction. Hence (i) and (ii) cannot hold simultaneously and in particular if (i) is true then (ii) is false. To complete the proof, it again suffices to show that if (i) is false then (ii) is true to complete the proof.

Suppose then that (i) is false. This implies that Ax + z = b has no solution for x ≥ 0_{n}, z ≥ 0_{m} as otherwise we would have that Ax = b−z ≤ b. Otherwise said, we have that

(A|I_{m})\binom{x}{z} = b and \binom{x}{z} ≥ 0_{n+m}

has no solution. Hence by Lemma 6.7 it follows that there exists a y ∈ ℝ^{m} such that y^{T}(A|I_{m}) ≥ 0_{n+m}^{T} and y · b < 0 as required. Note that the last but one inequality implies that y^{T}A ≥ 0_{n}^{T} and y^{T} ≥ 0_{m}^{T}.

Now we are ready to prove the Duality Theorem.

Proof of Theorem 6.4. Our strategy for proving this theorem is in two steps. Step 1 is to show that if (P) has no finite optimal feasible solution, then (D) has no finite optimal feasible solution. Since (P) is the dual of (D) it then follows by symmetry that (P) has (no) finite optimal feasible solution if and only if (D) has (no) finite optimal feasible solution. After that, step 2 will show that when both (P) and (D) have finite optimal solutions then their objective values must be equal.

Step 1. There are two ways that (P) may have no finite optimal feasible solution. Either because there is an unbounded optimal solution or because there is no feasible solution. If there is an unbounded optimal solution then by Corollary 6.2 it follows that (D) cannot have a feasible solution. Now suppose that (P) has no feasible solution. Then by Lemma 6.8 there exists a y_{∗} ≥ 0_{m} such that y_{∗}^{T}A ≥ 0_{n}^{T} and y_{∗} · b < 0. Suppose however that there exists a feasible solution to (D), say y, then it must be the case that for any λ ≥ 0 that y + λy_{∗} is also a feasible solution. However, as λ is unconstrained, it follows that we can make (y + λy_{∗}) · b arbitrarily large in magnitude showing that (D) has an unbounded feasible solution.

Step 2. From Step 1 it is now apparent that (P) has a bounded optimal solution if and only (D) has a bounded optimal solution. Assume now that both (P) and (D) have bounded optimal solutions. We know already from the Weak Duality Theorem that

max c · x ≤ min y · b.

It would therefore suffice to show that there exist feasible solutions x and y to (P) and (D) respectively such that c · x ≥ y · b. To this end, let us consider the existence solutions x ≥ 0_{n}, y ≥ 0_{m} to the linear system

Ax ≤ b ; y^{T}A ≥ c^{T} ; c · x ≥ y · b

or, in a more suitable notation for what follows,

( \begin{matrix}A & 0_{m×m} \\ 0_{n×n} & −A^{T} \\ −c^{T} & b^{T} \end{matrix} ) \binom{x}{y} ≤ ( \begin{matrix}b \\ −c \\ 0 \end{matrix} )  (9)

Now consider the linear system which requires one to find z ≥ 0_{m}, w ≥ 0_{n} and λ ≥ 0 such that

(z^{T} w^{T} λ) ( \begin{matrix}A & 0_{m×m} \\ 0_{n×n} & −A^{T} \\ −c^{T} & b^{T} \end{matrix} ) ≥ 0_{n+m}^{T}  (10)

and at the same time

( \begin{matrix}z \\ w \\ λ \end{matrix} ) · ( \begin{matrix}b \\ −c \\ 0 \end{matrix} ) < 0.  (11)

If we can prove that this system has no solution, then by Lemma 6.8 it follows that (9) must hold with x ≥ 0_{n}, y ≥ 0_{m}. In other words we want to show that

z^{T}A ≥ λc^{T} ; λb^{T} ≥ w^{T}A^{T} ; z · b < w · c.  (12)

has no solution to complete the proof.

If a solution to (12) were to exist and λ = 0 then let x_{∗} and y_{∗} be the optimal solutions to (P) and (D) (which are assumed to exist^{20}) and we have

z^{T}Ax_{∗} ≥ 0 and y_{∗}^{T}Aw ≤ 0.

But on the other hand, since Ax_{∗} ≤ b and y_{∗}^{T}A ≥ c^{T} then

c · w ≤ y_{∗}^{T}Aw ≤ 0 ≤ z^{T}Ax_{∗} ≤ z · b.

However, this leads to a contradiction of the third statement of (12).

Suppose instead that a solution to (12) were to exist and λ > 0. In that case we would see that

λz · b ≥ z^{T}Aw ≥ λc · w

which contradicts again the third statement of (12).

So in short, there can be no solution to (12), in other words no solution to (10) and (11), and hence there must be a solution to (9) and the proof is complete.

There is a very nice observation we can also make from the Duality Theorem.

Corollary 6.9 Suppose that there exists at least one feasible solution to each of (P) and (D) respectively. Then there exist optimal bounded optimal solutions to both (P) and (D) with equal objective values.

Proof. The result follows as a result of Corollary 6.2 and the Duality Theorem.

6.6 Complementary slackness

There is a another very nice relation between the solution of the primal and the dual which can be useful for verifying that a solution to the primal (and hence dual) is optimal.

20Remember that we are proving that IF bounded optimal feasible solutions to (P) and (D) exist then their objective values are equal.

Theorem 6.10 (Symmetric complementary Slackness) Consider the standard linear programming problem:

maximise z = c · x subject to: Ax ≤ b x ≥ 0_{n}

where A ∈ ℝ^{m×n}, b ∈ ℝ^{m}, c ∈ ℝ^{n} are given and x ∈ ℝ^{n}.

Then x ∈ ℝ^{n} is an optimal feasible solution for the primal and y ∈ ℝ^{m} is an optimal feasible solution for the dual if and only if they are feasible solutions and

\begin{matrix}y_{i}[(Ax)_{i} − b_{i}] = 0 ∀i \\ x_{i}[(y^{T}A)_{i} − c_{i}] = 0 ∀i. \end{matrix}

Note that the last two conditions can also be written equivalently as

\begin{matrix}y_{i} > 0 ⇒ (Ax)_{i} = b_{i} \\ x_{i} > 0 ⇒ (y^{T}A)_{i} = c_{i}. \end{matrix}

or equivalently again as

\begin{matrix}(Ax)_{i} < b_{i} ⇒ y_{i} = 0 \\ (y^{T}A)_{i} > c_{i} ⇒ x_{i} = 0. \end{matrix}

Proof. First we prove the ⇐ direction. To this end suppose that y_{i} > 0 ⇒ (Ax)_{i} = b_{i} and that x_{i} > 0 ⇒ (y^{T}A)_{i} = c_{i}. In that case

y^{T}Ax = ∑_{i} y_{i}(Ax)_{i} = ∑_{i:y_{i}>0} y_{i}(Ax)_{i} = ∑_{i:y_{i}>0} y_{i}b_{i} = ∑_{i} y_{i}b_{i}.

In other words

y^{T}Ax = y^{T}b = y · b.

On the other hand we also have

y^{T}Ax = ∑_{i}(y^{T}A)_{i}x_{i} = ∑_{i:x_{i}>0}(y^{T}A)_{i}x_{i} = ∑_{i} c_{i}x_{i}

so that we also have

y^{T}Ax = c · x

Hence c · x = y · b and thanks to the Weak Duality Theorem this shows that x and y are optimal for the primal and the dual respectively.

Now we deal with the ⇒ direction. Suppose that x and y are optimal for the primal and the dual respectively. Then by the Duality Theorem we have that c · x = y · b. Taking the latter into account together with the fact that Ax ≤ b and y^{T}A ≥ c^{T} we have

y^{T}Ax ≤ y^{T}b = c^{T}x ≤ y^{T}Ax.  (13)

In other words the inequalities in (13) are actually equalities and hence

y^{T}Ax − y^{T}b = ∑_{i} y_{i}{(Ax)_{i} − b_{i}} = 0.  (14)

As x is feasible it follows that (Ax)_{i} − b_{i} ≤ 0 and as y is feasible, it also follows that y_{i} ≥ 0. In consequence, re-examining (14) we see that

y_{i}{(Ax)_{i} − b_{i}} = 0

for each i = 1, ..., m. In other words y_{i} > 0 ⇒ (Ax_{i}) = b_{i} or equivalently Ax_{i} < b_{i} ⇒ y_{i} = 0.

To complete the proof note that (13) also implies that

∑_{i}{(y^{T}A)_{i} − c_{i}}x_{i} = 0

and similar reasoning to above shows that x_{i} > 0 ⇒ (y^{T}A)_{i} = c_{i} or equivalently (y^{T}A)_{i} > c_{i} ⇒ x_{i} = 0.

A similar result can be established for the asymmetric duality.

Theorem 6.11 (Asymmetric complementary slackness) Consider the canonical linear programming problem:

maximise z = c · x subject to: Ax = b x ≥ 0_{n}

where A ∈ R^{m×n}, b ∈ R^{m}, c ∈ ℝ^{n} are given and x ∈ ℝ^{n}.

Then x ∈ ℝ^{n} is an optimal feasible solution for the primal and y ∈ ℝ^{m} is an optimal feasible solution for the dual if and only if they are feasible solutions and

x_{i}[(y^{T}A)_{i} − c_{i}] = 0 ∀i.

Note that the last condition can also be written equivalently as

x_{i} > 0 ⇒ (y^{T}A)_{i} = c_{i}.

or equivalently again as

(y^{T}A)_{i} > c_{i} ⇒ x_{i} = 0.

Proof. The proof is somewhat simpler than the symmetric case. Suppose that x and y are optimal feasible solutions for the primal and dual respectively. Then by the Duality Theorem we have c · x = b · y. It follows that

y^{T}Ax = y^{T}b = c^{T}x

and hence

(y^{T}A − c^{T})x = ∑_{i}{(y^{T}A)_{i} − c_{i}}x_{i} = 0.  (15)

As x ≥ 0_{n} and y^{T}A ≥ c^{T} it follows that x_{i} > 0 ⇒ (y^{T}A)_{i} = c_{i}.

Now suppose that x_{i} > 0 ⇒ (y^{T}A)_{i} = c_{i} (or equivalently that (y^{T}A)_{i} > c_{i} ⇒ x_{i} = 0). In that case it follows that (15) holds. That is to say y^{T}Ax = c^{T}x. On the other hand, as Ax = b it follows that y^{T}Ax = y^{T}b. In conclusion c · x = b · y so that by the Weak Duality Theorem x and y are optimal.

Note that at the end of the above proof we have used the Weak Duality Theorem for the asymmetric dual. Strictly speaking the Weak Duality Theorem, as presented in this text, is a statement about the symmetric case. None the less it is still true for the asymmetric case. Can you prove it? (See Exercise sheets).

Here is an example of how to use complementary slackness. Consider the linear programming problem

minimise z = x_{1} + 6x_{2} + 2x_{3} − x_{4} + x_{5} − 3x_{6} subject to: x_{1} + 2x_{2} + x_{3} + 5x_{6} = 3 −3x_{2} + 2x_{3} + x_{4} + x_{6} = 1 5x_{2} + 3x_{3} + x_{5} − 2x_{6} = 2 x_{1}, ..., x_{6} ≥ 0.

Writing c^{T} = (1, 6, 2, −1, 1, −3) we can multiply the objective by −1 to convert the primal to a maximisation problem and then write down its (asymmetric) dual,

minimise z = 3z_{1} + z_{2} + 2z_{3} subject to: z^{T}A ≥ −c^{T}

Note in particular we have z ∈ ℝ. It is more convenient to multiply the constraints by −1 and since there is no restriction on the sign of the entries in z we can set −y = z. This leads us to the representation of the dual problem,

maximise z = 3y_{1} + y_{2} + 2y_{3} subject to: y^{T}A ≤ c^{T}

Suppose we are asked to verify that

x^{T} = (0, \frac{16}{29}, 0, \frac{66}{29}, 0, \frac{11}{29})

is optimal for the primal. If this were the case then the Asymmetric Complementary Slackness Theorem would hold. In particular we would have that as x_{2}, x_{4}, x_{6} are all strictly positive, then (y^{T}A)_{2} = c_{2}, (y^{T}A)_{4} = c_{4} and (y^{T}A)_{6} = c_{6}. Having written out carefully the matrix A we see that this is tantamount to

2y_{1} − 3y_{2} + 5y_{3} = 6
y_{2} = −1
5y_{1} + y_{2} − 2y_{3} = −3.

Solving these equations we find that

(y_{1}, y_{2}, y_{3}) = (−\frac{4}{29}, −1, \frac{19}{29})

and it is a straightforward exercise to check that y is feasible for the dual and that c · x = b · y = −3/29.

In conclusion, by assuming the given x is optimal, we have been able to use complementary slackness to deduce what the optimal y should be. However, having obtained this hypothetical value of y we discover that in fact x and y are necessarily optimal solutions to primal and dual respectively by by virtue of the fact that they have the same objective values and hence the Weak Duality Theorem comes into effect.

7 The transportation problem

7.1 Introduction

We are concerned with moving a commodity from various sources (factories) to various destinations (markets) at minimum cost. The standard problem has m sources of the commodity with respective (strictly positive) supplies s_{i}, i = 1, ..., m and n destinations for the commodity having respective (strictly positive) demands d_{j}, j = 1, ..., n. The transport cost per unit between supply i and demand j is c_{ij}. We wish to plan the transport from the supply to the demand at minimal total cost. The problem may be formulated as a linear programming problem with special structure as follows.

Let x_{ij} be the amount transported from source i to destination j. The problem is then:

minimise z = ∑_{i=1}^{m} ∑_{j=1}^{n} c_{ij}x_{ij}

subject to:

∑_{j=1}^{n} x_{ij} ≤ s_{i} i = 1, ..., m ∑_{i=1}^{m} x_{ij} ≥ d_{j} j = 1, ..., n x_{ij} ≥ 0 ∀i = 1, ..., m and j = 1, ..., n.

Remarks

1. Note that the general transportation problem is a linear programming problem with mn variables and m + n constraints.

2. In order for a feasible solution to exist we must necessarily have that total supply is greater or equal to demand. This is indeed the case since

∑_{i=1}^{m} s_{i} ≥ ∑_{i=1}^{m} ∑_{j=1}^{n} x_{ij} ≥ ∑_{j=1}^{n} d_{j}.  (16)

7.2 The canonical form

If it is the case that supply strictly outstrips demand, then we can introduce an additional ‘dump’ to the problem where additional supply, ∑_{i}^{m} s_{i} −∑_{j}^{n} d_{j} is consumed and cost of shipping to this ‘dump’ state as equal to zero. Note that the zero shipping cost means that objective function for the enlarged problem remains the same as the original one. In that case we may assume without loss of generality that supply is equal to demand. In other words the transportation problem is balanced and

∑_{i=1}^{m} s_{i} = ∑_{j=1}^{n} d_{j}.

In that case, the inequalities in (16) are in fact equalities and the transportation problem takes the canonical or balanced form as follows:

minimise z = ∑_{i=1}^{m} ∑_{j=1}^{n} c_{ij}x_{ij}

subject to:

∑_{j=1}^{n} x_{ij} = s_{i} i = 1, ..., m ∑_{i=1}^{m} x_{ij} = d_{j} j = 1, ..., n x_{ij} ≥ 0 ∀i = 1, ..., m and j = 1, ..., n

where c_{ij} ≥ 0, s_{i} > 0, d_{j} > 0 are given with ∑_{i} s_{i} = ∑_{j} d_{j}.

7.3 Dual of the transportation problem

It turns out that providing an effective algorithmic solution to the transportation problem involves making calculations which interplay between the primal and dual. We therefore devote a little time to understanding the relationship of the transportation problem to its dual.

Suppose we introduce the matrix

X = (x_{ij}) = ( x_{1}^{T} ; ⋮ ; x_{m}^{T} )

where for i = 1, ..., m

x_{i}^{T} = (x_{i1}, ..., x_{in}).

Then let x^{T} = (x_{1}^{T}, ..., x_{m}^{T}), in other words

x^{T} = (x_{11}, x_{12}, ..., x_{1n}, x_{21}, ..., x_{2n}, ..., x_{mn}).

Similarly write

C = (c_{ij}) = ( c_{1}^{T} ; ⋮ ; c_{m}^{T} )

where for i = 1, ..., m

c_{i}^{T} = (c_{i1}, ..., c_{in}).

Then let c^{T} = (c_{1}^{T}, ..., c_{m}^{T}), in other words

c^{T} = (c_{11}, c_{12}, ..., c_{1n}, c_{21}, ..., c_{2n}, ..., c_{mn}).

Now let

b^{T} = (s_{1}, ..., s_{m}, d_{1}, ...., d_{n})^{T}.

The primal problem can be written in the form

minimise z = c · x subject to: Ax = b x ≥ 0_{mn}

where A ∈ ℝ^{(m+n)×mn} and for 1 ≤ i ≤ m, 1 ≤ j ≤ mn we have^{21}

A_{ij} = δ_{j,{(i−1)n+1,....,in}}

and for m + 1 ≤ i ≤ m + n, 1 ≤ j ≤ mn

A_{ij} = δ_{j,{(i−m),(i−m)+n,...,(i−m)+(m−1)n}}.

Otherwise said, the main constraints can be seen in the following form

( \begin{matrix}1 & 1 ⋯ 1 & 1 \\ & & & 1 & 1 ⋯ 1 & 1 \\ & & & & & & . & . & . \\ & & & & & & & & & 1 & 1 ⋯ 1 & 1 \\ 1 & & & 1 & & & & & & 1 \\ & 1 & & & 1 & & & & & & 1 \\ & ... & & & ... & & & & & & ... \\ & & 1 & & & 1 & & & & & & 1 \end{matrix} ) ( x_{11} ; ⋮ ; x_{1n} ; x_{21} ; ⋮ ; x_{2n} ; ⋮ ; x_{m1} ; ⋮ ; x_{mn} ) = ( s_{1} ; ⋮ ; s_{m} ; d_{1} ; ⋮ ; d_{n} ) . (17)

It is important to note that the matrix A has only entries which are either zero or unity. Further, there are only two non-zero entries in each column.

21For scalar j and set S, we use the delta function δ_{j,S} which is equal to 1 if j ∈ S and zero otherwise.

This means that the matrix A^{T} has only two non-zero entries in each row. Another way of representing A is to characterise its columns. We know that A ∈ ℝ^{(m+n)×mn} and hence write

A = (c_{11}|c_{12}| ⋯ |c_{1n}|c_{21}| ⋯ |c_{mn})

where

c_{ij} = \binom{e_{i}^{(m)}}{e_{j}^{(n)}}

where e_{l}^{(k)} is the vector in ℝ^{k} whose entries all zero with the exception of the l-th entry which is unity.

Consider the asymmetric dual^{22} which requires one to

maximise z = b · y subject to: A^{T}y ≤ c y ∈ ℝ^{m+n}.

For convenience let us write

y = (u_{1}, ..., u_{m}, v_{1}, ..., v_{n})^{T}.

Then note that the condition A^{T}y ≤ c can be written as the system of inequalities

u_{i} + v_{j} ≤ c_{ij} for all i = 1, ..., m and j = 1, ..., n.

Further, the objective function b · y can be written as

∑_{j=1}^{n} v_{j}d_{j} + ∑_{i=1}^{m} u_{i}s_{i}.

To conclude, the dual problem to the balanced transportation problem can be summarised as the following linear programming problem

maximise z = ∑_{j=1}^{n} v_{j}d_{j} + ∑_{i=1}^{m} u_{i}s_{i}.

subject to:

u_{i} + v_{j} ≤ c_{ij} for all i = 1, ..., m and j = 1, ..., n. with no restriction on the sign of the u_{i}s and v_{j}s.

22In order to write down the asymmetric dual, it is useful to write the primal in the form: maximise z = (−c) · x subject to −Ax = −b and x ≥ 0_{mn}.

There is an intuitive way of looking at the formulation of the dual problem. The primal concerns a way of shipping goods from supply to demand at minimal cost. Suppose now there is a competing agent who wishes to offer the same service. This agent charges u_{i} for picking up goods at the i-th supply and charges v_{j} for delivering to the j-th demand^{23}. The charge for shipping from i → j is thus u_{i} + v_{j}. In order to remain competitive, this agent wishes to maximise his income, i.e. ∑_{j=1}^{n} v_{j}d_{j} + ∑_{i=1}^{m} u_{i}s_{i}, but at the same time offering a cheaper alternative to the operator in the primal problem, i.e. u_{i} + v_{i} ≤ c_{ij} on all routes i → j.

7.4 Properties of the solution to the transportation problem

Let us first examine existence of the solution to the transportation problem. Here The Duality Theorem will play an essential role.

Theorem 7.1 The transportation problem when in balance has a finite optimal solution.

Proof. Consider the solution to the primal

x_{ij} = \frac{s_{i}d_{j}}{∑_{i=1}^{m} s_{i}} .

It is easy to check that this solution satisfies ∑_{i} x_{ij} = d_{j}, ∑_{j} x_{ij} = s_{i}, x_{ij} ≥ 0. In other words it is a feasible solution.

On the other hand, if we look at the dual to the transportation problem, then it is also clear that v_{j} = u_{i} = 0 for all i = 1, ..., m and j = 1, ..., n is a feasible solution.

Hence by Corollary 6.9 a bounded optimal solution to both primal and dual must exist.

The following result is requires a proof which goes beyond the scope of this course and hence the proof is omitted.

23There is a slight problem with this interpretation in that the charges u_{i} and v_{j} may be negative valued! Can you think why this is not necessarily a problem?

Theorem 7.2 If the supplies s_{i} : i = 1, ..., m and the demands d_{j} : j = 1, ..., n are all integers, then any feasible solution to the transportation problem X = (x_{ij}) in balance is an integer valued matrix. That is to say X ∈ ℤ_{+}^{(m+n)×mn}.

Intuitively speaking, the reason why this result is true is due to the fact that the entries of the matrix A are either zero or unity.

Finally we conclude this section with a result concerning the rank of the matrix A which will be of use in the next section when we look at an algorithm for solving the transportation problem.

Theorem 7.3 Rank(A) = m + n − 1.

Proof. Recall that A ∈ ℝ^{(n+m)×mn}. Suppose that r_{1}^{T}, ..., r_{m+n}^{T} are the rows of A. Partitioning A between the m-th and (m + 1)-th row, in any one column we have just one unit entry on either side and the remaining terms zero. Hence

∑_{i=1}^{m} r_{i} = ∑_{i=m+1}^{n} r_{i} = ( 1 ⋮ 1 ) ∈ ℝ^{mn}.

It follows that the rows of A are linearly dependent. The main aim of this proof is to deduce from this fact that this means that rank of the matrix A is the number of rows minus one.

To this end, consider solutions (λ_{1}, ⋯, λ_{m+n−1}) to the equation

u(λ_{1}, ..., λ_{m+n−1}) = 0

where

u(λ_{1}, ..., λ_{m+n−1}) := ∑_{i=1}^{m+n−1} λ_{i}r_{i}.

As we do not include the bottom row of the matrix A it must be the case that

(∑_{i=1}^{m+n−1} λ_{i}r_{i})_{n} = λ_{1} = 0
(∑_{i=1}^{m+n−1} λ_{i}r_{i})_{2n} = λ_{2} = 0

⋮

(∑_{i=1}^{m+n−1} λ_{i}r_{i})_{mn} = λ_{m} = 0.

Note, it is worth looking at (17) to see why the above holds.

As a consequence of the above we also have

(∑_{i=1}^{m+n−1} λ_{i}r_{i})_{1} = λ_{1} + λ_{m+1} = λ_{m+1} = 0
(∑_{i=1}^{m+n−1} λ_{i}r_{i})_{2} = λ_{1} + λ_{m+2} = λ_{m+2} = 0
⋮
(∑_{i=1}^{m+n−1} λ_{i}r_{i})_{n−1} = λ_{1} + λ_{m+n−1} = λ_{m+n−1} = 0.

In conclusion, the unique solution to the equation u(λ_{1}, ..., λ_{m+n−1}) = 0 is λ_{1} = ⋯ λ_{m+n−1} = 0 showing that all but the last row of A are linearly independent and hence rank(A) = m + n − 1.

To some extent the conclusion of the last theorem is intuitively obvious since we have m equations of the form ∑_{j} x_{ij} = s_{i}, n equations of the form ∑_{i} x_{ij} = d_{j} with one constraint ∑_{i} s_{i} = ∑_{j} d_{j}.

7.5 Solving the transportation problem

The transportation problem is nothing more than a linear programming problem and hence there is no reason why not to solve it using the simplex algorithm. However, by virtue of the large number of variables involved this not the most practical option. Below we describe an algorithm for solving the transportation problem which is essentially based on the Asymmetric Complementary Slackness Theorem. Looking back to this theorem and putting it into the current context we have that, if x_{ij} and u_{i}, v_{j} are feasible, then they satisfy

x_{ij} > 0 ⇒ u_{i} + v_{j} = c_{ij}

and

u_{i} + v_{j} < c_{ij} ⇒ x_{ij} = 0

where i = 1, ..., m and j = 1, ..., n.

Here-now is the promised algorithm.

(a) Find a basic feasible solution to the primal. Since rank(A) = m + n − 1 it follows that that there should be at most m + n − 1 of the entries x_{ij} are non-zero.^{24} Assume that this solution is optimal.

24When there are less than m + n − 1 non-zero entries, we have a degenerate basic feasible solution.

(b) Write down the m + n − 1 equations which are of the form

u_{i} + v_{j} = c_{ij} for each i, j that x_{ij} > 0.

This will give the solution to the m + n variables u_{1}, ...., u_{m}, v_{1}, ..., v_{n} in terms of an unknown parameter. We can fix this parameter by choosing (for example) u_{1} = 0.

(c) If with this set values for u_{i}, v_{j}, i = 1, ..., m, j = 1, ..., n we also have

that

u_{i} + v_{j} ≤ c_{ij} for each i, j that x_{ij} = 0

then by the the Asymmetric Complementary Slackness Theorem it follows that the feasible solution (x_{ij}) must be optimal as assumed.

If however for some pair i, j we have that

u_{i} + v_{j} > c_{ij} and x_{ij} = 0

then we have violated the conditions of the Asymmetric Complementary Slackness Theorem (specifically the dual solution is not feasible) and the assumption that our feasible solution is optimal must be false. We can improve things by increasing the value of this particular x_{ij} (where i, j are the offending pair of indicies) to a non-zero value, adjusting the values of the other x_{ij}-s so that there are still m + n − 1 non-zero entries^{25} and returning to step (a).

This algorithm is not complete as at least the following two questions remain to be answered when utilising the suggested routine.

How does one choose an initial basic feasible solution?

How does one re-adjust the feasible solution in step (c)?

We address these two points in the next two sections.

25This is akin to moving to another basic feasible solution in the Simplex Method.

7.6 Initial feasible solutions

There are several methods, but here are two of the more popular.

The ‘North-West corner’ method: Start by putting x_{11} = min{s_{1}, d_{1}} and then move along the rows of x_{ij} using up either the supplies or filling the demands until all supplied are distributed and all demands are fulfilled. Below is an example.

7 11 15 ← d_{j}
14 7 7 0
13 0 4 9
6 0 0 6
s_{i} ↑

The ‘matrix’ (or ‘cheapest route’) method: We find min_{i,j} c_{ij} and then allocate as large an x_{ij} value as possible to route i → j (in other words satisfying either supply or demand, whichever is the smaller). We repeat with the next smallest^{26} c_{ij} and on until all factories are empty and all markets are full. This method is lengthier but gives a better starting point. Consider the previous example but with the transport costs given by the matrix

(c_{ij}) = ( \begin{matrix}14 & 13 & 6 \\ 15 & 14 & 8 \\ 9 & 11 & 2 \end{matrix} ) .

In this case we establish an initial feasible solution as

7 11 15 14 0 5 9 13 7 6 0 6 0 0 6

Note if the problem has been set up with an additional ‘dump’ then obviously it would not be advisable to allocate to the dump.

26If there is a tie between two values of c_{ij} then one picks one of them arbitrarily.

7.7 A worked example

We are given supply/demand/cost data in the following format

55 70 35 40 ← d_{j}
80 13 11 18 17
100 2 14 10 1
20 5 8 18 11
s_{i} ↑

Step 1: We start by constructing a feasible solution using the North-West corner method and obtain the following matrix for (x_{ij}) with the cost matrix (c_{ij}) embedded in.

55  70  35  40  u_{i} ↓
80 55 _{[13]} 25 _{[11]} 0 _{[18]} 0 _{[17]}  −3
100 0 _{[2]} 45 _{[14]} 35 _{[10]} 20 _{[1]}  0
20 0 _{[5]} 0 _{[8]} 0 _{[18]} 20 _{[11]}  10
v_{j} → 16  14  10  1 z = 2210

On the right hand column and bottom row, we have also included the values of (u_{1}, u_{2}, u_{3}, v_{1}, v_{2}.v_{3}) by solving the complementary slackness equations. Specifically

x_{11} > 0 ⇒ u_{1} + v_{1} = 13
x_{12} > 0 ⇒ u_{1} + v_{2} = 11
x_{22} > 0 ⇒ u_{2} + v_{2} = 14
x_{23} > 0 ⇒ u_{2} + v_{3} = 10
x_{24} > 0 ⇒ u_{2} + v_{4} = 1
x_{34} > 0 ⇒ u_{3} + v_{4} = 11.

We have solved this system by choosing arbitrarily u_{2} = 0.^{27} By inspection we can see that this solution is not optimal because

x_{31} = 0 and yet u_{3} + v_{1} = 26 > 5 = c_{31}.

in fact the entry in x_{31} is the ‘worst offender’ in this class so to speak in the sense that

min_{ij:x_{ij}=0}{c_{ij} − u_{i} − v_{j}} = c_{31} − u_{3} − v_{1} = −21

27Note the algorithm we have given suggests to take u_{1} = 0, however we have chosen u_{2} = 0 to show that it does not matter.

(which is necessarily a negative number when the solution is not optimal).

Step 2: We can improve our solution by introducing a non-zero entry to x_{31} by so that there are the same number of non-zero entries. Intuitively speaking, our initial solution was a basic solution and we now wish to change the basis to introduce x_{31}.

We introduce x_{31} with the following perturbation loop.

·  ·  ·  ·
· 55 − η 25 + η ·  ·  ·
·  ·  45 − η · 20 + η ·
· 0 + η  ·  · 20 − η ·
·  ·  ·  ·  ·

In this perturbation loop we add and subtract η from a loop of nodes in such a way that in any row or any column, for every time η has been added on to an entry, it has also be subtracted off from another entry. This has the effect that the perturbed solution is still feasible. Note that the maximum value of η we can work with is 20 in order to keep the maximum number of non-zero entries equal to m + n − 1 = 6. Reallocating resources with η = 20 and computing the u_{i}-s and v_{j}-s again using complementary slackness we have

55  70  35  40
80 35 _{[13]} 45 _{[11]} 0 _{[18]} 0 _{[17]}  −3
100 0 _{[2]} 25 _{[14]} 35 _{[10]} 40 _{[1]}  0
20 20 _{[5]} 0 _{[8]} 0 _{[18]} 0 _{[11]} −11
16  14  10  1 z = 1790

where we can see that again we have arbitrarily chosen u_{2} = 0.

Step 3: Note that

min_{ij:x_{ij}=0}{c_{ij} − u_{i} − v_{j}} = c_{21} − u_{2} − v_{1} = −14

and so we propose to introduce x_{21} with the following perturbation loop

·  ·  · ·
· 35 − η 45 + η ⋯
· 0 + η 25 − η ⋯
·  ·  ·  ⋯
·  ·  ⋯

The maximum value of η permitted is 25 and hence we obtain the new table

55  70  35  40
80 10 _{[13]} 70 _{[11]} 0 _{[18]} 0 _{[17]}  11
100 25 _{[2]} 0 _{[14]} 35 _{[10]} 40 _{[1]}  0
20 20 _{[5]} 0 _{[8]} 0 _{[18]} 0 _{[11]}  3
2  0  10  1 z = 1440

where the new values of u_{i} and v_{j} have been filled in with the arbitrary choice u_{2} = 0.

Step 4: We compute as usual

min_{ij:x_{ij}=0}{c_{ij} − u_{i} − v_{j}} = c_{13} − u_{1} − v_{3} = −3.

and work with the following perturbation loop

·  ·  ·  ·
· 10 − η · 0 + η · ·
· 25 + η · 35 + η · ·
·  ·  ·  ·  · ·
·  ·  ·  · ·

which permits η = 10 as its largest value. This results in the following table

55  70  35  40
80 0 _{[13]} 70 _{[11]} 10 _{[18]} 0 _{[17]}  8
100 35 _{[2]} 0 _{[14]} 25 _{[10]} 40 _{[1]}  0
20 20 _{[5]} 0 _{[8]} 0 _{[18]} 0 _{[11]}  3
2  3  10  1 z = 1410

By construction we have that x_{ij} > 0 ⇒ u_{i} + v_{j} = c_{ij} and a quick inspection shows that for the remaining entries satisfy u_{i} + v_{j} ≤ c_{ij} and x_{ij} = 0. Hence the conditions of the Asymmetric Complementary Slackness Theorem are fulfilled and this solution must be optimal.

7.8 Improvement at each iteration

Why did we work with the ‘worst offender’ min_{ij:x_{ij}=0}{c_{ij} − u_{i} − v_{j}} at each stage of the iteration in the example above?

If x_{kl} is the variable to be made basic, then we find a ‘pertubation loop’ of basic variables of the form

x_{kl} ↔ x_{kp} ↕ x_{rp} ↔ x_{r.}, ⋯, x_{s.} ↔ x_{sl} ↕ x_{kl}

These are alternately decreased and increased by the value η. The difference between the new value of the objective minus old value of the objective is therefore

η(c_{kl} − c_{kp} + c_{rp} − c_{r.} + ⋯ + c_{s.} − c_{sl})

= η[c_{kl} − (u_{k} + v_{p}) + (u_{r} + v_{p}) − (u_{r} + v_{.}) + ⋯ + (u_{s} + v_{.}) − (u_{s} + v_{l})]

= η(c_{kl} − u_{k} − v_{l})
< 0

whenever c_{kl} −u_{k} −v_{l} < 0. Note that the first equality follows becasue, apart from the pair k, l all other entries x_{ij} in the pertubation loop are positive and hence for such pairs i, j we have c_{i,j} = u_{i} + v_{j}. The second equality is the consequence of a large cancellation of terms. For a non-degenerate solution we will always have η > 0 and so there will be a definite reduction in the total transport cost. Choosing the ‘worst offender’, ie chooising x_{kl} such that

c_{kl} − u_{k} − v_{l} = min_{ij:x_{ij}=0}{c_{ij} − u_{i} − v_{j}} < 0,

is a good way of trying to maximise the reduction at each step^{28}. Of course, if one does not choose the ‘worst offender’ to include in the basis, it is not detrimental to the ultimate aim of minimising the objective, so long as one chooses an x_{kl} for which c_{kl} − u_{k} − v_{l} < 0.

7.9 Degeneracy

Recall that a basic feasible solution to a linear programming problem is degenerate if one of the variables in the basis is zero-valued. In the context of a transportation problem this would correspond to finding a solution x_{ij} for which there are strictly less than m + n − 1 non-zero entries. (Recall that rank(A) = m + n − 1).

28If one considers two different perturbation loops, then there is not guarantee that the quantity η is the same for each one and hence strictly speaking we cannot compare the relative reduction between loops at this level of abstraction even if the quantity c_{kl}−u_{k}−v_{l} is the smallest negative value possible.

It is not uncommon to obtain degenerate solutions to transportation problems (this can even happen by using, for example, the matrix method to construct an initial solution). One needs to be a little careful in the case that one has a degenerate solution. Consider for example the following example (taken from the 2004 exam).

In a land levelling project, soil must be taken from four zones and taken to five other zones. The amounts (in thousands of tons) to be removed from zones A,B,C and D are 25, 50, 30 and 40 respectively whilst the amounts needed at zones K, L, M, N and P are 30, 60, 20, 15 and 20. The average distances travelled (in meters) in moving soil between the various zones are shown below.

K L M N P
A 120 70 200 60 100
B 80 130 100 130 140
C 190 130 70 110 160
D 130 90 100 180 150

Solve the problem of moving the soil as required whilst minimising the effort required in terms of total tons times meters travelled.

This problem is a standard balanced transportation problem where the table above can be taken directly as the cost matrix. Applying the matrix method we can immediately produce the following initial solution.

30  60  20  15  20
25 0 _{[120]} 10 _{[70]} 0 _{[200]} 15 _{[60]} 0 _{[100]} ·
50 30 _{[80]} 0 _{[130]} 0 _{[100]} 0 _{[130]} 20 _{[140]} ·
30 0 _{[190]} 10 _{[130]} 20 _{[70]} 0 _{[110]} 0 _{[160]} ·
40 0 _{[130]} 40 _{[90]} 0 _{[100]} 0 _{[180]} 0 _{[150]} ·
·  ·  ·  ·  · ·

Note that m = 4, n = 5 and hence m + n − 1 = 8 and yet the solution above only has 7 non-zero entries. It is degenerate. This poses a problem if we try to compute the dual values u_{i}, v_{j} as, even choosing u_{1} = 0, there are not enough equations to obtain exact dual values. One way around this is to flag one of the zeros, say 0 = x_{22} as the degenerate part of the basic solution (we have done this by putting it in bold) and then including the equality u_{2} + v_{2} = c_{22} in the set of linear equations used to solve for the dual variables. We can now proceed to solve for the dual values (choosing u_{1} = 0) and obtain.

30  60  20  15  20
25 0 _{[120]} 10 + η _{[70]} 0 _{[200]} 15 − η _{[60]} 0 _{[100]} 0
50 30 _{[80]} 0 _{[130]} 0 _{[100]} 0 _{[130]} 20 _{[140]} 60
30 0 _{[190]} 10 − η _{[130]} 20 _{[70]} 0 + η _{[110]} 0 _{[160]} 60
40 0 _{[130]}  40 _{[90]} 0 _{[100]}  0 _{[180]} 0 _{[150]} 20
20  70  10  60  80 ·

On a practical note, there is a quick way of filling out the final column and row. Start by setting u_{1} = 0, move across the first row to non-zero entries, do a quick mental calculation and fill in the entries for those columns in the final row. Next by looking at the second row use these entries in the final row to obtain as many entries as possible in the final column and so on until all values of u_{i}, v_{j} have been filled out. Having filled out the final column and final row one may now quickly check for optimality by looking at all the zero entries and seeing if the corresponding values in the terminal row and column add up to something less than or equal to the corresponding cost in square brackets.

We have also included in the table above a perturbation loop around the ‘worst offender’ x_{34}.

Following through with the calculations (bearing in mind that we are still thinking of 0 = x_{22} as part of the basic solution), the next table takes the form

30 60 20 15 20
25 0 _{[120]} 20 _{[70]} 0 _{[200]} 5 _{[60]} 0 _{[100]} 0
50 30 _{[80]} 0 _{[130]} 0 _{[100]} 0 _{[130]} 20 _{[140]} 60
30 0 _{[190]} 0 _{[130]} 20 _{[70]} 10 _{[110]} 0 _{[160]} 50
40 0 _{[130]} 40 _{[90]} 0 _{[100]} 0 _{[180]} 0 _{[150]} 20
20  70  20  60  80 z = 13, 000

and the reader can easily check that we now have that u_{i} + v_{j} ≤ c_{ij} for all i, j such that x_{ij} = 0. Hence we have found the optimal solution.

7.10 Pricing out

Suppose that we wish to prevent that the solution to the transportation problem makes use of the route i → j. We can do this by setting c_{ij} = M where M is some arbitrarily large value (or indeed M = ∞). The transportation algorithm will naturally move away from solutions which use this route as it is a minimisation algorithm. We have thus ‘priced-out’ the route i → j.

Consider the example in the previous section. Suppose in addition, we want to solve this transportation problem with the additional constraint that at most 20,000 tons are taken from D to L. The way to handle this additional constraint is to spilt the zone D into two parts D_{1}, which corresponds to the first 20,000 tons of soil at D, and D_{2}, which corresponds to the remaining 20,000 tons of soil at D. We thus formulate the problem with the data in the table below.

30 60 20 15 20
25 120 70 200 60 100
50 80 130 100 130 140
30 190 130 70 110 160
D_{1} : 20 130 90 100 180 150
D_{2} : 20 130 M 100 180 150

This table shows in the last two rows that shipping any soil over and above 20,000 tons from D has identical costs to shipping any of the first 20,000 tons except to L which incurs an arbitrary large cost (and hence implicitly will be avoided by the algorithm). One may now proceed to solve this transportation problem in the usual way (note that the matrix method will produce a degenerate initial solution again).^{29} At any stage of the algorithm one should avoid introducing the variable x_{52} and when computing the values of the dual variables u_{i}, v_{j} note that it will always be the case that (when checking the zero 0 = x_{52}) u_{5} + v_{2} ≤ M since M is assumed to be arbitrarily large.

Below is another example of how to use pricing out (taken from the 2005 exam). Unlike the above example we use pricing out to build in a temporal feature of the given problem.

A factory makes a product for which there is a fluctuation but predictable demand. The following table shows the predicted demand for each of the next five months together with the production capacity and unit production costs for each month.

29Please refer to the model solutions provided on my webpage.

Month  1 2 3 4 5
Demand (‘000)  60 80 85 100 70
Production Capacity (‘000) 90 95 100 100 90
Unit production cost (pence) 50 55 60 60 65

Items produced in a month are available to meet demand in the same month but they may also be stored to meet demand 1,2,3 or 4 months later at a cost of 4,7,9 and 10 pence per item respectively. It is required to schedule production to meet the anticipated demand at minimum total cost.

Formulating this as a transportation problem is rather tricky as there is a temporal aspect that we need to build into the solution. The idea is to note that the supply at month n can in principle feed the demand at month n as well as all subsequent months. However the supply at month n cannot feed the demands at previous months.^{30} We may thus think of x_{ij} for i = 1, ..., 5 and j = i, .., 5 as the amount produced at time i which feeds the demand at time j. Further when computing c_{ij} we need to take account of the production cost at time i plus the storage costs for the additional months after the production month before the item is sold on. We think of the supply nodes as the production in the individual months (the second row of the above table) and the markets as the demand in each month (the first row of the above table). In that case we see that our transportation problem is not balanced and we need to introduce a ‘dump’. Hence we present the given data for use in a transportation problem as follows.

60 80 85 100 70 80
90 50 54 57 59 60 0
95 M 55 59 62 64 0
100 M M 60 64 67 0
100 M M M 60 64 0
90 M M M M 65 0

Using the matrix method we can introduce an initial solution given below. Note that it is degenerate and as before we have indicated in bold a zero which we wish to think of as an element of the basis and for which we shall force tightness in the corresponding dual inequality. We have also set as usual u_{1} = 0.

30We are making the assumption that time travel is not possible.

60  80  85  100  70  80
90 60 _{[50]} 30 − η _{[54]}  0 _{[57]} 0 _{[59]} 0 + η _{[60]}  0 _{[0]} 0
95 0 _{[M]} 50 + η _{[55]} 45 − η _{[59]} 0 _{[62]}  0 _{[64]}  0 _{[0]} 1
100 0 _{[M]}  0 _{[M]} 40 + η _{[60]} 0 _{[64]}  0 _{[67]} 60 − η _{[0]} 2
100 0 _{[M]}  0 _{[M]}  0 _{[M]} 100 _{[60]}  0 _{[64]}  0 _{[0]} 1
90 0 _{[M]}  0 _{[M]}  0 _{[M]} 0 _{[M]} 70 − η _{[65]} 20 + η _{[0]} 2
50  54  58  59  63  −2 ·

We have also indicated a perturbation loop and this yields the solution (setting again u_{1} = 0).

60  80  85  100  70 80
90 60 _{[50]} 0 _{[54]} 0 _{[57]} 0 _{[59]} 30 _{[60]} 0 _{[0]} 0
95 0 _{[M]} 80 _{[55]} 15 _{[59]} 0 _{[62]} 0 _{[64]} 0 _{[0]} 4
100 0 _{[M]} 0 _{[M]} 70 _{[60]} 0 _{[64]} 0 _{[67]} 30 _{[0]} 5
100 0 _{[M]} 0 _{[M]} 0 _{[M]} 100 _{[60]} 0 _{[64]} 0 _{[0]} 4
90 0 _{[M]} 0 _{[M]} 0 _{[M]} 0 _{[M]} 40 _{[65]} 50 _{[0]} 5
50  51  55  56  60 −5 z = 22, 885

A quick scan reveals that u_{i} + v_{j} ≤ c_{ij} for all i, j and hence we have reached an optimal solution.

8 Optimisation over Networks

8.1 Capacitated Networks

A capacitated network (N, γ) consists of a set of nodes N and a set of directed arcs^{31}, which can be regarded as a subset of of N × N. Associated with each arc is a non-negative real number γ which can be interpreted as the capacity of the arc. Formally γ is a a mapping γ : N ×N → [0, ∞) which satisfies

γ(i, j) ≥ 0 for all i, j ∈ N
γ(i, i) = 0 for all i ∈ N.

A network may be represented either

diagramatically, for example,

1  2 3  4 5  3 6 2  6

or equivalently by a matrix of values of γ

( \begin{matrix}0 & 5 & 0 & 0 \\ 0 & 0 & 6 & 3 \\ 0 & 2 & 0 & 6 \\ 0 & 0 & 0 & 0 \end{matrix} )

31To be clear, a directed arc simply means the direct path i → j where i, j ∈ {1, ..., N}.

8.2 Flows in Networks

Many problems are concerned with flows in networks where γ is interpreted as a capacity limiting the possible flow along an arc. One class of problems we shall be concerned with is finding the maximal flow through a network. In this we shall take some nodes, called sources, to be the origin of the flows and others, called sinks, to be where the flow is taken away. In general there can be many sources and sinks.

We need to introduce some new notation at this point.

1. Node s ∈ N is a source iff

γ(i, s) = 0 ∀i ∈ N and ∃i ∈ N s.t. γ(s, i) > 0

2. Node d is a sink iff

γ(d, i) = 0 ∀i ∈ N and ∃i ∈ N s.t. γ(i, d) > 0

3. A path is a sequence of arcs of the form:

(i_{1}, i_{2}), (i_{2}, i_{3}), ...., (i_{k−1}, i_{k})

where i_{1}, i_{2}, ..., i_{k} are distinct nodes.

4. A flow in a network is a function f defined on N × N satisfying:

(a) f(i, j) ≤ γ(i, j) ∀i, j ∈ N

(b) f(i, j) = −f(j, i)

(c) ∑_{j∈N} f(i, j) = ∑_{j∈N} f(j, i) ∀i ∈ N\(S ∪ D)

where S and D are the sets of sources and sinks respectively. Condition (c) simply requires conservation of flow at all nodes other than sources or sinks. Actually, both sides of the equation in (c) are both zero as a result of (b).^{32} The latter says that a flow in one direction is treated as equivalent to a negative flow in the opposite direction. This convention makes later computations somewhat more more convenient.

5. A path is saturated if it contains an arc (i, j) with f(i, j) = γ(i, j).

32The point being that from (b) we have that ∑_{j∈N} f(i, j) = − ∑_{j∈N} f(j, i) = ∑_{j∈N} f(j, i). Since the only number which is the negative of itself is zero, the claim follows.

8.3 Maximal flow problem and labelling algorithm

Henceforth we shall consider only networks with just one source and one sink.

For a network with single source and sink, the maximal flow problem is to determine arc flows f(i, j) in order to maximise the total flow from the source, A(f), which is defined by

A(f) = ∑_{i∈N} f(s, i).

In this section we shall describe the labelling algorithm.^{33} In outline, the algorithm is as follows. Suppose that a flow through the network has been established.

Step 1: Find an unsaturated path from source s to sink d. If none exists, the flow is maximal.

Step 2: Augment the existing flows to saturate the path. Return to Step1.

This looks easy, but the complicating feature is that the path may involve traversing arcs in the opposite direction to the present flow i.e. augmenting the flow may involve reducing the existing flow along one or more arcs. The most reliable way to find an unsaturated path is to adopt the labelling procedure now described.

Step 1: Suppose we have been able to assign a flow^{34} through the network, say g(i, j), then calculate the residual capacities, γ′(i, j), by the following

γ′(i, j) = γ(i, j) − g(i, j).

(a) Find all nodes i with the property

γ′(s, i) > 0

i.e all those nodes with an unsaturated arc from the source s. Label node i with (d_{i}, s) where d_{i} = γ′(s, i).

33Also known as the Ford-Fulkerson labelling algorithm.

34We mean specifically which respects the definition of a flow.

(b) Find nodes j with the property

γ′(i, j) > 0 where node i is labelled.

Label node j with (d_{j}, i) where

d_{j} = min(γ′(i, j), d_{i})

(c) Keep returning to Step 1 (b) until all nodes that can be labelled have been labelled.

Step 2: If the sink d is not labelled then all paths to the sink are saturated and the flow must be optimal. The optimal flow by the simple relation

g(i, j) = γ(i, j) − γ′(i, j)

where γ′ are the residual capacities in the current round of the algorithm.

If on the other hand the sink d is labelled then there is an unsaturated path from source to sink and the labels will show the nodes which make up the path as well as the minimal residual capacity along the path.

We now create a new set of residual capacities for the network, γ′′(ij) as follows

γ′′(i, j) = γ′(i, j) if i, j do not belong to the unsaturated path

and

γ′′(i, j) = γ′(i, j) − c if i, j belong to the unsaturated path

where c is the minimal residual flow along the unsaturated path.

Now return to Step 1 (a) replacing γ by γ′′.

This algorithm is to some extent incomplete as it does not tell us how to find an initial solution. Here is a rough description of how one may do this. Start at the source and find a path of arcs which connects the source to the sink. The maximum flow that can be allowed along this path is the minimum capacity of all the arcs along the path. Suppose this minimal capacity is c. Subtract this capacity from all the capacities along the the given path so that at least one arc has zero capacity. To get the best possible initial flow one can return to the source and proceed along another path with positive capacity at each arc and repeat the procedure and then continue in this way until there are no more paths possible from source to sink. Essentially this is going through the whole labelling procedure in one go and so it does not matter if not all paths are saturated on the initial flow as subsequent steps of the labelling algorithm will address these deficiencies.

Things are best visualised in an example!

8.4 A worked example

Consider the network whose structure and capacities are described by the matrix

( \begin{matrix}0 & 8 & 8 & 0 & 0 & 0 \\ 0 & 0 & 0 & 4 & 4 & 0 \\ 0 & 2 & 0 & 0 & 5 & 2 \\ 0 & 4 & 0 & 0 & 0 & 5 \\ 0 & 4 & 0 & 0 & 0 & 6 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{matrix} )

Which can otherwise be represented in diagrammatic form as follows.

1  2 3  45  6 0  8 2 0 8 0  00  5 4 4 4  4  0  5 2 0  6

We now construct an initial flow by looking at paths from the source to the sink and adding the minimal flow along these paths. In the diagram below we have inserted the residual capacities which follow from each included flow from source to sink. Note that final figures give us the initial values of γ′(i, j) to work with in step 1 (a) of the algorithm.

The paths we work with and their respective flows are

1 → 2 → 4 → 6 flow 4
1 → 3 → 5 → 6 flow 5
1 → 3 → 6 flow 2

1  2 3  45  6 0  8 2 0 8 0  00  54 4 4  4  0  5 2 0  6 3 4  4  8  0  4  1 1 5 0 5 2  0 5 7 1

The reader will note that strictly speaking we could in principle add more flow by adding a unit flow to the path 1 → 2 → 5 → 6. However, as we shall see below this is essentially what the first step of the labelling algorithm will do anyway.

Now following Step 1, (b) and (c) we insert the labels on the nodes where possible and observe that the sink has been labelled which means that saturation has not been reached. On the same diagram we have computed the new set of residual capacities.

1  2 3  4 5  6 0  4 4 2 7 4  4  8  0  4  1 0 1 5 0 5 2 1 (4,1) (1,1)  (4,2)  (1,5) 6  0 3 5 5  3

From the labels we see that a flow 1 → 2 → 5 → 6 of unit strength is still possible. We include this flow by augmenting the capacities again. With the new residual capacities we label the network as in the following diagram. Note that the sink is no longer labelled which means that the network is saturated. Using the relation that g(i, j) = γ(i, j) − γ′(i, j) where γ′ is the current residual capacities we find that

g(s, 1) = 8 − 3 = 5 and g(s, 2) = 8 − 1 = 7

making the maximal flow ∑_{j} g(s, j) = 12.

1  2 3  45  6 02 7  8  0  4  1 0 0 5 2 1 5  3  6  0 3 5 (3,1) (1,1)  (3,2)

Finally note that had we included the path 1 → 2 → 5 → 6 in the initial solution, the labelling algorithm would have concluded immediately that we have found the optimal solution.

8.5 Maximal flow and minimal cuts

The procedure described above is based upon ideas developed by Ford and Fulkerson in 1957. Central to the underlying theory is the idea of a cut, which is now described. This description again is limited to the case of netwoks having a single source and sink.

We need to introduce yet more notation.

1. A cut in a network is an ordered pair (A, B) of subsets of N such that

(a) s ∈ A

(b) d ∈ B (c) A ∪ B = N

(d) A ∩ B = ∅

2. The capacity of a cut γ(A, B) is defined as

γ(A, B) = ∑_{i∈A} ∑_{j∈B} γ(i, j)

This same notation will be used whenever A and B are subsets of N without necessarily being a cut.

3. If A and B are subsets of N, then the flow from A to B is defined as

f(A, B) = ∑_{i∈A} ∑_{j∈B} f(i, j)

Since for any i, j ∈ N, f(i, j) ≤ γ(i, j), double summation gives

f(A, B) ≤ γ(A, B)  (18)

for any cut (A, B). Lemma 8.1 For any cut (A, B) and any flow f, the amount of flow in the network equals the flow across (A, B). That is to say,

A(f) = f(A, B).

Proof. Note that by definition of a flow we ∑_{j∈N} f(i, j) = 0 for any i ∈ N\{s, d}. Now note that

A(f) = ∑_{j∈N} f(s, j) + 0
= ∑_{j∈N} f(s, j) + ∑_{i∈A,i≠s} ∑_{j∈N} f(i, j)
= ∑_{i∈A} ∑_{j∈N} f(i, j)
= ∑_{i∈A} ∑_{j∈A} f(i, j) + ∑_{i∈A} ∑_{j∈B} f(i, j)
= 0 + ∑_{i∈A} ∑_{j∈B} f(i, j)
= f(A, B)

where in the penultimate equality, the terms f(i, j) and f(j, i) = −f(j, i) appear in the double sum for i, j ∈ A which cancel one another out.

There is a nice intuition behind the above lemma. Think of arcs as pipes. Imagine that we are pumping fluid into the network at the source and the capacities represent the maximal flow of liquid that can pass through each pipe. A cut is a way of severing the network so that all liquid pumped in at the source leaks out of the system. Naturally, if this is the case then the total flow of the ‘leakage’ must match the flow pumped in at the source. This is what the lemma says.

Combining the conclusion of this lemma with this with the inequality (18), we see that for any cut (A, B),

A(f) ≤ γ(A, B).  (19)

In other words the amount of flow cannot exceed the capacity of any cut. This inequality allows us to deduce that if for some flow f we can find a cut (A′, B′) such that

A(f) = γ(A′, B′)

then it must be the case that f is an optimal flow. Indeed in that case then by (19) we necessarily have that

γ(A′, B′) = A(f) ≤ max_{g} A(g) ≤ min_{(A,B)} γ(A, B) ≤ γ(A′, B′)

where the maximum is taken over all flows g and the minimum is taken over call cuts (A, B). Hence the inequalities above must be equalities and we see that A(f) is optimal, (A′, B′) is a minimal cut in the sense that

γ(A′, B′) = min_{(A,B)} γ(A, B)

and A(f) = min_{(A,B)} γ(A, B).

Below, the celebrated maximal flow/minimal cut theorem, provides the converse to this conclusion.

Theorem 8.2 (Maximal flow/minimal cut theorem) Let f be a maximal flow, then there is a cut (A′, B′) such that

A(f) = γ(A′, B′) = min_{(A,B)} γ(A, B).

Thus the maximal flow at the source equals the minimal capacity across all cuts.

Proof. Let A′ be the set of nodes which have been labelled in the final step of the algorithm together with the source and let B′ be the set of unlabelled notes. Clearly A′ ∪ B′ = N and A′ ∩ B′ = ∅. We also claim that (A′, B′) is a cut. To verify this claim we need to show that d ∈ B′. However this is evident by virtue of the fact that the algorithm stops when the sink is not labelled.

Now recall from Lemma 8.1 that we must necessarily have A(f) = f(A′, B′). However, by virtue of the fact that the capacity of arcs from labeled nodes in A′ to unlabeled nodes in B′ must be saturated (this follows by definition of the terminal labeling), it follows that

f(A′, B′) = ∑_{i∈A′} ∑_{j∈B′} f(i, j) = ∑_{i∈A′} ∑_{j∈B′} γ(i, j) = γ(A′, B′).

Hence the discussion preceding this theorem implies the required result.

8.6 Example revisited

In the previously discussed example we computed the maximal flow to be equal to 12. In the diagram below we see all the possible cuts of the network and the associated capacities of the cuts.

1  2 3  45  6 0  8 2 0 8 0  00  54 4 4  4  0  5 2 0  6

Note that all the cuts are described as follows (we leave the reader to compute their capacites and verity that 12 is indeed the minimal flow across a cut).

{1}{2, 3, 4, 5, 6}
{1, 2}{3, 4, 5, 6}
{1, 3}{2, 4, 5, 6}
{1, 2, 3}{4, 5, 6}
{1, 2, 4}{3, 5, 6}
{1, 3, 5}{2, 4, 6}
{1, 2, 3, 4}{5, 6}
{1, 2, 3, 5}{4, 6}∗
{1, 2, 3, 4, 5}{6}

The cut marked with a ∗ is a minimal cut.

It is interesting to note that what is essentially going on here is Duality in disguise. Theorem 8.2 essentially says that if we think of the optimisation of flow over our given network as a linear programming problem, then when the the primal problem has a solution then so does the dual and their objectives are equal. The dual problem clearly manifests itself in cuts over the network. We leave the reader to puzzle over what precisely the dual problem is!

8.7 A second worked example

We are required to establish the optimal flow in the the following network

1  2  3 4  5 7 6  8 0  00  8 7  3 0  7  0  9 0  5 8  6 0 0 4 3 5 0 9  0 0 0 0  4

Note that this network has a source at node 1 and a sink at note 8. We identify the following paths and flows.

1 → 2 → 5 → 8 flow 5
1 → 3 → 6 → 8 flow 4
1 → 4 → 7 → 8 flow 3
1 → 2 → 6 → 8 flow 2

which results in the following residual capacities.

1  2  3 4  5 7 6  8 0  00  8 7  3 0  7 0  9 0  5 8 6 0 0 4 3 5 0 9  0 0 0 0 5  2  5  4  5  0 3 6  3  5  3 0 4  1  4  4  4  2 7  0  2  2  6  0 0

With a sharp eye, it should now be clear that we have saturated the network. However, to be sure, we can implement the labelling procedure to obtain the diagram below.

1  2  3 4  5 7 6  8 0  7 8 0 3  0 (1,1) (6,1)  (5,4) 7  0  5  0 0 6  0 3 0 4 4  1 5  4 2  2 3  5 3  6  (5,7) (5,5)  (2,2)

Since the sink at 8 is not labelled the network must be saturated. The proof of the maximum flow/minimal cut theorem tells us that the minimal cut corresponds to the the case that A is the source together with all labelled nodes and B is the remaining nodes. We have indicated this cut too in the diagram which corresponds to {1, 3, 4, 5, 6, 7}{8}. Note that the flow from the source is equal to 7 + 4 + 3 = 14. On the other hand the sum of the capacities across the indicated cut are equal to 5 + 6 + 3 = 14.