Linear Programming and Simplex Method Review - Prof. James C. Bremer, Study notes of Optimization Techniques in Engineering

Download Linear Programming and Simplex Method Review - Prof. James C. Bremer and more Study notes Optimization Techniques in Engineering in PDF only on Docsity! Linear Programming Redux Jim Bremer May 12, 2008 The purpose of these notes is to review the basics of linear programming and the simplex method in a clear, concise, and comprehensive way. The book contains all of this material, but it is unfortunately spread across several chapters and, in my opinion, confusing in part. These notes are a bit more demanding then the book; if you can read and thouroughly understand them then you are doing very well in the course. The difficulty that most students will encounter in these notes is that they assume a thorough knowledge of linear algebra. I have said several times in class and will repeat now: there is, in my experience, almost no subject with a better reward/effort ratio than linear algebra (statistics being perhaps more rewarding still per unit of effort). A thorough, complete knowledge of elementary linear algebra will serve anyone in a technical field well indeed. Note that this discussion is not comprehensive; in particular, I have omitted any discussion of numerical implementation and stability, as well as any discussion of the geometry underlying linear programs. 1. Underdetermined Systems We begin by reviewing underdetermined systems of linear equations. Let A be an m×n matrix with n > m. Then the system of equations (1) Ax = b is underdetermined — it has more variables than equations. Suppose we were to row reduce the augmented matrix ( A | b ) in order to solve the system (1). Because the matrix A has more columns than rows we will necessarily have free variables (if this isn’t clear to you, work out an example now). So we cannot expect a unique solution for this system. Of course, it is also possible that (1) has no solutions at all (if this isn’t clear, stop and find an example). By making an additional assumption on the matrix A, we can ensure that the system (1) has a solution for every b ∈ Rm. In particular, we will generally assume that the matrix A has rank m. Recall that the rank of a matrix is the dimension of its column and row spaces (these two dimensions are equal). So the assumption that A has rank m means that there are no redundant equations in the system (1), or, equivalently, that the column space of A must be all of Rm — in other words, for every b there is some solution of the system (1). It is also equivalent to saying that some set of m columns of the matrix A forms a basis for Rm. We will be particularly interested in a distinguished class of solutions of (1). We say that x ∈ Rn is a basic solution for the system (1) if there is a set of indices {i1, . . . , im} ⊂ {1, 2, . . . , n} such that: 1. xj = 0 for j /∈ {i1, . . . , im}, and 2. the corresponding columns Ai1 , . . . , Aim of the constraint matrix A form a basis for Rm. Note that the columns of the matrix corresponding to a basic solution form a submatrix B of A which is invertible; indeed, a set of m vectors v1, . . . , vm in Rm form a basis if and only if the matrix whose columns are v1, . . . , vm is invertible (check this!). Basic solutions for (1) are (relatively) easy to find: simply find a submatrix B consisting of the columns Ai1 , . . . , Aim of A which form a basis (perhaps via Gram-Schmidt Orthogonalization), solve the m × m system Bz = b for z, and form the basic solution x with entries xj = { zij if j ∈ {i1, . . . , im} 0 otherwise. Since the matrix B is invertible, there is one and only one basic solution x associated with the columns i1, . . . , im. We will call that solution the basic solution associated with the columns i1, . . . , im. Alternately, we will say that x is the basic solution corresponding to the invertible submatrix B. 1 These facts are sufficiently important that we will repeat them in Lemma form; we have proved the following: LEMMA 1.1. Suppose that A is an m× n matrix, n ≥ m, of rank m and further suppose that the columns i1, . . . , im of A form a basis for Rm. Then there is a unique vector x, which we will call the basic vector for the columns i1, . . . , im, such that 1. Ax = b 2. xj = 0 for all j /∈ {i1, . . . , im}. 2. Linear Programs A linear program is any optimization problem of the form max: ctx subject to: Ax = b(2) x ≥ 0, where A is an m× n matrix, n > m, of rank m, c is a given row vector of length n, b is a given row vector of length m, and x is a row vector of unknowns of length n. Remark 2.1. Note that the assumption that A has rank m ensures that there are no redundant constraints in (2). It differs from the usual definition of linear program only in that it excludes certain infeasible problems from consideration; e.g., the problem max: x1 + x2 subject to: x1 − x2 = −1 − x1 + x2 = −1 x1, x2 ≥ 0 is excluded. We will call a vector x which satisfies the constraints Ax = b and x ≥ 0 a feasible vector for the program (2). Moreover, a feasible vector x for which cx obtains a maximum among the set of all feasible vectors is called an optimal feasible vector for the program (2), or, more simply, a solution to (2). The linear function cx is called the objective function for the program (2). The first thing we should note about the problem (2) is that the constraint equations Ax = b are under- determined. Because of our assumption about the rank of A, we expect there to be an infinite number of solutions for this system of equations. That only stands to reason, since the problem is asking us to find from among the infinite number of x′s satisfying the constraints, a single x such that cx is maximized. Our second immediate observation is that, despite the fact that the system Ax = b has an infinite number of solutions, it is not necessarily true that there exists even a single x satisfying all of the constraints. For example, this is the case for linear program max: x1 + x2 + x3 subject to: x1 + x3 = −1(3) x2 − x3 = −1 x1, x2, x3 ≥ 0. We call a linear program for which there are no vectors x satisfying both Ax = b and x ≥ 0 infeasible. Clearly, infeasible linear programs do not admit a solution. There is another way in which a linear program can fail to have a solution. There might not be a maximum value of the objective function cx. In such cases, we say that the program (2) is unbounded. We close this section with a final observation about (2): if the set of x which satisfy the constraints is bounded (as a set in Rn), then the problem (2) cannot be unbounded. As with all sufficiently elementary facts, there are many different incantations we can invoke to see that this is so (e.g., a continuous function on a compact set in Rm obtains its maximum). This implies, for example, that a linear program of the We will call the form of the linear program (8) the tableau associated with the basic vector x̃. We say that (8) is a feasible tableau if x̃B ≥ 0 and it is an optimal tableau if x̃ is an optimal feasible vector. This form has several useful properties: 1. It explicitly exhibits the value of x̃ in the constants that appear on the right hand side of the constraint equation. 2. It is very easy to determine if x̃ is a feasible vector for the linear program: if x̃B ≥ 0, then x̃ is feasible, otherwise it is not. 3. Assuming x̃ is feasible, it is also easy to determine if x̃ is a solution to the LP. Although this form of the LP is associated with a particular vector x̃, it is, in fact, equivalent to the original LP (5). Thus, we can test the feasibility of any vector x and evaluate the objective function cx using (8). In particular, we examine the vector c̃. If none of the entries of c̃ are positive, then there can be no other basic feasible vector x for which cx takes on a value larger than ξ̃ (note the feasibility assumption is important here because it means the entries of xN we are testing are nonnegative). Remark 4.1. We are used to writing down the tableau form of a linear program in a table in the following manner: nonbasic basic B−1N Im x̃B −c̃ 0 ξ̃ 5. The Simplex Method: Overview In this section we introduce the simplex method, an improved algorithm for solving linear programs max: ctx subject to: Ax = b(9) x ≥ 0, where A is an m× n matrix, n > m, of rank m. Suppose that x̃ and ỹ are basic vectors for the LP (9), and further suppose that X̃ is associated with the columns {i1, . . . , im} of A and ỹ with the columns {j1, . . . , jm}. Then we say that x̃ and ỹ are adjacent if the sets {i1, . . . , im} and {j1, . . . , jm} have m− 1 elements in common. In other words, two basic vectors x̃ and ỹ are adjacent if their associated bases differ by one element. The idea behind the simplex method is quite simple: it is an iterative method, which starting with an initial basic feasible vector, moves from one basic feasible vector to another adjacent one in an effort to increase the value of the objective function. For each iteration j, we form the tableau max: ξ̃j + c̃jtxN subject to: ( I B−1N ) ( xB xN ) = x̃j(10) xB , xN ≥ 0, for the associated basic feasible vector x̃j (there is a slight abuse of notation here: we are identifying the basic part of x̃j with x̃j). If x̃j is a solution for the linear program, it is evident from the coefficients of c̃j (see the observations about the tableaus above). Otherwise, we can choose a nonbasic variable in xN whose coefficient is positive, called the entering variable, which we will move into the basis. Of course, we must swap this variable with a properly chosen basic variable from xB , called the leaving variable, in order to maintain a basis. In most cases, this results in an increase in the objective function. Remark 5.1. It is very important to understand what having an initial basic feasible vector entails. It means not only do we have a vector x̃ with no more than m nonzero entries which is feasible (x̃ ≥ 0 and Ax̃ = b), but also that the associated columns of the constraint matrix A form a basis. Not any old basis of columns of A will do — we need a basis consisting of columns of A such that the associated basic solution is nonnegative! 6. A Simplex Method Step In this section, we describe in detail what a single step of the simplex method entails. First, we assume that the underlying linear program is of the form max: ctx subject to: Ax = b(11) x ≥ 0, where A is an m × n matrix, n > m, of rank m. Next, we suppose that at the jth iteration we have the tableau max: ξ̃ + c̃txN subject to: ( I B−1N ) ( xB xN ) = x̃B(12) xB , xN ≥ 0, which is associated with the basic feasible vector x̃. Throughout this section, whenever we have a vector x ∈ Rn, we will denote by xB the part of x corresponding to the basic variables for the tableau (12) and by xN the part of x corresponding to the nonbasic variables for this tableau. Let us be explicit about the dimensions of the various blocks of the constraint matrix in (12): the matrix B is an m×m matrix, the idenity block in (12) is also m×m, and B−1N is m× (n−m). This, of course, means that xB represents the m basic variables and xN represents the n−m nonbasic variables. Our goal is to find an adjacent basic feasible vector x̂ for which the objective function is larger than x̂. Recall that a basic vector adjacent to x̃ is a basic vector which shares m− 1 of its basic columns from the matrix A in common with x̃. So, in other words, we are going to proceed by switching one basic variable with a nonbasic variable. We first choose the nonbasic variable from xN which will become basic in the next iteration. We do this by examining the coefficients of c̃ and looking for one which is positive. Increasing the corresponding nonbasic variable will increase the objective function (note that if no such nonbasic variable can be found, then x̃ is already optimal). We call the chosen nonbasic variable the entering variable since it will be entering the basis for the next iteration. Of course, when we increase the entering variable, we will need to adjust all of the basic variables in order to ensure that the constraint equation Ax̂ = b is still satisfied. We are proposing to change the nonbasic part of x̃ by letting x̂N = tek, where t is a suitably chosen positive real number and ek is the standard basis vector corresponding to the nonbasic variable we have chosen (i.e., it has a one in the correct position and zeros everywhere else). Since the equation (13) x̂B + B−1Nx̂N = xB must hold for any solution of the constraint equations, we have x̂B = xB − tB−1Nek. We define the vector ∆x by ∆x = B−1Nek. This is the direction of the change in x̂B . Notice that Nek is just the column of the original constraint matrix A which corresponds to the entering variable. We would like to increase t as much as possible — recall that the objective function will increase by c̃kt. The value of t is usually constrained, however, by the requirement that x̂B ≥ 0. In particular, when the value of coordinate i of ∆x is greater than 0, then the corresponding equation in (13) gives us a bound on t. So we let t be min i (xB)i (∆x)i , where the minimum ranges over all indices i for which (∆x)i is positive. If there are no such indices, then the problem is unbounded. That is, we can increase t to infinity without violating the constraints, and as we do so, the objective function will increase without bound. That is, we increase t until one or more of the variables in x̂B becomes zero. We now choose from among the set of x̂B a leaving variable. That is, a variable to move out of the basic set xB and into the nonbasic set xN at the next iteration. We have now computed the value of the basic feasible vector at the next iteration: x̂B and x̂N . It only remains to do bookkeeping: to update the list of variables that are basic and nonbasic for the next iteration. Remark 6.1. Note that it is never necessary while performing the simplex method algorithm to actually compute the form of the constraint matrix which appears in (12). Indeed, it is more convenient to leave the constraint matrix in its original form (11), and compute the direction ∆x by solving a system of equations. Remark 6.2. Assuming that we do not run into unboundedness, the new set of basic variables do indeed correspond to a set columns of A which forms a basis. To see this, note that the vector ∆x records the coefficients of the entering column with respect to the current basis. Because the entry of ∆x cooresponding to the leaving variable is nonzero (by definition), it means that the entering column is not in the span of the basis columns excluding the leaving column. This implies the the resulting set of columns is a basis. 7. Degeneracy 8. Initialization I There is an obvious unresolved difficulty with the simplex method: we need to find a basic feasible vector in order to start the simplex method (if we start with an infeasible vector, we might never find a solution to the LP). As was remarked upon above, this is a nontrivial task which entails more than finding an invertible submatrix of the constraint matrix A. We have been spoiled by the textbook’s practice of writing linear programs in the form: max: ctx subject to: Ax ≤ b(14) x ≥ 0. This form is particularly nice for initialization, because once m slack variables have been introduced (14) takes on the form max: ctx subject to: ( A I ) ( x w ) = b x, w ≥ 0. This makes it trivial to find an initial basis for the constraint matrix (we simply pick the identity submatrix I). Of course, the resulting basic vector is feasible only if b ≥ 0, so there is still work to be done even in this case, but nonetheless it is easier than the initialization of a general LP. Instead of considering this simple case, we will once again the consider the general LP max: ctx subject to: Ax = b(15) x ≥ 0, where A is an m × n matrix, n > m, of rank m. Without loss of generality, we may assume that b ≥ 0 (if bi < 0 for some i then we can multiply the corresponding constraint by −1). In order to initialize (15), we introduce the auxiliary linear program min: y1 + y2 + . . . + ym subject to: Ax + y = b(16) x, y ≥ 0. It should be clear that the original problem (15) is feasible if and only if (16) has a solution such that y1 = y2 = . . . = ym = 0. Moreover, it is easy to find an initial feasible vector for (16): we simply let y = b and x = 0. Then the corresponding submatrix of the constraint matrix is the identity (and so clearly invertible) and because b can be assumed to satisfy b ≥ 0, this vector is feasible. 10. The Dual Tableau The primary observation of this section — and the most important single fact in the theory of duality — is that the dual of a tableau for the primal problem is a itself a tableau for the dual problem. This fact and the form of the two corresponding tableaus imply all sorts of nice results (e.g., strong duality). To see that this is so, suppose that max: ξ̃ + c̃txN subject to: ( I B−1N ) ( xB xN ) = x̃(27) xB , xN ≥ 0, is a tableau associated with a basic feasible vector x̃ for a primal problem. Using what we learned in the last section (check this!), we find that the dual of (27) is min: ξ̃ + c̃ty subject to: ( I (B−1N)t ) y ≥ ( 0 c̃ ) (28) y free, Of course, we can rewrite this as min: ξ̃ + x̃ty subject to: (B−1N)ty ≥ c̃(29) y ≥ 0. If we introduce slack variables yB into (29), and rename the variables already present yN , then we arrive at the linear program min: ξ̃ + x̃tyB subject to: ( I −(B−1N)t ) ( yN yB ) = −c̃(30) yB , yN ≥ 0. We call this linear program the dual tableau corresponding to the primal tableau (27). We can immediately make the following observations: 1. There is a correspondence between the dual variables yB and the primal basic variables xB , and between the dual variables yN and the primal variables xN . 2. As written, the basic variables for the tableau (30) are the yN and not the yB ! So we have chosen to emphasize the relationship mentioned above rather than emphasize which variables are basic and which are nonbasic for the dual. 3. There is clearly a basic vector ỹ associated with the dual tableau: yB = 0 and yN = c̃. 4. The basic vector ỹ is not necessary feasible for the dual — we do not expect c̃ ≥ 0. We will call the basic vector ỹ associated with the dual tableau (30) the dual basic vector of x̃. This is a very important idea: we associate with each basic vector of the primal a basic vector of the dual. We close this section by proving the Strong Duality Theorem. Thanks to our discussion of tableaus and dual tableaus, the proof is trivial. THEOREM 10.1. (Strong Duality) Suppose that x̃ is a basic solution for the primal problem (e.g., a basic optimal feasible vector) max: ctx subject to: Ax ≤ b(31) x ≥ 0. Then the corresponding basic vector ỹ for the dual min: bty subject to: Ax ≥ c(32) y ≥ 0. of (31) is a basic solution (e.g., a basic optimal feasible vector for the dual). Moreover, (33) ctx̃ = btỹ. Proof: Write the tableau max: ξ̃ + c̃txN subject to: ( I B−1N ) ( xB xN ) = x̃B(34) xB , xN ≥ 0, for the basic vector x̃. That x̃ is feasible means that xB ≥ 0 and that it is optimal means that c̃ ≤ 0. Now the tableau min: ξ̃ + x̃BtyB subject to: ( I −(B−1N)t ) ( yN yB ) = −c̃(35) yB , yN ≥ 0. is the corresponding dual tableau — the one associated with the basic vector ỹ. That c̃ ≤ 0, means that ỹ is feasible and x̃B ≥ 0 implies that it is optimal (since the dual is a minization problem). So ỹ is a basic optimal vector for the dual. Equation (33) now follows because we have yB = 0 and xN = 0 for the pair of vectors x̃ and ỹ. QED. 11. Initialization II We now discuss the use of the dual problem to effect initialization. 12. Problems 1. Write down a feasible, bounded linear program and solve it. 2. What is the dual of the linear program max: ctx subject to: Ax = b l ≤ x ≤ u? Is it possible for the dual of this LP to be infeasible? What does that say about the primal, assuming it is feasible? 3. Why does the simplex method have to start at a basic feasible vector? Can we begin with a nonbasic feasible vector? How about an infeasible basic vector? 4. How could you go about making the first initialization procedure foolproof? 5. What is the dual of the linear program max: ctx subject to: Ax = b x free ? 6. Show that the proof in Section 9 of the fact that the dual of a program with an affine objective function has an affine objective function still holds if the equality constraint z = 1 is changed to z ≤ 1. 7. Consider the first algorithm for solving linear programs, which was introduced in Section 3. This algorithm can be easily modified to solve problems given in the form max: ctx subject to: Ax = b x free , where x is not constrained to be nonnegative. (It should be clear that any linear program can be written in this form — if not check it!!!) Indeed, the modified algorithm is clearly simpler than the algorithm as written. Why is it difficult to modify the simplex method to operate on programs written in this form?