Download A Simple Simplex Method Code – Programming Project Two | MAT 168 and more Study Guides, Projects, Research Optimization Techniques in Engineering in PDF only on Docsity!
Programming Project Two: A Simple Simplex Method Code
(Due Friday May 22, 2008)
Introduction
The object of this project is to implement your own simplex method code in Matlab. There are two big complications and one minor complication in “industrial strength” linear programming code: numerical stability and the preservation of sparsity are the major problems (and often contrary goals) and initialization is a minor problem. As novices, our first goal is to get a basic implementation working without worrying too much about these complications.
The Objective
You will be implementing three Matlab functions which together solve linear programs of the form
min: cx subject to: Ax = b
(1) x ≥ 0 ,
where A is a m × n matrix, x is a row vector of length n, and c is a column vector of length n. You may not assume that A has any special form; in particular, you should not assume that A is a constraint matrix arising from the introduction of slack variables of the form
A I
You will be implementing the “revised simplex method,” the algorithm described in Chapter 6 of the textbook, rather than the more familar tableau method. The two algorithms perform equivalent operations in the same order, but differ in what they store and in the nitty gritty details of the computations at each step. Although the tableau method is easier for human beings, there are numerous advantages to using the revised simplex method for a computer implementation.
A More General Initialization Procedure
The first initialization procedure we studied, which uses an auxillary linear program to find a basic feasible vector, is inadequate for the initialization of problems of the form (1). We could use the dual initialization method we studied, but that would require the implementation of the dual simplex method, which we would like to avoid. Instead, we will introduce another, slightly more general, initialization procedure. The downside will be that the new procedure will not always work.
Beginning with a problem of the form (1), we introduce m new variables y 1 ,... , ym and solve the auxillary problem
min: y 1 + y 2 +... + ym subject to: Ax + y = b x ≥ 0
(2) y ≥ 0.
It is clear that we can always find a basic feasible vector for the new auxillary problem (in particular, we choose the yj to be the basic variables, and the corresponding columns of the constraint matrix are the identity).
The original problem has an initial basic feasible vector if and only if the problem (2) has a solution such that y 1 = y 2 =... = yn = 0.
There is, however, one complication. If we arrive at a solution of (2) for which all of the basic variables are x′s, then we have a initial basic feasible vector for the initial problem. If, on the other hand, some of the y′ j s are basic feasible for our solution of (2), then we would have to perform further manipulations to get a basic feasible vector for the original problem. We won’t worry about this. Our initialization procedure will attempt to find the solution to the auxillary problem (2). If we find such a solution and the basic variables include the y′s, we will simply report that our initialization procedure failed. 1
Step One
The first step in the project is to implement a Matlab function called simplex step for executing a single step of the revised simplex method. We will keep track of the current basic feasible vector with three variables: iB, iN, and xB. The vector iB will hold the indices of the current set of basic variables, iN will hold the indices of the current set of nonbasic variables, and xB will hold the values of the basic variables.
The function simplex step should be placed in a file simplex step.m and it should have the calling sequence:
function [istatus,iB,iN,xB] = simplex_step(A,b,c,iB,iN,xB,irule) % % Take a single simplex method step for the linear program % % min: c*x % ST: Ax=b % x>=0, % % where A is an (m,n) matrix. % % That is, given a basic feasible vector vector described by the % variables iB,iN,xB return the values of iB,iN, and xB corresponding to % the adjacent basic feasible vector arrived at via a simplex method % step. % % Input Parameters: % % A - (n,m) constraint matrix % b - (m,1) POSITIVE vector appearing in the constraint equation above % c - (1,n) vector giving the coefficients of the objective function % % iB - (1,m) integer vector specifying the indices of the basic % variables at the beginning of the simplex step % iN - (1,n-m) integer vector specying the indices of the nonbasic % variables at the beginning of the simplex step % xB - (m,1) vector specifying the values of the basic % variables at the beginning of the simplex step % % irule - integer parameter speciying which pivot rule to use: % irule = 0 indicates that the smallest coefficient rule should be % used % irule = 1 indicates that Bland’s rule should be used % % Output Parameters: % % istatus - integer parameter reporting on the progress or lake thereof % made by this function % istatus = 0 indicates normal nondegenerate simplex method step % completed % istatus = 16 indicates the program is unbounded % istatus = -1 indicates an optimal feasible vector has been % found % % iB - integer vector specifying the m indices of the basic variables % after the simplex step
The final step in the project is the implementation of a Matlab function simplex method which used the preceding two functions in order to compute a solution to a linear program.
The calling sequence for the function simplex method, which should reside in the file simplex method.m, is as follows:
function [istatus,X,eta,iB,iN,xB] = simplex_method(A,b,c,irule) % % Find a basic optimal solution for the linear program % % min: c*x % ST: Ax=b % x>=0, % % where A is an (m,n) matrix. % % Input Parameters: % % A - (n,m) constraint matrix % b - (m,1) POSITIVE vector appearing in the constraint equation above % c - (1,n) vector giving the coefficients of the objective function % % irule - integer parameter speciying which pivot rule to use: % irule = 0 indicates that the smallest coefficient rule should be % used % irule = 1 indicates that Bland’s rule should be used % % Output Parameters: % % istatus - integer parameter reporting the results obtained by % this function % istatus = 0 indicates normal completeion (i.e., a solution % has been found and reported) % istatus = 4 indicates the program is infeasible % istatus = 16 indicates the program is feasible but our initialization % procedure has failed % istatus = 32 indicates that the program is unbounded % % X - vector of length n specifying the solution % eta - the minimum value of the objective function % iB - integer vector specifying the m indices of the basic variables % after the simplex step % iN - integer vector specifying the n-m indices of the nonbasic % variables after the simplex step % xB - vector of length m specifying the values of the basic % variables after the simplex step %
Extra Credit: Foolproof Initialization
Using any technique you would like, modify your function simplex init so that it will not fail unless the linear program is infeasible. Solving a modified dual of the input problem is one way to do this. Modifying the algorithm we used to initialize is another.
