Maximize or minimize a function with the Simplex method

There is a very common problem in linear programming which consists to find the values that make maximum or minimum the value of a linear function, given a number of restrictions on the values of their variables. For example, you may want to find a minimum cost or maximum production function. To do this, there is an algorithm developed in 1947 by George Dantzig, called Simplex, which allows perform such calculations in a simple and effective way.

The function for which we want to find the maximum or minimum value is called target function, and should be a linear function of n variables, in the form:

A1x1 + A2x2 + … + Anxn

Being An any positive or negative number.

The constraints are any number of conditions, in the form of inequations, which must accomplish the variables so that the solution is valid. If the target is to maximize the function, they must be written as follows:

B1x1 + B2x2 + … + Bnxn <= C

Whereas if you want to minimize it, they must be written in the form:

B1x1 + B2x2 + … + Bnxn >= C

Where B and C can be any number, positive, negative or 0.

In this link you can find a mathematical approach to the Simplex method, and in this other you can download the Simplex project, with the sample code accompanying this article, written in csharp with Visual Studio 2013.

The use of the application is very easy; just enter the equation and its restrictions, one per line:

In the first line you have to write the target function, and the constraints, in an unlimited number, in the following lines, without leaving any blank lines. As the inequation changes depending if the operation is for maximize or minimize, I have replaced the symbols <= and >= by the sign : that will be interpreted as one of the above, as appropriate.

The variables can have any name that starts with a letter, followed by any number of alphanumeric characters, numbers can have decimal places, and you can use the + or - sign to link the terms of the equations. If you are interested in this topic, in this link you can find a series of articles on the design of the grammar for an expressions analyser, along with a sample implementation

Once the equations are written, you have to press the Build button to compile them, and then, you can press the Maximize or Minimize button and the program will calculate the values of the variables that make the function maximum or minimum with the given constraints.

In the program, the class that is responsible for performing the calculations is SimplexCalculator, and the classes that implement the expression parser will be found in the Expressions folder. Of these, the most important ones are Expression, which implements the target function, and its derived class, ExConstraint, which implements the constraints. Those classes have a Variables property which provides a list of Variable class objects with all the variables of the expression.

The algorithm consists in construct a matrix with the coefficients of the target function and the constraints, in addition to the canonical basis vectors of a space whose dimension depends on the number of constraints. If, for example, the equation system is as follows:

2x1 + 4x2 + x3
x1 + x2 + 2x3:500
10x1 + 8x2 + 5x3:2000
2x1 + x2 :100

To maximize it, the matrix is constructed as follows:

1 1 2 1 0 0 500
10 8 5 0 1 0 2000
2 1 0 0 0 1 100
2 4 1 0 0 0 0

In the first three columns are the coefficients of the constraints, in the first rows, and, in the latter, those of the target function. In the following three columns are the canonical base vectors, in this case of three dimensions, since there are three constraints. In the last column are the constants that limit the values of the variables, and, in the last row, the final value of the function, which starts with the value 0.

The algorithm is defined to maximize function, so, if you want to minimize it, you have to rewrite it so that it becomes a maximization problem, transforming the target function Z = AX, with the BX >= C constraints, in the Z = CY function, with the constraints BY <= A.

In the case of the previous system, the matrix for minimizing must be written in the form:

1 10 2 1 0 0 2
1 8 1 0 1 0 4
2 5 0 0 0 1 1
500 2000 100 0 0 0 0

And you will proceed with the same algorithm as for maximizing.

In the case of SimplexCalculator class, the constructor receives as parameters the target function in the target parameter, the list of constraints in the constraints parameter, and the maximize parameter, which is true when you want to maximize the function, and false to minimize it.

The constructor calls the CreateArray function, which is responsible for building the matrix with the initial data from the expressions.

The working array is a two-dimensional one, of type double, the first dimension for the rows and the second for the columns, which is stored in the _simplex variable. The coefficients of the expressions are obtained with the method Coefficient of the Expression class, passing it the variable whose coefficient you want to extract as a parameter.

It also build an array in the _base variable that will contain the indices of the columns with the vectors of the base, which will be moving from its place as they are performed the calculations. In another array, in the _ccolumns variable, they are stored the indices of the columns we use to perform calculations, which will also be moving from its place on each step of the algorithm, when the basis change is performed.

The list of variables of the target function are stored in a list in the _lvars variable.

Once the object is built, to perform calculations simply call the Calculate function, that will return the maximum or minimum value and assigns the final value to the different variables of the target function.

public double Calculate()
{
while (PerformStep()) ;
if (_maximize)
{
for (int ib = 0; ib < _base.Length; ib++)
{
if (_base[ib] < _lvar.Count)
{
_lvar[_base[ib]].Value =
_simplex[ib, _simplex.GetLength(1) - 1];
}
}
}
else
{
for (int iv = 0; iv < _lvar.Count; iv++)
{
_lvar[iv].Value =
-_simplex[_simplex.GetLength(0) - 1, iv + _ccolumns.Length];
}
}
return -Math.Round(_simplex[_simplex.GetLength(0) - 1,
_simplex.GetLength(1) - 1], cprecision);
}

This function repeatedly calls another function, PerformStep, which performs each of the steps of the algorithm until the entire process is complete and it returns false. Then, the final values are assigned to the different variables. In the case of a maximization, they will be obtained from the last column of the _simplex array, the positions stored in the _base array will contain the indices of the variables to which must be assigned each.

For the minimization, the values of the variables will be in the last row of the _simplex array, at the initial positions of the basis vectors.

The maximum or minimum value of the function is obtained from the position corresponding to the last row and last column of the _simplex array.

The PerformStep function performs each of the steps of the algorithm:

private bool PerformStep()
{
Point pivot = SelectPivot();
if ((pivot.X < 0) && (pivot.Y < 0))
{
return false;
}
for (int col = 0; col < _simplex.GetLength(1); col++)
{
if (col != pivot.X)
{
_simplex[pivot.Y, col] =
Math.Round(_simplex[pivot.Y, col] /
_simplex[pivot.Y, pivot.X], cprecision);
}
}
_simplex[pivot.Y, pivot.X] = 1;
for (int row = 0; row < _simplex.GetLength(0); row++)
{
if (row != pivot.Y)
{
for (int col = 0; col < _simplex.GetLength(1); col++)
{
if (col != pivot.X)
{
_simplex[row, col] =
Math.Round(_simplex[row, col] -
(_simplex[pivot.Y, col] *
_simplex[row, pivot.X]), cprecision);
}
}
_simplex[row, pivot.X] = 0;
}
}
int tmp = _base[pivot.Y];
_base[pivot.Y] = pivot.X;
for (int ix = 0; ix < _ccolumns.Length; ix++)
{
if (_ccolumns[ix] == pivot.X)
{
_ccolumns[ix] = tmp;
break;
}
}

for (int col = 0; col < _simplex.GetLength(1) - 1; col++)
{
if (_simplex[_simplex.GetLength(0) - 1, col] > 0)
{
return true;
}
}
return false;
}

At each step, a column and a row of the _simplex matrix is selected, which element will be used as a pivot for a basis change. The position of this element is calculated using the SelectPivot function, which returns a Point structure with the column in X and the row in Y.

Then all values in that row are divided by the pivot element, and, for each of the elements on the rest of rows, we subtracts them the product of the element in the pivot row by this in the same column corresponding to their row. The elements in the pivot column are set to 0, except the pivot itself, which is set to 1.

The position of the base vector, in the _base array, corresponding to the selected row, is exchanged with the corresponding column position in the _ccolumns array.

Finally, it is checked whether all the values in the last row are all 0 or negative, in which case the calculations have been completed and the function returns false.

The SelectPivot function selects the row and column of the item to be used as a pivot:

private Point SelectPivot()
{
Point cmm = new Point(-1,-1);
double mmval = double.MinValue;
for (int col = 0; col < _ccolumns.Length; col++)
{
if (_simplex[_simplex.GetLength(0) - 1, _ccolumns[col]] > 0)
{
int row;
double tn =
Math.Round(ThetaN(_ccolumns[col], out row) *
_simplex[_simplex.GetLength(0) - 1,
_ccolumns[col]], cprecision);
if (tn > mmval)
{
mmval = tn;
cmm = new Point(_ccolumns[col], row);
}
}
}
return cmm;
}

This function processes all the columns contained in the _ccolumns array and select the one that maximizes the product of the value returned by the Thetan function, which also returns the row with the element candidate to pivot for that column, by the element of the last row corresponding to that column.

Finally, the Thetan function selects, for a given column, the row where the quotient of the value of the last column in that row and the corresponding value in the column, is minimal:

private double ThetaN(int col, out int row)
{
row = -1;
double min = double.MaxValue;
for (int r = 0; r < _simplex.GetLength(0) - 1; r++)
{
if (_simplex[r, col] > 0)
{
double m =
Math.Round(_simplex[r, _simplex.GetLength(1) - 1] /
_simplex[r, col], cprecision);
if (m < min)
{
min = m;
row = r;
}
}
}
return min;
}

The function returns this minimum quotient and the row where it has been found, that will be the row corresponding to the element used as a pivot.