American Journal of Computational and Applied Mathematics

p-ISSN: 2165-8935    e-ISSN: 2165-8943

2014;  4(3): 97-110

doi:10.5923/j.ajcam.20140403.05

Solving System of Linear Inequalities or Equalities on the Surface of the Unit Shell (LIS-III)

Paul T. R. Wang

WangPaul_Research

Correspondence to: Paul T. R. Wang, WangPaul_Research.

Email:

Copyright © 2014 Scientific & Academic Publishing. All Rights Reserved.

Abstract

The author presents an innovative concept of a generic algorithm that solves system of homogeneous linear inequalities for finite number of unknowns and constraints utilizing normalized unit vectors of length 1 on the surface of an n-dimensional hyper sphere with a radius of 1 coined as the unit shell and the concept of equal distanced points as sysmmetic points along the arc connecting two equal distanced points to a selected set of points on the unit shell with increasing ranks to locate the desired solution point or points of the given system of linear inequalities if such solution point or points exist. A direct application of this innovative technique applied to a linear program formulated as a system of self-dual homogenerous linear inequalities is illustrated to establish its validity. Furthermore such technique is also illustrated to extend its applicability to solve Differential Variation Inequalities (DVI) over the Unit-Shell. Such a new technique is shown to be extremely efficient, numerical stable, and quite suitable for very large system of linear inequalities with number of variables and/or constraints over millions. In addition, such a new approach does provide insight to solve linear inequalities that is compatible to Gaussian elimination solving linear equalities.

Keywords: Homogeneous Linear Inequalities, Unit-Shell, Liner Programs, Linear Projection and Distance to subspace, n-dimensional normed linear space, and differential variation inequalities

Cite this paper: Paul T. R. Wang, Solving System of Linear Inequalities or Equalities on the Surface of the Unit Shell (LIS-III), American Journal of Computational and Applied Mathematics , Vol. 4 No. 3, 2014, pp. 97-110. doi: 10.5923/j.ajcam.20140403.05.

1. Introduction

Many years back, the author was working on the hemispherical cover for unit vectors on n-dimensional unit shell as the surface of an unit sphere. If we normalize all row vectors of a linear inequalities system (LIS) in its feasibility form (), they are points on the unit shell of this n-dimensional hyper sphere. Note that has a solution implies that all the normalized row vectors of and the normalized solution, , has a dot product . Hence the hemispherical surface of the unit sphere has a center at with all points reachable within an arc distance of or less. In other words, we have a hemispherical cover covering all the points of . Such a solution of the LIS is also a point on this unit shell as the geometrical center (a point with its maximum arc distance to all the constraint points minimized) of this hemispherical cover. Note that LIS is generic for all linear systems equalities and inequalities are both included. In particular linear programs (LPs) is only a special case of LIS [ref. 17].
The pursuit to locate a unit vector for did not succeed then as the geometrical center of a fixed number of points (especially when we are dealing with millions of points!) on a unit shell is very troublesome to locate in n-dimensional hyper space. However, now an effective approach has been identified. Two new insights make such a pursuit feasible and verifiable. First, the concept of “equidistant” for a specific point as an equal distanced point (EDP) to a group of points on the unit shell can be maintained using orthogonal subspace projection. Second, equidistance can be maintained by moving an EDP on the n-dimensional great circle arc connecting two equidistant points with respect to a given set of points. By repeatedly applying this EDP tracking strategy to a set of linearly independent unit vectors with increasing rank as to increase the number of points that are equidistant to the center point as its geometrical center, the solution for the original LIS can easily be located with a minimum number of steps no more than the maximum possible rank of the defining matrix of the LIS. Such a technique is greedy, non-iterative, non-linear, and non-Gaussian, and use only the dot-product operator which is O(n) in operation. No basis vectors swapping is required to minimize the maximum arc distance from a given point to the geometrical center or centers of a set of points on the unit shell. Furthermore, most matrix operations carried out on the unit shell are over unit vectors; hence it is numerical stable avoiding ill-conditioned operations. Consequently, the unit shell approach can be used to resolve problems posted by any liner equalities or inequalities. As result of the perfectly symmetric nature of the unit shell, it can be demonstrated that it requires a minimum number of steps to reach a solution while maintaining numerical stability when it is compared to existing best known algorithms offered by the Simplex or interior point methods. The author illustrates this innovative new technique to solve a simple linear program (LP) system with 7 variables and 14 constraints. Excel is used to carry out such new approach step by step to illustrate and verify its correctness and efficiency. A twist of this approach is that it is capable of solving all the linear systems (equalities, Eigen systems, and inequalities alike) over the unit shell non-linearly! The normalization of vectors on the unit shell also offers numerical stability for ill-conditioned linear systems that are very troublesome to solve with existing techniques in dealing with very large linear systems with millions of variables or constraints. Such an approach may have other applications beyond linear systems such as physics as an alternative to the string theory to resolve the problem of renormalization and operations research NPC mystery (?).
The author presents such a geometrical center search strategy for the maximum possible rank with the same number of linearly independent rows in as a solution to . Such analysis shows that it is possible to solve LIS (linear inequalities system) with non-iterative and non-Gaussian alike algorithms. The author is putting this draft together to highlight this unit shell approach. This approach is coined as LIS-III or Unit Shell for Linear Systems or (US_LS). Note that LIS-II is a generalized Gaussian elimination (GGE) procedure solving system of linear inequalities. LIS-II is capable of determining the feasible interval of individual variable as a linear inequalities feasibility analyzer (LIFA). LIS-I is a set of three key algorithms that recursively reduce maximum infeasibility, sum of all infeasibility, and the number of constraints with the worst case infeasibility for any given linear system in feasibility form. US_LS has many advantages over existing interior point methods or the ellipsoid method as US_LS is symmetric and numerical stable for all the vectors are normalized to be unit vectors on the surface of a unit shell. It is a global, symmetric, stable, and one pass approach without the need of iterative application of some key algorithms.

2. Highlights

1. Any linear systems in feasibility form (LSF) [ref. 17]:
2. Linear Systems on the Unit Shell after normalization:
Note that each row, in U is a point on the surface of an n-dimensional unit sphere (with a radius of 1.0), namely the unit shell. The solution normalized solution, is also a point on the same unit shell. Moreover, the maximum arc distance from this point to all the points in U is minimized by the definition of a minimum cover on the unit shell for all points in U.
3. A great circle arc connecting two points on a unit shell
Let p and q are two points on the surface of an n-dimensional unit shell, then with with is a path on the unit shell connecting points and as part of the great circle that passing both points and . In other words, we have the equation for all the points on the great circle arc connecting p and q:
Note that is simply the normalization factor to project the vector onto the unit shell,
Consequently, the mid-point of is it is also the geometrical center of on the unit shell.
4. Equidistant Point (EDP) point of a set of points, on the surface of the unit shell is any point, p, on the unit shell such that for all i.
5. A Cover and the Geometrical Center of points on the unit shell
The disc centered at e with a radius of R is a cover for all the points with i.e., or simply.
Given a set of points, on the surface of the unit shell, a point, c, is a geometric center, GC, of if the maximum arc distance from c to all the points in is minimized. In other words,
Given a set of points, on the surface of the unit shell, let rank
Let be an EDP for on the surface of the unit shell, let. We have by definition of EDP; we also have for all i
Hence, we should have the solution for the coefficients from
If , from
then if is an interior, boundary, or exterior point of the cover, .
Swapping with we have .
Note: let and , if is outside the minimum disc cover of , the cover as a disc centered at an EDP, of (as the geometrical center of , namely, with a radius as does not cover , it may still be possible to find a disc cover centered at with a radius as an EDP for .
If rank = 3, then the normalized orthogonal projection of onto ,
i.e.,
is an EDP () of with radius.
The great circle that goes thru both and are all EDPs of . Note that an EDP, for is the normalized such that .
By induction, the same argument shows that exterior points of a disc cover on the unit shell may be covered by another disc as a boundary or interior point of another EDP at a higher rank. The minimum cover is one centered at an EDP with least number of boundary points and hence, smallest radius with the smallest rank. In summary, by moving to a higher rank, one can increase the number of boundary points centered at distinct EDP or EDPs with a larger radius. How to find a second EDP with increased rank (outside the subspace on the unit shell)?
Two methods are investigated:
(a) Using Liner Subspace Projection Operation, as described above.
(b) Using natural basis and subspace ranking
(a): Assume that we have set of unit vectors with , and another unit vector, , that is an EDP for . The following steps are needed to locate another for .
1. Compute the maximum arc distance from to any point in as
with .
2. Compute subspace projection of onto as
3. Identify the orthogonal component of to as
4. Normalize x to become an unit vector as
5. Move point for any along the great circle on the unit shell towards from y towards such that , i.e., we have
for
Alternative, for for
(b): Hint:
Let , where is the unit vector with 1 at its i-th position and 0 otherwise.
Note that in the above equality, and are n-dimensional vectors while is a scalar real number, in other words,
We have,
Solve
for
dimensional vector with its i-th component,
Namely, we treat as a row vector .
Note that has full rank n , if k<n, we can find at least one
from , such that we have
For to be an EDP distinct from for , we must have
for j=1,2,…..,k. Namely, ,
i.e.,
In other words, any unit vector as normalized with obtained from
is an EDP for
6. Four cases of relationship of with .
CASE I. Rank Rank and minimum covers.
we have also a minimum cover of .
CASE II. Rank Rank and minimum does not cover .
CASE III. Rank Rank and minimum covers .
It is possible to locate another EDP, , for such that covers
with . is obtained from the great circle arc connecting
where is the normalized for any full ranked .
CASE VI. Rank Rank and minimum does not cover.
Analysis: For CASE II,
Let , we have ; hence .
Move a point , from towards along the great circle arc , until we have
for at least one . Let and,
Let as the midpoint of , select from such that .
For i=2,3 ----, m we can select linearly independent subset
Such that we can find at least one such that , i.e.,
Let be any full ranked I by I sub-matrix of , we have
Solve for Normalized is an EDP of
and the great circle arc connecting and are all EDP for .
is selected from such that. When we have a new cover
for. For CASE IV, It is possible to locate another EDP for such that a hemispherical cover for is obtained. This is because is a hemispherical cover of
where
.
7. Arc connecting two EPDs on Unit Shell
If points p and q are distinct equidistant points (EDP) for on the surface of the unit shell, then, all the points on the great circle connecting p and q are all EDPs for .
Proof:
Let p and q be two distinct equidistant points for , we have
and .
Let be a normalized unit vector on the great circle connecting p and q. Then, it is clear that we have
Where and .
8. The US-LS algorithm, a unit shell algorithm to minimize the maximum infeasibility of a linear system in feasibility form, Excel is used to carry out this proposed algorithm step by step to solve a simple LP in LSF form with 7 variables and 14 constraints:
Step 1. Problem formulation: Let the LSF be normalized such that every row vector is an unit vector with it end point on the unit shell, we have the normalized solution of the LP also an unit vector on the unit shell and the problem becomes given an matrix with k unit vectors as
, locate a unit vector, such that we have .
Step 2. Select a pair of distinct points, from rows in such that compute the midpoint, , of as it is described in section 3.0. Note that is both an EDP (equal distanced point) and the geometrical center of .
Step 3. Compute the subspace projection of as .
Step 4. Find the first unit vector, from the unit vectors of such that the arc distance and that the normalized is not equal to . In other words, let such that and . Note that we have rank where
Step 5. Normalize vector to become an unit vector , then is also an unit vector with its end point on the surface of the unit shell. Note that both and are EDPs for points in . Consequently, the great circle are all EDP for . On this great circle, one can locate a unique point, as the shortest EPD for by the following two equations:
a. as
(note that we will always have
b. as is also an EDP for .
Solve for t, from a. and b., we have, note that we have , hence,
.
Substitute into equation a., we have the shortest EDP or the geometrical center, for .
Steps 2 to 5 can be repeated with increasing subscription for and for where m=rank (U).
Step 6. Once the highest rank for, m, for has reached, we have a geometrical center, of such that is an EDP for . For all rows (points) of , compute for maximum arc distance . If , then we have identified a minimum cover .
If , we have has rank m-1, the orthogonal projection normalized is an EDP for and has its arc distance to all points in as . In other words, we have .
If a hemispherical cover, for exists, it must be one of these m possible .
9. An Example: Consider the following sample LP pairs as primal and its dual [Refs 1 to 7].
In self-dual formulation as pure linear inequalities, we have:
Normalize (1) to have all row vectors as points, on the 7-dimensional unit shell, we have
Let , we have , , from we determine such that , we have , i.e.,
, the geometrical center of is the midpoint of , we have
Hence, .
is determined from such that , we have , i.e.,
, Note that is both EDP and the geometrical center of .
To locate the geometrical center which is also an EDP for, we can compute the orthogonal projection of of the subspace , i.e., is the normalized unit vector of where
. Note that this normalized unit vector
is orthogonal to points in and has an arc distance of to all the points in. In other words, we have identified another EDP for that is distinct from .
The great circle arc connecting and are all EDP for and the geometrical center for is also a point, , on this great circle arc as an EDP for with the minimum radius. There is only one such point as the normalized with t computed as shown in Step 5.
In this example, we have:
Using Excel, we have t=0.210480823 and
Hence, we have .
Continue this process for , and until we reach the maximum possible rank of 7 as:
with
Using Excel, we have t =0.460315637 and
We have
with
Using Excel, we have t=0.598264383 and
We have
with
Using Excel, we have t= 0.483043855 and
We have
From , we identified the unit vector orthogonal to from as
Using Excel, we have t= 0.468058988 and
Figure 1 illustrates the arc distance of and the arc distance of for where
. As for different value of t are all EDPs for , we have an EDP for at t= 0.468058988.
Figure 1. Computing EDP for as
Consequently, we have
With .
We have reached the geometrical center as an EDP for a full ranked as. The only possible hemispherical cover for is centered at the normalized and has solution if and only if such as hemispherical cover for exists.
It is trivial to check that the unit vector is the solution of for the original LP. In other words, we have identified the hemispherical cover for all the 14 points in with .
With Excel, we have:
To verify that is indeed the solution of the original LP in (1), we have
such that
Hence the solution to the LP in (1) is:
Figure 2. Solution of LP as the Center of a Hemispherical Cover, for
Figure 2 illustrates the arc distance from all points on U to the solution point, , such an EDP for as the geometrical center of . Note that such a solution point is a point with a minimum arc distance to all the constraint points. In other words, it has the minimum arc distance to all the constraint points. It is the best solution to any given set of linear inequalities that can be humanly resolved.
The solution obtained from the unit shell is identical to that obtained by either the simplex method, interior point method, LIS-I, or the GGE (LIS-II) approaches. Over the past few years, the author has developed three different methods for solving any linear systems and they are as efficient as the text book simplex and interior point methods with considerable simplification and numerical stability.
This paper details the theory and steps of such a generic algorithm to solve both linear equalities and inequalities on the Unit Shell with a minimum number of steps involving only the unit vectors and linear subspace projection operation . A simple LP solved with this new technique is demonstrated step by step with Excel tabulation and matrix operations. The support is attached for verification and validation of this new approach.
10. Solving Differential Variation Inequalities (DVI) on the Unit-Shell
We generalize the unit-shell algorithm over linear space to a normed Banach space with an inner-product operator with norm for solving variational inequalities (VI) or differential variational inequalities (DVI) as follows [refs 14 and 15]:
Let where B is an matrix as a nonempty convex compact polyhedron in
Let be a continuously differentiable function from into with Jacobian .
The variational inequality problem (VIP) associated with and is to locate a solution in satisfying the variational inequality (VI): in. Note that in , we have
Let the gap function associated with a VIP be defined for in as:
While the dual gap function associated with a VIP is defined as
Using Newton’s first order Taylor linear approximation around a point in , a linearized VIP as LVIP can be computed iteratively for as:
Consider the following nonconvex, nonlinear constrained mathematical program:
subject to
Note that optimality occurs at:
subject to:
Consequently, we have the following homogeneous linear feasible system of inequalities
Note that can be solved effectively over the unit shell using only linear projection and the concept of equi-distanced points to selected set of normalized unit vectors derived from row vectors of L as described and demonstrated in this draft paper from Section 1.0 to Section 8.0.
11. Other Examples of Linear Programs Solved by the Unit Shell Approach
First, Any non-homogeneous linear system with equalities, inequalities, or mixed as constraints may be converted to homogeneous linear system by increased its number of columns as variables by 1 shifting the right hand side vector to the left hand side so that the right hand side is always a vector with only zero coefficients. The row vectors of such a homogeneous linear system can always be normalized into vector of unit length as unit vector on the surface of a n-dimensional unit sphere as the unit shell. The normalized solution of such a homogeneous linear system is also a unit vector on the same unit-shell. In this section, the author provides both theoretical analysis and examples for four possible cases of solving homogeneous linear systems on the unit shell, namely, with unique solution, no solution, infinite solutions, and solutions that are unbounded. Both unbounded and infeasible LPs are examined and analyzed to provide necessary and sufficient conditions for linear systems with null or infinite many solutions such that the solvability of any system of linear equalities or inequalities are resolved with a computational complexity that is compatible to that of the traditional Gaussian Elimination for system of linear equalities, namely, where k is the rank of the linear system (equalities and inequalities included) with m constraints and n variables.
11.1. Necessary, Sufficient Conditions, and Examples for the Four Cases of Linear Systems with Unique Solution, No Solution, Infinite Many Solution, and Unbounded Linear Solution
Example in Section 8 clearly illustrated the fact that when the rank of the homogenous linear inequalities equals to the number of unknowns plus 1, and all the unit vectors associated with the homogenous linear inequalities are reachable (or covered) by the hemispherical cover of the unit shell centered at:
If all the unit vectors associated with a given homogenous linear system can not be covered by such a hemispherical cover, it does not have a solution and the conflicting inequalities can be easily identified. On the other hand, if the rank of hemispherical cover is less than the number of unknowns, we will have an infinite number of EDPs for all the constraint inequalities, hence, it will have infinite many solutions or as unbounded cases.

3. Conclusions and Future Work

In conclusion, the author presents an innovative approach that relates algebraic relation in vectors to geometrical relation as arc distance on the unit shell, i.e., the surface of a n-dimensional hyper sphere.
The solution or solutions to any given linear system (both equalities and inequalities included) is shown to be the geometrical center or centers that are points on the unit shell with its maximum arc distance to all the constraints point minimized. This technique is illustrated and verified with an example of solving a linear program (LP) with 7 variables and 14 constraints. Only the dot product and linear subspace projection operation are used to locate the geometrical center or centers with a minimum number of steps. As dot product of two unit vectors is a direct measurement of correlation between two sets of numbers, the geometrical center or centers do reveal maximum correlation with minimum discrepancy. When we can compute efficiently the geometrical centers of a set of constraint points, it is applicable to a very wide range of scientific and engineering applications in aviation, finance, economics, operations research (OR), data mining, signal processing, and pattern recognition … etc. Additional benefit of dealing with only unit vectors is its numerical stability free from ill-conditioned operations that other approaches encountered very frequently. A future plan is to implement the unit shell algorithm in java as a generic optimization tool for general purpose scientific and engineering use.
In addition, the author is working on recursive projection and distance functions for unit vectors onto its k-dimensional subspaces and its direct application in solving homogeneous linear inequalities. This approach does provide effective and numerical stable algorithms that can handle very large system of linear inequalities with a computational complexity compatible to the classic Gaussian for solving system of linear equalities. This approach is extremely handy in dealing with linear inequalities with millions of unknowns and/or constraints. A drafted paper is currently under peer review by both mathematician and engineers.
Figure 3 illustrates the possibility of using recursively projection and distance functions in trigonometry as Sin and Cos functions to replace the classic linear space projector, , and its orthogonal operators: .
In Figure 3, is a linear subspace with dimension and rank = k and an orthogonal basis.
Note that the computational complexity of locating such a geometry center for the set of normalized homogenous linear inequalities can be shown to be compatible to that of the classic Gaussian Elimination solving systems of linear equalities.
Support from NSF funding or private foundations source in innovative mathematical computation for further clarification and elaboration will be pursued.
Figure 3. Recursive Linear Subspace Projection and Distance Functions as Trigonometry Identities

ACKNOWLEDGEMENTS

The author acknowledges valuable suggestions, comments, and initial peer review of the concept of the Unit Shell by his colleagues, in particular, David Hammrick, Dr. Leone C. Monticone, and Dr. William P. Niedringhaus during past 10 years at the MITRE Corporation. Latest discussion and comments from Oliver C. Wang and Brian M. Jones have contributed to the clarity of the role of linear projection and orthogonal distance functions relevant to solving linear inequalities. Peer reviews and editorial enhancements from the referees and editors at Scientific & Academic Publishing (SAP), USA is also essential for the correctness and readability of this paper.

References

[1]  Strang, Gilbert, “Introduction to Applied Mathematics”, John Wiley & Sons Inc., New York, 1979.
[2]  Strang, Gilbert, “Karmarkar’s algorithm and its place in applied mathematics”, The Msthematical Intelligencer 9(2): pp. 4-10, New York: Springer, 1987.
[3]  Dantzig, G. G. “Maximization of a llinear function of variables subject to linear inequalities”, 1947, Published pp. 339-347, in T.C. Koopmans (ed.): Activity Analysis of Production and Allocation, Wiley & Chapman-Hall, New York-Lodon, 1951.
[4]  Dantzig, G. B. “Linear Programming and Extensions”. Princeton, NJ: Princeton University Press, 1963.
[5]  Fukuda, Komei and Terlaky, Tamas, “Crisis-cross methods: A fresh view on pivot algorithms”, Mathematical Programming: Series B, No. 79, Papers from 16th International Symposium on Mathematical Programming, Lausanne, 1997.
[6]  Khachiyan, L. G., “Polynomial algorithms in linear programming”, U.S.S.R., Computational Mathematical and Mathematical Physics 20 (1980) pp. 53-72.
[7]  Karmarkar, N., “A New Polynomial Time Algorithm for Linear Programming”, AT&T Bell Laboratories, Murray Hill, New Jersey, September, 1984.
[8]  Gondzio, Jackek and Terlaky, Tamas, A computational view of interior point method”, Advances in linear and integer programming, Oxford Lecture Series in Mathematics and its Applications, 4, New York, Osford University Press. pp. 103-144, MR1438311, 1996.
[9]  Nocedal, Jorge and Wright, Stephen J, : “ Numerical Optimization”, Springer Science+Business Media, Inc., 1999.
[10]  Michael. R. Garey and David. S. Johnson, COMPUTERS AND INTRACTABILITY, A Guide to the Theory of NP-Completeness, Bell Laboratories, Murray Hill, New Jersey, 1979.
[11]  Nemirovsky, A. and Yudin, N. “Interior-Point Polynomial Methods in Convex Programming”, Philadelphia, PA: SIAM, 1994.
[12]  Alexander Schrijver, “Theory of Linear and Integer Programming”, Department of Econometrics”, Tilburg University, A Wiley-Interscience Publication, New York, 1979
[13]  Niedringhaus, W., “Stream Option Manager (SOM): Automated Integration of Aircraft Separation, Merging, Stream Management, and Other Air Traffic Control Problems”, IEEE Transactions Systems. Man & Cybernetics, Vol. 25 No. 9, Sept. 1995.
[14]  Niedringhaus, W., “Maneuver Option Manager (MOM): Automated Simplification of Complex Air Traffic Control Problems”, IEEE Transactions Systems. Man & Cybernetics, May 1992.
[15]  Marcotte, Patrice, A New Algorithm for Solving Variational Inequalties with Application to the Traffic Assignment Problem”, Centre de Recherche sur les Transports, University de Montreal, Canada, Mathematical Programming 33 (1985) pp. 339-351, North-Holland
[16]  Sun, Min, “A New Alternating Direction Method for Co-corecive Variational Inequality Problems with Linear Equality and Inequality Constraints”, pp. 161-176, Advanced Modeling and Optimization, Vol. 12, Number 2, 2010
[17]  Wang, Paul T. R., “Solving Linear Programming Problems in Self-dual Form with the Principle of Minmax”, MITRE MP-89W00023, The MITRE Corporation, July, 1989.
[18]  Wang, Paul T. R., Niedringhaus, William P., and McMahon, Matthew T., “A Generic Linear Inequalities Solver (LIS) with an Application for Automated Air Traffic Control”. America Journal of Computational and Applied Mathematics, pp 195-206, Volume 3, Number 4, August 2013.
[19]  EE236-A, Lecture 15,” Self-dual formulations”, University of California, Department of Electrical Engineering, 2007-08.
[20]  Self-Dual Form: http://www.ee.ucla.edu/ee236a/lectures/hsd.pdf
[21]  ILOG, “Introduction to ILOG CPLEX”, 2007, http://www.ilog.com/products/optimization/qa.cfm?presentation=3
[22]  ILOG, “CPLEX Barrier Optimizer”, 2008, http://www.ilog.com/products/cplex/product/barrier.cfm
[23]  Steven Skiena, “LP_SOLVE: Linear Programming Code”, Stony Brook University, Dept. of Computer Science”, 2008 http://www.cs.sunysb.edu/~algorith/implement/lpsolve/implement.shtml