by Alan Jones (alan.jones@qcs.org)
This article describes the Linear Programming problem. LP software is being developed for the C64 to emphasize that the C64 is capable of solving HEAVY MATH problems. It describes some of the C64 program design issues and gives readers an opportunity it influence the development of the software.
Linear Programming, LP, is the simplest type of constrained numerical optimization problem. It's simplest (standard) form is:
max P = c'x
s.t. Ax = b
x >= 0
That is, we want to find the solution vector x which maximizes the function P, which is a linear combination of elements of x, while satisfying the linear equation Ax = b, and with each element of x >= 0. A is a matrix with typically more columns than rows. The maximization problem itself is easy, except for the inequalities and determining which elements of x are fixed at a bound and which are free.
Assume x is a vector of 20 variables and we have 10 equality constraints making A a 10 by 20 matrix. Solving the problem involves choosing 10 of the variables to be bound (at zero), eliminating those 10 columns from A and solving the reduced linear system, say B*xb = b, for the other 10 elements of x, xb. Then sort through the solutions to find the feasible solutions that satisfy all of the constraints, x >= 0, and chose from them the solution that will maximize the function P. Taking 20 columns of A 10 at a time we will have 184,756 solutions to solve and evaluate!
LP is an important problem that requires a computer and clever algorithms to solve non-trivial problems. The development of LP algorithms parallels the early development of computers. George Dantzig published his first paper on his simplex method in 1947. LP problems have long been studied in mathematics, economics, and business fields, but the simplex method and the computer were the big breakthroughs. Just imagine the problem of allocating and distributing limited resources to millions of soldiers in WW II. There are many types of problems that can be solved as LP problems, and some have specialized algorithms for their efficient solution.
The inequality constraint occurs naturally. Many processes are irreversible. For example you can burn fuel in rocket at a positive fuel flow rate to produce thrust, but you can't unburn the exhaust to refill the tanks. A table leg can push on the floor but not pull (unless bolted). A rope, cable, or chain can pull but not push. You can't buy or sell negative quantities of a new item. You can't start assembling a machine before the parts arrive.
LP problems can also be written in more general forms. Perhaps the most general is:
max P = c'x
s.t. bl <= Ax <= bu
xl <= x <= xu
Where bl an xl are vectors of lower bounds and bu and xu are vectors of upper bounds.
The LP forms can be manipulated with simple algebraic operations and by adding constraints and "slack" variables. Most numerical optimization problems are set up to minimize a cost function. Historically, LP problems are set up as maximization problems. If your function P is a total cost to be minimized, simply multiply the vector c by -1 and maximize instead. Multiplying the i'th constraint equation by -1 will change the sign of all the elements of the i'th row of A, bl, bu, and change the direction of the inequalities, but not change x(i). The i'th constraint above is actually two equations and can be written (neglecting the subscripts):
Ax <= bu
Ax >= bl
To convert them to equality constraints:
Ax + y1 = bu
Ax + -y2 = bl
y1 and y2 are additional variables both >= 0 added to take up the "slack" in the inequality equations. When the constraint is free (i.e. could be ignored) the slack variable is >0. When the constraint is active, the slack variable is constrained to zero. The slack variables are simply appended to the vector x and treated as normal variables.
Lower bounds on x can be removed with a change of variable, x = y - xl, without adding more rows or variables.
xl <= x <= xu, becomes 0 <= y <= xu - xl
Upper bounds on x can be handled by adding slack variables. The expanded problem could look like:
Ax = b
x + xs = xu
x >= 0, xs >= 0
This can double the length of the x vector and add an equal number of equality constraints. A fair amount of algebraic manipulation is possible to get a LP into the form required by the solution software. This can significantly increase the problem size. Good software will use the more general forms of the LP program and use a more sophisticated solution logic instead of increasing the problem size.
Mixture problems are probably the easiest to describe for illustrating the importance of LP. Suppose you are producing cans of mixed nuts. You want to find the percentage of each nut to mix and package at the lowest cost, or maximum profit. You check the commodities prices of the various nuts available at different times and places. c(i) could be the price per pound of each nut, i, and x(i) would be the percent of that nut to be mixed. One constraint is that the percentages total 100. The Marketing Department may demand no more than 25% peanuts and no less than 5% of Cashews, pecans, and almonds. Each nut may be scored for crunchiness, and other aspects of taste to meet other constraints. Of course no x(i) can be less than zero. You can do the same thing with "real fruit juice" drinks. Mix the cheapest available fruit juices while balancing sweetness, acidity, color, taste, and other factors.
Another type of problem is project scheduling. Suppose you are going to design and build a bridge, or spacecraft. You identify all of the tasks that must be accomplished, and then estimate their completion time and constraints. E.g. task 43 can't begin until tasks 16, 17, and 41 have been completed. You could then set up a LP problem to schedule each task so as to complete the project in the minimum time. Various other performance criteria can be embedded in a LP problem.
LP is also a special case of Non Linear Programming, NLP. The NLP problem is similar to the LP problem except that the objective function P and constraint equations can be nonlinear functions. NLP is a much more expensive problem to solve. An approach to solving a NLP is to linearize the problem at a point to get a LP subproblem. Solving the LP gives a search direction for the NLP problem. You then advance in the search direction until you have made sufficient improvement in the objective function or diverged too far from the constraints. You then correct back to the constraints and linearize again. The sequence is repeated until the solution is found.
Computers were not created by secretaries and typesetters. Although computers can do a lot of things, they were created to do Heavy Math. The first electronic computer, the Attanasoff-Berry Computer, ABC, was built to solve systems of up to 29 linear equations, although it could be adapted to solve other problems as well. The second and third, ENIAC and EDVAC, were intended to solve ODEs (Ordinary Differential Equations, i.e. computing artillery range tables), as well as more general usage.(ENIAC's first operational use was solving a 25 by 25 set of linear equations from Los Alamos to see if an atomic bomb might be feasible. The ABC could have solved this problem, but it had a read/write reliability problem and Attanasoff was called up for other war related work instead of perfecting the ABC.) And I have already discussed the LP problem.
Software for the C64 has been readily available and published in the C64 related literature for solving systems of linear equations and numerically integrating differential equations (4th order Runge-Kutta being typical). However nothing is available or published for solving LP problems on a C64. I do recall seeing a one inch ad in some Commodore magazines trying to sell software to solve small LP problems on a C64, maybe 15-20 equations, but I never read any reviews of it and I always suspected it to be poor software. I intend to plug this hole by writing a good LP program for the C64.
This is Part 1 of at least a two part series on LP for the C64. Part 2 will include the presentation of the C64 LP program. I usually hate multi-part articles that could have been published as one large article. However, in this case Part 2, and indeed the C64 LP computer program, has not been written yet. Part 2 will be dependent on reader feedback! No feedback means no reader or user interest and I won't bother to submit Part 2. I have not written a LP program before, or even used one much on other systems. There is a very good chance that some readers will have experience with LP software and can offer practical suggestions for writing C64 LP software. I would also like to hear from any reader who may have some type of LP problem that they would like to be able to solve, perhaps something related to a hobby or home computing application. Also, someone may be able to point me to an ideal LP reference book or even provide suitable source code to examine.
Choosing a LP solution method for the C64 is not an easy task. Appropriate reference books and papers are hard to come by. There is nothing in the CBM literature that I am aware of. Many books have been published on the LP subject. However, most are old books written for business students or for training clerks to solve small problems by hand, and they use outdated methods. These are easily found in local libraries and at small colleges. LP is still an active research area and many new journal articles are still being written, as well as new books. However, the focus in on solving VERY large LP problems using sparse matrix algebra techniques, and newer methods such as Karmarker or interior-point methods. There is even a LP FAQ available. Two of the best references that I have found so far are: "Computer Solution of Linear Programs", J. L. Nazareth (1987), and "Advanced Linear Programming", B. A. Murtagh (1981). The LP source codes that I have looked at were either too crude, suitable only for small problems, too poorly documented, or suited only for large problems on large computers. This gives you some idea of where I am at, should you want to provide helpful feedback.
Our C64 (and 128) does arithmetic slowly, but reliably, and has limited memory. We would not want to solve a large LP problem with thousands of constraints and variables on a C64 even if we could. Still, we want to be able solve problems as large as practical using efficient methods. At the heart of the solution algorithm we must be able to solve an n by n system of equations, where n is the number of constraint equations. For, say n = 70, we will need 24,500 bytes to store the full matrix, leaving the rest for the OS/language/LP program. The simplex type solutions typically require computation time = K * n~3, so we will not likely want to solve any LP problems on the C64 with more than 70 constraints.
The desire to solve huge LP problems also motivated the development of sparse matrix algebra techniques. Matrices associated with LP tend be very sparse. The idea is to be able to store only the non-zero elements, to store them in data structures that can be used effectively by linear equation solvers, and to use larger but often slower external memory efficiently. These codes have a higher overhead in terms of code size and computation/FP multiply. However, these codes can compute the same results with substantially fewer multiplies than simpler code for dense or full matrices. They often perform faster overall, at least on scalar computers. Our 6502 can zip through complicated data structures and ML code fairly fast compared to the cost of a FP multiply. Using sparse matrix techniques is possible with the C64, but not required. I think the complications of using sparse matrix methods is not warranted in this case. Users wanting to solve large sparse problems are free to disagree.
There are several LP solution methods: (primal) simplex, revised simplex, dual simplex, primal-dual simplex, Karmarker, and others. There are also many variations among these. I have chosen the revised simplex method for the C64 program. Some other choices might be better for some types of LP problems, and I may also write a primal-dual simplex program.
In the revised simplex method, the constraint matrix A can be stored outside of main memory (on disk or REU) and brought in a column at a time. The matrix A can be stored in a sparse matrix form. This is probably a good idea for our slow disk storage, but unnecessary for our REU storage. A square matrix B will be based on a set of n columns of A. At each iteration a new column of A is chosen to enter B and an old column is re moved. Each column represents an element of x and the objective is to get the unbound variables selected into B with the remaining variables set to their limits. At each iteration we have to solve many linear systems using B with different right hand side vectors and then update B to account for the column addition and removal.
There are several ways to work with the matrix B. B can be explicitly inverted. The inverse can then be explicitly updated using the Sherman- Morrison-Woodburry formula, or updated in product form by storing a sparse matrix factor that accounts for the update. The later is preferred for sparse implementations, but the storage used grows with each iteration and eventually a new explicit inverse of B will have to be computed. Both methods are numerically unstable and can have excessive growth of roundoff errors, necessitating error checks and occasional reinversion.
B can also be efficiently used in an LU factored form. B = LU where L is lower triangular and U is upper triangular. The updating is difficult, especially when exploiting sparcity and external storage. Bartels and Golub (1969) developed a numerically stable way of updating the L and U factors that enabled L to be stored externally, with U stored in main memory. Forrest and Tomlin (1972) developed a method that allows both L and U to use external storage. This is the method used in most commercial codes for solving large LP problems. However, it is not numerically stable and needs to be monitored and refactored. Sanders (1976) developed a variation of the Bartels-Golub method that is stable and allows most of U to be stored externally. Fletcher and Matthews (1984) developed a stable method for updating the LU factors explicitly in main memory.
In the Fletcher-Matthews update we have: PBQ = LU. B is the unfactored matrix, P is a row permutation matrix, Q is a column permutation matrix, L is a unit lower triangular matrix, and U is an upper triangular matrix. L and U are stored explicitly (not in a factored form) in a square matrix B. L and U are first formed using a normal LU factorization (Gaussian elimination). P is a row permutation matrix that represents the partial pivoting normally used to stabilize Gaussian elimination. Q is an optional column permutation that represents the full pivoting that could be done with Gaussian elimination to get better numerical stability. P and Q are completely defined by integer arrays. In the update, one column of B is removed, the rest of the vectors are shifted left to fill the gap, and the entering vector placed in the right most column of B. L, U, and P are then updated to a stable factorization of the new B matrix. Pivoting is still required for stability, but only two adjacent rows can be exchanged. This is a weaker form of stability than Gaussian elimination with partial pivoting. It may be advisable to use some more extreme measures to assure greater accuracy, such as preliminary scaling of the problem (equilibration), full pivoting in the initial factorization, or doing a final refactorization. Q is not changed during the update (no column pivoting for stability), but we will have to keep track of where each column of A is located in B (or the LU factorization of B). The update is also cheaper when the column of B removed is farther to the right.
I have chosen to use the Fletcher-Matthews update. This choice is ideal for dense LP problems, but less so for sparse problems. It limits the number of constraints that we can handle. It also has consequences for other program design choices that must be made for the C64 LP program. Using the simple standard form keeps the program logic simple, but converting problems with inequality constraints and upper and lower bounds to standard form will make the problem size bigger. This will be more expensive to solve when not taking advantage of sparcity, and make our size limit more significant. I am leaning toward solving LP problems in the form:
Max P = c'x
s.t. Ax = b
xl <= x <= xu
There are many more program design choices to make. There are different strategies for selecting the entering and leaving variables (columns). There are different ways to pick a starting B matrix. There are different ways to perform "Phase 1" of the revised simplex algorithm. Many of the internal tests of floating point results involve tolerances. These tolerances do affect the performance of the program. However, I have only seen arbitrary suggested tolerances. These tolerances should be determined by the floating point precision used, the problems size, matrix condition number, and perhaps some problem specific parameters. I have not seen any numerical analyses indicating how to set the tolerances.
I have completed the subroutines to perform the Fletcher-Matthews update, matrix factorization, and solving linear systems. Your interest and feedback can influence other major parts of the program, or even convince me to use a different updating method. Or perhaps everyone would rather not see the sequel? You can reach me via e-mail at: alan.jones@qcs.org
Since I am writing this software it will be developed in the COMAL 2.0 language, which is fast and very readable (and thus easily translated to other languages). A typical commercial LP code represents about 10 man years of development. The C64 software will represent perhaps 10 days, plus past research, or 10 weeks of part time effort, and be submitted for the next issue of C=Hacking. It won't be the best software, but it will be good. Over a period of perhaps a year, we can incorporate and publish changes (hopefully not actual bugs). Then it can be translated be to ACE assembly or some other language to make it available to more C64/128 users.
The C64 LP software is being developed simply to fill an empty slot in C64 software. It is also to emphasize that the C64 can indeed be used for HEAVY MATH, subject only to constraints of limited speed, memory, and software. There is little demand for solving LP problems on a C64, otherwise it would have been done many times already. This article will not interest many C64/128 users. It does not even come close to describing how to write a working LP program. It does describe what LP is and some of the issues involved in designing a LP program for the C64. I do thank those that read this far. Don't forget to write.
C= Hacking Home | Issue 14 Contents
Copyright © 1992 - 1997 Commodore Hacking
Commodore, CBM, its respective computer system names, and the CBM logo are either registered trademarks or trademarks of ESCOM GmbH or VISCorp in the United States and/or other countries. Neither ESCOM nor VISCorp endorse or are affiliated with Commodore Hacking.
Commodore Hacking is published by:
Last Updated:
1997-03-31 by Jim Brain