����������� ������ 54 (2002), 163–169

UDK 512.644

���������� ������ ��� research paper

A MODIFICATION OF REVISED SIMPLEX METHOD

Marko D. Petković, Predrag S. Stanimirović and Nebojša V. Stojković

Abstract. We introduce a modification of the revised simplex method, which accelerates the process of finding the first basic feasible solution in the two phases revised simplex method. We report computational results on numerical examples from Netlib test set.

1. Introduction

Consider the linear program in the general form

Maximize  f(x) = f((x_{N})_{1},…, (x_{N})_{n_{1}}) = ∑_{i=1}^{n_{1}} c_{i}(x_{N})_{i} + d

subject to  N_{i}^{(1)} : ∑_{j=1}^{n_{1}} a_{ij}(x_{N})_{j} ≤ b_{i}, i = 1,…, r

N_{i}^{(2)} : ∑_{j=1}^{n_{1}} a_{ij}(x_{N})_{j} ≥ b_{i}, i = r + 1,…, s

J_{i} : ∑_{j=1}^{n_{1}} a_{ij}(x_{N})_{j} = b_{i}, i = s + 1,…, m

(x_{N})_{j} ≥ 0, j = 1,…, n_{1}.

(1.1)

We change each inequality of the form N_{i}^{(1)} (LE constraint) into an equality by adding a slack variable x_{B,i}:

N_{i}^{(1)} : ∑_{j=1}^{n_{1}} a_{ij}(x_{N})_{j} + (x_{B})_{i} = b_{i}, i = 1,…, r.

Also, we transform each inequality of the form N_{i}^{(2)} (GE constraint) into an equality by subtracting a surplus variable (x_{B})_{i}:

N_{i}^{(2)} : ∑_{j=1}^{n_{1}} a_{ij}(x_{N})_{j} − (x_{B})_{i} = b_{i}, i = r + 1,…, s.

AMS Subject Classification: 90 C 05

Keywords and phrases: Linear programming, revised simplex method, MATHEMATICA

Communicated at the 5th International Symposium on Mathematical Analysis and its Applications, Niška banja, Yugoslavia, October, 2–6, 2002.

After these transformations we get the equivalent linear program in the standard form

Maximize c_{1}x_{1} + ⋯ + c_{n}x_{n} + d = c^{T} x + d, c_{i} = 0, i > n_{1}
subject to Ax = b,
b = (b_{1},…, b_{m}),
x=((x_{N})_{1},…, (x_{N})_{n}, (x_{B})_{1},…, (x_{B})_{s}),
(x_{N})_{j} ≥ 0, j = 1,…, n, (x_{B})_{i} ≥ 0, i = 1,…, s,

(1.2)

where the real matrix A is of order m × n, for n = n_{1} + s.

The paper is organized as follows. In the second section we describe two-phases revised simplex method which is analogous with corresponding two-phases simplex method from [3], [4], [6] and [7]. In the third section an algorithm is described for finding the first basic feasible solution. In the last section we report some numerical experience with several Netlib test problems.

2. The revised simplex method

By A_{i} we denote i-th column of A, i = 1,…, n. We suppose that matrix A is of full rank (rank(A) = m), so there exists a base {A_{i_{k}}, k = 1,…, m} in ℝ^{m}. The matrix B = [A_{i_{1}},…, A_{i_{m}}] containing corresponding columns from A is called the basic matrix. Variables x_{i_{1}},… x_{i_{m}} are called basic variables, and a vector containing basic variables is denoted by x_{B}. The m × (n − m) matrix containing columns from A corresponding to nonbasic variables is denoted by N, and the (n − m)-vector containing nonbasic variables is denoted by x_{N}. Columns T_{i} of Tucker table T are representations of vector A_{i} in the base B, i.e. T = B^{−1}N. Here the objective function is represented as equality c_{1}x_{1} + ⋯ + c_{n}x_{n} + d = −x_{c} where x_{c} is a basic variable representing a value of the objective function. After corresponding modifications of the corresponding algorithm from [3], [4], [6], [7], we state the following two phases revised simplex algorithm. In this algorithm the notation A^{i} means i-th row of the matrix A.

Algorithm 1. (Revised simplex method for basic feasible solution) Step S1A. Reconstruct the vector b_{T} = B^{−1}b, Tucker table T = B^{−1}N and the objective function c_{T} = T^{n+1}. Step S1B. If (c_{T})_{1},…, (c_{T})_{n} ≤ 0, then the basic solution is an optimal solution. Step S1C. Choose an arbitrary (c_{T})_{j} > 0. Step S1D. If t_{1j},…, t_{mj} ≤ 0, stop the algorithm. Maximum is +∞. Otherwise, go to the next step. Step S1E. Compute

min_{1≤i≤m} \{\frac{(b_{T})_{i}}{t_{ij}} | a_{ij} > 0\} = \frac{(b_{T})_{p}}{t_{pj}} ,

replace the nonbasic and basic variables (x_{N})_{j} and (x_{B})_{p}, respectively, and go to Step S1A.

If the condition (b_{T})_{1},…, (b_{T})_{m} ≥ 0 is not satisfied, we use the following algorithm to search for the first basic feasible solution. Algorithm in this version is not found in the literature.

Algorithm 2. (Find the first basic feasible solution). Step S1. Reconstruct the vector b_{T} = B^{−1}b and Tucker table T = B^{−1}N. Step S2. Select the last (b_{T})_{i} < 0. Step S3. If t_{i1},…, t_{in} ≥ 0 then STOP. Linear program cannot be solved. Step S4. Otherwise, find t_{ij} < 0, compute

min_{k>i} (\{\frac{(b_{T})_{i}}{t_{ij}} \} ∪ \{\frac{(b_{T})_{k}}{t_{kj}} | t_{kj} > 0\}) = \frac{(b_{T})_{p}}{t_{pj}} ,

replace the nonbasic variable (x_{N})_{j} and the basic variable (x_{B})_{p} and go to Step S1. We use the last t_{ij} < 0.

3. Modification of algorithm for finding first basic solution

The problem of the replacement of a basic and a nonbasic variable in the general revised simplex method is contained in Step S4. We observed two drawbacks of Step S4.

1. If p = i and if there exists an index y < i = p such that

\frac{(b_{T})_{y}}{t_{yj}} < \frac{(b_{T})_{p}}{t_{pj}} , (b_{T})_{y} > 0, t_{yj} > 0,

in the next iteration (x_{B})_{y} becomes negative:

(x_{B})_{y}^{1} = (b_{T})_{y}^{1} = (b_{T})_{y} − \frac{(b_{T})_{p}}{t_{pj}} t_{yj} < (b_{T})_{y} − \frac{(b_{T})_{y}}{t_{yj}} t_{yj} = 0.

2. If p > i, in the next iteration (b_{T})_{i}^{1} is negative:

\frac{(b_{T})_{p}}{t_{pj}} < \frac{(b_{T})_{i}}{t_{ij}} ⇒ (b_{T})_{i}^{1} = (b_{T})_{i} − \frac{(b_{T})_{p}}{t_{pj}} t_{ij} < 0.

We can find solution of this problem in two ways.

For fixed (b_{T})_{i} try to find t_{ij} < 0 such that

min_{k} (\{\frac{(b_{T})_{i}}{t_{ij}} | t_{ij} < 0\} ∪ \{\frac{(b_{T})_{k}}{t_{ij}} | t_{kj} > 0, (b_{T})_{k} > 0\}) = \frac{(b_{T})_{i}}{t_{ij}} .  (3.1)

In this case, it is possible to choose t_{ij} for the pivot element and obtain a new solution with a smaller number of negative coordinates (b_{T})_{i} in B_{T}.

This way is suitable for the basic simplex but not for the revised simplex because we need to reconstruct columns corresponding to all t_{ij} < 0.

In the second one, the main idea is to find a minimum

\frac{(b_{T})_{p}}{t_{pj}} = min \{\frac{(b_{T})_{k}}{t_{kj}} | (b_{T})_{k} < 0, t_{kj} < 0, k = 1,…, m\}

and if relation (3.1) is valid then choosing (b_{T})_{p} > 0 as a pivot element we obtain a new solution with a smaller number of negative coordinates (b_{T})_{i}.

In order to do that, we propose a modification of Step S4. This modification follows from the following lemma.

Lemma. Let the problem (1.2) be feasible and let x be the basic infeasible solution with q negative coordinates. Then there exists t_{ij} < 0. Also, in the following two cases:

a) q = m,

b) q < m and

Neg = \{\frac{(b_{T})_{k}}{t_{kj}} | (b_{T})_{k} < 0, t_{kj} < 0, k = 1,…, m\} ,

P os = \{\frac{(b_{T})_{k}}{t_{kj}} | (b_{T})_{k} > 0, t_{kj} > 0, k = 1,…, m\} ,

min (Neg ∪ P os) = \frac{(b_{T})_{r}}{t_{rj}} ∈ Neg  (3.2)

it is possible to produce the new basic solution with at most q−1 negative coordinates in only one iterative step of the revised simplex method, if we choose t_{rj} for the pivot element, i.e. if we replace nonbasic variable (x_{N})_{j} by the basic variable (x_{B})_{r}.

Proof. Consider an arbitrary (b_{T})_{i} < 0. If the condition t_{ij} > 0 is satisfied for all j = 1,…, n, then

−(x_{B})_{i} = t_{i1}(x_{N})_{1} + ⋯ + t_{im}(x_{N})_{m} − (b_{T})_{i} > 0.

Hence, for (x_{N})_{1} > 0,…, (x_{N})_{n} > 0 we have (x_{B})_{i} < 0. That is a contradiction because problem (1.2) is feasible.

We now prove the second part of theorem.

a) If q = m, for an arbitrary pivot element t_{ij} < 0 we get a new solution with at least one coordinate positive:

(x_{B})_{i}^{1} = (b_{T})_{i}^{1} = \frac{(b_{T})_{i}}{t_{ij}} > 0.

b) Assume now that the conditions q < m and (3.2) are satisfied. Choose t_{rj} for the pivot element.

For (b_{T})_{k} > 0 and t_{kj} < 0 it is obvious that

(b_{T})_{k}^{1} = (b_{T})_{k} − \frac{(b_{T})_{r}}{t_{rj}} t_{kj} > (b_{T})_{k} ≥ 0.

For (b_{T})_{k} > 0 and t_{kj} > 0, using \frac{(b_{T})_{k}}{t_{kj}} ≥ \frac{(b_{T})_{r}}{t_{rj}}, we conclude immediately

(b_{T})_{k}^{1} = (b_{T})_{k} − \frac{(b_{T})_{r}}{t_{rj}} t_{kj} ≥ 0.

Hence, all positive (b_{T})_{k}, remain positive. Moreover, for (b_{T})_{r} < 0 we get

(b_{T})_{r}^{1} = \frac{(b_{T})_{r}}{t_{rj}} ≥ 0,

which completes the proof.

Remark 3.1. Let the condition (3.1) is not satisfied, min (Neg ∪ P os) = \frac{(b_{T})_{r}}{t_{rj}} ∈ Neg. If we choose t_{rj} for the pivot element we will also have

(b_{T})_{k}^{1} = (b_{T})_{k} − \frac{(b_{T})_{r}}{t_{kj}} t_{rj} ≥ 0,

for (b_{T})_{k} > 0 and t_{rj} both positive or negative. But for negative (b_{T})_{k} > 0 we can prove similarly that

(b_{T})_{k}^{1} = (b_{T})_{k} − \frac{(b_{T})_{r}}{t_{kj}} t_{rj} < 0,

for both t_{rj} positive or negative. Finally new solution has the same number of negative coordinates q.

Algorithm 3. (Modification of Algorithm 2 ). Step SM1. Let B is a basic and N is a non-basic matrix. Reconstruct vector b_{T} = B^{−1}b. Step SM2. Construct the set

B = {(b_{i})_{1},…, (b_{i})_{q}} = {(b_{i})_{k} | (b_{i})_{k} < 0, k = 1,…, q}.

Step SM3. Select (b_{i})_{s} < 0. Step SM4. Reconstruct i_{s}-th row. If t_{i_{s},1},…, t_{i_{s},n} ≥ 0 then STOP. Linear program is not solvable.

Otherwise, choose first t_{i_{s},j} ≤ 0. Step SM5. Reconstruct j-th column T_{j}. Compute

min_{1≤i≤n} \{\frac{(b_{T})_{i}}{t_{ij}} | (b_{T})_{i}t_{ij} > 0, i = 1,…, m\} = \frac{(b_{T})_{p}}{t_{pj}} ,

replace nonbasic and basic variables (x_{N})_{j} and (x_{B})_{p} and go to Step SM2.

4. Implementation and numerical experience

In the previous algorithms the reconstruction of the Tucker table is based on the matrix inversion. For large numerical examples this method is too slow and numerically unstable. Moreover, the matrix inversion should be avoided, since we require only few rows and columns of Tucker table. In the sequel we will show how to reconstruct a simple row or column of Tucker table using linear equation solver. It is well known that it is much more efficient to solve a linear system than to invert a matrix [2]. Reconstruction of i-th column T_{i} is based on the application of the MATHEMATICA function LinearSolve to solve the system BT_{i} = A_{i} with respect to unknown variables in the vector T_{i}. To find a row T^{i} we need to know (B^{−1})^{i}, so:

T^{i} = (B^{−1})^{i}N.

Vector (B^{−1})^{i} can be found from B^{−1}B = I, which implies (B^{−1})^{i}B = I^{i}. If we apply the transpose to the left- and right-hand side we get: B^{T} ((B^{−1})^{i})^{T} = (I^{i})^{T}. This method is incorporated in the following algorithm :

Algorithm 4. (Reconstruction of row T^{i} of Tucker table). Step 1. Let B be a basic, N a non-basic, and I the identity matrix. Solve the system

B^{T} ((B^{−1})^{i})^{T} = (I^{i})^{T}

with respect to unknown variables contained into the vector ((B^{−1})^{i})^{T}. Step 2. Compute and return T^{i} = (B^{−1})^{i}N.

The revised simplex method and the modification are implemented in the program called RevMarPlex, written in programming language MATHEMATICA [8], and tested on some Netlib test problems. Results are compared with robust LP solver PCx.

From the columns two and three in Table 1 we can see that our results are in accordance with the results obtained by PCx. In the columns four and five of Table 1 we give the number of iterations needed for the construction of the first basis feasible solution, using Algorithm 3 and Algorithm 2, respectively.

Table 1.

After observing these columns, it is easy to see that Algorithm 3 is faster than Algorithm 2 in all cases. In some of them, there were no iterative steps before the first basic feasible solution. Table cells corresponding to these cases are marked by a streak.

Also, RevMarPlex solved a class of extremely ill-conditioned test examples, taken from [1], [5], which PCx is unable to solve (it ends with UNKNOWN status). These test problems are avaiable on the web page www.psmtmath.s5.com. Results are presented in Table 2. Let us observe that numbers of iterations of Algorithm 2 and Algorithm 3 were the same in all cases.

Table 2.

Program RevMarP lex (by means of MATHEMATICA function LinearSolve) indicates that matrices of these examples are ill conditioned (with the condition number between 10^{15} and 10^{20}). Nevertheless RevMarP lex solved these examples giving solutions close to optimal.

REFERENCES

[1] M.D. Ašić and V.V. Kovačević-Vujčić, Ill-conditionedness and interior-point methods, Univ. Beograd Publ. Elektrotehn. Fak., 11 (2000), 53–58.

[2] M.A. Bhatti, Practical optimization with MATHEMATICA applications, Springer Verlag Telos, 2000.

[3] B.D. Bounday, Basic Linear Programming, Edvard Arnold, Baltimore, 1984.

[4] J.P. Ignizio, Linear Programming in Single-multiple-objective Systems, Englewood Cliffs: Prentice Hall, 1982.

[5] V.V. Kovačević-Vujčić, and M.D. Ašić, Stabilization of interior-point methods for linear programming, Computational Optimization and Applications, 14 (1999), 1–16.

[6] M. Sakarovitch, Linear Programming, Springer-Verlag, New York, 1983.

[7] S. Vukadinović, S. Cvejić, Mathematical Programming, Univerzitet u Prištini, Priština, 1996. (In Serbian)

[8] S. Wolfram, Mathematica Book, Version 3.0, Wolfram Media and Cambridge University Press, 1996.

(received 12.12.2002)

M. D. Petković, P. S. Stanimirović, University of Niš, Department of Mathematics, Faculty of Science, Višegradska 33, 18000 Niš, Yugoslavia N. V. Stojković, University of Niš, Faculty of Economics, Trg Kralja Aleksandra 11, 18000 Niš, Yugoslavia E-mail: dexter of nis@inmail.sk, peckois@ptt.yu, nebojsas@orion.eknfak.ni.ac.yu