MANAGING COST UNCERTAINTIES IN TRANSPORTATION AND ASSIGNMENT PROBLEMS

In a fast changing global market, a manager is concerned with cost uncertainties of the cost matrix in transportation problems (TP) and assignment problems (AP).A time lag between the development and application of the model could cause cost parameters to assume different values when an optimal assignment is implemented. The manager might wish to determine the responsiveness of the current optimal solution to such uncertainties. A desirable tool is to construct a perturbation set (PS) of cost coeffcients which ensures the stability of an optimal solution under such uncertainties.


Introduction
The widely-used methods of solving transportation problems (TP) and assignment problems (AP) are the stepping-stone (SS) method and the Hungarian method, respectively.Managerial applications are many and go beyond these prototype problems to include job scheduling, production inventory, production distribution, allocation problems, and investment analysis, among others.As a dual-simplex algorithm, the SS proved successful for solving a TP and became the standard technique for over 50 years.In practice, however, the SS algorithm encounters major obstacles as a solution procedure.It has di culties identifying an initial basic feasible solution, resolving SS degeneracy, and enumerating SS paths.A recent urry of activity has improved the SS algorithm to an extent, enabling it to over-come some of these de ciencies.
(a) Finding an initial basic feasible solution: A prerequisite for SS is a basic feasible solution with a certain number of positive e n tries.The best known method to achieve this is Vogel's, which is not applicable to an unbalanced problem.Goyal 39] modi ed Vogel's method to cover unbalanced cases, and Ramakrishnan 68] improved Goyal's approach.Sultan 76] suggested a heuristic to nd an initial feasible solution that may produce an SS degenerate basis.Kirka and Satir 50] present another heuristic method with minimal improvement o ver Vogel's method for the unbalanced case.
(b) Resolving SS degeneracy: Shafaat and Goyal's 72] is the rst paper with a systematic approach for handling SS degeneracy (i.e., failure of the optimality test.)Their algorithm, however, is limited to a special case of degeneracy.
(c) Enumerating SS paths: Finding SS paths can be very tedious, particularly for larger problems.Intrator and Szwarc 43] extended an existing decomposition method.Their approach m a k es it easier to nd SS paths but requires solving several sub-problems.Later, Wilsdon 82] provided a systematic method for nding the SS path that must be applied in each iteration.
The AP is a special case of TP and is traditionally solved using the Hungarian method.Even though the method is more e cient than the SS method to solve an AP, to "draw the minimum number of lines to cover all the zeros in the reduced cost matrix may not be an easy task " 36].Moreover, the Hungarian method is limited 55] to 2-dimensional problems (e.g., people to projects, jobs to machines).
Although we h a ve gained new insights, these above-mentioned improvements do not lend themselves to a uni ed approach for an e cient solution algorithm.The network-based approaches have similar di culties 25], 6 5 ] .
There are strong motivations for performing perturbation analysis (PA) to deal with a collection of managerial questions related to so-called "what-if" problems.Managers are concerned with unforeseen changes in the input cost parameters.They are likely to be unsure of their current v alues and even more uncertain about their future values at the time when the solution is to be implemented.Uncertainty is the prime reason why P A is helpful in making decisions.We c a n u s e P A t o g i v e us information like: the robustness of an optimal solution critical values, thresholds, and breaking-even values where the optimal strategy changes sensitivity of important v ariables sub-optimal solutions exible recommendations the values of simple and complex decision strategies and the "riskiness" of a strategy or scenario.This information provides systematic guidelines for allocating scarce organizational resources to data collection, data re nement, and prediction activities of the cost parameters, and could be considered as a model validation activity.Examining the e ects of an increase or decrease in unit costs of transportation would help, for example, to negotiate rates with trucking companies and in reacting to changed conditions caused by rate changes.
Although it has long been known that TP and AP can be modeled as linear programs, this is generally not done, due to the relative ine ciency and complexity of the simplex methods (primal, dual, and other variations).These problems are usually treated by one of over 20 specialized algorithms for each (see, e.g., 4], 9], 8], 10], 13], 14], 12], 19], 20], 29], 30], 31], 32], 37], 42], 45], 48], 49], 57], 58], 59], 6 1 ] , 65], 66], 74], 78], 83]).This leads to several di culties.The solution algorithms are not uni ed as each algorithm uses a di erent strategy to exploit the special structure of a speci c problem.Furthermore, a small variation in the problem, such as the introduction of side-constraints, destroys the special structure and requires a new solution algorithm.These algorithms obtain solution e ciency at the expense of managerial insight, as the nal solutions from these algorithms do not easily provide su cient information to perform post optimality analysis for TP and AP.
Another approach is to adapt the simplex network optimization through network simplex.This provides uni cation of the various problems but maintains all ineciencies of the simplex, such a s p i v otal degeneracy, a s w ell as in exibility t o h a n d l e side-constraints.Even ordinary (one-change-at-a-time) sensitivity analysis (OSA) limits, long available in the simplex, are not easily available in network simplex.
Most linear programming books, including management science and operations research books, have extensive discussion of linear programming (LP) sensitivity analysis (SA) but remain silent about the SA and side-constraints for TP and AP.Both methods, the SS and the Hungarian, lack t h e ability to test the validity of the current optimal solution with changes in cost parameters without resolving the problem.Some textbooks encourage the reader to formulate the problem as an LP and solve with a software package.The available computer software packages also have not proven to be e ective in dealing with the SA of TP and AP as an LP.
Netsolve, Net o, and Genos are three network-based software packages that have a subroutine for TP and AP.Netsolve p r o vides OSA for problems with nondegenerate nal tableau.Net o and Genos have no SA capability 4 7 ] .Cplex 22], the current state-of-the-art code for solving network problems, incorporates nothing special for TP or AP's SA.It uses LP-type SA which can be misleading because of a degenerate optimal solution, especially for AP.Another popular software package, QSB, has a module for TP that automatically deletes the last constraint from the input to an LP formulation.(See 63] for details of these packages).Lindo, a general LP package, yields one change at a time for cost coe cient.However, the results are misleading for AP, and could be so for TP with degenerate optimal tableau.To the best of our knowledge, the literature still lacks managing cost uncertainties for a degenerate nal tableau case, which is common in these problems.
We propose a single uni ed algorithm that solves both TP and AP, and also provides useful information to perform cost-sensitivity analysis to a decision maker.Similar to the simplex algorithm and its many v ariants, the proposed solution algorithm is pivotal.The algorithm initiates the solution with a warm-start and does not require any s l a c k/surplus or arti cial variables.Unlike the Hungarian method, the algorithm can solve higher than 2-dimensional AP.T h e proposed solution algorithm also facilitates incorporation of side-constraints, which are frequently encountered in real-life applications.This algorithm makes available the full power of LP's SA extended to handle optimal degenerate solution.The essential calculations of proposed PA i n volve the same data and the same manipulation as the solution algorithm.With little extra computational e ort we can obtain the perturbed solution, which provides the necessary information to construct the "critical" region, that is, the largest set of cost values for which the current solution remains optimal.In contrast to OSA, the proposed PA p r o vides ranges for which the current optimal basis remains optimal, for simultaneous dependent or independent c hanges of the cost coe cients from their nominal values.The preliminary computational results demonstrate that the algorithm is more e cient than the simplex method in terms of number of iterations and size of tableaux.The proposed approach is easy to understand, easy to implement, and thus can serve as e ective tools for solving TP, A P , and related problems in all phases, namely design, analysis, and operation.
The next section presents the solution algorithm which includes an illustrative T P numerical example.This is followed by a discussion on geometric representation and theoretical properties of the proposed solution algorithm in Section 3. Section 4 is devoted to the solution of a 3-dimensional AP.Computational behavior and computer implementation issues are presented in Section 5. General perturbation analysis of cost coe cients is covered in Section 6, followed by special cases of sensitivity analysis in Section 7. Handling the side-constraints and PA of the optimal degenerate case problems are given in Sections 8 and 9 respectively.The last section presents some concluding remarks.The proofs are given in the appendix.

The Solution Algorithm
A shipper having m warehouses with supply S i of goods at his i th warehouse must ship goods to n geographically dispersed retail centers, each with a given customer demand D j which m ust be met.The objective is to determine the minimum possible transportation cost given that the unit cost of transportation between the i th warehouse and the j th retail center is C ij .This problem is known as the TP and has the following standard LP formulation:

Minimize
with the balanced condition and X ij D j S i 0. Note that any u n balanced problem can be converted to a balanced one by adding a row o r c o l u m n .Clearly any o f t h e a b o ve m + n constraints which are redundant can be deleted.When S i and D j = 1, this problem refers to the assignment problem that involves determining the most e cient assignment o f people to projects, jobs to machines, salespeople to territories, contracts to bidders, a n d s o o n .The objective is to minimize total costs or total time for performing the tasks at hand.In a 2-dimensional problem, the optimal solution of AP would be an integer, i.e., only one job or worker is assigned to one machine or project.Also note that in an AP, n = m and in any optimal solution exactly n out of the 2n ; 1 basic variables have v alue 1 and (n;1) 2 non-basic variables along with (n;1) basic variables have v alue 0. Thus, AP as an LP is extremely (primal) degenerate.The PA using any software packages provide unrealistic results since all these packages assume a problem to be non-degenerate.
Although the classical TP also requires that supply and demand (S i and D j ) b e integers (traditionally the number of units), we can relax this condition without loss of generality.Also, we relax the condition of an integer solution.We maintain, however, the condition in an AP that all workers and jobs etc., should be fully allocated, i.e., S i and D j = 1 .The integrality condition of variables in these problems is super uous, imposed by researchers only.The existing solution algorithms associated with TP and AP, such as the SS, Hungarian and network algorithms, can provide a solution only under such limiting condition.In real life situations, however, as in a higher dimensional AP, it is feasible and may be optimal to rotate workers among jobs, or salespeople among territories, and a TP could involve non-discrete demand and supply, s u c h a s w eights and volume.When S i and D j are integers, our algorithm does provide an integer solution.Our approach, however, does not restrict these variables to being discrete.This is a signi cant di erence between existing methods and the approach presented below.One can, however, introduce Gomory's cutting planes 6], 16] of LP to ensure an integer solution, if so desired.

A Pivotal Algorithm
We present steps of a simplex-type pivotal algorithm for the solution of general TP and AP, hereafter referred to as the Push-and-Pull algorithm, to re ect the two main strategies of the approach.The following additional notations of the simplex method are used: The algorithm consists of cost adjustments and initialization followed by t wo phases.
Step 1 : Reduced c ost matrix 1.1 : Row-column reduction From each r o w subtract the smallest cost.
Then subtract the smallest cost in each column.
Accumulate the e ect of row and column reductions into the base cost.
1.2 : Eliminate redundant constraint Identify the row/s and/or column/s with the most zeros in the reduced cost matrix.
Eliminate the corresponding constraint/s (in case of an AP).
In this step we reduce the cost in the matrix cells by a base amount.This helps identify lowest cost cells with zeros for variables to enter into basis.Dropping the row with zeros facilitates a good start in the initialization phase, i.e., the next step.
Step 2 : Initialization phase 2.1 : Set up the simplex tableau Use a row for each constraint and a column for each v ariable.Enter reduced C ij 's in the cost row with a negative e n try for the base cost under the RHS.

: Identify BVs
For each unit-vector column, label the row containing the "one" with the name of the variable for the column.Label the remaining rows as open rows.

: Delete BV columns
This phase tries to set up the initial tableau with as many BVs as possible in it without any iterations.These are the variables with only single +1 in its column.Since we drop the constraint with the most zeros, the phase gives many basic variables immediately.At the end of this phase, the tableau has entries of 1 and 0 only.
Step 3 : Push phase The purpose of this phase is to make the current solution feasible while maintaining the optimality condition.This phase is similar to the dual simplex in that it starts from a vertex which is infeasible but close to the optimal solution (a warm-start).This is in contrast to the usual dual simplex which often starts far from the optimal vertex.

Numerical Example 1
We select a 3x3 TP used by Shih 73] to illustrate the Push-and-Pull algorithm.Step 2 : Initialization Phase: The following tableau is obtained after steps 2.1 through 2.3.Step 4 : Pull phase: The RHS < 0, hence Pull phase is needed.Variable X 32 enters to replace variable X 21 .After GJP, the following optimal tableau is obtained with non-negative R H S .The optimal solution is: X 11 = 1 5 , X 22 = 8 0 , X 31 = 5 5 , X 32 = 2 0 , a n d X 13 = 40, with a total cost of $ 3320.

Comments on the Push-and-Pull Algorithm
i) The entries in the body of the tableau are 1, -1, and 0 only.ii) Cost row a l w ays has non-negative e n tries.
iii) RHS entries are all non-negative at optimality.iv) There are only m + n ; 1 basic variables in any solution.This is consistent with the theory that in a TP, no more than m + n ; 1 cells should be occupied.
v) At t h e e n d o f S t e p 3 w e h a ve an initial solution.The entries under non-basic variables would have a count of one more 1 than of -1.In fact, these entries are analogous to corner cells of a path of the SS method.When a NB variable comes into basis, we add to the BV cell with -1 and subtract from the cell with 1.For example in Table 3, to add a unit to NB variable X 12 (an unoccupied cell), we h a ve to take a way a unit from X 11 and X 32 each w i t h 1 e n try and add a unit to variable X 31 with -1 entry.This indeed is a closed SS path from X 12 in this solution.
vi) The entries in the cost row for non-basic variables represent the amount by which the current total cost would increase if a unit of that variable is to be brought i n to the current basis.
vii) The algorithm avoids the use of the extra variables (slack a n d surplus) that are needed for the equality constraints, each of which must be converted to two inequalities, in the dual simplex 11] .The Push-and-Pull algorithm is also arti cial-free, as compared with the simplex method.This reduces computational complexity considerably.

Alternate Solutions
The proposed algorithm provides a clear indication of the presence of alternate, optimal solutions upon termination.Clearly, di erent alternate solutions give the same cost.The decision maker, however, has the option of deciding which optimal solution to implement on the basis of all the other factors involved.
Note that generating all alternative solutions can be computationally burdensome and may not be necessary.Moreover, the occurrence of alternative solutions is rare and can easily be avoided by a slight perturbation of the cost parameters.A decision maker, however, can easily generate alternate optimal solutions, if necessary.Given a nal solution, check the entries in the cost row to nd if there are any alternate solutions.An entry of 0 in the cost row represents a non-basic variable which can come into the basis to provide an alternate solution.Discard any that result in the same optimal strategy, that is, the basis changes but optimal shipment routes do not change.This can happen if the nal solution is primal degenerate.From the distinct optimal strategies, pick one to implement.Failure to identify alternate routes deprives a manager of the exibility regarding decisions to reallocate capacity from one route to another while maintaining the same cost.The optimal solution for example 1 given in Table 3 indicates that there are no alternate solutions in this case.

Properties of The Push-and-Pull Algorithm
The Push-and-Pull algorithm operates in the space of the original variables and has a geometric interpretation of its strategic process.We p r o vide this geometric interpretation of the algorithm by comparing it to the geometry behind the conventional simplex method.
The simplex method is a vertex-searching method.It starts at the origin, which is far away from the optimal solution for the TP and the AP.I t t h e n m o ves along the intersection of the boundary hyper-planes of the constraints, hopping from one vertex to neighboring vertex until an optimal vertex is reached in two phases.It requires adding arti cial variables since it lacks feasibility a t t h e o r i g i n .In the rst phase, starting at the origin, the simplex hops from one vertex to the next vertex to reach a feasible one.Upon reaching a feasible vertex, that is, upon removal of all arti cial variables from the basis, the simplex moves along the edge of the feasible region to reach an optimal vertex, improving the objective value in the process.Hence the rst phase of simplex tries to reach feasibility, and the second phase strives for optimality.The simplex works in the space of m n+ ( m+n;1) dimensions, because there are m n X ij s a n d m + n ; 1 arti cial variables, where m is the number of supply points and n is the number of demand points for the TP or AP.
In contrast, the Push-and-Pull algorithm strives to create a full BVS, that is, the intersection of m + n ; 1 constraint h yper-planes of the TP or AP that provides a v ertex.The initialization phase provides the starting segment of a few intersecting hyper-planes and yields an initial BVS with some open rows.The algorithmic strategic process is to arrive at the feasible part of the boundary of the feasible region, initial BVS (never empty for TP or AP), which m a y c o n tain an optimal vertex or a vertex that is close to it.In the Push phase the algorithm pushes toward an optimal vertex, unlike the simplex, which o n l y s t r i v es for a feasible vertex.Occupying an open row means arriving on the face of the hyper-plane of that constraint.Each successive iteration in the Push phase augments the BVS by including another hyper-plane in the current i n tersection.By restricting incoming variables to open rows only, this phase ensures movement in the space of intersection of hyper-planes selected in the initialization phase only until we hit another hyper-plane.Recall that no replacement of variables is done in this phase.At each iteration we reduce the dimensionality of the working region until we ll up the BVS, indicating a v ertex.This phase is free from pivotal degeneracy.The selection of an incoming variable as the one having the smallest C ij helps push toward an optimal vertex.As a result, the next phase starts with a vertex in the neighborhood of an optimal vertex.
At the end of the Push phase the BVS is complete, indicating a vertex.If feasible, this is an optimal solution.If this basic solution is not feasible, it indicates that we h a ve pushed too far.Note that, in contrast to the rst phase of the simplex, this infeasible vertex is to the other side of the optimal vertex.Like the dual simplex, now the Pull phase moves from vertex to vertex to retrieve feasibility while maintaining optimality, and it is free from pivotal degeneracy since it removes any negative (not zero) RHS elements.The space of the Push-and-Pull algorithm is m + n ; 1 dimensions in the Push phase and m n dimensions in the Pull phase.Note that m+n;1 i s t h e n umber of constraints and m n is the numberofvariables.
The importance of the Push-and-Pull algorithm is recognized by the fact that whereas 95 percent of the pivots in the simplex method are degenerate, which m a y cause cycling 62], in solving AP 28], the proposed algorithm is completely free of pivotal degeneracy.
Theorem 1: The Push-and-Pull algorithm is free from pivotal degeneracy which may cause cycling.
Theorem 2: The Push-and-Pull algorithm terminates successfully in a nite number of iterations.
The proofs of Theorem 1 and Theorem 2 are provided in the appendix.

Numerical Example 2
The optimal solution of AP is X 112 = X 211 = X 222 = X 121 = 0:5, and all other X ijk = 0.This implies that Adams should do job 1 on machine 2 for 50 percent of the time and job 2 o n machine 1 for the other 50 percent o f the time.Brown, likewise, should divide his time equally betwe e n j o b 1 o n m a c hine 1 and job 2 o n m a c hine 2. So each person spends 50 percent o f his time on each job using both machines at the total cost of 9. Note that this problem does not have a totally unimodular matrix to ensure an integer solution.Solving this problem for an integer solution using Gomory's cutting planes 6] of LP yields solution X 111 = X 222 = 1 , with inferior optimal value of 12.
The theoretical basis for the Push-and-Pull algorithm rests largely upon the total unimodularity of the coe cient matrix in the simplex tableau for TP, that is, the values of the coe cients are 0, -1, or 1 (see the appendix for a detailed proof).The nominal problem always has a unimodular constraint matrix.However, as observed in numerical example 2, this may not be the case when side-constraints are present.The violation of the total unimodularity does not a ect the solution procedure.This is a signi cant di erence between existing methods dealing with side-constraints and our approach.If the solution is not integral, one may i n troduce cutting planes to ensure an integral solution, if so desired.

Computational Behaviour and Implementation
As mentioned earlier, there are well over 40 solution algorithms for the TP and AP.Many di erent algorithms are presented and coded in various sources (see, e.g., 19], 24], 4 0 ] , 51]).The network simplex algorithm is consistently superior to the ordinary simplex and out-of-kilter algorithms.Most of this computational testing has been done on random problems generated by t h e w ell-known computer program, Netgen.To perform a meaningful comprehensive computational comparison, one needs to code all existing algorithms using similar data structure and processing.The computerized version of the algorithm must utilize the same compiler and the same platform.Then one can perform a thorough computational comparison either in terms of CPU, storage, or both.This is a major undertaking.Because of a lack of manpower resources at the present time, we are not able to present t h i s t ype of comprehensive computational comparison.
The Push-and-Pull algorithm is pivotal like simplex-based algorithms and is more e cient t h a n a n y other pivotal algorithm, such as the dual simplex and many other variants of simplex.A preliminary experiment to study computational behavior of the algorithm supports this assertion.Table 7 presents results comparing the proposed algorithm to the closely-related pivotal solution algorithm, namely the simplex algorithm, via the widely-used package Lindo on the VAX VMS 5.5-2.This computational comparison includes several small-to medium-sized published problems, which include some of the largest speci c (non-random) problems of which we are aware.Note that Lindo is equipped with many helpful features such a s a n tidegeneracy devices, steepest-edge selection 38], 41], 84], crashing start, and so on.In spite of all these features in Lindo, the Push-and-Pull algorithm is more e cient in terms of number of iterations.As demonstrated by the results in Table 7, the proposed algorithm reduces the numberby almost 50 percent as compared to Lindo.

Computer Implementation Issue
Up to now, the description we h a ve g i v en of the Push-and-Pull algorithm aimed at clarifying the underlying strategies and concepts.A second important aspect that we consider now is e cient computer implementation and data structures similar to those already developed by computer scientists for network-based solution algorithms 40].Practical TP and AP may h a ve thousands of constraints and even more variables.These problems have a large percentage of zero elements in the cost matrix (i.e., a sparse matrix).The proposed algorithm as described is not e cient for computer solution of large-scale problems since it updates all the elements of the tableau at each iteration.Some useful modi cations of the proposed algorithm as well as its e cient implementation will reduce the numberofelementary operations and the amount of computer memory needed.This section describes how an ecient implementation of the Push-and-Pull algorithm capable of solving large-scale problems could be developed.Adaptation of sparse technology is motivated by the fact that, for a large sparse problem, an iteration of the revised simplex method takes less time than an iteration of the standard simplex GJP.Our implementation is a straightforward adaptation of the proposed algorithm through appropriate use of data structure, sparse technology, pivot strategy rules, and triangularity of the coe cient m a t r i x .Computer implementation of all phases of the algorithm is discussed.We are concerned with computational e ciency, storage, accuracy, and ease of data entry.

Sparsity and Data Structure
The density of a matrix is the ratio of the number of non-zero entries in the matrix to the total number of entries in the matrix.A matrix that h a s a l o w d e n s i t y is said to be sparse.The sparse matrix technology is based on the following principles: -only non-zero elements are stored, -only those computations that lead to changes are performed, -the number of ll-in (i.e., non-zero) elements is kept small, -any u n used locations (zero elements) can be used for storing of ll-ins (non-zeros produced later).
Sparse technology 56] enables us to reduce the total amount o f w ork that is done in each iteration.Many schemes are available to avoid storing the zero entries of matrices.Since our calculations operate on tableaux by columns, we discuss a storage structure, called a linked list, that makes it easy to access data by columns.A piece of computer memory must contain three entries: one to hold an entry from the tableau, one to hold a row n umber, and one for the pointer to hold the next memory address.For example, the initial tableau of our numerical example 1 and its linked list for holding the data is as follows: Table 8.Initial tableau of the numerical example 1 RHS 1 1 0 0 0 0 0 55 0 0 1 1 1 0 0 80 0 0 0 0 0 1 1 75 1 0 0 1 0 1 0 100 0 1 0 0 1 0 1 40 Cost 25 0 2 0 5 10 1 -3120 The empty spaces in Table 9 are unused locations.By having a matrix generator program in order to generate constraints and the objective function automatically at our disposal, we can store the non-zero values in a linked-list form.This scheme has minimal storage requirements and has proven to be very convenient for several important operations, such a s p e r m utation, addition, and multiplication of sparse matrices, in our solution algorithm and PA.For example, it can be shown that if b is a 1 by m] matrix and A is an m by n] matrix that is stored as a linked list, then bA can be evaluated in d*m*n multiplications, where d is the density o f A. An advanced variant o f t h i s scheme is Sherman's compression, which is useful in storing triangularized matrices.Such an approach can be applied in all phases of the Push-and-Pull solution algorithm and the needed PA, since the constrained coe cients matrix can be triangularized.

Matrix Triangularization
Working with the triangularization of matrix A rather than A itself has been shown to reduce the growth of new non-zero elements signi cantly.It is well known that matrix A is reducible to a unique triangular structure.The very structure of matrix A a l l o ws for rearranging it in a lower triangular form.The ll-in is zero when A is triangularized.The general idea is to permute the rows and columns of the basis to make the rearranged matrix lower triangular.There are two phases as follows: Front phase: Among the unassigned rows, nd the one with the minimum count of non-zeros.If this count is 1, assign the row and its (only) unassigned column to the front and repeat otherwise, exit.
Rear phase: Among the unassigned columns, nd the one with the minimum count of non-zeros.If this count is 1, assign the column and its (only) unassigned row to the rear and repeat otherwise, exit.Table 10 shows a given matrix at a given iteration: 1 Upon entering the Front phase, we nd that row 2 is a singleton, so it is assigned to the front with column 1.Then, row 3 is assigned to the (new) front with column 3. Next the Rear phase is entered, and the singleton column 4 is assigned to the rear with row 4 .Finally, the singleton column 2 is assigned to the rear with row 1 .
The result is shown in Table 11

Push phase
In the Push phase one can construct a lower (block) triangular matrix by permuting matrix A with interchanging rows and columns.In order to preserve sparsity, an alternative approach for this reordering is block triangular diagonalization.This approach is attractive since under GJP, i t m a i n tains sparsity.In other words, this triangular structure is invariant under row or column interchanges to achieve sparsity within each block.
We were able to develop a Push phase employing the strategy of avoiding the use of arti cial variables with large positive coe cients 15], 18], 2 6 ] , 2 7 ] .In this phase, we start the initial tableau with a partially lled basic variable set.Therefore, the basic variable set is being developed.This particular crashing technique pushes toward a "vertex close to the optimal."Most commercial codes incorporate some other form of crashing techniques to provide the initial basis.

Pull phase
Using the matrix A, any iteration in this phase consists of solving A]x B = b to determine the direction of movement (say, k) and updating the pivot column y k by solving A]y k = C k .This basic solution is used in the Pull phase to determine which current basic variable must leave the basic variable set (if needed) and to move t o the new basic solution by performing the pivot operation (i.e., updating the A matrix).In the computer implementation, the inverse of the basis matrix is not needed 81].If A is a triangular matrix, then the foregoing systems of equations can be solved easily by one forward and two backward substitutions, respectively 60].
We can also achieve e ciency in performing GJ operations in this phase from another point of view known as A = LU factorization of A (see, e.g., 69], 75]).Matrices L and U are lower and upper triangular matrices for which t h e inverses can be computed easily and then solved for the above systems of equations by i nversion.In this approach L ;1 is stored as a sequence of "elementary" matrices and U is stored as a permuted upper triangular matrix.
In 33], 79], 80] some ingenious variant of the Bartles-Golub triangular factored updating procedure to exploit sparsity is proposed.This approach c o u l d lead to a decrease in accuracy through the introduction of round-o error.Working with integers (not real numbers), however, and not using any large numbers such as big-M, the computations would be precise and reliable.Since the A] elements are 0, 1 or -1, the column ratio (C/R) division operations can be replaced by a simple comparison operation.
The steepest-edge pivoting rule, also referred to as the largest-decrease rule, chooses the candidate whose entrance into the basis brings about the largest decrease in the objective function.This usually reduces the number of iterations, but may increase total computing time.A large fraction of the time needed to perform an iteration can be spent in retrieving the data.However, this approach i s best suited to the Push-and-Pull algorithm, given the sparsity o f TP and AP models after the row-column reduction.
Our discussion of computer implementation is by necessity and in no way is comprehensive.There are many interesting variants of what we mentioned that deserve further study to reduce the total number of iterations (the sum of all the iterations in all phases) and the amount o f w ork that is done in each iteration.

Managing The Cost Uncertainties
The cost vector C, is composed of elements C ij appearing row-wise, i.e., C = fC 11 , C 12 , . . . ,C 1n , C 21 , C 22 , . . . ,C 2n , C m1 , C m2 , . . . ,C mn g.The cost parameters are more or less uncertain.We are likely to be unsure of their current v alues and even more uncertain about their future values at the time when the solution is to be implemented.A time lag between the development and application of the model could create such uncertainty.In a volatile global market, an ever-changing currency exchange rate could a ect the cost of stationing executives in various countries in an AP.Therefore, managers are concerned with the stability o f the optimal solution under uncertainty of the estimated cost parameters.
Uncertainty is the prime reason why P A is helpful in making decisions.We can use PA t o g i v e u s i n f o r m a t i o n l i k e: the robustness of an optimal solution critical values, thresholds, and break-even values where the optimal strategy changes sensitivity of important v ariables sub-optimal solutions exible recommendations the values of simple and complex decision strategies and the "riskiness" of a strategy or scenario.For a more detailed discussion of PA in practice, see 64].

Current approaches to deal with cost uncertainties include:
Scenario Analysis -In this approach one assumes scenarios (a combination of possible values of uncertain costs) and solves the problem for each.By solving the problem repeatedly for di erent scenarios and studying the solutions obtained, a manager observes sensitivity and heuristically decides on a subjective a p p r o ximation.
Worst-Case Analysis -This technique attempts to account for putting safety margins into the problem in the planning stage.
Monte-Carlo Approach -In this approach, stochastic models assume that the uncertainty i n c o s t i s k n o wn by its distribution (see, e.g., 54]).
Reoptimization -The combinatorial and the network-simplex approaches in networkbased SA can only handle integer changes in C ij .The combinatorial approach requires the solution of a completely di erent maximum ow network problem for each unit change in any cost coe cient C ij , until infeasibility is reached.The network-simplex approach can handle any number of unit changes 3].Neither of these two approaches produce any managerially-useful prescriptive ranges for sensitivity analysis.A prerequisite for these approaches is to have a n a n ticipated direction of change.From a manager's point o f v i e w , a n ticipation of possible scenarios may not be an easy task.Application of a transformation scheme to attain the integrality condition for the nominal problem in the network-based algorithm makes PA too complicated to interpret.
PA deals with a collection of questions related to so-called "what-if" problems in preservation of the current optimal strategy generated by the proposed solution algorithm.PA starts as soon as one obtains the optimal solution to a given nominal problem.There are strong motivations for a manager to perform PA for the parameters C ij to: -help determine the responsiveness of the solution to changes or errors in cost values, -adapt a model to a new environment with an adjustment in these parameters, -provide systematic guidelines for allocating scarce organizational resources to data collection and data re nement activities by using the sensitivity information, -determine a cost perturbed region in which the current strategic decision is still valid.Although all aspects of PA are readily available for LP models, even OSA is rarely performed on TP or AP.While the SS method, the Hungarian method, or other traditional network-based solution algorithms may be e cient i n s o l v i n g a T P a n d AP, they are not designed to perform PA.The nal solutions generated by t h e S S method for TP or the Hungarian method for AP do not contain enough information to perform PA.If the needed sensitivity information were readily available in these algorithms, the operations research and management science textbooks would have c o vered the SA of the TP and AP.Considerable additional work is involved in obtaining the needed information.Moreover, one must be careful in using the existing LP computer packages for performing SA for AP as the optimal solutions are degenerate.
The basis inverse matrix is an essential prerequisite to SA.The SS and Hungarian methods do not provide this matrix directly.One may suggest using the optimal solution obtained by these methods, or any other traditional algorithm, to construct the basis matrix.This basis matrix can then be inverted to obtain the basis inverse matrix and other needed information.However, in addition to these extra computational costs, there is the problem of recognizing the basic variables contained in the optimal solution since the optimal solution may be degenerate.This is always the case for the AP since the number of positive v ariables is much smaller than the size of the BVS, due to a degenerate optimal solution.Also, note that it is this degeneracy of the optimal solution which causes the regular simplex to provide very misleading SA for the AP (and probably the TP).
Two recent d e v elopments in network-based SA, the combinatorial approach and the network approach (including the SS and the Hungarian methods), can handle only integer changes 3].The combinatorial approach requires the solution of a completely di erent problem, until infeasibility i s r e a c hed. Speci cally, it requires a maximum ow problem for each unit change in a cost coe cient.More complex changes such as simultaneous SA must be dealt with a sequence of this simple change.The network algorithm approach to SA can handle any number of unit changes.The change, however, is restricted to one change at a time as in the combinatorial approach.These limitations are caused by the inability of these solution algorithms (SS, Hungarian, network-based, and others) to handle a larger scope of SA.The network community h a s put much e ort on developing e cient solution algorithms at the expense of SA.
Neither of these approaches produces any managerially useful prescriptive ranges on costs for the SA.Additionally, a prerequisite for some of these approaches is to have a n a n ticipated direction of changes.From a managerial point of view, anticipation of possible scenarios may not be an easy task.
De ne the perturbed TP to be Minimize X X (C ij + C 0 ij )X ij (4) with the same constraints as before and the admissibility condition C 0 ij ;C ij .
The PA starts as soon as we obtain the nal tableau of the Push-and-Pull algorithm.At this stage we need to introduce some more notations as shown in Table 12.The initial tableau is partitioned into B, the BV coe cients, and N, the NB coe cients, as they appear in the nal tableau, and the RHS column.Note that, as it progresses, the solution algorithm removes the nonbasic columns.
Table 12.Necessary components for cost perturbation analysis from the nal tableau As mentioned earlier, the occurrence of alternative solutions is rare and can easily be avoided by a slight perturbation of the cost parameters.Note, however, that if a decision maker wants alternate optimal solutions and related SA, he needs to generate many ( a c o m binatorial problem) optimal solutions to get basis matrices, B, if using any other solution methods.Then each B h a s t o b e c o n verted and multiplied by the N matrix.However, A] = B ;1 N is provided by the Push-and-Pull solution algorithm as a by-product.We can easily generate other distinct A]'s by using this one to generate all other solutions, if needed.
Having obtained the necessary information from the nal tableau, cost PA can be started.To nd the critical region, that is, the largest set of perturbed costs for which the current solution remains optimal, we m ust nd the allowable changes in the cost coe cients.The necessary and su cient condition to maintain the current optimal strategy is that C N 0, where 0 stands for a zero vector with the appropriate dimension.Let denote the set of perturbed costs C 0 ij to maintain optimality: = fC 0 ij jC N 0 and C 0 ij ; C ij g: The set is non-empty since it contains the origin C 0 ij = 0, for all i, j.It is convex with linear boundary functions, by virtue of GJP operations.This set can be used to check whether the given perturbed (dependent or independent) cost coe cients have the same optimal solution as the nominal problem.To help clarify what we have done up to now, consider the numerical example 1.From the nal tableau in table 3, we h a ve: A] = 2 6 6 6 6 4 The perturbed cost vector is: The set can be used to check whether a speci c scenario has the same optimal basis as the nominal problem.For example, assume that the perturbed cost vector for a given scenario is: C + C 0 = f9 22 15 13 16 15 17 28 26g with C 0 = f4 ;8 3 ;7 ;2 ;15 2 3 3g.By substituting the values of C 0 ij , one can easily check that the current solution is still optimal.
For our numerical example 2 of AP, the critical region is:

Special Cases of Sensitivity Analysis
The PA construction given in the previous section is the most general one that can handle the simultaneous and independent c hanges in the cost parameters.In this section we derive some special cases of the sensitivity analysis in a direct method.
The PA set would generate the same results.

Parametric Perturbation Analysis
Parametric sensitivity analysis is of particular interest whenever there is dependency among the cost parameters.This analysis can be considered as simultaneous changes in a given direction.De ne a perturbation vector P specifying a perturbed direction of the cost coe cients.Introducing a scalar parameter 0 ( t h us, perturbation analysis,) we w ould like to nd out how far we can move in the direction of P, being the step size, while still maintaining optimality of the current assignment.
Step 1 : Identify P N and P B , s u b -v ectors of P corresponding to the NB and BVS, respectively, as they appear in the nal tableau of the nominal problem.
Step 6 : New C N = Old C N + (P N ; P B A]) for any 2 0 0 ].

Ordinary Sensitivity Analysis
From a managerial point o f v i e w , a n ticipation of possible scenarios or directions of change may not be possible.In this sub-section, we nd the range for any particular arc cost, holding all other costs unchanged.Ordinary sensitivity analysis, OSA, is very popular and readily available in LP.One cannot, however, use existing LP computer packages for performing SA for these problems when the optimal solution is degenerate.The other references present in current literature require signi cant additional computation beyond the solution algorithm.
OSA is a special case of parametric analysis where we w ould like t o n d o u t h o w far we can move in the direction of any one of the axes in the parametric space C 0 ij .Here, P is a unit row v ector or its negative, depending on whether we w ant to compute an upper or lower limit.The step size is the amount of increase or decrease in that direction.Alternately, note that we can also nd the allowable changes for any particular cost C 0 ij by setting all other costs to zero in the set , equation (6).The results are as follows: Lower Limit Upper Limit Similarly, for numerical example 2 for AP, these limits are obtained from the corresponding set , equation (7), by letting all C 0 ijk = 0, except the one for which we are calculating the bounds.

The 100% Rule
The above analysis is for one-change-at-a-time. Suppose we w ant t o n d t h e s i m ultaneous allowable increases in all cost coe cients.Bradley et al 17] discuss the use of the 100 percent rule for simultaneous increase or decrease of all costs.This rule is based on the ordinary sensitivity limits.The 100 percent rule says that optimal basis will be preserved if where the sum is over all i and j and the denominators (C ij ) are the allowable increases from the ordinary sensitivity analysis.That is, as long as the sum of all of the percentages based on the allowable changes for each cost coe cient is less than 100 percent, the current optimal basis remains unchanged.For the above example, the current routings are optimal as long as cost increases are such t h a t C 0 11 =15 + C 0 13 + C 0 22 =12 + C 0 31 + C 0 32 =15 1.This condition is su cient and not necessary.Similarly, the application of the 100 percent rule when decreasing all cost coe cients provides P P C 0 ij = ;C ij 1, where the sum is over all i and j and the denominators are the allowable decreases from the ordinary sensitivity analysis with a similar interpretation.Note that the upper limit and lower limit must be rounded down and up respectively.

Models with Side-Constraints
In real-life situations, it is common for a few side-constraints to evolve during the time period of development and implementation of an optimal solution 2], 46], 71], 73].For example, there could be a limit on the shipment along a particular route or combination of routes.Alternately, there could be a constraint on one route in relation to other routes.A general side-constraint is to determine the amount of shipment from source i to destination j under the conditions that L ij P a ij X ij U ij where X ij denotes the number of units shipped from i to j and L ij , U ij and a ij are constants.Without loss of generality, we assume the L ij to be zero.If an L kt > 0, then a change in variable X kt reduces the lower bound to zero.We provide a method to accommodate these constraints through the Push-and-Pull algorithm.

Method Capacitated
Step 1 : Ignore the upper bounds.Solve this nominal problem by Push-and-Pull algorithm.
Step 2 : Derive the nal tableau of the nominal TP.
Step 3 : Check if all conditions P a ij X ij U ij are satis ed.If yes, go to step 7.
Step 4 : Pick the basic variable X kt for which U ij ; P a ij X ij is the largest.
Add an open row with constraint P a ij X ij = U kt .
Step 5 : Pivot on X kt column so X kt remains basic.
By construction, RHS is infeasible with a negative e n try in the open row.
Step 6 : Perform the Pull phase of the Push-and-Pull algorithm.Go to step 3.
Step 7 : This is the optimal tableau.STOP.Obtain the solution.
Consider the example of Table 1 with constraint X 31 2X 32 .This constraint i s not satis ed in the optimal solution of the nominal TP, T able 3.After performing steps 4 and 5 of the Method Capacitated, we o b t a i n t h e f o l l o wing table: Note that due to the unrestricted nature of the given constraints, the non-zero entries of a tableau in Method Capacitated could deviate from the regular pattern of being 1 and -1 only.This deviation from unimodularity does not a ect the solution algorithm.After performing the Pull phase to bring X 33 into BVS, we obtain the following optimal tableau: Step 1: Perform the row and column cost reductions.Solve the TP by using Push-and-Pull algorithm.
Step 2: Obtain the nal tableau as presented in Table 3 in the previous section.
Step 4: Open a row with X 22 = 70, corresponding to the most violated constraint X 22 70.

Degenerate Optimal Solution
The e ect of degenerate optimal solution on SA has been studied by many research workers (see, e.g., 34], 53], 5 2 ] ) and is well known.However, to the best of our knowledge, no algorithm or LP package exists for cost SA in case of an optimal degenerate problem.The researchers have completely overlooked the degeneracy phenomenon, a common occurrence in a TP, in cost SA.Lindo and Netsolve g i v e wrong cost sensitivity results if the optimal solution of a TP turns out to be degenerate, which i s a l w ays the case with an AP.Note that the set given by equation ( 5 ) i s v alid for a non-degenerate problem only, that is, when RHS > 0. For a degenerate optimal solution, we modify the procedure for solving the perturbed problem (4) as follows: Step 1 : L e t Z = c o u n t of zeros in the RHS of the optimal tableau.Set I = Z+1, i = 0 .Step 2 : Let i = i+1.Determine i for the current optimal tableau.
Step 3 : If i = I, go to Step 8.
Step 5 : PC = pick the column with the largest negative row ratio, i.e., cost row/PR.Step 6 : Generate the next distinct nal tableau.
Note that as before, we c a n u s e for any scenario analysis.The set can also be used for parametric, or ordinary cost SA.Alternately, w e n d 0 for parametric PA, and lower and upper limits for OSA, for each distinct nal tableau and then use the corresponding intersection for analysis.We illustrate using the AP from Anderson et al. 5].

Numerical Example 3
Cost matrix of the AP: 10 15 9 9 18 5 6 14 3 We solve this AP using the Push-and-Pull algorithm.The following nal tableau is obtained.
There are no alternative optimal strategies.As expected, however, the optimal tableau is degenerate.The other optimal tableaux are as follows: Table 19.Second degenerate nal tableau BVS X 33 X 22 X 13 X 32 RHS X 31 1 Table 20.Third degenerate nal tableau BVS X 21 X 22 X 13 X 11 RHS X 31 1 The critical regions, i 's, based on these nal tableaux are: The 100% rule limits for C 0 ij 0 and C 0 ij 0 can be obtained directly from these above limits.Similarly, to carry out parametric analysis given any perturbation vector P, w e rst nd 0 1 , 0 2 , 0 3 .The current basis are optimal for any step size 2 0 0 ] where 0 = m i n f 0 1 0 2 0 3 g.

Concluding Remarks
We have proposed a general-purpose uni ed algorithm to solve the classical TP and AP.The algorithm provides one treatment for both problems, in place of the SS and Hungarian methods.The proposed algorithm solves these problems with a w arm-start and uses Gauss-Jordan pivoting only.It is computationally practical in the sense that it does not require any slack/surplus variables (as in simplex) or any arti cial variables (as in dual simplex).In addition, the Push-and-Pull algorithm and the proposed PA have t h e advantage of being easy for managers to understand and implement.The results empower the manager to assess and monitor various types of cost uncertainties encountered in real-life situations.The nal solution tableau provides the information required to handle any special cases of side-constraints.For such a case, it is not necessary to restart the algorithm.If the current solution does not satisfy some of the side-constraints, the method brings the most violated constraint i n to the nal tableau and uses the "catch-up" operations of the Pull phase to generate the updated nal tableau.This catch-up capability i s also useful in coping with structural changes, such as the deletion of a route that becomes inaccessible due to a storm, construction, or bad road conditions.Finally, the algorithm is free from pivotal degeneracy (cycling) and can solve AP of higher dimensions than the traditional 2-dimensional AP.
The proposed solution algorithm deals e ectively with change and chance and is compatible with other solution algorithms.The PA in this paper is e cient since it uses the same data and the same manipulation as the solution algorithm.All the information needed to carry out PA of cost coe cients is readily provided by the nal tableau.In contrast to OSA, our algorithm allows for simultaneous, independent, or dependent c hange of the cost coe cients from the estimated values.Unlike the existing network-based SA, the PA provides prescriptive uncertainty ranges.As long as the estimated parameters remain within ranges speci ed by the PA, the current optimal strategy remains optimal.In cases where the optimal solution is not unique, one must be careful in using any part of the PA results, because the results may not be correct for an alternate optimal solution.We demonstrate the applications of PA when an optimal solution is degenerate and provide modi ed steps 44].
From a managerial point o f v i e w , t h e P A results provide an assessment and analysis of the stability of the optimal strategy by monitoring the admissible range that preserves it.This permits evaluation of the impact of uncertainty in the data.PA g i v es a manager more leverage in allocating scarce resources.Monitoring the admissible ranges of PA that preserve t h e current optimal strategy aids in determining how resources should be applied.Information about departure from these limits allows a manager to anticipate the consequences and to pre-determine backup strategies.The results empower a manager to assess, analyze, monitor, and manage various types of cost uncertainties in all phases of the model, namely design, solution, and implementation.
The proposed approach can also provide a rich modeling environment for strategic management of complex transportation and assignment problems using Monte Carlo simulation experiments.In this treatment, one may select a p r i o r i cost coefcients for each route or assignment with the support domain produced by the PA.The stability of optimal shipments with respect to changes in supply and demand in TP is given in 1].
Computational results comparing the proposed algorithm to the closely-related pivotal solution algorithm, the simplex, via the widely-used package Lindo, indicate that the Push-and-Pull algorithm reduces the number of iterations by almost 50 percent as compared to the simplex.Future work should develop e cient codes to implement this approach for computational studies, design and experiment to compare the Push-and-Pull algorithm to existing algorithms.An immediate extension would apply our approach to a study of the more-for-less paradox, that is, given an optimal solution, is it possible to nd a less (or equivalent) cost solution by shipping more total goods under the restriction that the same amountofgoods are shipped from each origin and to each destination and all shipment costs are non-negative 70].Another possible direction for future research w ould be design of an intelligent computer-assisted system of various parts of PA for the decision maker.
through the Push and Pull phases does not contain any l o o p s , it su ces to show that each phase terminates successfully.
The Initialization phase uses the structure of the problem to ll-up the BVS as much as possible without requiring GJP iterations.The Push phase completes the BVS.The number of iterations is nite since the size of the BVS is nite.At the end of the Push phase we h a ve a basic solution that may not be feasible.The Pull phase terminates successfully by the well-known theory of dual simplex for these problems.

Call for Papers
Important advances in mathematics, physics, biology, and engineering science have shown the importance of the analysis of instabilities and strongly coupled dynamical behavior.New investigation tools enable us to better understand the dynamical behavior of more complex structures.However, the increasing interest in mechanical structures with extreme performances has propelled the scientific community toward the search for solution of hard problems exhibiting strong nonlinearities.As a consequence, there is an increasing demand both for nonlinear structural components and for advanced multidisciplinary and multiscale mathematical models and methods.In dealing with the phenomena involving a great number of coupled oscillators, the classical linear dynamic methods have to be replaced by new specific mathematical tools.
This special issue aims to assess the current state of nonlinear structural models in vibration analysis, to review and improve the already known methods for analysis of nonlinear and oscillating systems at a macroscopic scale, and to highlight also some of the new techniques which have been applied to complex structures.
We are looking for original high-quality research papers on topics of interest related to specific mathematical models and methods for nonlinear and strongly coupled (correlated) oscillating systems and for distributed-parameter structures that include but are not limited to the following main topics: • Vibration analysis of distributed-parameter and multibody systems, parametric models • Global methods, wavelet methods, and fractal analysis for spatially and temporally coupled oscillators • Nonlinear time series methods for dynamic systems • Control of nonlinear vibrations and bifurcations, control of chaos in vibrating systems.Transient chaos.Chaotic oscillators.Bifurcations • Micro-and nanovibrating structural systems Before submission authors should carefully read over the journal's Author Guidelines, which are located at http://www .hindawi.com/journals/mpe/guidelines.html.Prospective authors should submit an electronic copy of their complete manuscript through the journal Manuscript Tracking Sys-tem at http://mts.hindawi.com/according to the following timetable: The base cost is 5(55) + 18(80) + 15(75) + 7(40) = 3120 dollars.Delete the rst demand constraint w i t h t wo zeros from the initial simplex tableau.

Table 7 .
Number of Iterations using the Push-and-Pull algorithm versus Lindo Problem Type Simplex Algorithm Problem Source and

Table 1 .
Cost matrix of Shih's problem.

Table 5 .
Initial tableau using the Push-and-Pull algorithm BVS X 111 X 112 X 121 X 122 X 211 X 212 X 221 X 222 RHS

Table 9 .
The linked list for the coe cient matrix of example 1

Table 13 .
Tableau after adding an open row

Table 14 .
Optimal tableau of TP with side-constraints Numerical Example for Fixed Upper Bounds in TPShih 73]discusses the special case of L ij X ij U ij where X ij are bound by constants only and o ers a solution by modifying the SS method.Though conceptually the modi cation Shih o ers is simple, creating new tables and nding new alternative paths repeatedly, turns out to be cumbersome.We illustrate the steps of Method Capacitated further by walking through the example of Table1with xed upper bound constraints X 22 70, X 31 50, X 21 8.

Table 15 .
Modi ed tableau after adjusting column X 22 BVS X 12 X 21 X 23 X 33

Table 17 .
Optimal tableau with xed upper bounds