Chapter 4

LINEAR PROGRAMMING

 Formulating the Linear Programming Problem-A Simple Refinery top

To this point in the discussion of linear programming, the emphasis has been on the solution of problems by the Simplex Method. In this section procedures will be presented for the formulation of the linear programming problem for a plant or process. This will include developing the objective function from the cost or profit of the process or plant and the constraint equations from the availability of raw materials, the demand for products, and equipment capacity limitations and conversion capabilities. A simple petroleum refinery will be used as an example to illustrate these procedures. Also an optimal solution will be obtained using a large linear programming code to illustrate the use of one of these types of programs available on a large computer. In the following section the optimal solution of the general linear programming problem will be extended to a sensitivity analysis, and these results will be illustrated using the information computed from the large linear programming code for the simple refinery example.

Figure 4-7 shows the flow diagram for the simple petroleum refinery, and Table 4-2 the definition is given for the name of each of the process streams. There are only three process units in this refinery, and these are a crude oil atmospheric distillation column, a catalytic cracking unit, and a catalytic reformer. The crude oil distillation column separates the crude oil into five streams: are fuel gas, straight-run gasoline, straight-run naphtha, straight-run distillate, and straight run fuel oil. Part of the straight run naphtha is processed through the catalytic reformer to improve its quality, i.e., increase the octane number. Also parts of the straight run distillate and straight-run fuel oil are processed through the catalytic cracking unit to improve their quality so they can be blended into gasoline. The refinery produces four products: premium gasoline, regular gasoline, diesel fuel and fuel oil. Even for this simple refinery there are 33 flow rates for which the optimal values have to be determined. This small problem points out one of the difficulties of large linear programming problems. The formulation of the problem is quite straight-forward. However, there is a major accounting problem in keeping track of a large number of variables, and the collection of reliable data to go with these variables is usually very time consuming (9).

Table 4-3 lists the capacities, operating costs, process stream, mass yields, and volumetric yields for the three process units in the refinery. These are typical of a medium size refinery in the Gulf coast area. The mass yields were taken from those reported by Aronfsky, Dutton, and Tayyaabkhan (10) and were converted to volumetric yields by using API gravity data. The operating costs were furnished by the technical division of a major oil company which has refineries on the Gulf Coast.

The quality specification and physical properties are given in Table 4-4 for the process streams, and the crude oil cost and the product sales prices are given in Table 4-5. The data in Table 4-4 was reported by Aronfsky et.al. (29), and the cost and prices in Table 4-5 were obtained from the Oil and Gas Journal (11). The information given in Table 4-3, 4-4, and 4-5 is required to construct the objective function and the constraint equations for the linear programming model of the refinery.

It is standard practice to present the linear programming problem for the refinery in matrix form, as shown in Figure 4-8. In the first row the coefficients of the terms in the objective function are listed under their corresponding variables. The sales prices are shown as positive, and the cost are shown as negative, so the problem is formulated to maximize the profit. These numbers were taken from Table 4-5, and it was convenient to combine the crude cost (\$32.00/barrel) with the operating cost of the crude oil atmospheric distillation column (\$1.00/barrel) to show a total cost of \$33.00 per barrel of crude oil processed in Figure 4-7. Consequently, the first row of Figure 4-8 represents the objective function given below:

-33.0 CRUDE  +  0.01965 FGAD - 2.50 SRNRF  +  0.01965 FGRF - 2.20 SRDSCC

-2.20 SRFOCC  +  0.01965 FGCC  +   45.36 PG  +  43.68 RG  +  40.32 DF  +  13.14 FO

The constraint equations begin with the second row in Figure 4-8. They are grouped in terms of quality and quantity constraints on the crude and products, in terms of the performance of the process unit using the volumetric yields, and in terms of the stream splits among the process units and blending into the products.

The second row is the crude availability constraint limiting the refinery to 110,000 barrels/day. This is followed by the four quantity and quality constraints associated with each product. These are the daily production and blending requirements and two quality constraints. These have been extracted from Figure 4-8 and are shown in Table 4-6 for the four products. The minimum production constraint states that the refinery must produce at least 10,000 bbl/day of premium gasoline to meet the company's marketing division's requirements. The blending constraints state that the sum of the streams going to produce premium gasoline, must equal the daily production of premium gasoline. The quality constraints use linear blending, and the sum of each component weighted by its quality must meet or exceed the quality of the product. This is illustrated with premium gasoline octane rating blending constraint, which is written as the following, using the information from the matrix:

78.5 SRGPG  +  104.0 RFGPG  +   65.0 SRNPG  +  93.7 CCGPG - 93.0 PG  >  0                                     (4-11)

Here the premium gasoline must have an octane number of at least 93.0. Corresponding inequality constraints are specified in Table 4-6 using the same procedure for premium gasoline vapor pressure, regular gasoline octane number and vapor pressure, diesel fuel density and sulfur content, and fuel oil density and sulfur content.

The next set of information, given in the constraint equation matrix, Figure 4-8, is the description of the operation of the process unit using the volumetric yield shown in Table 4-3. This section of the matrix has been extracted and is shown in Table 4-7 for the three process units. Referring to the volumetric yields for the crude oil distillation column, this data states that 35.42 times the volumetric flow rate of crude produces the flow rate of fuel gas from the distillation column, FGAD i.e.:

35.42 CRUDE - FGAD = 0                                      (4-12)

Corresponding yields of the other products from the crude oil distillation are determined the same way. For the catalytic reformer the yield of the Table 4-6 Quantity and Quality Constraints for the Refinery Products fuel gas (FGRF) and the reformer gasoline (RFG) are given by the following equations:

158.7 SRNRF - FGRF = 0                                      (4-13)

0.928 SRNRF - RFG = 0                                        (4-14)

Similar equations are used in the matrix, Figure 4-8, and are summarized in Table 4-7 for the catalytic cracking unit.

The use of volumetric yields to give linear equations to describe the performance of the process units is required for linear programming. The results will be satisfactory as long as the volumetric yields precisely describe the performance of these process units. These volumetric yields are a function of the operating conditions of the unit, e.g., temperature, feed flow rate, catalyst activity, etc. Consequently, to have an optimal solution these volumetric yields must represent the best performance of the individual process units. To account for changes in volumetric yields with operating conditions, sometimes a separate simulation program is coupled to the linear programming code to furnish best values of the volumetric yields. Then an iterative procedure is used to converge optimal linear programming flow rates with corresponding values of volumetric yields from the simulation program. (See Figure 6-5.)

The last group of terms in Figure 4-8 gives the material balance around points where streams split among process units and blend into products. The stream to be divided is given a coefficient of 1, and the resulting streams have a coefficient -1. For example, the straight-run naphtha from the crude oil distillation is split into four streams. One is sent to the catalytic reformer, and the other three are used in blending premium gasoline, regular gasoline, and diesel fuel. The equation for this split is:

SRN - SRNRF - SRNPG - SRNRG - SRNDF = 0                                      (4-15)

There are a total of seven stream splits as shown in Figure 4-8.

This information is now available to determine the optimum operating conditions of the refinery. The optimal solution was obtained using the Mathematical Programming System Extended (MPSX) program run on the IBM 4341. The format used by this linear programming code has become an industry standard according to Murtagh (12) and is not restricted to the MPS series of codes developed originally for IBM computers. Consequently, we will also describe the input procedure for the code because of its more general nature. Also, we will use these refinery results to illustrate the additional information that can be obtained from sensitivity analysis.

 Solving the Linear Programming Problem for the Simple Refinery top

Having constructed the linear programming problem matrix, we are now ready to solve the problem using a large linear programming computer program. The input and output for these programs has become relatively standard (12), making the study of one beneficial in the use of any of the others. The solution of the simple refinery has been obtained using the IBM Mathematical Programming System Extended (MPSX). The detailed documentation is given in IBM manuals (15,16) and by Murtagh(12) on the use of the program, and the following outlines its use for the refinery problem. The MPSX control program used to solve the problem is given in Table 4-8. The first two commands, PROGRAM and INITIALZ, define the beginning of the program and set up standard default values for many of the optional program parameters. TITLE writes the character string between the quotation marks at the top of every page of output. The four MOVE commands give user specified names to the input data (XDATA), internal machine code version of the problem (XPBNAME), objective function (XOBJ), and right-hand-side vector (XRHS). Next, CONVERT calls a routine to convert the input data from binary coded decimal (BCD), or communications format, into machine code for use by the program, and BCDOUT has the input data printed. The next three commands, SETUP, CRASH, and PRIMAL, indicate that the objective function is to be maximized, a starting basis is created, and the primal method is to be used to solve the problem. Output from PRIMAL is in machine code so SOLUTION is called to produce BCD output of the solution. The RANGE command is used in the sensitivity analysis to determine the range over which the variables, right-hand-sides, and the coefficients may vary without changing the basis. The last two statements, EXIT and PEND, signal the end of the control program and return control over to the computer's operating system.

Input to the MPSX program is divided into four sections: NAME, ROWS, COLUMNS, and RHS. The first two are shown in Table 4-9. The NAME section is a single line containing the identifier, NAME, and the user-defined name for the block of input data that follows. (MPSX has provisions for keeping track of several problems during execution of the control program). When the program is run it looks for input data with the same name as that stored in the internal variable XDATA. The ROWS section contains the name of every row in the model, preceded by a letter indicating whether it is a non-constrained row (N), the objective function, a less-than-or-equal-to constraint (L), a greater-than-or-equal-to constraint (G), or an equality constraint (E).

The COLUMNS section of the input data is shown in Table 4-10. It is a listing of the non-zero elements in each column of the problem matrix (Figure 4-8). Each line contains a column name followed by up to two row names and their corresponding coefficients from Figure 4-8.

The last input section is shown in Table 4-11. Here, the right-hand-side coefficients are entered in the same way that the coefficients for each column were entered in the COLUMS section, i.e., only the non-zero elements. The end of the data block is followed by an ENDATA card.

The solution to the refinery problem is presented in Table 4-12 (a) and (b), as listed in the printout from the MPSX program. It is divided into two sections, the first providing information about the constraints (rows), and the second giving information about the refinery stream variables (columns).

In the ROWS section there are eight columns of output. The first is the internal identification number given to each row by the program. The second column is the name given to the rows in the input data. Next is the AT column which contains a pair of code letters to indicate the status of each row in the optimal solution. Constraint rows in the basis have the code BS, non-basis inequality constraints that have reached their upper or lower limits have the code UL or LL, and equality constraints have the status code EQ. The fourth column is the row activity, as defined by the equation:

This is the optimal value of the left-hand side of the constraint equations. However, it is computed by subtracting the slack variable from the right hand side. The column labeled SLACK ACTIVITY contains the value of the slack variable for each row. The next three columns are associated with sensitivity analysis. The sixth and seventh columns show the lower and upper limits placed on the row activities. The final column, DUAL ACTIVITY, gives Lagrange Multipliers, which are also called the Simplex multipliers, shadow prices and implicit prices. As will be seen subsequently in sensitivity analysis, they will relate changes in the activity to changes in the objective function. Also, the dot in the table means zero, the same convention used in the Simplex Tableau.

Examination of this section of output shows that the activity (or value) of the objective function (row 1, OBJ) is 701,823.4, i.e., the maximum profit for the refinery is \$701,823.40 per day. Checking the rows which are at their lower limits, LL, for production constraints one finds that only row 15, FOMIN, is at its lower limit of 10,000 barrels/day, indicating that only the minimum required amount of fuel oil should be produced. However, row 3, PGMIN, row 7, RGMIN, and row 11, DFMIN, are all above their lower limits with values of 47,113 barrels/day for premium gasoline, 22,520 barrels/day for regular gasoline, and 12,491 barrels/day for diesel fuel. More will be said about the information in this table when sensitivity analysis is discussed.

The COLUMNS section of the solution also has eight columns. The first three are analogous to the first three in the ROWS section, i.e., an interval identification number, name of the column, and whether the variable is in the basis BS or is at its upper or lower limit, UL or LL. The fourth column, ACTIVITY, contains the optimal value for each variable. The objective function cost coefficients are listed in the column INPUT COST. REDUCED COST is the amount by which the objective function will be increased per unit increase in each non-basis variable and is part of the sensitivity analysis. It is given by cj' of equation(4-29).

For this simple refinery model, there were 33 variables whose optimal value were determined, and 38 constraint equations were satisfied. For an actual refinery there would be thousands of constraint equations, but they would be developed in the same fashion as described here. As can be seen, the model (constraint equations) was simple, and only one set of operating conditions was considered for the catalytic cracking unit, catalytic reformer, and the crude distillation column.

If the optimal flow rates do not match the corresponding values for volumetric yields, a search can be performed by repeating the problem to obtain a match of the optimal flow rates and volumetric yields. This has to be performed using a separate simulation program which generates volumetric yields from flow rate through the process units. (See Figure 6-5)   Thus, the linear model of the plant can be made to account for nonlinear process operations. Another procedure, successive linear programming, uses linear programming iteratively, also, and it will be discussed in Chapter 6. The state of industrial practice using both linear programming and successive linear programming is described by Smith and Bonner(13) for configuration of new refineries and chemical plants, plant expansions, economic evaluation of investment alternatives, assessment of new technology, operating plans for existing plants, variation in feeds, costing and distribution of products, evaluation of processing and exchange agreements, forecasting of industry trends, and economic impact of regulatory changes.

top