Features

Solving Linear Equations Using C

Matt Weisfeld

Matt Weisfeld is currently employed by the Allen-Bradley Company in Highland Heights, Ohio. He is responsible for the design and development of test software on VAX/VMS, UNIX, DOS and other platforms. Matt is currently working on a book entitled Building and Testing Portable Libraries in C to be published by QED later this year. He can be reached on Compuserve at [71620,2171].

Resource allocation, a fundamental issue in any discipline, involves maximizing value and minimizing costs. As early as Junior High School, students solve simultaneous equations to find where these two variables break even. However, solving simultaneous equations by hand works only when the number of unknowns is small. When the problem involves dozens of equations and variables, pencil and paper just aren't enough. Linear Programming (LP) is used to solve these large resource-allocation problems. This article presents a program for solving linear equations using C.

Sample Application

To illustrate how Linear Programming works in a real life application, consider the following example. The Acme soft drink company markets two different lines of cola: diet and regular. Operating at maximum daily capacity, Acme can produce four tons of diet and six tons of regular. Most ingredients (such as water and sweetener) are obtained upon demand, so they do not affect potential output. However, the secret ingredient that makes Acme cola taste so special is limited to 24 pounds a day. Each ton of diet needs three pounds of the secret ingredient while each ton of regular needs four pounds. Finally, each ton of diet sells for $200,000, while each ton of regular sells for $600,000. The Acme company must maximize income by producing the optimal amount of diet and regular cola.

The standard LP problem consists of variables, an objective, and constraints. The variables in this example represent the two types of cola, diet (x1) and regular (x2). The goal is to allocate these variables in a fashion that maximizes the company's income (z). The objective function can be written

z = 2x1 + 6x2 (in units of $100,000)
To complete the LP model three constraints must be factored in:

x1 <= 4 (max diet in tons)
x2 <= 6 (max regular in tons)
3x1 + 4x2 <= 24 (in pounds)
The final problem is:

Maximize

z = 2x1 + 6x2 (objective function)
such that

x1 <= 4 (constraint 1)
x2 <= 6 (constraint 2)
3x1 + 4x2 <= 24 (constraint 3)
x1, x2 >=0 (it is impossible to have a negative amount of cola)

The Graphical Solution

Since this problem has only two variables, you can depict the solution graphically. Understanding the problem in terms of a graph makes the transition to the algorithmic solution much easier. However, once a third variable (and thus a third dimension) enters the picture, graphing becomes quite a challenge.

Figure 1 shows the graph of the Acme problem with only the constraints plotted. Notice the area contained within the constraint lines and the x1 and x2 axis. This area is called the Feasible Solution Space. Any point within this area will yield a combination of diet and regular cola that can be produced and yield some income. However, the goal is to find the best solution, not simply a feasible one. Notice also that there are five edges that surround the Feasible Solution Space. An edge is where two lines meet to form a single point. The optimal solution will be found at one of these edges.

To find the optimal solution, you need the objective function. Figure 2 shows two possible solutions based on a different value for z. Since the ratio of x1 to x2 is known, the slope of the line can be calculated. By varying the position of the line (the slope is always the same so the lines must move in a parallel manner), different solutions to this problem can be explored. The values of x1 and x2 are obtained by examining the points where the objective function intersects the lines bounding the feasible solution space. Line 1 intersects the feasible solution space at edge b (point x1=0, x2=6: z=36). Notice that as the line moves away from the origin, the value of z increases. Also remember that the objective line must intersect the feasible solution space for a valid solution to exist. Looking at the graph, it is apparent that the largest value of z that intersects the Feasible Solution Space is at edge c (point x1=1.5, x2=6, z=39). An optimal solution not falling on an edge represents a situation where the objective function is parallel to one of the constraints and more than one (actually infinite) optimal solution exists.

The Standard Form of the LP Model

In order to have a mathematical algorithm that can solve general LP models as well as lend itself to a computer implementation, all LP problems must follow the Standard Form of the LP model. The Standard Form includes these features:

To satisfy the first rule the constraints must be converted to equations. This may seem confusing, but in their present form they are not now equations. The operator used in all the constraints is <=. To make the constraints conform to the Standard Form, equal signs must separate the left-hand-side from the right-hand-side. However, simply replacing the <= with an = is inappropriate since it actually alters the constraint (eg: i<=1 is obviously not the same as i=1). To properly convert the constraint to an equation, a slack variable must be added. The slack variable accounts for the fact that the <= represents a wide range of values and permits the use of the equality.

Even though the slack variables are necessary in the conversion of the constraints, they do not contribute to the value of the objective function. Only x1 and x2 affect the income. The slack variables are created to assist in the computations, they do not represent real unknowns. To illustrate this fact, the objective function multiplier for each slack variable is 0. Thus the Acme Cola problem, represented in the Standard Form, is:

Maximize

z = 2x1 + 6x2 + 0s1 + 0s2 + 0s3 (objective function)
such that

x1 + s1 = 4 (constraint 1)
x2 +s2 = 6 (constraint 2)
3x1 + 4x2 s3 = 24 (constraint 3)
x1, x2 >= 0 (it is impossible to have a negative amount of cola)
In this case, the slack variables are all positive because all the constraints are of the <= type. If the constraints were >=, the slack variables would be negative. Figure 3 shows how the slack variables are represented graphically. There is a direct relationship between the number of graph sides (5) and the number of variables (x1, x2, s1, s2, s3).

The Starting Solution

To begin solving the Acme problem you must identify an initial, feasible starting solution. Since the optimal solution must be an edge, this problem has five possible candidates (as can be seen on the graph). Normally, you would begin at the origin when a problem is bounded by the x1 and x2 axis (where x1, x2 > 0), since this is both a feasible solution and an edge. The concept of the algorithm (presented later) is to move from the current edge to the next adjacent edge on a path to the optimal solution. The direction to move is determined by the algorithm. Thus, from the origin (edge a in Figure 2) , the algorithm moves to edge b and then stops at edge c, the optimal solution. The algorithm must follow adjacent edges (that is, the solution cannot move directly from edge a to edge c).

The Initial Table

A table format is used to group all the information needed to solve the problem. The C data structures mimic this table. Table 1 contains the initial table. The top row of the table is simply a header. The second row represents the objective function. The remaining number of rows are variable and depend on the number of constraints.

The first column identifies the variables that are currently in the basic solution. As mentioned before, the starting edge is at the origin (x1 and x2 are zero). This leaves the slack variables (s1, s2, s3) as the variables in the starting solution, which is called the basis. Notice that the rows representing the slack variables form an identity matrix (ones down the horizontal). Note, if some of the constraints had been >=, the corresponding slack variables would be negative and thus not a candidate for entry into the basis. In this instance, the algorithm presented would need to be modified to construct a proper basis. Finally, the column at the far right represents the current values for the variables in the basis.

The C Data Structures

Since the user interface is not a topic of this discussion, all data structures are defined internally. However, in any practical application the input data should not be hard-coded. Listing 1 contains all definitions. The C program includes three arrays that hold all the data necessary to represent a table. The arrays basis and objective hold the names (labels) that represent the variables. The actual computations are performed on the array called table. The dimensions of the table are always rows by columns. In the example table, the dimensions are four rows by seven columns. The only other information needed are the number of variables and constraints, which in this example are two and three respectively.

The Simplex Method

The algorithm used to solve this problem is called the simplex method. With this method you determine if there are any non-basic variables that, if in the basic solution, would enhance the final result. Based on the initial table, the basis variables are s1, s2, and s3. This is a feasible solution, however, since the variables x1 and x2 are not in the basis, the result is to produce nothing (z=0). To identify any variables that should enter the basis, inspect the objective function. The objective function is

 z  x1  x2  s1  s2  s3  sol
---------------------------
 1  -2  -6  0   0   0   0
Any variable that is negative (in a maximization problem) will increase the value of z, and so should enter the basis. If there is more than one negative value, the most negative one is chosen. Thus, there are two candidates for the entering value (--2 and --6). In this case, x2 should enter the basis.

It is obvious that if one variable enters the basis, another currently in the basis must leave. Since s1, s2, and s3 are in the basis, one must leave if x2 is to enter. Determining the leaving variable is done by inspecting the numbers in the column under entering variable, x2. They are

       x2
     -----
s1    0.0
s2    1.0
s3    3.0
All numbers that are not positive are not considered, thus s1 cannot leave the basis at this time. To determine which candidate actually leaves the basis, the ratios between the values in the entering column and the values in the solution are taken. For example, since there are two possible candidates to leave the basis, the ratios considered are as follows:

    sol     x2   ratio
   -----  -----  -----
s2  :  6    / 1   =  6
s3  : 24    / 3   =  8
The lowest value is chosen to leave the basis, thus s2 receives the honor. This makes the value 1.00, at the intersection of column x1 and row s2, the pivot element.

Now that the entering and leaving variables have been identified, the next iteration of the table is calculated by using the Gauss-Jordon method. Looking at the table, it is evident that for x2 to become part of the basis, it must replace s2 as part of the identity matrix. To do this the pivot element must be a 1. Thus, all values in the leaving equation are divided by the pivot element. In this instance, the value of the pivot element is already a 1 so the division does not change anything. The row that was s2 and is now replaced by x2 is called the pivot equation.

As mentioned above, placing x2 in the basis means that it must replace s2 as part of the identity matrix. Since the pivot element is now a 1, all the other values in the x2 column must be 0. To do this the transformation

new equation = old equation - (entering column coefficient 
   new pivot equation)
is applied to all elements of all the rows (including the objective function) other than the pivot equation.

Thus the objective function transformation is

old  :1      -2       -6       0         0        0        0

     -(-6)0  -(6)(0)  -(6)(1)  -(6-(0)  -(6)(1)  -(6)(0)  -(6)(6)
-------------------------------------------------------------------

new  :1      -2       0        0         6        0        36
The s1 equation is unchanged due to the fact that the entering column coefficient is already zero. The s3 equation transformation is as follows:

old  :0     4        3        0        0        1        24

     -(3)0  -(3)(0)  -(3)(1)  -(3)(0)  -(3)(1)  -(3)(0)  -(3)(6)
------------------------------------------------------------------

new  :0     4        0        0        -3       1        6
The first iteration is complete and the second table is presented in Table 2. The important information that can be gathered from inspecting this table is in the solution column. The 36 in the solution column of the objective function indicates that with the current solution the ACME Co. can make $36.00 (in units of 100,000). To determine the product mix, examine the variables in the basis. In this case only x2 is in the basis. Its solution is 6 (if six units of x2 are produced, ACME will make $36.00). However, since the objective function still has a negative value, the optimum solution has not yet been reached.

The variable x1, which has a value of -2, is the next to enter the basis. The ratios are:

   s1: 4/1 = 4
   s3: 6/4 = 1.5
Thus s3 is the leaving variable. In this case the pivot element is not one, so all the elements of the pivot equation are divided by the pivot element, which is four. The new pivot equation is as follows:

   0   1   0   0   -0.75   0.25   1.5
Applying the Gauss-Jordon technique yields the table in Table 3. The first thing to notice is that there are no negative values in the objective function. This means that the optimal solution has been found (z=39). Also notice that both variables, x1 and x2, are in the basis. Thus, to maximize profits, ACME should produce 1.5 units of diet and 6 units of regular.

The C Routines

The amount of C code necessary to perform the table operations is surprisingly small. Listing 2 contains the code for the mainline which loops until there are no more entering variables, printing each iteration of the table. The entering variable (enter_pos) and leaving variable (leave_pos) are determined by the code in Listing 3. If enter_pos is 0 then the optimal solution has already been obtained and the loop terminates. The location of the pivot element is placed in the variable pivot_element.

The routines new_pivot and new_equation are presented in Listing 4. new_pivot calculates the next pivot equation by dividing all elements in the pivot row by the pivot element, while new_equation performs the Gauss-Jordon algorithm on the entire table and completes the current iteration.

To help illustrate the program's performance, and to aid in debugging, two additional routines are included. The first, build_basis, labels all the rows and columns for better readability when the table is printed. The second, print_table, simply prints the contents of the table so that the information it contains can be inspected. These routines are presented in Listing 5. All of the listing files can be combined into one program file for compilation and execution.

Conclusion

Linear Programming is a powerful tool when attempting to allocate resources. The example used in this article presents a solution for one specific case, a maximization problem with all constraints of the <= variety. As stated before, there are a number of variations to the standard LP model. In practice, there can be a mixture of constraints (>=, <=, or = ) and the final goal may be to minimize the result (as in costs). When the constraints are mixed, extensions to the simplex algorithm are required. Furthermore, the results obtained from the final table are anything but final. All the numbers in the table represent useful information. (The act of gathering this information is called sensitivity analysis.)

References

Hamdy A. Taha. 1982. Operations Research, An Introduction, 3rd ed. New York, NY: MacMillan Publishing Company, Inc.

Hillier, Frederick S. and Gerald J. Lieberman. 1980. Introduction to Operations Research, 3rd ed. San Francisco, CA: Holden-Day Inc.