Optimization in Practice with MATLAB for Engineering Students and Professionals (2015)

---

11

---

Linear Programming

11.1 Overview

Linear programming (LP) is a technique for optimizing a linear objective function, subject to linear equality and linear inequality constraints. Linear programming is an important field of optimization for several reasons (Refs. [1, 2]). Many practical problems in engineering and operations research can be expressed as linear programming problems. Linear programming problems can be solved in an easy, fast, and reliable manner.

Linear programming was first used in the field of economics, and still remains popular in the areas of management, economics, finance, and engineering. During World War II, George Dantzig of the U.S. Air Force used LP techniques for planning problems. He invented the Simplex method, which is one of the most popular methods for solving LP problems.

Linear programming is a well developed subject in operations research, and several references that focus on linear programming alone are available [3, 4, 5, 6]. The reader is encouraged to consult these references for a more detailed discussion of linear programming.

In this chapter, the basics of linear programming problems will be studied. First, an introduction of the basic terminology and important concepts in LP problems will be presented in Sec. 11.2. In Sec. 11.3, the graphical approach for solving LP problems will be discussed, and four types of possible solutions for an LP problem will be introduced. Details on the MATLAB linear programming solver are presented in Sec. 11.4. Sections 11.5 and 11.6 discuss thevarious concepts involved in the Simplex method. In Sec. 11.7, the notion of duality and interior point methods are briefly discussed. Section 11.8 concludes the chapter with a summary.

11.2 Basics of Linear Programming

In this section, some basic terminology and definitions in linear programming will be discussed.

A generic linear programming problem consisting of linear equality and linear inequality constraints is given as

(11.1)

subject to

	(11.2)
	(11.3)
	(11.4)
	(11.5)
	(11.6)
	(11.7)

where z is the linear objective function to be minimized; c represents the vector of coefficients for each of the n design variables; A and b represent the matrix and vector of coefficients for m inequality constraints, respectively; A_eqand b_eq represent the matrix and vector of coefficients for p equality constraints, respectively; and x_i-lb and x_i-ub are the lower and upper bounds on the i-th design variable, respectively. Note that the above problem is not defined inthe so-called standard form, which is commonly used in some LP solvers. The standard form will be discussed later.

In a matrix notation, the above formulation can be written as

(11.8)

subject to

	(11.9)
	(11.10)
	(11.11)

where x_lb and x_ub are vectors of the lower and upper bounds on the design variables, respectively. Note that in some references, the objective function, z, is called the cost function. The coefficients c₁,...,c_n are also known as cost coefficients. Some references pose the LP problem as a maximization problem, which is different from the convention in the present chapter, a minimization problem.

11.3 Graphical Solution Approach: Types of LP Solutions

The graphical approach is a simple and easy technique to solve LP problems. The procedure involves plotting the contours of the objective function and the constraints. The feasible region of the constraints and the optimal solution is then identified graphically. Note that the graphical approach to solving LP problems may not be possible as the number of variables increases.

There are four possible types of solutions in a generic LP problem: (1) unique solution, (2) segment solution, (3) no solution, and (4) solution at infinity. These four types of solutions will be studied using examples.

11.3.1 The Unique Solution

Consider the following LP problem.

(11.12)

subject to

	(11.13)
	(11.14)
	(11.15)

Plot the constraints and the objective function contours, as shown in Fig. 11.1(a). There is a unique solution for this problem, where the objective function contour has the least value while remaining in the feasible region.

Figure 11.1. Types of LP Solutions

Here is an important point to keep in mind. If a constant value is added to the objective function in the above problem, for example, x₁ - 2x₂ + 5, and the constraints remain unchanged, how would it affect the optimal solution?The optimal values of the design variables will not change because of the added constant. This is because the additional constant in the objective function does not change the slope of the function contours. However, the optimalvalue of the objective function will change when compared to the original objective function.

11.3.2 The Segment Solution

Consider the LP problem with infinitely many solutions in the form of a segment solution. Examine the following example:

(11.16)

subject to

	(11.17)
	(11.18)
	(11.19)

From Figure 11.1(b), the slope of the objective function and the slope of the constraint function, -0.5x₁ + x₂ ≤ 10, are the same. Upon optimization, the objective function coincides with the constraint function, resulting in infinitely many solutions along the segment shown in Fig. 11.1(b).

A discussion of the third type of LP solutions is presented next.

11.3.3 No Solution

In this case, we will study LP problems that do not have a feasible solution. Consider the following example.

(11.20)

subject to

	(11.21)
	(11.22)
	(11.23)

Observe the feasible region of this problem from Fig. 11.1(c). The respective feasible regions of the inequality constraints do not intersect. There is no solution that satisfies all the constraints. Thus, there is no solution to this problem.

11.3.4 The Solution at Infinity

Consider the following problem.

(11.24)

subject to

	(11.25)
	(11.26)

The above problem leads to a solution at infinity. The feasible region, as illustrated in Fig. 11.1(d), is not bounded to yield a finite optimum value. The solution for this problem lies at infinity.

An unconstrained linear programming problem has a solution at infinity. These problems are rarely encountered in practice. On the contrary, unconstrained nonlinear programming problems are fairly common in practice. Nonlinear optimization problems do not require the presence of a constraint in order to yield a finite optimum.

Thus far, four possible types of LP solutions have been described. The graphical approach to solving LP problems is feasible only for small scale problems. For large scale problems, software implementations of LP solution algorithms are used. However, note that graphical visualization means can be used to examine the solution obtained for a large scale LP problem.

11.4 Solving LP Problems Using MATLAB

In this section, the command linprog in MATLAB is revisited. The MATLAB LP solvers use variations of the numerical methods that will be learned later in this chapter. The command linprog employs different solution strategies, such as the Simplex method and the interior point methods based on the size of the problem. The command allows for linear equality and linear inequality constraints and bounds on the design variables.

Different software codes follow their own respective formulations. Before using the solver, the problem is transformed into the formulation specified by the solver. The default problem formulation for linprog is given as

	(11.27)

subject to

	(11.28)
	(11.29)
	(11.30)

In the above formulation, f is the vector of coefficients of the objective function; A and A_eq represent the matrices of the left-hand-side coefficients of the linear inequality and linear equality constraints, respectively; b and b_eqrepresent the vectors of the right-hand-side values of the linear inequality and equality constraints, respectively; and x_lb and x_ub are the lower and upper bounds on the design variables, respectively.

Example: Consider the following example to demonstrate the use of the linprog command.

(11.31)

subject to

	(11.32)
	(11.33)
	(11.34)

The constraints can be rewritten as -4x₁ + 6x₂ ≤ 9 and x₁ + x₂ ≤ 4, according to the MATLAB standard formulation. The MATLAB code that solves the above problem is summarized below.

f = [1;-2] % Defining Objective
A= [-4 6; 1 1]; % LHS inequalities
b = [9; 4] % RHS inequalities
Aeq = []; % No equalities
beq = []; % No equalities
lb = [0;0] % Lower bounds
ub = [] % No upper bounds
x0 = [1;1] % Initial guess

x = linprog(f,A,B,Aeq,beq,lb,ub,x0)

The output generated by MATLAB is given below.

Warning: Large scale (interior point) method uses a
built-in starting point;
ignoring user-supplied X0.
> In linprog at 235
Optimization terminated.

x =

1.5000

2.5000

The warning displayed above informs the reader that the interior point algorithm used by the MATLAB solver does not need a starting point. The starting point provided has been ignored by the solver. These messages are warningsthat often help the user understand the solver algorithm and should not be mistaken for error messages.

Thus far, we have studied the graphical approach to solving LP problems, followed by the use of the MATLAB software options. Next, we present a popular numerical technique for solving LP problems that forms the basis of many commercial solvers.

11.5 Simplex Method Basics

This section introduces the Simplex method. Before discussing the algorithm of the Simplex method, some terminology and basic definitions associated with the Simplex method are presented.

11.5.1 The Standard Form

In order to apply the Simplex method, the problem must be posed in standard form. The standard form of an LP problem for the Simplex method is given below.

(11.35)

subject to

	(11.36)
	(11.37)
	(11.38)
	(11.39)
	(11.40)

The standard formulation does not contain inequality constraints. In a matrix notation, the standard formulation can be written as follows.

(11.41)

subject to

	(11.42)
	(11.43)

where c = [c₁,...,c_n] is the vector of the cost coefficients for the objective function; A is an m × n matrix of the coefficients for the linear equality constraints; and b is the m× 1 vector of the right-hand-side values.

The feasible region of the standard LP problem is a convex polygon. The optimal solution of the LP problem lies at one of the vertices of the polygon. In the Simplex method, the solution process moves from one vertex of the polygon to the next along the boundary of the feasible region.

11.5.2 Transforming into Standard Form

In the standard definition of LP problems, there are no inequality constraints. All design variables have non-negativity constraints. If a given problem is not in this standard form, certain operations are performed to transform the problem into the Standard Form. How the operations are performed for inequality constraints and for unbounded variables is discussed below.

Inequality Constraints

If the given problem formulation contains inequality constraints of the form g(x) ≤ 0, they are transformed into equality constraints by using slack variables. For constraints of the form g(x) ≤ 0, add a non-negative slack variable, s₁, to the left-hand-side of the constraint, yielding g(x) + s₁ = 0. The variable s₁ is called a slack variable because it represents the slack between the left-hand-side and the right-hand-side of the inequality.

For inequalities of the form g(x) ≥ 0, they are transformed into equality constraints by using surplus variables. Subtract a non-negative surplus variable, s₂, from the left-hand-side of the constraint, yielding g(x) - s₂ = 0. The variable s₂ represents the surplus between the left-hand-side and the right-hand-side of the inequality.

The slack/surplus variables are unknowns, and will be determined as part of the LP solution process.

Example: Consider the following linear programming problem.

	(11.44)

subject to

	(11.45)
	(11.46)
	(11.47)

The standard form for the above formulation can be given as

	(11.48)

subject to

	(11.49)
	(11.50)
	(11.51)

The standard form for the LP formulation requires the design variables to be non-negative.

Unbounded Design Variables

In the standard form, the design variables should be non-negative, x ≥ 0. However, in practice, the bounds on the design variable could also be of the form x ≤ 0. In some problems, the design variables may be indefinite (i.e., no bounds may be specified). In these cases, how are the design variables put into a standard form?

If a design variable, x_i, does not have bounds imposed in the problem, the following technique is used. Let x_i = s₁ - s₂, where s₁,s₂ ≥ 0. In the standard form, the variable x_i is then replaced by s₁ - s₂, and the additional constraintss₁,s₂ ≥ 0 are added to the problem. In other words, the unbounded design variable is rewritten as the difference between two non-negative additional variables. The additional variables will be determined as part of the LP solution process.

Example: Consider the following formulation.

	(11.52)

subject to

	(11.53)
	(11.54)
	(11.55)

Note that the variable x₃ is unbounded in the above formulation. Assume that x₃ = s₁ - s₂ and s₁,s₂ ≥ 0. The standard formulation can be written as follows.

(11.56)

subject to

	(11.57)
	(11.58)
	(11.59)

Next, we introduce the Gauss Jordan method, which forms an important component of the Simplex method.

11.5.3 Gauss Jordan Elimination

Note that the number of variables (including the n design variables and the p slack variables) is not necessarily equal to the number of equations, m. If the number of variables is equal to the number of equality constraints, then the solution is uniquely defined. In most LP problems, there exists more variables than equations. This results in an under-determined system of equations, resulting in infinitely many feasible solutions for the equality constraint set. The optimization problem then lies in determining which feasible solution(s) results in the minimization of the objective function, while satisfying the non-negativity constraints for the design variables.

Example: Consider the following under-determined system of equations:

	(11.60)
	(11.61)
	(11.62)

The above set of equations have more variables than equations. Therefore, there are infinitely many solutions for this case.

To efficiently work with the constraint set of the LP problem, reduce the constraint set into a special form. The set of equations in the special form is said to be in a canonical form. The original constraint set and the canonical form are equivalent (i.e., they have the same set of solutions). By transforming the LP constraint set into a canonical form (which is easier to solve), the solutions can be found more efficiently. To solve LP problems, use a canonical form known as the reduced row echelon form. This approach of using a reduced row echelon form to solve a set of linear equations is known as the Gauss Jordan elimination.

How can a given constraint be reduced into a row echelon form? A canonical form is usually defined with respect to a set of dependent or basic variables, which are defined in terms of a set of independent or non-basic variables.The choice of basic variables is arbitrary. A system of m equations and n variables is said to be in a reduced row echelon form with respect to a set of basic variables, x₁,...,x_m, if all the basic variables have a coefficient of one in only one equation, and have a zero coefficient in all the other equations. A generic matrix based representation of a set of equations in a reduced row echelon form with m basic variables and p non-basic variables is given as

(11.63)

where d is the matrix of coefficients for the non-basic variables; x_b is the set of basic variables; and x_nb is the set of non-basic variables. The number of basic variables is equal to the number of equations.

The reduced row echelon form is illustrated with the help of the following example.

Example: The following equations are in a reduced row echelon form with respect to the variables x₁,x₂,x₃, and x₄.

	(11.64)
	(11.65)
	(11.66)
	(11.67)

Writing the above equation set in a matrix form allows us to obtain the following.

(11.68)

Notice the matrix of coefficients on the left-hand-side. Each row of the matrix represents one constraint. We note that the coefficients corresponding to the basic variables x₁,x₂,x₃, and x₄ are equal to one in only one equation, and zero elsewhere.

The given set of constraints of an LP problem is not usually given in the row echelon form. A series of row operations must be performed to transform it into the standard form. The details are discussed next.

11.5.4 Reducing to a Row Echelon Form

A pivot operation consists of a series of elementary row operations to make a particular variable a basic variable (i.e., reduce the coefficient of the variable to unity in only one equation, and to zeroes in the other equations). A particular variable to be made basic, x_bi, could exist in some or in all of the m equations. We need to decide which coefficient of the basic variable should be made one (i.e., in which equation). By definition, each equation is to have only one basic variable with unit coefficient. The choice as to which equation corresponds to which basic variable is arbitrary, or is based on algebraic convenience for the reduced row echelon form. There will be further restrictions on this issue when the Simplex method is discussed.

There are two types of elementary row operations that can be performed to reduce a set of equations into a reduced row echelon form: (1) Multiply both sides of an equation with the same non-zero number. (2) Replace one equation by a linear combination of other equations.

Example: Reduce the following set of equations into a row echelon form.

	(11.69)
	(11.70)
	(11.71)

Choose x₁, x₂, and x₃ as basic variables. In order to obtain a row echelon form, perform operations such that the variables x₁, x₂, and x₃ appear in only one equation with a unit coefficient, and do not appear in the other equations. Use the following representation, where the first four columns represent the coefficients of each variable, and the last column represents the right-hand-side of the equation.

(11.72)

First, make x₁ a basic variable. Choose x₁ to have a unit coefficient in R₁, and zero coefficients in R₂ and R₃. Replace R₂ by R₂+R₁ and R₃ by R₃- R₁ to obtain

(11.73)

With the above transformations, x₁ is made a basic variable (appears only in R₁). Now make x₂ a basic variable. That is, make its coefficient one in R₂ and zeroes in the other rows. First, divide R₂ by 4 to obtain a unit coefficient in R₂.

(11.74)

Now replace R₁ by R₁ - R₂ and R₃ by R₃ - R₂ to make the coefficients of x₂ zeroes in the other equations to obtain

(11.75)

Next, make x₃ a basic variable by making its coefficient unity in R₃ and zeroes in the other equations. Divide R₃ by 3 to obtain the following.

(11.76)

Replace R₂ by R₂ + 2R₃ and R₁ by R₁ - R₃ to make the coefficients of x₃ zeroes in the other equations.

(11.77)

The above set of equations are in the reduced row echelon form with respect to the basic variables x₁,x₂, and x₃. Note that the particular choice of making x₁, x₂, and x₃ the basic variables is arbitrary.

11.5.5 The Basic Solution

Thus far, our discussion has dealt with how to reduce a given set of equations into a canonical form. The next step is to find the solution for the set of equations in the canonical form. A basic solution is obtained from the canonicalform by setting the non-basic (or independent) variables to zero. A basic feasible solution is a basic solution in which the values of the basic variables are non-negative. In other words, the basic feasible solution is feasible for the standard LP formulation, whereas a basic solution is not always feasible.

Example: Consider the following canonical form derived previously.

(11.78)

Set the non-basic variable, x₄, to zero. Then, obtain the basic solution by solving for x₁, x₂, and x₃. The basic solution can then be written as x₁ = -1,x₂ = 4,x₃ = 1, and x₄ = 0. Note that this basic solution is not a basic feasible solution for the standard LP formulation since x₁ < 0.

The choice of the set of basic variables is not unique, and can be decided according to computational convenience. In the example considered above, x₂,x₃, and x₄ could have been chosen as basic variables instead of x₁, x₂, and x₃. For a generic problem, any set of m variables (recall that there are m equations) from the possible n variables can be chosen as basic variables. This implies that the number of basic solutions for a generic standard LP problem with mconstraints and n variables is given as

(11.79)

Example: Consider the following set of equations from the previous example.

	(11.80)
	(11.81)
	(11.82)

Here, m = 3 and n = 4. Therefore, C₃⁴ = = 4 basic solutions.

As mentioned earlier, the feasible region of the standard LP problem is a convex polygon. Each vertex of this convex polygon corresponds to a basic feasible solution of the constraint set. The optimal solution of the LP problem,which is the basic feasible solution with the minimum objective function value, lies at one of the vertices of the polygon. In the Simplex method, the solution process moves from one vertex of the polygon to the next along the boundary of the feasible region. Now that the supporting mathematics of the Simplex method has been presented, we may proceed to discuss the Simplex algorithm.

11.6 Simplex Algorithm

11.6.1 Basic Algorithm

The optimal solution of the LP problem lies at one of the vertices of the feasible convex polygon. A cumbersome approach to solve an LP problem would be to list all the basic feasible solutions from which the optimal solution could be chosen. In the Simplex method, the solution process efficiently moves from one basic feasible solution to the next, while ensuring objective function reduction at each iteration (Ref. [7]).

Before discussing the details of the Simplex algorithm, please note the following. The standard LP problem was posed as a minimization problem. The rules of the following algorithm apply to minimization problems only. Some other textbooks may treat the standard LP problem as a maximization problem, and the rules of the algorithm would differ accordingly.

The Simplex method algorithm is presented in 6 generic steps:

1.Transformation into the Standard LP Problem

2.Formation of the Simplex Tableau

3.Choice of the Variable that Enters the Basis - Identify the Pivotal Column

4.Application of the Minimum Ratio Rule - Identify the Pivotal Row

5.Reduction to Canonical Form

6.Checking for Optimality

The presentation of these steps follows. The algorithm of the Simplex method is summarized in Fig. 11.2.

Figure 11.2. Simplex Algorithm

1.Transform into Standard LP Problem: Transform the given problem into the standard LP formulation by adding slack/surplus variables. Consider a generic formulation with n design variables, m inequality constraints, and mslack variables. The standard LP formulation is given as

(11.83)

subject to

	(11.84)
	(11.85)
	(11.86)
	(11.87)

Example: Consider the following LP problem.

(11.88)

subject to

	(11.89)
	(11.90)
	(11.91)

The standard form for the above problem is given below.

	(11.92)

subject to

	(11.93)
	(11.94)
	(11.95)

2.Form the Initial Simplex Tableau: List the constraints and the objective function coefficients in the form of a table, known as the Simplex tableau, shown below.

(11.96)

Example: The Simplex tableau for the example is shown in Table Reproduce.

The basic feasible solution for the initial Simplex tableau with s₁ and s₂ as basic variables is x₁ = 0,x₂ = 0,s₁ = 4, and s₂ = 3, and the function value is f = 0 (see Fig. 11.3).

Figure 11.3. Simplex Method Example

3.Choose the Variable that Enters the Basis - Identify the Pivotal Column: The above set of equations is in the reduced row echelon form with respect to the slack variables, s₁,...s_p. The Simplex algorithm begins with this initial basic feasible solution. By observing the coefficients of the objective function row in the initial Simplex tableau, the algorithm moves to the adjacent basic feasible solution that reduces the objective function. The other adjacent basic solution(s) with objective function value(s) higher than the current solution are ignored.

In order to move to the next basic feasible solution, the algorithm proceeds by making an existing basic variable a non-basic variable. In addition, an existing non-basic variable is made into a basic variable. How do we choose which non-basic variable becomes a basic variable? Similarly, how do we choose which will be the next basic variable?

The non-basic variable with the highest negative coefficient in the objective function is selected to become the basic variable in the next iteration. This choice is driven by our interest in minimizing the objective function value. The variable with the highest negative coefficient has the potential to reduce the objective function value to the maximum extent, when compared to the other variables. In other words, this choice can incur the maximum improvement to the objective value per unit of increase of the variable. In contrast, the coefficient of a non-basic variable would not have the opportunity to play an active role in the minimization process, since the value of allnon-basic variables are set to zero - in determining the basic feasible solutions. As the non-basic variable enters the basis (in the next iteration), it will directly impact the function minimization process. Similarly, in a functionmaximization process, the non-basic variable with the highest positive coefficient would be chosen to enter the basis in the next iteration.

According to the above approach, each iteration of the Simplex method makes a non-basic variable a basic variable. The corresponding variable is then said to enter the basis. When all the coefficients of the objective function are positive at any iteration, then the corresponding basic solution is the optimal solution. Therefore, the Simplex algorithm steps are repeated until all the coefficients in the objective function row are positive.

Example: Consider the initial Simplex tableau given in Table 11.1. Currently, the basic variables are s₁ = 4 and s₂ = 3 since they appear with unit coefficient in only one equation and have zero coefficients in the other equations. The non-basic variables are x₁ = 0 and x₂ = 0, and the corresponding function value is f = 0. Observing the entries of the objective function row (R₃), the coefficient of x₂ is negative (see bold-faced entries inTable 11.2). The current basic variables, s₁ and s₂, are not part of the objective function (see the corresponding coefficients in R₃ in Table 11.2). If x₂ is made a basic variable, the objective function value can be reduced from its current value. Therefore, the variable x₂ enters the basis, as per Table 11.2. (Note that enters the basis means will enter the basis in the next iteration.

Table 11.1. Initial Simplex Tableau