Apparatus for optimizing multiunit processing systems

Abstract

A processing system comprising a plurality of individual units is optimized by first estimating the yields of the individual units at a standard set of operating conditions and then establishing optimum flow rates using a linear programming model or similar mathematical tehniques. Individual units are then controlled and locally optimized consistent with a sensitivity analysis which is performed by treating a proposed change in operating conditions of an individual unit as a disturbance in the unit yield column of the linear programming model. The overall system may then be optimized for the changes in operating conditions by determining and establishing new flow rates. The steps of changing operating conditions and establishing new flow rates may be repeated until the sensitivity analysis reveals that any further change in operating conditions will not further improve the profit of the overall system.

Parent Case Text

RELATED APPLICATION

This is a continuation-in-part of copending application Ser. No. 246,445, filed Apr. 21, 1972, now abandoned which is incorporated herein by reference.

Description

BACKGROUND OF THE INVENTION

This invention relates to a method for optimizing large, complex processing systems comprising a plurality of individual units.

A number of techniques are available for optimization of such large processing systems. Linear programming is one of the most widely used techniques, and modern refinery complexes are optimized by linear programming almost without exception. However, a linear programming model is at best an approximation of a real physical system. For example, it is commonly assumed that the yield information incorporated into the LP (linear programming model) is relatively fixed and the coefficients of the constraint equations which represent these yields are also fixed. In reality, however, these yield coefficients are dependent upon the operating conditions of the units and any change in the operating conditions for a single unit will of necessity affect the operation of other units.

In many instances, the operating conditions may be closely controlled. Local unit managers have appropriate tools (i.e., mathematical models, optimizers, and process control computers) for the local optimization and control of their individual units. As a result, large amounts of detailed information about the processes are continuously generated which are valuable for accurate adjustment of operating conditions. Each process unit can thus be locally optimized and controlled continuously.

However, any change in operating conditions of a unit which locally optimizes the unit will not necessarily optimize the overall system. In fact, a change in operating conditions which optimizes the local unit may have an adverse effect on the overall system. As a result, changes in operating conditions of the individual units have been discouraged or avoided for fear of the adverse effect on the overall system.

Changes in the operating conditions of individual units have been avoided for another reason. Any attempt to change the operating conditions at an individual unit would necessarily render the LP and its solution obsolete since the yield coefficients in the LP for a particular unit would change. Such a change would therefore require the LP to be solved again and prior art computer techniques for solving LP's are extremely cumbersome and complex. In the case of a digital computer, a great deal of computer time is required to solve an LP, particularly where the LP described a complex processing system.

Because of the foregoing difficulties associated with making the changes in operating conditions in the prior art, individual local objective functions consistent with the overall objective function of the system have been assigned and every effort has been made to maintain those operating conditions at the individual units which will satisfy the local objective functions rather than change to better operating conditions. In other words, no effort is made to deliberately change the operating conditions at the units to optimize the overall system.

Sensitivity analyses such as that disclosed by C. S. Bightler and D. J. Wilde, Hydrocarbon Processing, 44, No. 2, 111 (1965) have been proposed to determine the effect of changes in LP constraints on the operation of a system. However, no specific method has been proposed to take advantage of the sensitivity analysis for changes in the yield coefficients of an LP which optimizes the overall system, i.e., making changes in the operating conditions which will optimize the overall system as indicated by the sensitivity analysis.

SUMMARY OF THE INVENTION

It is an object of this invention to provide an improved method of and apparatus for controlling the operation of a processing system including a plurality of individual processing units supplied by a plurality of feed streams of materials being transformed by the units into a plurality of product streams flowing from the units where the marginal product yield and the marginal product cost for each unit are dependent upon the operating conditions.

It is a more specific object of this invention to provide an improved method of and apparatus for controlling the operation of the processing system in a manner so as to encourage and implement changes in operating conditions for the individual units even though those changes do require a change in the linear program model of the processing system.

In accordance with these objects, a preferred embodiment of the invention comprises feed stream computer means for generating initial feed stream signals representing the initial feed stream flow rates for a given yield and cost at each of the units under a given set of operating conditions so as to maximize the profitability of the overall system. Feed stream control means are coupled to the feed stream computer means for controlling the feed stream flow rates in response to the initial feed stream signals. Unit yield and cost computer means generate new yield and cost signals representing new yield and cost for one of the units corresponding to an increase in the profitability of the overall system. Unit control means are coupled to the unit yield and cost computer means for controlling the operating conditions at the one unit in response to the new yield and cost signals. Another feed stream computer means is coupled to the unit yield and cost computer means for repeatedly and continuously generating new feed stream signals representing new feed stream flow rates in response to the new yield and cost signals. The other feed stream computer means is coupled to the feed stream control means for controlling the feed stream flow rates in response to the new feed stream signals.

It is also a specific object of this invention to provide an improved method of and apparatus for operating and utilizing a linear program model.

In accordance with this specific object, the feed stream computer means includes means for computing the feed stream flow rates x.sub.i from the linear program model having an overall system objective function of maximizing the profitability.

Z = c.sub.1 x.sub.1 + . . . . c.sub.n x.sub.n

subject to

p.sub.1 x.sub.1 + . . . . p.sub.n x.sub.n = Q

x.sub.i .gtoreq. O for all i

where p.sub.j = a.sub.ij , Q = b.sub.1 , . . . . . . . . . . . . . . . . . . . . a.sub.mj b.sub.m

c.sub.j is the marginal profit coefficient of the jth unit, x.sub.i is the flow rate of the ith stream, p.sub.j is the yield column of the jth unit and a function of operating conditions, Q is the demand constraint column, a.sub.ij is the yield coefficient of the ith feed stream producing the jth product stream unit, and b.sub.i is the demand coefficient of the ith product stream.

It is another specific object of this invention to assure that all changes in operating conditions at the units improve the overall profitability of the system.

In accordance with this specific object, the unit yield and cost computer means includes means for computing a local objective function

f.sub.m = (.DELTA.c.sub.m - .pi..degree..DELTA.p.sub.m)> 0

where .DELTA.c.sub.m is the change in the marginal profit coefficient for the mth unit, .DELTA.p.sub.m is the change in the mth yield column for the mth unit and .pi..degree. are the simplex multipliers [c.sub.i, . . . c.sub.n ].beta..degree. at the previous operating conditions where .beta..degree. is the inverse of the basis matrix in [p.sub.1, . . . p.sub.n ].

It is a further specific object of this invention to provide an improved method of and apparatus for correcting the linear program model after permitting a change in operating conditions at one of the operating units.

In accordacne with this specific object, the other feed stream computer means includes means for computing the feed stream flow rates ##EQU1## where

s.sub.i = .beta..degree..DELTA.p.sub.m. . . . . . s.sub.n

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram of a multiunit processing system which may be operating in accordance with this invention;

FIG. 2 is a block diagram of the method of this invention;

FIG. 2a is a block diagram showing a slight modification of the method shown in FIG. 2;

FIG. 3 is a block diagram of a modified method embodying the invention;

FIG. 4 is a simplified refinery system operated in accordance with the method of this invention;

FIG. 5 is a graphical solution of a local optimization problem;

FIG. 6 is a block diagram of a processing system operated in accordance with the method depicted in FIG. 2;

FIGS. 7 (a-d) are schematic circuit diagrams of a local computer and process control shown in block form in FIG. 6 where FIG. 7a shows circuitry for computing unit yield and cost coefficients, FIG. 7b shows circuitry for computing the local constraints and the local objective function, and FIGS. 7(c and d) shows circuitry for correcting the linear program model;

FIGS. 8 (a- c) are schematic circuit diagrams of the system computer and control shown in block form in FIG. 6 where FIG. 8a shows circuitry for computing the solution to the linear program model, FIG. 8b shows circuitry for determining the bases in the linear program model and FIG. 8c shows circuitry for controlling the feed stream flow rates to the units;

FIG. 9 is a flow diagram for a programmed digital computer which performs the same correction of the linear program performed by the circuitry of FIGS. 7(c-d);

FIG. 10 is a block diagram of a typical refinery system which is controlled in accordance with this invention;

FIG. 11 is a block diagram of a fluid catalytic cracking unit of the system shown in FIG. 6; and

FIG. 12 is a block diagram depicting the relation and interconnections between the circuit diagrams of FIGS. 7(a- d) and FIGS. 8(a-c).

TABLE OF CONTENTS FOR THE DETAILED DESCRIPTION

I. complex Processing System Employing Invention

Ii. method of Operating Complex Processing System

Iii. general Description of a Simple Processing System

Iv. detailed Description of the Simple Processing System Including Analog Computer Control

A. computing Yield and Cost Coefficients at the Local Computer

B. computing an LP Solution at the System Computer

C. computing the Bases of the LP Solution at the System Computer

D. setting the Flow Rate Controls

E. computing and Storing the Inverse Basis Matrix

F. computing the Simplex Multipliers

G. checking Local Optimization

H. checking Local Constraints

I. checking the Local Objective Function

j. Correcting the LP Solution

1. Computing Sensitivity Coefficients

2. Computing New X's

3. Computing New Basis Matrix

4. Computing New Simplex Multipliers

5. Computing New Cost Coefficients

V. control of the Simple Processing System with a Digital Computer

A. correction of the LP

B. computing the Local Objective Function

Vi. numerical Examples

Vii. method of Operating a Complex Refinery System.

Viii. modified Method of Operating a Complex System.

Ix. another Numerical Example.

DETAILED DESCRIPTION OF A PREFERRED EMBODIMENT

I. Complex Processing System Employing Invention

FIG. 1 is a schematic diagram of a complex processing system including a plurality of individual operating units 10.sub.1-7 which are interconnected by a plurality of streams of material 12.sub.1-15. The units 10.sub.1-7 have local computers 14.sub.1-7 which are utilized to determine operating conditions which optimize the operation of the associated local units 10.sub.1-7. An overall planning computer 16 which is in communication with the various computer 10.sub.1-7 through communications links 18.sub.1-7 is also provided. The computer 16 is utilized in planning the optimum operation of the entire system. The operation of the system of FIG. 1 in accordance with the method of this invention will now be discussed with reference to FIG. 2.

II. Method of Operating Complex Processing System

Initially, the yields at standard or design operating conditions and the corresponding costs for various units 14.sub.1-7 are estimated as indicated at block 20. These yields and costs may be estimated utilizing the local unit tools available to the local unit managers including mathematical models, optimizers and process control computers.

After the yields and costs for the individual units have been estimated, data representing the yields and costs is transmitted to the overall system planning computer 16 through data communication links 18.sub.1-7 as indicated at block 22. The overall system computer 16 solves the LP utilizing the transmitted yields as indicated at block 24, and the LP solution is then transmitted over the data links 18.sub.1-7 to each of the local units 10.sub.1-7. In addition, the local objective function based on the LP solution is also transmitted over the data links 18.sub.1-& as indicated by block 26.

In accordance with this invention, the local manager of each unit and/or his local computer is now in a position to perform a sensitivity analysis and thereby determine if a change in operating conditions at his particular unit which will result in optimization of his unit, e.g., improved product quality or quantity, will also satisfy the local objective function transmitted from the overall system computer 16. Using a proposed change indicated by his mathematical model or optimizer, he tests the local objective function on his local computer 14 as indicated by block 28. If the change satisfies the local objective function, the change in operating conditions may be made and the resulting change in yield is transmitted through the overall system computer 16 through the appropriate data link 18 as indicated at block 30.

In accordance with another important aspect of this invention, the overall system computer 16 may now correct the LP solution as indicated at block 32 without again solving the entire LP. The steps represented by blocks 26 and 28, 30 and 32 may now be repeated until any further change in yield of a unit will no longer satisfy the local objective function.

It will therefore be seen that the foregoing method advantageously permits local unit optimization which is consistent with overall planning of the system by providing communication between the local unit level and the overall system planning level as by data links 18 while permitting those changes in yield of individual units which will satisfy an overall objective function which is consistent with the optimization of the overall system. Thus in accordance with another aspect of the invention, the two levels of decision making, the overall system level and the individual unit level, are cordinated to establish overall optimum solution. The relationships between overall system optimization, local unit optimization, the overall system objective function and the local objective function will now be described in terms of a small and somewhat simple processing system.

In describing the invention as depicted in FIG. 2, a particular nomenclature will be utilized. Although this nomenclature is more or less standard in linear program modeling, a glossary providing definitions of the nomenclature will now be set forth for convenience and clarity.

glossary ______________________________________ b.sub.i = demand coefficient for ith stream B = base matrix b = .beta..degree. Q c.sub.j = marginal profit coefficient of jth unit c.sub.j = relative cost coefficient of jth nonbasic column d.sub.i = local decision vector of ith unit ##EQU2##

g.sub.m (p.sub.m, c.sub.m) = set of convex constraints for mth unit L (x,p' s, d.sub.i 's) = equality constraints of demand and supply P.sub.j = vector function of product yields of jth unit p.sub.j = yield column of jth unit p.sub.ij = yield coefficient for ith stream of jth unit Q = demand constraint column R.sub.m = convex region of feasible operation of mth unit s = .beta..degree..DELTA.p.sub.m = sensitivity coefficients T.sub.1, T.sub.2 = temperatures X.sub.i = flow rates for the ith stream Z = overall objective function for the system Greek .beta. = inverse of B .gamma. = row vector of profit coefficients of base columns .DELTA. = disturbances .pi. = simplex multipliers ______________________________________ Superscript .degree. represents original values Superscript 1 represents disturbed values Superscripts -, + represent the lower and the upper limits, respectively Underline represents a vector Underlines represent a matrix

Let x.sub.i be the ith stream (one of the streams 12 in FIG. 1) and d.sub.i is a decision vector of operating conditions for the ith unit (one of the units 10 in FIG. 1). Loosely speaking, the overall system objective is to maximize an overall profit function of the entire system Z; i.e.,

Maximize Z = Maximize Z (x, d.sub.1, . . . d.sub.N)

{x, d.sub.1, . . . d.sub.N } (1)

where N is the number of process units in the system, subject to:

(i) Yield equations

p.sub.i = p.sub.i (x, d.sub.i), for all i, (2)

where p.sub.i is a column vector representing the product yields of ith unit. P.sub.i is generally a vector of nonlinear functions which define product yields as functions of its operating conditions (d.sub.i) and feedstock (x.sub.i);

(ii) Non-negativity constraints

x.sub.i .gtoreq. 0 for all i; (3)

(iii) Demand-supply constraints

L (x, p.sub.i 's, d.sub.i 's) = b; (4)

Inequality constraints are assumed to be converted to equality constraints in the form of Equation (4) by adding slack variables; and

(iv) Limits on decision vectors

d.sub.i .sup.- .ltoreq. d.sub.i .ltoreq. d.sub.i .sup.+ for all i. (5)

Equations (1)-(5) define the problem symbolically. More specifically, we are concerned with the following special case.

______________________________________ N Maximize Z = Max .SIGMA. c.sub.i (d.sub.i) x.sub.i (6) i=1 ______________________________________

{x, d.sub.i, . . . d.sub.N }

where c .sub.i is the marginal cost or profit coefficient for the ith unit subject to

p (d.sub.1) x.sub.1 + p.sub.2 (d.sub.2) x.sub.2 + . . . + p.sub.N (d.sub.N) x.sub.N = Q (7)

and Equations (2) and (3) where p.sub.j is the yield column for the jth unit and Q is the demand constraint column. Notice that if a standard set of operating conditions are assumed, d.sub.1 .degree., d.sub.2 .degree., . . . , d.sub.N .degree., the above problem just reduces to a regular LP problem for determining an optimal schedule for the production flow rates (X*). As one of the d.sub.i 's changes, the yields from ith unit change, this in turn requires changes in another part of the system either in x.sub.i 's or d.sub.i 's. Equations (2), (6) and (7) are generally highly nonlinear and complex. If the dimensionality of the problem is reasonable, an optimum seeking method may be used to directly search the optimum. However, for large systems such as refineries, thousands of equations are commonly involved, precluding direct search methods, and the simplex type of computational algorithms must be relied upon. The Simplex method and simplex algorithms are described in Linear Programming and Extensions, by George B. Dantzig, 1963, pages 94-100 which is incorporated herein by reference.

Assuming a set of standard operating conditions (designated by adding superscript o), the problems of overall system optimization becomes

N Maximize Z = .SIGMA. c.sub.i .degree.x.sub.i (8) i=1

subject to

p.sub.1 .degree. x.sub.1 + p.sub.2 .degree. x.sub.2 +. . . + p.sub.N .degree. x.sub.n = Q (9)

and x .gtoreq. o

where c.sub.i .degree.'s and p.sub.i .degree.'s are calculated at the standard operating conditions, d.sub.i .degree.'s. Q represents the demands of the products, and is given by ##EQU3## b.sub.i equals the demand coefficient of the ith stream. This is nothing but a regular LP problem. The optimum operating conditions of individual units d.sub.i *'s), i.e., for a given schedule of production (x = x.degree.), are

Maximize f.sub.m (d.sub.m, x.sub.m .degree.) with respect to d.sub.m (10)

subject to yield equations

p.sub.m = p.sub.m (x.sub.m .degree., d.sub.m) (11)

where max f.sub.m is the local objective function for the mth unit, and defined as

f.sub.m = c.sub.m - .pi..degree.p.sub.m (x.sub.m .degree., d.sub.m) (12)

Here, .pi..degree. represents the simplex multipliers of the basic feasible solution x.degree..

The optimization of the overall system in accordance with Equations (6) and (7), the optimization of the local units in accordance with Equations (10) and (11) and the use of local objective function (12) in accordance with this invention will now be discussed in detail mathematically.

As indicated at block 20, the yields and costs of the various units must be estimated at standard operating conditions. This is accomplished through the use of a mathematical model of each unit usually in the form of a computer program which predicts the yield product quantity and quality as function of standard operating conditions. The estimation which may be performed on the local computers 14 generates the yield column ##EQU4## in Equation (9) where a.sub.ij equals the yield coefficient of the ith stream for the jth unit.

After this initial step, the values in the columns p.sub.j .degree. and the value for each c.sub.i .degree. are transmitted to the overall system planning computer 16. With the transmission of the values for the yield columns p.sub.j .degree. and the c.sub.i .degree. as indicated at block 22, Z = c.sub.l .degree. x.sub.l +. . . c.sub.n .degree. x.sub.n subject to p.sub.l .degree. x.sub.l +. . . + p.sub.n .degree. x.sub.n = Q. Since Q, all p.sub.j .degree. and all c.sub.j .degree. are known, flow rates x.sub.i may be determined to maximize Z for the standard operating conditions as indicated at block 24. The solution to the LP along with the local objective function is now transmitted to each unit as indicated at block 26 where the sensitivity analysis is performed.

For simplicity, only one of the p.sub.j columns is assumed free to be chosen subject to unit constraints. This will be called the mth unit. For the time being, it will be assumed that a local objective function f.sub.m is known and the local optimum solution is obtained from the mathematical model of that unit and the sensitivity analysis so as to be consistent with the objective function f.sub.m. This local optimization has generated new operating conditions resulting in a change in product yield coefficient p.sub.m and a new cost coefficient c.sub.m. To show that the overall objective function actually increases, the planning problem is analyzed as follows.

Since the effect of local optimization on the planning model appears as some perturbation in the mth yield column p.sub.m, Equation (9) becomes

p.sub.l .degree. x.sub.l + . . . + p.sub.m .degree. x.sub.m + . . . + p.sub.n .degree. x.sub.n = Q - .DELTA.p.sub.m x.sub.m (13)

and the overall system objective function Z becomes

Z = c.sub.l .degree. x.sub.l + . . . + c.sub.m .degree. x.sub.m + . . . c.sub.n .degree.n + .DELTA.c.sub.m x.sub.m (14)

where

.DELTA.p.sub.m = p.sub.m - p.sub.m .degree.

.DELTA.c.sub.m = c.sub.m - c.sub.m .degree..

The mth column could be either basic or nonbasic, where a basic column has positive relative cost coefficients and a nonbasic column has negative relative cost coefficients. (A basis is any collection of the variables which contain exactly as many variables as there are constraints. See also the aforesaid book by Dantzig, page 81.) If nonbasic originally, the mth column can come into the bases in the first iteration only if its relative cost factor becomes positive. Otherwise, the problem is already solved. Since the object is to improve the standard operation of the overall system by simultaneously considering the local strategies and the overall system problem, the units of concern should generally be in the bases. If the first m columns are the bases for the original optimal solution, Equations (13) and (14) may be rewritten as

1 x.sub.1 + 0 x.sub.2 +...+ 0 x.sub.m + .beta..degree.p.sub.m.sub.+1 x.sub.m.sub.+1 +...+ .beta..degree.p.sub.n x.sub.n = .beta..degree.b - .beta..degree..DELTA..sub. 1 x.sub.m 0 1 . . . . . . (15) . . . 0 0 1

ox.sub.1 + ox.sub.2 + . . . + ox.sub.m + c.sub.m.sub.+ 1 x.sub.m.sub.+ 1 + . . . + c.sub.n x.sub.n = z-z.degree. + (.pi..degree..gamma.P.sub.m -c.sub.m)x.sub.m (16)

where .beta..degree. is the inverse of the basis matrix .beta..degree.,

.uparw. .uparw. .uparw. B.degree.= p.sub.1 .degree. p.sub.2 .degree. ... p.sub.m .degree. (17) .dwnarw. .dwnarw. .dwnarw.

and .pi..degree. = .gamma..degree..beta..degree. where .gamma..degree. = [c.sub.l, . . . c.sub.n ]. (17')

A sensitivity analysis may now be performed which allows the local manager to determine the effect of the change in .DELTA.p.sub.m and .DELTA.c.sub.m without solving the LP again. From Equation (14), the objective function z is, if x.sub.m.sub.+ 1, . . . , x.sub.n equal to zero.

Z = Z.sub.o + (.DELTA.c.sub.m - .pi..degree. .DELTA.p.sub.m) x.sub.m (18)

Equation (18) is verbally interpreted as follows: If the quantity in the parentheses becomes positive, the overall objective function increases, since x.sub.m must always be positive. Notice that, for given x.degree., only the last term of Equation (18) is dependent upon the local optimization. x.sub.m in Equation (18) can be represented in terms of the original value x.sub.m .degree. and the sensitivity coefficient s.sub.m, and hence, the local objective function f.sub.m is given by ##EQU5## It is important to realize that if f.sub.m can be made positive by changing d.sub.m , the overall objective function can be increased further. No deliberate changes in local operating conditions should be made, unless the local objective function defined by Equation (19) is rendered positive.

The local manager's next question concerns how far he can move his unit in the positive direction of his objective function. If his local optimization is subject only to yield constraints of Equation (11), then he can be allowed to move until a maximum value of f.sub.m is reached. This point should be an extreme corner of the space defined by yield equations, p.sub.m = P.sub.m (x.sub.m .degree., d.sub.m). In reality, however, if the local changes are too large, the overall scheduling problem may easily become unfeasible (or degenerate) at the new conditions. Therefore the local manager can only be allowed to move within the valid limits of the overall planning problem. Additional constraints from the overall problem thus must be satisfied always at the local level. To assure that the overall solution does not become degenerate, the constraints given the following Equations are required for the local optimization of the mth unit: ##EQU6## By multiplying the mth row by s.sub.i /1 + s.sub.m and substracting it from the ith row of Equation (15), it is easily shown that the new solution x.sup.1 is simply, ##EQU7## where b.sub.i = the ith element of [.beta..degree.Q]. The new basis matrix, B.sup.1, with the perturbed column is

.uparw. .uparw. .uparw. .uparw. .uparw. B.sup.1 = B ... (p.sub.m + .DELTA.p.sub.m) = B.degree. + 0 ... 0 .DELTA.p.sub.m (24) .dwnarw. .dwnarw. .dwnarw. .dwnarw. .dwnarw.

The inverse of B.sup.1 is obtained from Equations (22) and (23) after considerable algebraic manipulation. ##EQU8## The last term of Equation (25), denoted by .DELTA..beta., is an m .times. m matrix where the ith row is the mth row of .beta..degree. multiplied by a scalar factor + s.sub.i /(1+s.sub.m). In other words, the new inverse of basis matrix is the sum of the original inverse matrix and .DELTA..beta. defined above.

The changes in simplex multipliers are also obtained similarly. The new simplex multipliers .pi..sup.1 is, by definition,

.pi..sup.1 = .gamma..sup.1 .beta..sup.1 = (C.sub.1, . . . , C.sub.m + .DELTA.C.sub.m) .beta..sup.1 = [.gamma..degree. + (0, . . . , 0,.DELTA.C.sub.m)] (.beta..degree. + .DELTA..beta.) (26)

Substitutions of Equations (17') and (25) into Equation (26) result in the following Equation: ##EQU9##

Changes in simplex multipliers are given by the mth row of .beta..degree. multiplied by a scalar constant. It is important to notice that the numerator of f.sub.m is the same as the quantity in the parenthesis in Equation (18).

Now that .pi..sup.1 is known, the relative cost coefficients are readily obtained from the definition.

C.sub.j = C.sub.j - .pi..sup.1 .sup.. p.sub.j = C.sub.j - (.pi..degree.+.DELTA..pi.) p.sub.j. (29) = C.sub.j .degree. - .DELTA..pi. (30) b.j.

Changes in the overall solution will occur only in the rows for which p.sub.m has non-zero elements. Therefore, although the total number of rows of the planning model may reach to several hundred, the number of the additional constraints for local optimization defined by Equation (21) is generally reasonable, so as to be included into the local optimization packages. Equations (11), (19), (21) and (22) define the local optimization problem.

From the results of the local optimization, only the yield column p.sub.m.sup.l and the cost coefficient c.sub.m.sup.l are transmitted to the overall system computer 16 as indicated at block 28. These replace the old values p.sub.m .degree. and c.sub.m .degree.. The overall objective function increases monotonically, since we have already assured that f.sub.m is positive at each interation. The overall solution is then adjusted as indicated at block 32 to respond to the local changes just implemented.

Tests are now made to see whether the corrected solution is optimal: If any of the relative cost coefficients c.sub.j of the nonbasic columns became positive in the last iteration, Z can be further increased by forcing this column into the bases. This involves changes in the bases. Otherwise, the system is at an optimum with the latest yield information. This procedure is iterated until f.sub.m cannot be made positive over the entire range of d.sub.m.

The object has been to improve the optimal overall solution (x.degree., Z.degree.) obtained based on the standard yields and the standard operating conditions by considering the effect of a local optimization on the overall problem. The overall objective function (Z) can be improved only in two ways: First, when the local objective function (f.sub.m) becomes positive, Z will be increased by adjusting the original solution (x.degree.) to accommodate the local changes. In this portion of the improvement, no bases change is required. Moreover, should this adjustment render sign changes in at least one relative cost coefficients of the nonbasic columns, Z can be improved move by forcing one of the nonbasic columns with positive c.sub.j into the base. Second, although f.sub.m cannot be positive, it is possible to increase Z if at least one of the nonbasic c.sub.j becomes positive at some value of f.sub.m (even if negative). It is not, however, guaranteed that Z increases for the second case. Since f.sub.m is negative, the adjustment in x to accommodate the local changes will first decrease Z by f.sub.m x.sub.m .degree., and then Z will be increased when the nonbasic column with a positive c.sub.j enters into the base. Therefore, the net change in Z could be either positive or negative.

The foregoing iterative procedure has neglected the second possibility described above. To make the procedure complete, it is preferred that the range of negative f.sub.m be examined at least once during the entire procedure, preferably in the first iteration. Fortunately, it can be easily proved that if no sign changes in c.sub.j of nonbasic columns occur at both the maximum and the minimum values of f.sub.m, then no such change is possible over the entire feasible range of f.sub.m. Therefore, it is sufficient to test the optimality of the overall solution at only these two values.

The foregoing method is particularly well suited for the optimization of large refinery complexes. Almost without exception, the operation of refineries is planned by linear programming. Many of the constraint equations are reasonably linear and yields from some units are almost fixed except for a few key process units, such as catalytic cracking, hydrocracking, and catalytic reforming units. Detailed mathematical models for these key units will become the main tool for the local optimization and control of individual units. The key process units of the refinery may be individually optimized and controlled continuously with small process control computers at the plant site, and the overall planning program is solved once every day (or more often if required) with a large computer located at the computing center of the refinery complex. Shortterm scheduling for production targets is derived from the overall planning program. To meet this target under varying situations in the refinery, day-to-day operation should be adjusted frequently, if not continuously.

The local managers must have detailed information about their individual units for local optimization and control. This should include, for instance, unit operating conditions, the state of the catalyst, performance of auxiliary equipment, and the like. In addition to these, the above described method allows the local managers to consider the inventory problem (between processes) as a local decision. They can thus adjust their production rate based on the tankage availability of their feedstock and products. Most of the process information remains in the process level and only a fraction of it passes to the planning level as iteration proceeds.

The method of this invention takes advantage of the current practice as much as possible. It allows use of the current LP model for the overall planning and the detailed process models for the local optimization tools. It provides the communication link between these two levels for adjusting the overall planning as well as the local strategies to keep the refinery operation always close to the optimum.

III. Description of a Simple Processing System

As shown in FIG. 6, a simple processing system comprises four processing units 110.sub.1-7 which are respectively supplied by four raw material feed streams 112.sub.1-4 to produce two product streams 112.sub.5 and 112.sub.6. As in the complex case shown in FIG. 1, local computers 114.sub.1-4 are associated with the respective operating units 110.sub.1-4 to control the operating conditions of these units. An overall system computer 116 is in communication with the local computers 114.sub.1-4 through data communication links 118A.sub.1-4 and 118B.sub.1-4. The data communication links 118A are utilized to transmit signals representing unit operating conditions to the system computer 116 while data links 118B are utilized to transmit signals representing the local objective functions from the system computer 116 to the local units 114.sub.1-4. In addition, control links 120.sub.1-4 are provided between the system computer 116 and valve 122.sub.1-4 for controlling the flow of raw materials through the feed streams 112.sub.1-4.

The units 110.sub.1-4 may comprise any of a number of processing units such as those found in a refinery, e.g., a catalytic cracker, a hydrocracker, a catalytic reformer, a gas plant, etc. For purposes of this embodiment of the invention, the units 110.sub.1-4 will be considered as processing units wherein temperature is a specific operating condition sensed by sensors 124.sub.1-4, which temperature is a function of the setting of a valve 126.sub.1-4 as controlled by the local computers 114.sub.1-4. The individual operating conditions, i.e., the operating temperature, in turn controls or affects the individual yields of the individual units 110.sub.1-4 as well as the profit or cost of the units. It will of course be understood that the other operating conditions may be sensed including the quality of the products produced by the units.

In the system comprising the units 110.sub.1-4 as shown in FIG. 6 streams 112.sub.1-4 and the product streams 112.sub.5 and 112.sub.6 may be defined in terms of a LP model having the overall system objective of maximizing the profit.

Z = C.sub.l x.sub.1 + c.sub.2 x.sub.2 + c.sub.3 x.sub.3 + c.sub.4 x.sub.4 (31)

subject to ##EQU10## where the c's are the cost coefficients for the four operating units, the p's are the yield coefficient for the two product streams of the four units and the b's are the demand constraints for the two feed streams. The LP model also includes a certain constraints which require that the x's of all feed streams be greater than or equal to zero. In addition, the individual units may be subjected to local constraints which require that, for instance, for Unit 4,

p.sub.14 0, p.sub.24 0, c.sub.4 0,

and p.sub.14 + 2 p.sub.24 + 3 c.sub.4 = 2.

The system of FIG. 6 including the local computers is operated in accordance with the method described in the block diagram of FIG. 2a which differs only slightly from that shown in FIG. 2. In the block diagram of FIG. 2, the correction of the LP is performed at the system planning level on the system computer as indicated at block 32. In the block diagram of FIG. 2a, the correction of the LP solution is performed at the local computer as indicated at block 30a and the corrected LP is transmitted to the system computer as indicated at block 32a.

IV. Detailed Description of the Simple Processing System Including Analog Computer Control

The simple processing system of FIG. 6 under analog computer control will now be described. Because of the complexity of the circuitry involved, it may be helpful from time to time to refer to the drawing of FIG. 12 which illustrates the relationship of the circuits shown in FIGS. 7(a-d) and 8(a-c) with interconnections identified by the appropriate reference characters.

A. computing Yield and Cost Coefficients At the Local Computer

FIG. 7a, which represents a portion of the analog computer circuitry of the local computer 114.sub.4, includes circuitry for estimating the cost and yield coefficients for the unit 110.sub.4 at standard operating conditions as depicted by block 20 in FIG. 2a. As shown in FIG. 7a, the estimating circuitry includes a differential amplifier 150 for comparing a signal sensed by the sensor 124.sub.4, which may comprise a thermocouple, and a setpoint signal taken from a movable contact 152 associated with a slide wire 154 supplied by a suitable DC source designated as +V. With the sliding contact 152 set at a position corresponding to standard operating or design conditions for the unit 110.sub.4, the output from the differential amplifier 150 will drive a motor 156 which in turn controls the position of the valve 126.sub.4 through a mechanical linkage 158 until the temperature as measured by the sensor 124.sub.4 equals the setpoint determined by the position of the sliding contact 152.

In order to estimate the yield coefficient p.sub.- .degree. and the cost coefficient c.sub.4 .degree. for standard operating conditions, the temperature signal from the sensor 124.sub.4 is applied through resistors 162, 164 and 166, designated respectively as kp.sub.14 .degree., kp.sub.24 .degree. and kc.sub.4 .degree.. These resistors 162, 164 and 166 are appropriately chosen such that the temperature signal applied therethrough will produce signals representing the yield coefficients p.sub.14 .degree. and p.sub.24 .degree. and the cost coefficient c.sub.4 .degree.. Thus the resistors 162, 164 and 166 actually represent a model of the unit 110.sub.4 relating the standard operating temperature of the unit to the yield and cost coefficients.

Although not shown in any drawings, similar yield and cost coefficient estimating circuitry is provided in the other local computers 114.sub.1-3. Since the units 110.sub.1-3 are similar to the unit 110.sub.4, the circuitry may be identical except for the choice of the resistors corresponding to resistors 162, 164, and 166 which, of course, must represent a model of the relationship between the standard operating temperature and the yield and cost coefficients for the other units.

After the yield and cost coefficients for all of the units 110.sub.1-4 have been estimated, signals representing the yield coefficients and the cost coefficients are transmitted to the central computer 116 as depicted by block 22 of FIG. 2. As shown in FIG. 7a, the signals p.sub.14 .degree., p.sub.24 .degree. and c.sub.4 .degree. are transmitted to the system computer 116 through lines 168, 170 and 172. Similar lines are provided from the other local computers 114.sub.1-3 to the central computer 116.

B. computing an LP Solution at the System Computer

Once the yield coefficient and cost coefficient signals have been transmitted to the system computer 116, the system computer 116 may solve the LP defined by the equations (31) and (32) as depicted by block 24 of FIG. 2a. The specific circuitry at the system computer 116 which is capable of solving the LP will now be described with reference to FIG. 8a.

As shown in FIG. 8a, the system for solving the LP comprises a number of variable resistors 174 which have resistance values which represent and are proportional to the various yield coefficients of the various units. In order to solve the LP for standard operating conditions, it is necessary that the variable resistors 174 be set so as to correspond with the yield coefficients which are transmitted over lines 168 and 170 from the local computer 114.sub.4 and similar lines from other local computers 110.sub.1-3. For purposes of illustration, the manner in which the resistors 174 are set so as to represent the yield coefficient p.sub.14 .degree. will now be described.

A switch 176 in the line 168 as shown in FIG. 8a is closed when control of the overall system is initiated, either manually or by automatic means not shown, so as to apply the signal p.sub.14 .degree. to the non-inverting terminal of a differential amplifier 178 while a signal from a movable contact 180 associated with a slidewire or variable resistor 182 connected across a DC power supply +V is applied to the inverting terminal. The output of the differential amplifier 178 is connected to a motor 184 which is in turn mechanically coupled to the slidewire contact 180 through a mechanical linkage 186 so as to move the contact 180 until the output of the differential amplifier 178 is nulled. Another mechanical linkage 188 is connected from the motor 184 to the variable resistor 174 corresponding with the yield coefficient p.sub.14 .degree. so as to control the magnitude of the resistance such that it equals or is at least proportional to the input signal p.sub.14 .degree. to the differential amplifier 178. Similar circuitry bearing the same reference characters is associated with the line 170 from the local computer 114.sub.4 for controlling the magnitude of the resistor 174 representing the yield coefficient p.sub.24 .degree.. It will of course be appreciated that the other resistors 174 corresponding to the other yield coefficients for the other units 110.sub.1-3 will be set in an identical manner by circuitry not shown.

The circuitry for solving the LP also includes a number of variable resistors 190 which represent or are proportional to the cost coefficients for the various units 110.sub.1-4. In the case of the resistor 190 representing the cost coefficient c.sub.4 .degree., the magnitude of the resistor will be set in response to the signal c.sub.4 .degree. applied over the line 172 from the local computer 110.sub.4 to the system computer 116 utilizing circuitry which is substantially identical to that for setting the yield coefficient including the differential amplifiers 178 and 184. The magnitude of the resistors for the cost coefficients of the other units 110.sub.1-3 will be set in an identical manner.

Having now set the resistances in the circuitry for solving the LP in accordance with the estimated yield and cost coefficients transmitted from the local units 110.sub.1-4, the circuitry shown in FIG. 8a is now ready to solve the LP. Generally speaking, the circuitry of FIG. 8a, represents well known analog circuitry for solving an LP. For example, substantially identical circuitry is disclosed in a book entitled "Analog Computation," Albert S. Jackson, McGraw Hill Book Company, Inc., 1960, pages 353-363. However, for the sake of completeness of disclosure, the circuitry in FIG. 8a will now be described in considerable detail.

The LP solution circuitry includes a resistance 192 having a magnitude equal to .epsilon. which represents the speed of movement of the solution point for the LP through various coordinates of four-dimensional space defined by the variables x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree. which are the feed stream flow rates to the various units 110.sub.1-4. The resistor 192 is located in a line 194 connecting a +1 volt DC power supply to each of the cost coefficient resistors 190. The output of the resistors 190 are then applied to integrating amplifiers 196.sub.1-4. The other inputs to the integrating amplifiers 196.sub.1-4 are provided by feedback signals which are operated upon by yield coefficient resistors 174, which feedback signals will be described subsequently. Signals representing the value for the feed stream flow rates x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree. which maximize the overall objective function Z, are generated at the output of the integrating amplifiers 196.sub.1-4 respectively.

In order to generate the previously referred to feedback signals, the various feed stream flow rate signals x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree. are applied through various yield coefficient resistors 174 to two summing amplifiers 198.sub.1 and 198.sub.2. In particular, the signals x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree. are applied, through resistors 174 representing the yield coefficients p.sub.11 .degree., p.sub.12 .degree., p.sub.13 .degree. and p.sub.14 .degree. respectively, to the four different inputs of the summing amplifier 198.sub.1. In addition, a -1 volt signal is applied through a resistor 200.sub.1 representing the demand constraint b.sub.1 .degree. to another input of the summing amplifier 198.sub.1. The sum signal at the output of the amplifier 198.sub.1 is applied to a multiplying amplifier 202.sub.1 which generates an output signal equal in magnitude to the sum signal multiplied by any arbitrary, large, real number K. The output signal from the multiplying amplifier 202.sub.1 is then applied as one of the feedback signals to the yield coefficient resistors 174 representing the yield coefficients p.sub.11 .degree., p.sub.12 .degree., p.sub.13 .degree. and p.sub.14 .degree. at the input to the integrating amplifiers 196.sub.1-4 respectively.

Similarly, the feed stream flow rate signals x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree. may be applied to the inputs of the summing amplifier 198.sub.2 through yield coefficient resistors 174 representing p.sub.21 .degree., p.sub.22 .degree., p.sub.23 .degree. and p.sub.24 .degree.. A -1 volt signal is also applied through resistor 200.sub.2 representing the demand constraint b.sub.2 .degree., to another input of the summing amplifier 198.sub.2, and the output from the summing amplifier 198.sub.2 is applied to a multiplying amplifier 202.sub.2 having a multiplication factor of K. The output from the multiplying amplifier 202.sub.2 is then applied as a feedback signal to the resistors 174 at the input of the integrating amplifiers 196.sub.1-4 where the resistors 174 represent the yield coefficients p.sub.21 .degree., p.sub.22 .degree., p.sub.23 .degree. and p.sub.24 .degree. respectively. Note that the integrating amplifiers 196.sub.1 and the summing amplifiers 198.sub.1 include diode limiting circuits 204 and 206.

In order to solve for the X.degree., which represents the maximum Z at standard operating conditions for the units 110.sub.1-4, a summing amplifier 208 is provided. Signals representing x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree. are applied to the inputs of the summing amplifier 208 through resistances 190 representing the cost coefficient c.sub.1 .degree., c.sub.2 .degree., c.sub.3 .degree. and c.sub.4 .degree. respectively.

By way of general explanation, the above described circuitry of FIG. 8a solves the LP by searching out the point in four-dimensional space which represents the maximum value of Z at the standard set of operating conditions subject to particular demand constraints. As mentioned previously, .epsilon. represents the movement of the solution point through four-dimensional space where the four dimensions are x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree.. In the circuit shown in FIG. 7a, the rate of movement .epsilon. of the solution point along the x.sub.1 .degree. direction, the x.sub.2 .degree. direction, the x.sub.3 .degree. direction and the x.sub.4 .degree. direction have been equalized. Thus the signal applied to the resistors 190 representing the cost coefficients c.sub.1 .degree., c.sub.2 .degree., c.sub.3 .degree. and c.sub.4 .degree. at the inputs to the integrating amplifiers 196.sub.1-4 represent the rates of change with respect to the time of the variables x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree. or x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree.. The quantities x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree., after multiplication by the cost coefficients c.sub.1 .degree., c.sub.2 .degree., c.sub.3 .degree., and c.sub.4 .degree. are integrated with respect to time by the integrating amplifiers 196.sub.4-1. By choosing a +1 volt source connected to the resistor 192 having a value .epsilon. , the output signal from the resistor 192 is positive such that the solution point is constrained to move in a direction which increases the value of Z.degree.. In order to assure that the solution point is constrained to move within boundaries determined by the demand constraints b.sub.1 .degree. and b.sub.2 .degree., the feedback signals from the multiplying amplifiers 202.sub.1 and 202.sub.2 are utilized to cancel out any velocity component normal to the boundaries, while permitting movement of the solution point along lines parallel to the boundaries. As long as the signal at the output of the resistor 192 is positive, the solution point will always have a component of movement along the gradient of Z.degree. until the maximum of Z.degree. is reached. When the maximum of Z.degree. is reached, the last component of movement of the solution point in the direction of the gradient of Z.degree. will have been cancelled by the feedback signals to the integrating amplifiers 196.sub.1-4 and the solution point will come to rest. The diodes 204 of the diode limiting circuits are necessary to assure that the values of the variables x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree. always remain positive. Diodes 206 of the amplifiers 198.sub.1 and 198.sub.2 are required to assure that the outputs from the summing amplifiers 198.sub.1 and 198.sub.2 are always negative, i.e.,

b.sub.1 .degree. .gtoreq. p.sub.11 .degree.x.sub.1 .degree. + p.sub.12 .degree.x.sub.2 .degree. + p.sub.13 .degree.x.sub.3 .degree. + p.sub.14 .degree.x.sub.4 .degree.

b.sub.2 .degree. .gtoreq. p.sub.21 .degree.x.sub.1 .degree. + p.sub.22 .degree.x.sub.2 .degree. + p.sub.23 .degree.x.sub.3 .degree. + p.sub.24 .degree.x.sub.4 .degree.

C. computing the Bases of the LP Solution at the System computer

Signals representing the solution of the LP, that is signals representing the value of x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree., are now applied over lines 210 to circuitry shown in FIG. 8b for determining which of the feed stream flow rates are the bases, i.e., which of the feed stream flow rates have negative cost coefficients associated therewith. This circuitry, which is located at the system computer, is responsive to cost coefficient signals c.sub.1 .degree., c.sub.2 .degree., c.sub.3 .degree. and c.sub.4 .degree. which are transmitted from the various local units 110.sub.1-4 to the system computer 116 over lines 172.

Two checking circuits 214 determine which of the various cost coefficients c.sub.1 .degree., c.sub.2 .degree., c.sub.3 .degree. and c.sub.4 .degree. are greater than or equal to zero. As a result of the check performed by the circuits 214, various switches 216.sub.1-4, which are associated with the various checking circuits 214, are positioned so as to place the movable contact a in contact with stationary contact b or stationary contact c. If the cost coefficient checked by one of the circuits 214 is greater than or equal to zero, the movable contact a, which is moved by a mechanical linkage 218 associated with a relay of the checking circuit 214 will move so as to be in contact with stationary contact b. If the cost coefficient is less than zero, the movable contact a will be moved or remain in contact with the stationary contact c. Additional switches 220.sub.1 and 220.sub.2 having a movable contact d and stationary contacts e and f are actuated by the linkage 218. In switches 220.sub.1 and 220.sub.2, the movable contact d is in contact with the stationary contact f if the cost coefficient c.sub.1 .degree. and c.sub.3 .degree. are greater than or equal to zero. If these cost coefficients are less than zero, movable contact d is in contact with stationary contact e. Additional switches 222.sub.1 and 222.sub.2 are provided having movable contacts g and stationary contacts h and i. When the cost coefficient c.sub.1 .degree. and c.sub.3 .degree. are less than zero, the linkage 218 associated with the c.sub.1 .degree. and c.sub.3 .degree. checking circuits 214 will place the movable contacts g in contact with the stationary contacts i. If the cost coefficient c.sub.1 .degree. and c.sub.3 .degree. are greater than or equal to zero, the movable contact g will be placed in contact with the stationary contacts h.

In the particular LP which has just been solved by the circuitry shown in FIG. 8a, the cost coefficients c.sub.1 .degree. and c.sub.3 .degree. have been assumed to be less than zero and the cost coefficients c.sub.2 .degree. and c.sub.4 .degree. have been assumed to be greater than zero. As a result, the movable contacts a, d and g of the switches 216.sub.1-4, 220.sub.182 and 222.sub.182 respectively have assumed the position shown in FIG. 8b. As a result, the solutions to the LP x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree. which are applied over lines 224 to the movable contacts a of the various switches 216.sub.1-4 are sorted out such that the signals representing the bases x.sub.4 .degree. and x.sub.2 .degree. appear on output lines 226 and 228 and signals representing the non-bases x.sub.1 .degree. and x.sub.3 .degree. appear on output lines 230 and 232. Had the polarity of the cost coefficients c.sub.1 .degree., c.sub.2 .degree., c.sub.3 .degree. and c.sub.4 .degree. differed from that assumed, the signals appearing on the output lines 226, 228, 230 and 232 would of course differ since the bases and non-bases would differ accordingly.

Similar switching circuitry 300 is provided for sorting out the basic cost and yield coefficient signals so as to obtain basic yield coefficient signals p.sub.14 .degree., p.sub.12 .degree., p.sub.24 .degree. and p.sub.22 .degree. on lines 302, non-basic yield coefficient signals p.sub.11 .degree., p.sub.13 .degree., p.sub.21 .degree. and p.sub.23 .degree. on lines 304, basic cost coefficient signals c.sub.4 .degree. and c.sub.2 .degree. on lines 306 and non-basic cost coefficient signals on lines 308. The entire LP solution is then transmitted over lines 226, 228, 230, 232, 302, 304, 306 and 308 to the local computer 114.sub.1-4 as indicated in block 26 of FIG. 2

D. setting the Flow Rate Controls

Now that the LP has been solved, the feed stream flow rate x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 .degree. may be established by appropriately setting the valves 112.sub.1-4 shown in FIG. 8c. In order to set the valves 112.sub.1-4, signals x.sub.1 .degree., x.sub.2 .degree., x.sub.3 .degree. and x.sub.4 20 are applied to the local computer 114.sub.4 as shown in FIG. 7a over output lines 226, 228, 230 and 232.

The setting of the valves 122.sub.4 and 122.sub.2 is done in response to the signals x.sub.4 .degree. and x.sub.2 .degree. on lines 226 and 228 having switches 234 and 236 therein as shown in FIG. 7c. The switches 234 and 236 include a movable contact a-1 and stationary contacts b-1 which are connected to lines 226 and 228 and stationary contacts c-1. The movable contacts a-1 are momentarily placed in contact with the stationary contacts c-1 by a linkage not shown when the switch 176 in the LP solution circuitry of FIG. 8a is closed. This results in the charging of capacitors 238 and 240 of FIG. 7c so as to store signals representing the flow rates x.sub.4 20 and x.sub.2 .degree.. The x.sub.4 .degree. and x.sub.2 .degree. signals are then transmitted back to system valve control circuitry shown in FIG. 8c over lines 239 and 241.

The lines 239 and 241 are connected to the inverting input terminals of differential amplifiers 242 and 244 of the valve control circuitry shown in FIG. 8c. The inverting terminals of the differential amplifiers 242 and 244 are connected to the sliding contacts 246 and 248 of slidewires or potentiometers 250 and 252 which are connected between a DC power supply +V and ground.

The output from the differential amplifiers 242 and 244 are utilized to drive valve control motors 256 and 258 which are mechanically linked to the slidewire contacts 246 and 248 as well as the feed stream flow rate control valves 122.sub.4 and 122.sub.2 respectively. Thus the motors 256 and 258 reposition the valves 122.sub.4 and 122.sub.2 until the output of the differential amplifiers 242 and 244 are nulled when the voltages on the slidewire contacts 246 and 248 equal the voltage across the capacitors 238 and 240 which represent or are proportional to the feed stream flow rates x.sub.4 .degree. and x.sub.2 .degree..

Similar motor control circuitry is provided for positioning the valves 122.sub.1 and 122.sub.3. This circuitry includes switches 260 and 262 in the lines 230 and 232 as shown in FIG. 8c which are actuated by and linked with the switches 176 in the LP solution circuitry of FIG. 8a so as to apply the feed stream flow rate signals x.sub.1 .degree. and x.sub.3 .degree. to capacitors 264 and 266 to the non-inverting terminals of differential amplifiers 268 and 270 and ground. The inverting terminals of the differential amplifiers 268 and 270 are connected to sliding contacts 272 and 274 of slidewires 276 and 278. Motors 280 and 282 which are connected to the outputs of the amplifiers 268 and 270 reposition the valves 122.sub.1 and 122.sub.3 while moving the sliding contacts 272 and 274 to null the output of the differential amplifiers 268 and 270. Since the switches 176 shown in FIG. 8b are closed infrequently, only when the LP solution circuitry of FIG. 8a is utilized switches 260 and 262 are of the self-locking type so as to remain closed since the charge on the capacitors 264 and 266 would leak off and this would in turn erroneously reposition the valves 122.sub.1 and 122.sub.3.

E. computing and Storing the Inverse Basis Matrix

The p.sub.14 .degree., p.sub.12 .degree., p.sub.24 .degree. and p.sub.22 .degree. signals which are transmitted over lines 302, 304, 306 and 308 and form the basis matrix in the LP; i.e., p.sub.14 .degree. = B.sub.14 .degree., p.sub.12 .degree. = B.sub.12 .degree., p.sub.24 .degree. = B.sub.24 .degree. and p.sub.22 .degree. = B.sub.22 .degree., are applied to a circuit 310 for computing the inverse of the base matrix .beta..degree.. The output of the inverse base matrix circuit 310 is applied through switches 312.sub.1, 312.sub.2, 312.sub.3 and 312.sub.4 to capacitors 314.sub.1, 314.sub.2, 314.sub.3 and 314.sub.4. The signals .beta..sub.11 .degree., .beta..sub.12 .degree., .beta..sub.21 .degree., and .beta..sub.22 .degree. are applied to the capacitors 314.sub.1-4 when movable contacts a-2 are in contact with stationary contact b-2 of the switches 312.sub.1-4. The movable contacts a-2 are placed in this position by means of a linkage with the switches 176 of the LP solution circuitry shown in FIG. 8a but only for a sufficient length of time so as to charge the storage capacitors 314.sub.1-4 to voltages representing the values of .beta..sub.11 .degree., .beta..sub.12 .degree., .beta..sub.21 .degree. and .beta..sub.22 .degree..

F. computing and Storing the Simplex Multipliers

The signals representing .beta..degree. are now applied over lines 316.sub.1-4 to circuitry shown in FIG. 7d for computing simplex multipliers .pi..sub.1 .degree. and .pi..sub.2 .degree.. The simplex multiplier circuitry comprises multiplying amplifiers 318.sub.1-4 which combine the signals .beta..sub.11 .degree., .beta..sub.21 .degree., .beta..sub.22 .degree. and .beta..sub.12 .degree. with signals representing c.sub.2 .degree. and c.sub.4 .degree. on lines 306 to obtain the product signals representing .beta..sub.21 .degree. c.sub.4 .degree., .beta..sub.22 .degree. c.sub.2 .degree., .beta..sub.11 .degree. c.sub.4 .degree. and .beta..sub.12 .degree. c.sub.2 .degree.. The signals c.sub.2 .degree. and c.sub.4 .degree. are obtained from lines 306. Outputs from the multiplying amplifiers 318.sub.1-4 are applied to a pair of summing amplifiers 322.sub.1 and 322.sub.2 to obtain the simplex multipliers .pi..sub.2 .degree. and .pi..sub.1 .degree. on lines 324.sub.182.

G. checking Local Optimization

Unlike the prior art practice of striving to maintain the operating conditions at the local units 110.sub.1-4 on which the LP solution is based, the operating units 110.sub.1-4 in accordance with this invention are free to vary or change the operating conditions so as to change the LP solution as long as the local objective functions for the individual units 110.sub.1-4 are satisfied. The process of changing the operating conditions to change the yield and cost coefficients in the LP begins when, as shown in FIG. 7a, switches 328.sub.1 and 328.sub.2 are actuated so as to place movable contacts a-3 in the contact with stationary contacts b-3 so as to charge capacitors 330.sub.1 and 330.sub.2 to a voltage proportional to the signals p.sub.24 .degree. and c.sub.4 .degree.. Movement of the switches 328.sub.1 and 328.sub.2 to this position is controlled automatically by a linkage with the switch 176 of the LP solution circuitry shown in FIG. 8a. This in turn applies signals p.sub.24 .degree. and c.sub.4 .degree. to summing amplifiers 332.sub.2 and 332.sub.3 while a signal p.sub.14 .degree. from the resistor 162 is applied to a summing amplifier 332.sub.1. Simultaneously, a switch 334 which connects a trigger pulse generator with a pulse generator for generating pulses of random amplitude is closed by a linkage connected to the switch 176 in FIG. 7b. The trigger pulse generator which generates pulses at a constant rate triggers the random pulse generators 338.sub.1-3 which apply the pulses of random amplitude to the summing amplifiers 332.sub.1-3 so as to produce new yield coefficient signals p.sub.14.sup.1 and p.sub.24.sup.1 and a new cost coefficient signal c.sub.4.sup.1 after a polarity inversion by polarity changing (PC) amplifiers 340.sub.1-3.

Since the pulses from the generators 338.sub.1-3 are of purely random amplitude, there is no certainty that the change in operating conditions represented by these pulses will be consistent with the optimization of the local unit as well as local constraints let alone the local objective function which is consistent with the overall optimization of the system. Accordingly, a first test is made to determine if the change represented by the random pulses will optimize the local unit. This is accomplished by checking if the yield coefficient for the most valuable product stream which flows through line 112.sub.5 more than offsets the change in yield coefficient for the less valuable product stream flowing through the line 112.sub.6 as shown in FIG. 6.

The check is performed by first applying signal p.sub.14 .degree. and p.sub.24 .degree. to the polarity changing amplifiers 342.sub.1 and 342.sub.2 to the summing amplifiers 344.sub.1 and 344.sub.2. The resulting signals -p.sub.14 .degree. and -p.sub.24 .degree. are applied to summing amplifiers 344.sub.1 and 344.sub.2 along with the new yield coefficient signals p.sub.14.sup.1 and p.sub.24.sup.1. If the resulting change .DELTA.p.sub.14 in the p.sub.14 yield coefficient is more than four times greater than the resulting change .DELTA.p.sub.14 in the yield coefficient, then the unit 110.sub.4 will be locally optimized since the yield in stream 112 is considered to be four times as valuable as the yield in stream 112. This check is performed by applying a signal proportional to .DELTA.p.sub.14 to the non-inverting terminal of a differential amplifier 346 and the signal .DELTA.p.sub.24 to the inverting terminal of the amplifier 346 after multiplication by a factor of 4 due to a resistance 348 in the input circuit. If the output is greater than zero as determined by a checking circuit 348, the relay of the comparison circuit 348 will be actuated to close switches 345.sub.1-3 by means of a mechanical linkage 356 so as to apply signals p.sub.14.sup.1, p.sub.24.sup.1 and c.sub.4.sup.1 to the circuit for checking local constraints. However, if .DELTA.p.sub.14 - 4.DELTA.p.sub.24 is negative, the switches 354.sub.1-3 will remain open awaiting the testing of the next set of random pulses from the generator 338.sub.1-3 as controlled by the trigger pulse generator 336. Note that the initiation of a change in operating conditions is entirely automatic in this embodiment and not a result of a managerial decision.

H. checking Local Constraints

Assuming now that the signal .DELTA.p.sub.14 - 4.DELTA.p.sub.24 is positive indicating local optimization so as to close the switches 354.sub.1-3, the new yield coefficients p.sub.14.sup.1 and p.sub.24.sup.1 and the new cost coefficient c.sub.4.sup.1 are applied over lines 356.sub.1-4 to a summing amplifier 358 shown in FIG. 7b. As pointed out previously, the local constraints on the unit 110.sub.4 require that p.sub.14.sup.1 + 2p.sub.24.sup.1 + 3c.sub.4.sup.1 = 2 or p.sub.14.sup.1 + 2p.sub.24.sup.1 + 3c.sub.4.sup.1 -2 = 0. In order to satisfy this equation representing the local constraints, the line 356.sub.2 includes a resistor 360.sub.1 having a value representing the integer 2 and the line 356.sub.3 includes a resistor 360.sub.2 representing the integer 3. Also, a DC source generating a -1 volt is connected to a fourth input terminal of the summing amplifier 358 through a resistor 360.sub.3 having a value representing the integer 2.

If the local constraints are satisfied, the error at the output of the summing amplifier 358, after polarity inversion by an amplifier 362 and integration by a capacitor 364 should equal zero. In order to determine when the charge on the capacitor 364 does equal zero, the zero voltage across a capacitor 366 which is connected to ground is compared with voltage across the capacitor 364 by comparison circuit 368. If the error does equal zero, switches 370 will be closed by a relay linkage 371 so as to apply the signals p.sub.14.sup.1, p.sub.24.sup.1 and c.sub.4.sup.1 to the circuitry for checking the local objective function. If the error from the amplifier 358 is more or less than zero, the switches 370 will remain open awaiting the next set of random pulses from the generator 338.sub.1 which do pass the local optimization check.

I. checking the Local Objective Function

Assuming that the local constraints are satisfied and the switches 370 are closed by the relay of the comparison circuit 368, the new yield coefficient signals p.sub.14 ' and p.sub.24 ' and the new cost coefficient signal c.sub.4 ' are applied as inputs to summing amplifiers 382.sub.1-3. Signals representing -p.sub.14 .degree., -p.sub.24 .degree. and -c.sub.4 .degree. are also applied to the summing amplifiers 382.sub.1-3 having been stored on capacitors 374.sub.1-3. Note that the charge on the capacitors 374.sub.1-3 is a result of the application of signals p.sub.14 .degree., p.sub.24 .degree. and c.sub.4 .degree. over lines 168, 170 and 172 to polarity changing amplifiers 378 through switches 380 which are opened by partially shown linkages with the switches 176 in FIG. 8b and the switch 334 in FIG. 7a when these switches are closed.

After polarity changes performed by amplifiers 384.sub.1-3 connected to the outputs of the summing amplifier 382.sub.1-3 signals representing .DELTA.p.sub.14, .DELTA.p.sub.24 and .DELTA.c.sub.4 are generated. The .DELTA.p.sub.14 and .DELTA.p.sub.24 signals are applied to multiplying circuits 386.sub.1 and 386.sub.2 where they are combined with the simplex multiplier signals .pi..sub.1 .degree. and .pi..sub.2 .degree. which represent the voltage across capacitors 388.sub.1 and 388.sub.2 charged by the signals passed through closed switches 378 in the lines 324.sub.1-2 which are opened by means of partially shown linkages when the switches 176 in FIG. 8a and the switch 334 in FIG. 7a are closed. The output from the multiplying circuits 386.sub.1 and 386.sub.2 are signals representing the product of .pi..sub.1 .degree..DELTA.p.sub.14 and .pi..sub.2 .degree..DELTA.p.sub.24 which are applied to a summing amplifier 392. The output from the summing amplifier 392 is then summed with the signal .DELTA.c.sub.4 at summing amplifier 394 to produce a signal, after a polarity change at the amplifier 396, representing .DELTA.c.sub.4 - .pi..sub.1 .degree..DELTA.p.sub.14 - .pi..sub.2 .degree..DELTA.p.sub. 24. This signal is then multiplied at a multiplication circuit 398 by a signal x.sub.4 .degree. which is applied over line 326 from the circuitry shown in FIG. 7c to charge capacitor 400 through a switch 402 which is closed by means of a partially shown linkage until the switches 176 shown in FIG. 8a and the switch 334 shown in FIG. 7a are closed. The product at the output of the multiplying circuit 398 is a signal f.sub.4 which is the local objective function (.DELTA.c.sub.4 - .pi..sub.O.DELTA. p.sub.4) x.sub.4 .degree..

Since the overall system profitability will be maximized only when f.sub.4 is maximized, thee new f.sub.4 signal which is applied across a capacitor 404 is compared with the previous f.sub.4 signal which was applied through a closed switch 406 across a capacitor 408. For performing this comparison, a comparison circuit 410 is provided. If the comparison circuit 410 determines that the voltage across the capacitor 404 is greater than the voltage across the capacitor 408, the comparison 410 will actuate its relay to close the switch 406 momentarily. The relay will also close switch 412 in the line 326, and switches 414 connected between lines 324.sub.1 and 324.sub.2 and the capacitors 388.sub.1 and 388.sub.2 to provide a new x.sub.4 .degree. signal (x.sub.4.sup.1) and new .pi..sub.1 .degree. and .pi..sub.2 .degree. signals (.pi..sub.1.sup.1 and .pi..sub.2.sup.1) to the local objective function circuitry in preparation for the next check. The relay of the comparison circuit 410 will also close switches 416 and 418 to apply signals .DELTA.p.sub.14 and .DELTA.p.sub.24 over lines 420 and the signal c.sub.4.sup.1 over lines 422 to the circuitry of FIGS. 7c and d. In addition, closure of the switches 418 momentarily will charge the capacitors 474.sub.1-3 to values representing c.sub.4.sup.1, p.sub.14.sup.1 and p.sub.24.sup.1 in preparation for the next check of the local objective function.

Should the local objective function f.sub.4 fail to be further maximized by the change in operating conditions the relay in the comparison circuit 410 will not be actuated and the switches 412, 414, 416 and 418 will remain open in the position shown awaiting the next set of pulses from the random pulse generator 338.sub.1-3.

J. correcting the LP Solution

1. Computing sensitivity Coefficients

Assuming that the change in operating conditions does maximize the local objective function, the LP solution is corrected in accordance with block 30a of FIG. 2a by circuitry shown in FIGS. 7c and 7d. In this connection, the signals .DELTA.p.sub.14, and .DELTA.p.sub.24 are applied to multiplying circuits 424.sub.1-4 along with signals representing the coefficients of .beta..degree. to generate product signals representing .beta..sub.14 .degree..DELTA.p.sub.14, .beta..sub.12 .degree..DELTA.p.sub.24, .beta..sub.24 .degree..DELTA.p.sub.14 and .beta..sub.22 .degree..DELTA.p.sub.24. The signals .beta..sub.14 .degree..DELTA.p.sub.14 and .beta..sub.12 .degree..DELTA.p.sub.24 are summed at a summing amplifier 426.sub.1 to generate signals representing the negative of the sensitivity coefficient s.sub.1. At the same time, signals .beta..sub.24 .degree..DELTA.p.sub.14 and .beta..sub.22 .degree..DELTA.p.sub.24 are summed at the summing amplifier 426.sub.2 to generate the negative of the sensitivity coefficient s.sub.2. To obtain the sensitivity coefficients s.sub.1 and s.sub.2, the signals -s.sub.1 and -s.sub.2 are applied to polarity changing amplifiers 428.sub.1 and 428.sub.2.

2. Computing New x's

Sensitivity coefficient signals s.sub.1 and s.sub.2 may now be utilized to compute the new values of the feed stream flow rates x.sub.4.sup.1 and x.sub.2.sup.1 to determine if they satisfy the constraints requiring that x.sub.4.sup.1 be greater than or equal to zero and x.sub.2.sup.1 be greater than or equal to zero. This check is performed by summing the signal -s.sub.1 with a signal equal to -1 at a summing amplifier 430. The sum signal s.sub.1 + 1 is applied to a dividing circuit 432.sub.1 along with the sensitivity coefficient signal s.sub.1 to obtain a signal representing ##EQU11## Similarly, the sum signal s.sub.1 + 1 is applied to a summing circuit 432.sub.2 along with the sensitivity coefficient signal s .sub.2 to obtain a signal representing ##EQU12## These signals are then multiplied by x.sub.4 .degree. appearing on lines 326 at multiplying circuits 434.sub.1 and 434.sub.2 to obtain a signal ##EQU13## These signals which appear at the output of the multipliers 434.sub.1 and 434.sub.2 are then applied to summing amplifiers 436.sub.1 and 436.sub.2 and undergo a polarity change at amplifiers 438.sub.1 and 438.sub.2 to generate signals representing the new feed stream flow rates x.sub.4.sup.1 and x.sub.2.sup.1. The signals x.sub.4.sup.1 and x.sub.2.sup.1 are applied to checking circuits 440.sub.1 and 440.sub.2 to determine if x.sub.4.sup.1 is greater than or equal to zero and x.sub.2.sup.1 is greater than or equal to zero so as to satisfy the constraints on the feed stream flow rates.

If this last check indicates that x.sub.4.sup.1 and x.sub.2.sup.1 are greater than zero, the LP solution is corrected by applying the signals x.sub.4.sup.1 and x.sub.2.sup.1 which appear across capacitors 442.sub.1 and 442.sub.2 at the output of the checking circuits 440.sub.1 and 440.sub.2 to the capacitors 238 and 240 as the relays associated with the checking circuits 440.sub.1 and 440.sub.2 move the movable contacts a-1 of switches 234 and 236 to a position of contact with the stationary contacts c-1. This in turn applies the x.sub.4.sup.1 and x.sub.2.sup.1 signals to the differential amplifiers 244 and 242 over lines 239 and 241 respectively so as to reposition the valves 122.sub.2 and 122.sub.4. Of course, if the check at the checking circuits 440.sub.1 and 440.sub.2 had indicated that x.sub.4.sup.1 and x.sub.2.sup.1 were less than zero, the LP solution would not be corrected but rather wait for the next set of random pulses from the generators 338.sub.1-3 representing another change in operating conditions at the unit 110.sub.4.

3. Computing New Basis Matrix

Correction of the LP further requires that the signals representing .beta..sub.24.sup.1, .beta..sub.22.sup.1, .beta..sub.14.sup.1 and .beta..sub.12.sup.1 be applied across the capacitors 314.sub.1-4 to update the inverse of the basis matrix to reflect the changes in operating conditions at the unit 110.sub.4. The signals .beta..sub.24.sup.1, .beta..sub.22.sup.1, .beta..sub.14.sup.1 and .beta..sub.12.sup.1 are generated by the application of the signals .beta..sub.24 .degree., .beta..sub.14 .degree. and .beta..sub.12 .degree. to multiplying circuits 446.sub.1-4 along with the signals ##EQU14## to generate signals representing .DELTA..beta..sub.24, .DELTA..beta..sub.22, .DELTA..beta..sub.14 and .DELTA..beta..sub.12. The signals at the output of the multipliers 446.sub.1-4 are then added to the signals .beta..sub.24 .degree., .beta..sub.14 .degree. and .beta..sub.12 .degree. with summing amplifiers 448.sub.1 having outputs connected to polarity changing amplifiers 450.sub.1-4. The signals .beta..sub.21.sup.1, .beta..sub.22.sup.1, .beta..sub.14.sup.1 and .beta..sub.12.sup.1 may then be applied to the capacitors 314.sub.1-4 for computation of the base matrix in response to closure of the switches 444 in lines 452 connected to the stationary contacts c-4 of switches 312.sub.1. Thus when the checking circuits 440.sub.1 and 440.sub.2 indicate that the constraints on x.sub.4.sup.1 and x.sub.2.sup.1 are satisfied, the relays of the checking circuits 440.sub.1 and 440.sub.2 will automatically close the switches 444 and place the movable contacts a-2 of the switches 312.sub.1 in contact with the stationary contacts c-4 so as to charge the capacitor 314.sub.1-4 representing .beta..sub.24.sup.1, .beta..sub.22.sup.1, .beta..sub.14.sup.1 and .beta..sub.12.sup.1.

4. Computing New Simplex Multipliers

The .beta..sub.21.sup.1, .beta..sub.22.sup.1, .beta..sub.14.sup.1 and .beta..sub.12.sup.1 signals are also utilized to generate new simplex multipliers .pi..sub.2.sup.1 and .pi..sub.1.sup.1. Referring now to FIG. 7d, the signals .beta..sub.24.sup.1, .beta..sub.22.sup.1, .beta..sub.14.sup.1 and .beta..sub.12.sup.1 are applied to multipliers 452.sub.1-4 along with the new cost coefficient signal c.sub.4.sup.1 on line 422 which corresponds with the change in operating conditiions at unit 110.sub.4 and the cost coefficient signal c.sub.2 .degree. on line 306. The signals .beta..sub.24.sup.1 c.sub.4.sup.1 and .beta..sub.22.sup.1 c.sub.2 .degree. are summed by the summing amplifier 454.sub.1 having an output coupled to a polarity changing amplifier 456.sub.1 to obtain the simplex multiplier signal .pi..sub.2.sup.1. Similarly, signals .beta..sub.14.sup.1 c.sub.4.sup.1 and .DELTA..sub.12.sup.1 c.sub.2 .degree. are summed at a summing amplifier 454.sub.2 having an output coupled to a polarity changing amplifier 456.sub.2 to obtain a simplex multiplier signal .pi..sub.1.sup.1. The signals .pi..sub.2.sup.1 and .pi..sub.1.sup.1 are applied through momentarily closed switches 458 in line 460 to the local objective function circuitry shown in FIG. 7b. Closure of the switches 458 is accomplished by the relays of the checking circuits 440.sub.1 and 440.sub.2 through a linkage 462. In this connection, it will be understood that momentary closing of the switches 458 in FIG. 7c and the switches 414 in FIG. 7b are sufficiently coincident to allow the signals .pi..sub.2.sup.1 and .pi..sub.1.sup.1 to be applied to and stored on the capacitors 388.sub.1 and 388.sub.2.

Conditions at the local unit 110.sub.4 must of course be changed such that the operating conditions actually correspond with those represented by the corrected LP. This is actually accomplished simultaneously with the correction of the LP by the closure of a switch 464 as shown in FIG. 7a so as to apply the signal p.sub.14.sup.1 at the output of the amplifier 340.sub.1 through a resistor 466 which is connected to the inverting input terminal of a differential amplifier 468. Closure of the switch 464 accomplished by the linkage 462 associated with the checking circuits 440.sub.1 and 440.sub.2 shown in FIG. 7c remains closed for a sufficient length of time to charge the capacitor 466 to a voltage corresponding with the new yield coefficient p.sub.14.sup.1.

Differential amplifier 468 which has its non-inverting terminal connected to the sliding contact 152 of the slidewire 154 generates an output which is utilized to drive a motor 470 which in turn repositions the contact 152 so as to null the output from the amplifier 468. This in turn will produce unbalance at the amplifier 150 so as to drive the motor 156 and reposition the valve 126.sub.4 until the change in temperature of the unit 110.sub.4 as sensed by the sensor 124.sub.4 again produces a null at the output of the amplifier 150. This control of the operating conditions assumes a linear relationship between the coefficient p.sub.14 and the temperature of the operating unit. This linearity permits resistance 162 having a value kp.sub.14 .degree. to convert temperature into yield while permitting the resistor 466 having a value 1/kp.sub.14 .degree. to convert yield into temperature.

The linkage 462 to the switch 464 also produces momentary closure of movable contacts a-3 with stationary contact c-3 of switches 328.sub.1-3 shown in FIG. 7a. This allows capacitors 330.sub.1 and 330.sub.2 to charge to values representing p.sub.24.sup.1 and c.sub.4.sup.1 in preparation for the next proposed change in operating conditions at the unit 110.sub.4. In other words, the circuitry shown in FIG. 7a as well as that shown in FIGS. 7 (b-d) is now ready to process the next proposed change in operating conditions in response to the next trigger pulse from the trigger pulse generator 336 which will result in the application of pulses of random amplitude to the summing amplifier 332.sub.1-3. the entire process may be repeated for the local unit 110.sub.4.

Although the local computer circuits 114.sub.1-3 have not been shown, it should be understood that similar changes in operating conditions may be proposed and checked at the various local units 110.sub.1-3. This of course requires any change in operating conditions at the unit 110.sub.4 to be transmitted to the system computer 116 in accordance with block 32a and then transmitted to the other units 110.sub.1-3 so that each unit or its local computer 114.sub.1-3 is aware of any correction in the LP and only makes changes in operating conditions with respect to the corrected LP. (Of course, all changes in the solution of x's are transmitted over lines 239 and 241.) In order to assure that changes in operating conditions are not implemented at the same time, a single trigger pulse generator 336 may be utilized for all of the local computers 114.sub.1-4 such that only one out of four trigger pulses from the trigger pulse generator 336 is applied to the pulse generators 338.sub.1-3 at any one local computer 114.

5. Computing New Cost Coefficients

The correction of the LP may proceed indefinitely without utilizing the LP solution circuitry shown in FIG. 8a as long as the nonbasic cost coefficients c.sub.1.sup.1 and c.sub.3.sup.1 remain negative. However, if the cost coefficients c.sub.1.sup.1 and c.sub.3.sup.1 become positive, this indicates that the basis of the LP has changed requiring that the LP be solved again and the new bases determined. This is done automatically with each correction of the LP by the linkage 462 associated with the checking circuits 440.sub.1 and 440.sub.2 in FIG. 7c which closes switches 472 so as to apply the .pi..sub.2.sup.1 and .pi..sub.1.sup.1 signals to lines 474.sub.1 and 474.sub.2 to summing amplifiers 476.sub.1 and 476.sub.2 along with the signals -.pi..sub.2 .degree. and -.pi..sub.1 .degree.. The output of the summing amplifiers 476.sub.1 and 476.sub.2 are coupled to polarity changing amplifiers 478.sub.1 and 478.sub.2 to generate signals .DELTA..pi..sub.1 and .DELTA..pi..sub.2. The signal .DELTA..pi..sub.2 is multiplied by signals p.sub.23 .degree. and p.sub.21 .degree. appearing on lines 304 at multipliers 480.sub.1 and 480.sub.2 while the signal .DELTA..pi..sub.1 is multiplied by signals .pi..sub.13 .degree. and .pi..sub.11 .degree. on lines 304 at multipliers 480.sub.3 and 480.sub.4. Summing amplifiers 482.sub.1 and 482.sub.2 at the output of the multipliers 480.sub.1-4 have outputs connected to summing amplifiers 484.sub.1 and 484.sub.2 which produce signals, after polarity changes by amplifiers 486.sub. 1 and 486.sub.2, which represent c.sub.1.sup.1 = c.sub.1 .degree. - .DELTA..pi..sub.1 p.sub.11 .degree. - .DELTA..pi..sub.2 p.sub.21 .degree. and c.sub.1.sup.1 = c.sub.1 .degree. - .DELTA..pi..sub.1 p.sub.13 .degree. - .DELTA..pi..sub.2 p.sub.23 .degree..

The signals c.sub.3.sup.1 and c.sub.1.sup.1 are then applied to checking circuits 488.sub.1 and 488.sub.2 to determine if c.sub.1.sup.1 is greater than zero or c.sub.3.sup.1 is greater than zero. If either c.sub.1.sup.1 or c.sub.3.sup.1 is greater than zero, the relay of the checking circuit 488.sub.1 or the relay of the checking circuit 488.sub.2 will open the switch 334 connected to the output of the trigger pulse generator 336 as shown in FIG. 7a. Simultaneously, switches 490 and 492 on line 494 and 496 will close in the LP solution circuit shown in FIG. 8a so as to apply the signals p.sub.14.sup.1 and p.sub.24.sup.1 generated in the circuit of FIG. 7a to the differential amplifiers 178 shown in FIG. 8a while switches 176 open to reset the value of the yield of the yield coefficient resistors 174 corresponding to the yield coefficients p.sub.14 and p.sub.24. Similarly, the signal c.sub.4.sup.1 on line 498 in FIG. 7a will be applied to the LP solution circuit to reset the value of the cost coefficient resistor 190 corresponding to c.sub.4. The other resistors in the LP solution circuit will also be reset in response to signals from the local computers 114.sub.1-3 corresponding to the then existing yield and cost coefficients.

Once the LP has been solved again utilizing the circuit shown in FIG. 8a and the bases determined utilizing the circuitry shown in FIG. 8b, the new LP solution may be transmitted to the local computers 114.sub.1-4. Changes in the operating conditions at the local operating units may then again be implemented and the LP solution corrected accordingly utilizing the circuitry of FIGS. 7(a-d).

V. Control of the Simple Processing System With a Digital Computer

In the foregoing embodiment of the invention, an analog computer has been utilized for the system computer 116 and analog computers have also been utilized for the local computers 114.sub.1-4. It is however possible to replace all or any portion of the analog computers with a properly programmed digital computer. In this connection, a program for performing the LP correction provided by the analog circuitry shown in FIG. 7c and 7d will now be described with reference to the block diagram of FIG. 9 and the computer program listing attached hereto as Appendix A.

A. correction of the LP

As shown in FIG. 9, the original solution to the LP is stored as indicated by block 500. In particular, the various values of the last LP solution are read into register means within a digital computer such as the IBM 1800 in accordance with instructions 2-22 of the program listing including the storage of B.degree., .beta..degree., .gamma..degree., .pi..degree., c.degree. and b.

In the next step as depicted by block 502 of FIG. 9, the change in yield and cost coefficients for a particular column in the LP are stored in register means within the digital computer. The step depicted in the block 502 is performed by instructions 23-26 of the program listing.

The sensitivity coefficients for the particular change in yield are now computed where the sensitivity coefficients s are computed in accordance with the equation shown in block 504 of FIG. 9. The step depicted in block 504 is performed in response to instructions 27-32 of the program listing. When utilizing the program listing to correct the LP of Equations 31 and 32, the sensitivity coefficients s.sub.1 and s.sub.2 will be computed.

Now that the sensitivity coefficients have been computed, the new values of x may be computed utilizing the equation shown in block 506 of FIG. 9 which is the same equation which was utilized to compute the new values of x in the analog embodiment of FIGS. 7c and d where b.sub.j .degree. and b.sub.m are equal to x.sub.j and x.sub.m. The computation of the new values of x, i.e., x.sub.4.sup.1 and x.sub.2.sup.1, is performed in response to instructions 33 and 34 of the program listing.

A new inverse base matrix is computed utilizing the equation depicted in block 508 of FIG. 9 which is the same equation which was utilized in the analog embodiment to compute the new inverse basis matrix .beta..sup.1. Instructions 35-37 of the listing provide for the computation of the inverse basis matrix.

Block 508 corresponds to the computation of new simplex multipliers in accordance with the indicated equation, which computation is provided by instructions 38-43 of the program listing. Finally, block 512 provides for the computation of new nonbasic cost coefficients c.sub.j.sup.1 utilizing the indicated equation which is identical to that utilized in computing these cost coefficients in the analog circuitry of FIG. 7d.

B. computing the Local Objective Function

The local objective function may be computed by a digital computer such as the IBM 1800 which has been programmed with the nonlinear optimization program of "Least Squares Estimation of Nonlinear Parameters," by D. W. Marquardt, Engineering Dept., E. I. duPont de Nemours and Co., Inc., Wilmington Del., March, 1964, commercially available under Distribution No. 3094 as described in the manual attached hereto as Appendix B-1. Listings for subroutines SUBZ, F code and P code referred to on page 1 of the Marquardt manual for computing the local objective functions for the system shown in FIG. 6 are attached hereto as Appendices B-2, B-3 and B-4 respectively.

In the very short listing of Appendix B-2, instruction 1 reads in the data of the LP parameters and solution and calls the Marquardt program.

In the program listing of Appendix B-3, instructions 1-4 define the parameters and the dimensions of the system. Instruction 4 and 6 define penalty functions which are utilized in computing the objective function in instruction 23. Instructions 7-10 define the x's and the .pi.'s, and instructions 12 and 13 define the p.sup.1 's.

New p's and the .DELTA.c are computed in response to instructions 13-15 and sensitivity coefficients are computed in response to instructions 16 and 17.

The DO loop of instructions 18 and 19 computes a function EQUAL (i) for i = 1, 2. If as computed at instruction 21, the function Equal (i) is greater than zero, a penalty 1 equal to zero is stored. If a computed at instruction 22, the function Equal (i) is greater than zero, a penalty 2 equal to zero is stored..

A local objective function is then computed in response to instruction 3 where the penalties as defined in instructions 5 and 6 are added to the local objective function if Equal (i) for i = 1, 2 is less than or equal to zero.

The P code of Appendix B-4 merely indicates that no analytic derivatives are provided and those of the Marquardt program must be utilized as described on page 1 of the Marquardt manual.

VI. A Numerical Example

A simple numerical example based on the system of FIG. 6 will now be solved to illustrate the foregoing method. The problem is:

Maximize Z = -6x.sub.1 - 4x.sub.2 - x.sub.3 + c.sub.4 x.sub.4 (33) ##EQU15## and

x.sub.i .gtoreq. 0 for i = 1, ...., 4 (35) p.sub.i4 .gtoreq. 0 i = 1, 2 c.sub.4 .gtoreq. 0

(35a)

furthermore,

p.sub.14 + 2p.sub.24 + 3c.sub.4 = 2 (35b)

where p.sub.1, p.sub.2, p.sub.3, c.sub.1,c.sub.2 and c.sub.3 have been estimated. Notice that the 4th column of Equation (34) is not specified, and free to be chosen such that constraints (35a) and (35b) are satisfied.

Step 1: Estimate p.sub.j and c.sub.j values for all units including p.sub.4 and c.sub.4 using the circuitry of FIG. 7a such that

p.sub.14 2 p.sub.24 .degree. = 0 (36) c.sub.4 0

Notice that Equation (36) satisfies local constraints (35a) and (35b).

Step 2: Solve Equations (33), (34) and (36) using the circuitry of FIG. 8a. This corresponds to the overall planning problem. The optimial solution is: ##EQU16##

Bases are x.sub.4 and x.sub.2. ##EQU17## Step 3: Local optimization for the 4th column has the local objective function:

Max f.sub.4 = Max (.DELTA.c.sub.4 - .pi..degree..DELTA.p.sub.4) x.sub.4 (40)

The local constraints are Equations (35a), (35b) and the following equations: ##EQU18## This local optimization may be solved using any of the standard non-linear optimization program such as IBM's POP-II (refer to "A Process Optimization Program for Nonlinear Systems: POP-II," by H. V. Smith, Mar. 18, 1965, IBM Program No. 7090 H9 IBM 0021) or aforesaid "Least Squares Estimation of Nonlinear Parameters" by D. W. Marquardt of Appendix B-1. It may also be solved on the analog computer circuitry of FIG. 7b. The solution to this local optimization is given as ##EQU19##

The local objective function, Equation (40), is made positive which indicates that if the new p.sub.4 replaces the old column p.sub.4 .degree., the overall profit, Z, increases. This optimum solution is illustrated graphically in FIG. 5 where p.sub.14 is the abscissa and p.sub.24 is the ordinate. The constraints (41) and (42) are depicted by vertical lines and various values of the cost coefficients are depicted by lines representing c.sub.4 = -1/2, c.sub.4 = 0 and c.sub.4 = 1/2 and the area within the triangle BOA represents the local constraints (35a). The original column p.sub.4 .degree. is located at point A(2,0,0) and the local objective function indicates that f.sub.4 is increasing along line ACB so as to be maximized at point C where .DELTA.p.sub.14 = 4/7, .DELTA.p.sub.24 = 2/7 and .DELTA.c.sub.4 = 0.

Step 4: The overall solution based on the new p.sub.4 yield is

x.sub.1 1 = 0 p.sub.14 1 = 10/7 x.sub.2 0 Z.sup.1 = O, p.sub.24 (47) x.sub.3 0 c.sub.4 2/7 x.sub.4 7/2 0

which may be obtained using a digital computer and the program of Appendix A or the analog computer circuitry of FIGS. 7c and d. The improved overall solution is obtained by replacing the 4th column of Equations (33) and (34). ##EQU20## Equation (48) is transformed to Equation (49) at the optimal solution ##EQU21## and

0.sup.. x.sub.4 + 0.sup.. x.sub.2 - 12x.sub.1 - 2x.sub.3 = Z - 0 (50)

the bases are still x.sub.4 and x.sub.2. The new basis matrix B.sup.1 is ##EQU22## The inverse .beta..sup.1 of B.sup.1 becomes ##EQU23## The new simplex multipliers, .pi..sup.1 is

.pi..sup.1 = .gamma..sup.1 .beta..sup.1 = (1, -5)

These new quantities (denoted by a superscript 1) have been easily obtained from the original values (denoted by a superscript o) by Equations (21) through (30) using the computer program in Appendix A.

Step 5: Going back to step 3 and repeating the local optimization, using any of the standard computer programs described previously, the local objective function now is: ##EQU24## where s.sub.1, s.sub.2 are defined by ##EQU25## The local constraints are now Equations (35a), (35b), and the following equations: ##EQU26## From Equations (53) and (54), if .DELTA.c.sub.4 = 0, ##EQU27## Since .DELTA.p.sub.24 .ltoreq. 0 is the limiting constraint, it is clear that f.sub.4 cannot be positive.

Step 6: It is now assured that Z cannot be increased further without changing the bases. In this case, Z cannot be increased further even if it is allowed to change bases, because for the whole range of f.sub.4, - 4.ltoreq.f.sub.4 .ltoreq. 0, the optimality holds: That is, the relative cost coefficients of nonbasic columns remain negative always.

Although very simple, this example clearly has illustrated the procedure of this method.

The solution and correction of the LP which describes the system shown in FIG. 6 may, as indicated in the foregoing, be performed on an analog computer or a properly programmed digital computer. However, where the system becomes extremely complex as in the case of a modern day refinery such as that shown in FIG. 10, the use of an analog computer becomes impractical to either solve or correct the LP. A system requiring a digital computer to solve and correct the Lp will now be described.

VII. Method of Operating a Complex Refinery System

As shown in FIG. 10, the refinery comprises nine operating units in all including an atmospheric plus vacuum distillation tower 600 with a crude feed stream and eight product streams which form the feed streams for the other units including a gas plant 602, an alkylation unit 603, a treating and blending unit 604, a catalytic reformed 606 and a hydrodesultfurization unit 603. In addition, the distillation tower 600 includes a product stream which provides a feed stream for a catalytic cracking unit 610 and a delayed coker 612. The system further comprises a hydrocracker unit including a hydrogen manufacturing capability 614. A system computer 616 such as the IBM 7500 is linked with each of a plurality of feed stream valves 617 to control the flow of feed streams to the various units. In addition, the system computer is linked with the various units through lines 618 for the transmission of flow rate control signals and unit yield and cost signals.

It will of course be appreciated that individual units within the typical refinery system shown in FIG. 10 are extremely complex. Consider the catalytic cracker 610 which is shown in schematic form in FIG. 11.

As shown in FIG. 11, the catalytic cracker is supplied with feed which enters a feed heater 620 with recycled heavy cycle oil. The temperature of the feed at the outlet of the heater 620 is controlled in response to a temperature sensor 622 such as a thermocouple which transmits a signal to a local digital process control computer 624 such as the IBM 1800 after conversion by an analog-to-digital converter 626. The local computer 624 then resets a valve 628 in the fuel gas line to the heater 620 by means of a linkage with the computer 624.

Reslurry recycle which flows through a valve 630 controlled by the local computer 624 and heavy cycle oil recycle which flows through a valve 632 controlled by the local computer 624 join the feed at the outlet of the heater 620. The combined feed is then applied to a reactor 634. The temperature of the reactor is controlled by a temperature sensor such as a thermocouple 636 which is coupled to the computer 624 through an analog-to-digital converter 638. The computer 624 in turn controls the catalyst circulation rate from a regenerator 640 through a valve 642. The catalyst bed level within the reactor is controlled by a level sensor 644 which is coupled to the local computer 624 through an analog-to-digital converter 646. The computer 624 in turn controls the level by repositioning a valve 648 to control the circulation rate of the catalyst to the regenerator 640.

In order to control the temperature of the regenerator 640, a thermocouple 650 is provided which is coupled to the local computer 624 through an analog-to-digital converter 652. The local computer 624 adjusts a valve 654 in the steam flow line accordingly. The flow of air to the regenerator 640 is controlled by a valve 656 associated with an air blower 658.

The output from the reactor 634 flows to a fractionator 660. A slurry settler 662 is connected to the bottom of the fractionator 660 through a valve 664 which is controlled by the local computer 624. Gasoline is obtained from an overhead receiver 666 as well as gas where the flow of gas from the receiver 666 is controlled by the computer 624 through a valve 668 associated with a blower 670. A valve 672 is provided in the heavy cycle oil line under the control of the local computer 624.

Having now described a typical refinery system and a typical unit within the refinery system, the control of the system and the unit in accordance with this invention will now be described.

In accordance with block 20 of FIG. 2a, the cost and yield coefficient at standard or design operating condition are estimated. This is performed at the local computers of the individual units such as the local computer 624 of the catalytic cracker 610. Various computer programs are known to those of skill in the art for estimating the yields of various units. For purposes of illustration, a listing for a computer program for estimating the yields of the catalytic cracker 610 is provided in Appendix C for use with the local computer 624. It will of course be appreciated that the local computor 624 sets the various operating conditions of the catalytic cracking unit 610 by appropriately positioning the various valves so as to correspond with the estimated cost and yield coefficients.

In accordance with block 22 of FIG. 2a, the estimated yield and cost coefficients are now transmitted from the local computer 624 of the unit 610 to the system computer 616 where the LP is solved in accordance with block 24. Although various programs are commercially available in the art, one particularly suitable program listing for solving the LP is attached hereto as Appendix D-2. The resulting LP solution may then be transmitted to the local computers associated with the various units including the local computer 624 of the catalytic cracker 610.

In accordance with block 28 of FIG. 2, the proposed change in one of the units such as the catalytic cracker 610 is now checked utilizing the local objective function to see if the proposed change will be consistent with overall system optimization. This check of the local objective function may be performed utilizing a standard linear optimization program such as the aforesaid "Least Squares Estimation of Nonlinear Parameters" commercially available as SHARE Distribution No. 3094.

If the local objective function for the proposed change in the unit is satisfied, the local computer such as the local computer 624 may actually make the change in operating conditions by adjusting the valves of the catalytic cracker 610 or one of the other units. In addition, the local computer 616 may actually correct the LP as depicted by block 30a of FIG. 2a, using the program described in the listing attached hereto as Appendix A and previously referred to in conjunction with the flow diagram of FIG. 9. This change in local operating conditions and correction of the LP may be reiterated so as to optimize the local units consistent with the overall optimization of the system.

In the analog computer embodiment of FIGS. 7(a-d) and 8(a and b) and the digital computer embodiment for the system of FIGS. 9 and 10, the correction of the LP has actually been performed at the local computer where it has been assumed that the local computer has sufficient capacity to perform this function. However, the block diagram of FIG. 2 actually indicates that correction of the LP can be performed at the system computer where the local computer has insufficient capacity to perform this function. In reality, this seldom is the case since the computer program for correcting the LP is relatively simple, at least as compared with the program for initially solving the LP. This simplicity of the LP correction program which may readily be seen by reference to the program listing of Appendix A and comparison with the program listing for solving the LP of Appendix D-2, not only permits the correction to be performed at the local computer but also saves considerable computer time by eliminating the necessity for solving the LP again utilizing the rather involved LP solution program such as that represented by the listing of Appendix D-2.

XIII. Modified Method of Operating a Complex Processing System

In the foregoing, a sensitivity analysis has been performed on a perturbation of the yield coefficients in an LP. This sensitivity analysis has been performed at the local operating units to permit the local unit managers to optimize their units consistent with the overall optimization of the system. However, this sensitivity analysis need not be performed at the local unit level and it need not be performed on a perturbation which is consistent with local optimization. It may for example be advantageously performed at the overall system planning level where it becomes desirable to determine how a slight change in one or more yield columns p.sub.j will affect the overall optimization of the system regardless of the effect on the local optimization. A method utilizing this type of sensitivity analysis will now be described with reference to FIG. 4.

Assuming that the yields for the various units have been estimated, the LP solution for the system is solved at the overall system planning level using the estimated yields as indicated at block 34. A check is then made at the planning level to determine if a change in yields for a particular unit will satisfy the local objective function of that unit as indicated at block 36 so as to further optimize the overall system. If the change in yield will satisfy the overall objective function as indicated, the change in yield is transmitted to the local unit concerned as indicated at block 38. Then, using a mathematical model, the local unit determines if a change in the local operating conditions may be made as indicated at block 40. If and when the local unit makes the change, an indication to that effect is transmitted to the overall system planning as indicated at block 42. In block 44, the LP solution is corrected for the new yield by overall system planning. This process may also be iterated.

IX. Another Numerical Example

This type of sensitivity analysis will now be described with reference to the simplified refinery of Bightler and Wilde shown in FIG. 4. The refinery comprises a fuel processing unit 46 including a local processing computer 48 and a lube processing unit 50 including a local process computer 52. The scheduling for the units of this system is planned at the overall system planning level through the use of an overall system computer 54. Four grades of crude oil are fed to this refinery with all four going to the fuel processing unit 46 and only crude 4 going to the lube processing unit 50. An overall system or planning level computer 54 is in communiction with the units 46 and 50 through data links 56.

The products produced by this refinery are gasoline, oil, jet fuel and lube oil. A set of constraints on the availability of four different crudes and product demands are shown on the following Table 1.

Table 1 __________________________________________________________________________ Product Yields and Constraints Crude Identity Product 1 2 3 4 on order, Fuel Lube bbl/wk __________________________________________________________________________ Product Yield, vol % Gasoline 0.6 0.5 0.3 0.4 0.4 .ltoreq.170,000 Heating Oil 0.2 0.2 0.3 0.3 0.1 .ltoreq. 85,000 Lube Oil 0 0 0 0 0.2 .ltoreq. 20,000 Jet Fuel 0.1 0.2 0.3 0.2 0.2 .ltoreq. 85,000 Loss 0.1 0.1 0.1 0.1 0.1 -- Crude Available, bbl/wk 100,000 100,000 100,000 --200,000-- -- Profit, $/M bbl Crude 100 200 70 150 250 -- __________________________________________________________________________

The first column in Table 1 implies that 1 bbl of crude 1 is converted to 0.6 bbl of gasoline, 0.2 of heating oil, no lube oil, 0.1 jet fuel, and 0.1 bbl loss during the process for the profit of $100/M bbl of crude. The sixth element of this column represents the 100,000 bbl of crude available per week. The first row of Table 1 represents the constraint on gasoline demand. The total gasoline production rate must be lower than 170,000bbl/wk. The objective is to find the set of feasible crude production rates which maximize the profit by the equation

Z = 100x.sub.1 + 100x.sub.3 + 70x.sub.3 + 150x.sub.4 + 250x.sub.5 (58)

where x.sub.i represents the processing rate of the ith crude. Using the notations utilized by Beightler and Wilde, the Simplex Tableau of this example is given by Table 2 where u.sub.1-4, and u.sub.g-j are slack variables.

Table 2 __________________________________________________________________________ Simplex Tableau x.sub.1 x.sub.2 x.sub.3 x.sub.4 x.sub.5 u.sub.1 u.sub.2 u.sub.3 u.sub.4 u.sub.g u.sub.h u.sub.i u.sub.j 1 1 100 1 1 100 1 1 100 1 1 1 200 0.6 0.5 0.3 0.4 0.4 1 170 0.2 0.2 0.3 0.3 0.1 1 85 0.2 1 20 0.1 0.2 0.3 0.2 0.2 1 85 100 200 70 150 250 0 0 0 0 0 0 0 0 Z __________________________________________________________________________

The following optimal solution may be obtained:

x.sub.1 = 37.5 M bbl x.sub.2 = 100 M bbl x.sub.3 = 58.33 M bbl with z = $67,833/week x.sub.4 = 100 M bbl x.sub.5 = 100 M bbl

The optimal table can now be modified as shown in Table 3.

Table 3 __________________________________________________________________________ Simplex Tableau of Optimal Solution __________________________________________________________________________ u.sub.1 x.sub.2 u.sub.3 x.sub.4 x.sub.1 x.sub.3 x.sub.5 u.sub.j u.sub.2 u.sub.4 u.sub.g u.sub.h u.sub.i b 1 0 0 0 0 0 0 0 .75 .25 -.25 2.5 2.5 62.5 0 1 0 0 0 0 0 0 1 0 0 0 0 100. 0 0 1 0 0 0 0 0 0.167 0.833 1.667 -5 -5 41.67 0 0 0 1 0 0 0 0 0 1 0 0 -5 100. 0 0 0 0 1 0 0 0 -.75 -.25 2.5 -2.5 -2.5 37.5 0 0 0 0 0 1 0 0 -0.167 -0.833 -1.167 5 5 58.33 0 0 0 0 0 0 1 0 0 0 0 0 5 100. 0 0 0 0 0 0 0 1 -0.075 0.075 .25 1.25 -1.25 3.75 -133.3 - 66.7 -113.3 -100 .600 -67.833 __________________________________________________________________________

The base matrix consists of columns of u.sub.1, x.sub.2, u.sub.3, x.sub.4, x.sub.1, x.sub.3, x.sub.5, u.sub.j and

B.degree. = u.sub.1 x.sub.2 u.sub.3 x.sub.4 x.sub.1 x.sub.3 x.sub.5 u.sub.j (59) 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0.5 0 0.4 0.6 0.3 0.4 0 0 0.2 0 0.3 0.2 0.3 0.1 0 0 0 0 0 0 0 0.2 0 0 0.2 0 0.2 0.1 0.3 0.2 1 ______________________________________

The inverse of this is

.beta..degree. = (60) 1.0 .75 0 .25 -2.5 2.5 2.5 0 0 1 0 0 0 0 0 0 0 0.167 1 0.833 1.667 -5 -5 0 0 0 0 1 0 0 -5 0 0 -.75 0 -.25 2.5 -2.5 -2.5 0 0 -0.167 0 -0.833 -1.667 5 5 0 0 0 0 0 0 0 5 0 0 -0.075 0 0.075 0.25 -1.25 -1.25 1 __________________________________________________________________________

The row vector of the cost coefficients of basic columns is given by the following Equation:

.gamma..degree. = (0, - 200, 0, - 150, - 100, - 70, - 250, 0) (61)

The simplex multipliers then become

.pi..degree. = .gamma..degree..beta..degree. = (0, - 113.33, 0, - 66.67, - 113.33, - 100, - 600, 0) (62)

Having now solved the LP as indicated at block 34, the effect of a change in yield at a particular unit is determined as indicated at block 36. More particularly, a determination is made to see if change in yield at a particular unit will improve the overall objective function.

To demonstrate this type of sensitivity, it is assumed that the product yield distribution of the lube oil part of the refinery is slightly changed.

Assume that the yield of gasoline can be reduced by 5% to 0.35, heating oil increased by 5 % to 0.15, lube oil increased by 5% to 0.25, and finally, jet fuel decreased by 5% to 0.15. As the distribution of the product yield changes, the marginal profit should also change. Thus, let the profit from the lube plant with the new product yields by $200/M bbl, down by $50/M bbl from the original case. One of the immediate questions is whether the lube plant should be moved to the new operating conditions or not. An intuitive answer is no because local profit would decrease since the lube plant makes less gasoline, generally the most valuable product, thus losing some of the marginal profit. It will however be shown otherwise using the sensitivity analysis of this invention.

As indicated at block 38, it is important to know if the objective function is satisfied, i.e., whether the overall profit increases. If the local objective function (.DELTA.c.sub.x.sbsb.5 -.pi..degree. .DELTA.p.sub.x.sbsb.5) is greater than zero, the overall profit of the refinery will increase. It may be seen, substituting the suitable values into this equation, that ##EQU28## Therefore, the overall refinery profit is assured to increase if the manager can move the lube plant to the new point. Notice that this is contrary to intuition. Having thus satisfied the local objective function, the manager of the local unit may be notified that he can make the change. When the manager determines if he can make the change as indicated at block 40, the overall system planning is notified and the LP solution is corrected as indicated at block 42.

From Equation (22), the sensitivity column is ##EQU29##

From Equation (23), the new solutions are also directly obtained. For instance, ##EQU30##

Similarily, x.sub.2 = 100 - 0 = 100 u.sub.3 = 41.67 - 0.5833(80) = 88.33 x.sub.4 = 100 + .25(80) = 120 (65) x.sub.1 = 37.5 + 0.375(80) = 67.5 x.sub.3 = 58.33 - 0.5833(80) = 11.67 x.sub.5 = 100 - 0.25(80) = 80 u.sub.j = 3.75 + (0.1875)(80) = 18.75

Since none of these become negative, the new solution is feasible. It remains to show if the new solution is optimal. The row of .beta..degree. corresponding to x.sub.5 is from Equation (60),

.beta..sub.x.sbsb.5 .degree. = (0, 0, 0, 0, 0, 0, 5, 0) (66)

According to Equation (25), the inverse of the new base matrix, .beta..sup.1, is exactly the same as .beta..degree. except for the seventh column for ##EQU31## The changes in Simplex multipliers are by ##EQU32## Thus, .pi..sup.1 = (0, - 113.33, 0, - 66.67, - 133.33, - 100, -913.33, 0)

The optimality is confirmed for

C.sub.j = C.sub.j .degree.- .DELTA..pi.p.sub.j .ltoreq. 0 for j = all nonbasic (69)

columns

Although not described in detail, it will of course be appreicated that three or four or more levels of optimization may be achieved utilizing this invention. For example, it is possible to correct LP describing the operation of two or more plants utilizing this invention while also correcting an LP describing the operation of various units in each plant. Such a system would involve a system computer, a plant computer and local computers for the units in the plant.

Although specific embodiments have been described, it will be understood that the various modifications may be made without departing from the spirit and scope of the invention as set forth in the appended claims. ##SPC1##

APPENDIX "B-1"

march 1964

LEAST-SQUARES ESTIMATION OF NONLINEAR PARAMETERS

by

Donald W. Marquardt

This program is a revision of Share Distribution No. 1428 by T. Baumeister, III and Donald W. Marquardt. In addition to Mr. Baumeister, the following programmers have contributed significantly to the present version: J. Ann Sheldon, Ruby M. Stanley.

ENGINEERING DEPARTMENT

E.I. DU PONT DE NEMOURS AND CO., INC.

Wilmington, Delaware

LEAST-SQUARES ESTIMATION OF NONLINEAR PARAMETERS

I, purpose

Given a model

Y.sub.- = 1(X.sub.11, X.sub.12. . . . . X.sub.1m ; b.sub.1, b.sub.2. . . . . b.sub.11)

which predicts the value Y, of a dependent variable Y, where the model 1 contains m independent variables X.sub.11 and k parameters b.sub.1, and given n observations

(Y.sub.1, X.sub.11, X.sub.12. . . . . X.sub.im) i = 1,2. . . . . n,

this program will compute the least-squares estimates b.sub.1. That is, the program will adjust the b.sub.1 to minimize ##EQU33## The program provides options to use analytic or estimated derivatives .differential.f/.differential.b.sub.1, to control printing, to omit parameters (i.e., fix their values and obtain a constrained minimization of .PHI.), and to plot the observed and predicted values. Confidence limits for the estimated parameters are also computed. The program is written in FORTRAN IV, employing the maximum neighborhood algorithm and other computational procedures described in Reference 1 and the Appendix of this writeup. The program makes no use of machine-code inserts.

Ii. procedure

The user must supply the following items:

a. A FORTRAN subroutine named F CODE to evaluate Y for a specific combination of X.sub.11, and b.sub.1. F coding evaluates the function Y.sub.1 for one data point each time it is executed.

b. A FORTRAN subroutine named P CODE to evaluate .differential.Y.sub.1 /.differential.b.sub.1 fr h.sub.1, ... , b.sub.k. P coding evaluates the analytic derivatives of the function Y, for the same data point just processed by F coding.

c. A FORTRAN subroutine named SUBZ which is to be performed once as soon as the data are read in, e.g., to compute constants needed by F coding and P coding.

d. The data (Y.sub.1, X.sub.11, . . . , X.sub.1m) = 1, . . . , n.

During execution of the program, F coding will sometimes be used alone and will sometimes be followed immediately by P coding. If analytic derivatives (in P Coding) are not supplied, or if it is desired not to use them, the program will calculate estimates of .differential.Y.sub.1 /.differential.b.sub.1 using finite-difference approximations. It should be recognized that use of estimated derivatives will usually increase the computing time relative to the use of analytic derivatives, since the finite-difference derivative expressions for k parameters require (k = 1) evaluations of Y. The initial guesses for the b.sub.1 must be different from zero if estimated derivatives are used, since the estimated derivatives are calculated by the form: ##EQU34## where the multiplier .DELTA. = 10.sup.-.sup.5 is used unless overridden by input option. If the model is insensitive to some of its parameters, then use of a larger value of .DELTA. may be required to prevent zero increments for these parameters due to the finite-difference derivative rounding to zero.

Iii. initialize Coding

The FORTRAN coding to read in case constants, perform preliminary calculations required by a case, or any other one-time operations desired, is called for by CALL SUBZ (Y, X, B, PRNT, NPRNT, N). The dimensioned variables are Y(500),X(500,10),B(50),PRNT(5). Variables used in more than one subroutine must be linked through COMMON statements. As described under Output (SECTION IX), as many as five words of auxiliary information (e.g., values of independent variables, etc.) may be printed alongside the Y.sub.1, Y.sub.1, (Y.sub.1 - Y.sub.1) values. NPRNT is the number of such auxiliary words (NPRNT .ltoreq.5). The program assumes NPRNT = 0 unless SUBZ sets NPRNT .noteq.0.

Iv f coding

The FORTRAN coding to evaluate the function F is called for by CALL F CODE (Y, X, B, PRNT, F, I) with the dimensioned variables being Y(500), X(500,10), B(50), PRNT(5).

______________________________________ Mathematical Symbol FORTRAN Variable ______________________________________ X.sub.11 X(I,L) Y.sub.1 Y(I) Y.sub.1 F b.sub.1 B(J) ______________________________________

If NPRNT .noteq.0, then F coding must set up PRNT(JJ), JJ = 1, . . . , NPRNT, to contain the desired output information.

V. p coding

The FORTRAN coding to evaluate .differential.Y.sub.1 /.differential.b.sub.1, j = 1, . . . , k is called for by CALL P CODE (P, X,B, PRNT, F, I). The dimensioned variables are P(50), X(500,10), B(50), PRNT(5). The functions .differential.Y.sub.1 /.differential.b.sub.1 are labeled P(J). This subroutine must be included even if it is empty.

Vi. input Preparation

The input cards for a particular case are assembled in the following order, all items being required always except where indicated otherwise:

Input Mathematical FORTRAN Card Item No. Symbol Label Format Columns Comments __________________________________________________________________________ 1 n N 13 1-3 No. of data points k K 13 4-6 Total No. of parameters p IP 13 7-9 No. of held count parameters (IP-50) m M 13 10-12 No. of independent variables IFP 13 13-15 IFP = 000 to tabulate Y.sub.1, Y.sub.17, (Y.sub.1 = Y.sub.1), PRNT(1),...,PRNT(5) IFP = 001 to plot Y.sub.1, Y.sub.1 IWS 1 13 1-3 Doesn't apply IWS 2 13 4-6 = 0 Analytic Derivatives = 1 Estimated Derivatives IWS 3 13 7-9 = 0 Abbreviated Printout = 1 Detail Printout IWS 4 13 10-12 = 0 No Force Off = No Forced Branch to Confidence Region Calculation after No. Iterations IWS 5 13 13-15 = 0 Sense Switches Not Interrogated = 1 Sense Switches Interrogated IWS 6 13 16-18 = 0 Nonlinear Confidence Limits Desired = 1 Omit Nonlinear Confidence Limits YMN F10.0 1-10 Left side of plot SPRD F10.0 11-20 Spread of plot __________________________________________________________________________

(Item 3 is required only if IFP = 001; Format statement for item 3 is 930. Card for item 3 must be omited if IFP = 000.)

______________________________________ IB(1) 13 1-3 Subscripts of omitted b.sub.1 's. IB(2) 13 4-6 . . . . . . . . . IB(IP) . . ______________________________________

(item 4 is required only if IP > 0. Omit item if IP = 0. Format statement is 900; note that two cards are required if IP > 25. A zero or blank subscript will given an error mesage and cause program to go to next case.)

__________________________________________________________________________ F.sub.1.sub.-.sub..alpha. (k-p, n-k+p) FF F10.0 1-10 Variance ratio statistic t.sub.1.sub.-.sub..alpha. (n-k+p) T F10.0 11-20 Student's t .epsilon. E F10.0 21-30 Convergence criterion .tau. TAU F10.0 31-40 Convergence criterion .lambda..sup.(o) XL F10.0 41-50 Program parameter .gamma..sup.o GAMCR F10.0 51-60 Critical angle .DELTA. DEL F10.0 61-70 Used in finite-difference derivatives Z ZETA F10.0 71-80 Singularity criterion for matrix __________________________________________________________________________ inversion

Any or all of the quantities in Item 5 may be left blank on the card. If this is done, the program will supply the following values as being reasonable for most situations:

FF = 4.0, T = 2.0, E = 5 .times. 10.sup.-.sup.8, TAU = 10.sup.-.sup.8, XL = 10.sup.-.sup.2, GAMCR = .pi./4, DEL = 10.sup.-.sup.5, ZETA = 10.sup.-.sup.21

(Format statement is No. 931.)

______________________________________ b.sub.1.sup.(0) B(1) F10.0 1-10 Initial guess for b.sub.2.sup.(O) B(2) F10.0 11-20 parameters, 7 per card etc. etc. F10.0 thru 61-70 ______________________________________

Input Item No. Comments ______________________________________ 7 This item is a single card containing the format statement according to which the data (Item 8) are to be read. Column 1 contains the opening parenthesis of the format statement; the statement must end at or before column 60. 8 This item consists of n sub-items. Each sub-item is the input data for one observation, punched according to the format statement in item 7. The sequence of the variables in a sub-item must be; Y(I) X(I,1) X(I,2) through X(I,M). Each Y(I) must begin on a new card. As supplied, the program has the dimensions N = 500, K = 50, M = 10. These limits can readily be altered by changing the dimension statements. 9 Any case data read in from subroutine SUBZ should go here. 10 Sequential cases may be stacked by repeating Items 1 through 9. A blank card after the last case will cause a normal stop without error message. Without the blank card a systems stop indicating end of file on data tape will occur. ______________________________________

Vii. operating Procedure

Assemble the main program deck (FORTRAN source deck or binary deck), as well as subroutines SUBZ, F CODE, P CODE according to IBSYS procedure. The subroutines must all be included even if they are empty (i.e., contain only a DIMENSION statement and a RETURN statement). Data for as many cases as desired can follow.

Note: If k > 25,the program automatically uses FORTRAN logical Tape 3 for intermediate storage. Thus if k > 25 a scratch tape must be mounted on Tape 3 for use during execution of the object program. (See "Comment," Section XI.) VIII. Sense switches

If IWS 5 = 0, no sense switches will be interrogated; if the user is not present at run time he can thereby protect against operator error.

If IWS = 1, the program will stop on PAUSE 5. The Sense Switches can now be set, and will be interrogated.

______________________________________ Switch ON OFF ______________________________________ 1 Detailed Output No Detailed Output on On-Line Printer on On-Line Printer 2 Estimated Derivatives Analytic Derivatives 3 Detailed Printout Abbreviated Printout on Output Unit on Output Unit 4 Forced Branch to Confidence Region Calculations 5 Forced Branch to Next Case ______________________________________

Clearly, these switch settings imply certain assignments of output units in the IBSYS system. The program refers to file 6 for the system output unit and to file 12 for the on-line printer.

Ix. output

The program initially prints

N= K= P= M= IFP= GAMMA CRIT= DEL= FF= T= E= TAU= XL= ZETA=

a. With IWS 3 = 1, or IWS 5 = 1 and Sense Switch 3 ON, and IFP = 000, detailed printout is obtained at each iteration in the following format:

PARAMETERS B(1) B(K) OBS PRED DIFF PRNT(1) PRNT(2) PRNT(3) n of these; format PRNT(4) PRNT(5) statement is 925 PHI SE LAMBDA ANALYTIC PARTIALS USED INCREMENTS DB(1) ... DB(K) .phi.(.gamma.) PHI LAMBDA GAMMA LENGTH INCREMENTS DB(1) ... DB(K) .phi.(.gamma./10) PHI LAMBDA GAMMA LENGTH

(Similar printouts for .PHI. (10.lambda.), etc., if needed.) (LENGTH is the length of the DB(J) vector.)

If IFP = 001, then the lines containing Y(I), F, etc. are replaced by a plot showing observed (O) and predicted (P) values for visual monitoring of the fit. See Section X for details.

b. With Sense Switch 3 OFF, or IWS 3 = 0, the abbreviated printout at each iteration is in the format:

PARAMETERS B(1) ... B(K) PHI SE LENGTH GAMMA LAMBDA ANALYTIC PARTIALS USED

c. The program may branch to the confidence-region calculations by any of four routes. See Exhibit A for details. One of the following messages will appear on the printout to identify the route taken:

Message Remarks ______________________________________ EPSILON TEST Standard convergence route GAMMA LAMBDA TEST Alternate convergence route GAMMA EPSILON TEST Alternate convergence route FORCE OFF Via IWS 4 or Sense Switch 4 ______________________________________

After this message, the following output is printed:

N= K= P= M= FF= T= E= TAU=

One complete iteration of detailed printout is printed, either with tabulation of Y(1), F, etc. (IFP = 0 ) only, or with plotting (IFP > 0) as well.

Ptp inverse

(listed by row, 5 columns at a time. Error message will occur if any diagonal element is negative.)

PARAMETER CORRELATION MATRIX (Listed by row, 10 columns at a time.)

One-parameter and support-plane linear confidence intervals, along with the standard errors of the b.sub.1, in the format:

STD ONE PARAMETER SUPPORT PLANE B ERROR LOWER UPPER LOWER UPPER NONLINEAR CONFIDENCE LIMITS PHI CRITICAL = PARA LOWER B LOWER PHI UPPER B UPPER PHI

X. plotting Option

The plotting option is obtained whenever IFP .gtoreq. 001 in Input Item 1. The plotting option replaces the tabulation of Y(I), F, (Y(I) - F), PRNT (1), etc., with a plot showing for ech data point (each line of output) the observed and predicted values. If the data points as supplied to the computer are ordered with respect to an independent variable, X.sub.1, and are equally spaced with respect to that independent variable, then the plot option will provide a good visual portrayal of the curve of Y.sub.1 and Y.sub.1 versus X.sub.1. If the ordering of the data points is otherwise, the plotting option will still provide a qualitative feel for the goodness of fit, since the dependent variable is always accurately represented. An example of plotting-option output is shown in FIG. 1.

The plots are accurate to 1 part in 100. Specifically, the user must supply

Ymn = y value at left end of plotting area

Sprd = range of Y values to be covered in plot

The plotting area is 101 character positions wide beginning at position 2. A point, PT, is plotted from the (II + 2) position, according to the formula ##EQU35## The symbols which may be plotted are

Symbol Meaning ______________________________________ O Observed value of Y P Predicted value of Y Y Both observed and predicted values occupy the same plotting position X At left edge of plotting area indicates that II<0 At right edge of plotting area indicates that - II ______________________________________ >100

At the top of the plot a + symbol is printed at the left (II =0) and the right (II = 100) edges of the plotting area. Just above these symbols are printed YMN and (YMN + SPRD) respectively as a record of the scale used for plotting.

Xi. program Listing

Attached is a complete listing of the FORTRAN statements which make up the program.

Comment: If it is necessary to change program dimensions to accommodate a particular machine or problem, the criterion dimension (Cards 111, 339, 463, 572, 636) for writing matrix A on logical tape 3 must always be .ltoreq.1/2 (Number of rows given in dimension statement for matrix A).

REFERENCE

1. Marquardt, D. W., "An Algorithm for Least-Squares Estimation of Nonlinear Parameters," J. Soc. Indust. and Appl. Math., 11, No. 2, (1963) 431-441.

TABLE I __________________________________________________________________________ LABELS USED BY PROGRAM FORTRAN Label Mathematical Symbol* Meaning __________________________________________________________________________ BA(J) Working storage for b.sub.1 's BC Working storage B(J) b.sub.1 Value of jth parameter BL Lower nonlinear confidence limit BS(J) Working storage for b.sub.1 's BU Upper nonlinear confidence limit CGAM cos.gamma. D D Used in nonlinear confidence DEL .DELTA. Multiplier used in finite-difference derivatives DD Working storage DB(J) K.delta..sub.1 Increment to b.sub.1 DBW Working storage for b.sub.1 in est. partials DTD .SIGMA..sub.1 .delta..sub.1.sup.2 DTG .SIGMA..sub.1 .delta.g.sub.1 E .epsilon. Convergence criterion F .gamma..sub.1 Predicted value of dependent variable at ith data point FF F1-.alpha.(k-p,n-k+p) Variance Ratio Statistic FWS Working storage for Y in est. partials GAMCR .gamma..sub.o Critical angle G(I) g.sub.1 Right-hand side of normal equations GTG .SIGMA..sub.1 G.sub.1.sup.2 HJTD Half-length of support-plane confidence interval I Counter; during F-coding and P-coding I=data point number __________________________________________________________________________ *See Reference 1.?

TABLE 1-2 __________________________________________________________________________ FORTRAN Label Mathematical Symbol Meaning __________________________________________________________________________ IBCH Symbol for a blank column on card or printer IB(J) Subscript of an omitted b.sub.1 IBKA Return switch for main program IBKM Return switch for matrix inversion IBKN Switch showing upper or lower nonlinear conf. limits IBKP Indicator showing which nonlinear conf. limits have not been found IBKT Switch showing use of A or Tape 3 IBK1 Return switch from DB(J) Calculation IBK2 Return switch from PHI Calculation IBOUT Flag indicating bad data IFP Plotting option indicator II Counter III Counter IM1 I --1 IOCH The symbol O (alphabetic) IP .rho. No. of omitted parameters IPCH The symbol P IPP No. of blanks preceding character P in plot IP1 First symbol plotted IP2 Second symbol plotted ITCT Switch showing if first iteration - or not IWHER A logical indicator to show calculations to be performed and where to return to in convergence program IWS Working storage IXCH The symbol X IYCH The symbol Y IIMI No. of blanks before first symbol plotted __________________________________________________________________________

TABLE 1-3 __________________________________________________________________________ FORTRAN Label Mathematical Symbol Meaning __________________________________________________________________________ IIM2 No. of blanks between plotted symbols JJ A counter JP1 j + 1 K k No. of parameters KIP K-IP L A counter M m No. of independent variables MSING Indicator for singular matrix N n No. of data points NPRNT No. of words of auxiliary info. to be printed with Y.sub.1, Y.sub.1, (Y.sub.1 -Y.sub.1) OPL One parameter lower limit OPU One parameter upper limit P(J) Y.sub.1 /.differential.b.sub.1 PC .phi..sub.c Phi critical PHI .phi. Sum of squares of residuals PHID .phi.(D) PHIL .phi.(.gamma.) PHIL 10 .phi.(.gamma./10) PHIT 10 .phi.(10.gamma.) PHIZ .phi.(O) PKN Working storage = (K-IP)/ (N - K + IP) PL Phi for lower nonlinear confidence limit PRNT(1) Auxiliary info. for printing with . Y.sub.1, Y.sub.1, (Y.sub.1 = Y.sub.1) . PRNT(5) PU Phi for upper nonlinear confidence limit SA(I) Working storage __________________________________________________________________________

TABLE 1-4 __________________________________________________________________________ FORTRAN Label Mathematical Symbol Meaning __________________________________________________________________________ SE s.sub.e Sp'd error of estimate SPL Support plane lower limit SPU Support plane upper limit SPRD Spread of plot STE St'd error of b.sub.1 T t.sub.1.sub.-a a(n-k+p) Student's TAU .tau. Constant used in convergence test TWS Working storage WS Working storage WS1 Working storage WS2 Working storage X(I,L) X.sub.11 Value of lth independent variable at ith data point* XKDB K Multiplier of .delta. XK1 K.sub.1 Used in nonlinear confidence interval calculation XK2 K.sub.2 " XK3 K.sub.2 " XL .gamma. Lambda XLL Length of DB(J) vector XLS Working storage for .gamma. Y(I) Y.sub.1 Observed value of dependent variable at ith data point YMN Left side of plot ZETA Singularity criterion for matrix inversion __________________________________________________________________________ *Note: The order of the i and j subscripts in this program is interchange relative to the notation of Reference 1.

REFERENCES

B1. booth, G. W., G. E. P. Box, M. E. Muller, and T. I. Peterson, "Forecasting by Generalized Regression Methods, Nonlinear Estimation (Princeton-IBM)" Feb. 1959, International Business Machines Corp. Mimeo (IBM Share Program No. 687 WL NLl).

B2. box, G. E. P. and G. A. Coutie, Proc. Inst. Elec. Engrs. 103 Part B, Suppl. No. 1 (1956) Paper No. 2138 M.

B3. marquardt, D. W., R. G. Bennett, and E. J. Burrell, "Least Squares Analysis of Electron Paramagnetic Resonance Spectra," Jour. Molec. Spectroscopy 7, 269-279 (1961).

B4. scheffe, H., "The Analysis of Variance," John Wiley and Sons, New York, 1959.

B5. stone, H., Discussion on paper by E. M. L. Beale, J. Roy, Statist. Soc., Ser. B., 22 (1960), p. 41 ff. See espec. pp. 84,85. ##SPC2## ##SPC3## ##SPC4## ##SPC5## ##SPC6##

Claims

What is claimed is:

1. Apparatus for controlling the operation of a processing system including a plurality of individual processing units supplied by a plurality of feed streams of materials being transformed by the units into a plurality of product streams flowing from the units, the product yield and the product cost or profit for each unit being dependent upon the operating conditions of the unit, said apparatus comprising:

feed stream computer means for generating initial feed stream signals representing the initial feed stream flow rates for a given yield and cost at each of the units under a given set of operating conditions so as to maximize the profitability of the overall system;

feed stream control means coupled to said feed stream computer means for controlling the feed stream flow rates in response to said initial feed stream signals;

unit yield and cost computer means for repeatedly and continuously generating new yield and cost signals representing new yields and costs at one of the units corresponding to an increase in the profitability of the overall system;

unit control means coupled to said unit yield and cost computer means for controlling the operating conditions at said one of said units in response to said new yield and cost signals; and

another feed stream computer means coupled to said unit yield and cost computer means for repeatedly and continuously generating new feed stream signals representing new feed stream flow rates in response to said new yield and cost signals;

said other feed stream computer means being coupled to said feed stream control means for controlling the feed stream flow rates in response to said new feed stream signals.

2. The apparatus of claim 1 wherein said feed stream computer means and said other feed stream computer means comprise an overall system computer and said yield and cost computer means comprises a local computer located at said one unit and remote from said system computer, said system including means for transmitting said new yield and cost signals from said local computer to said system computer.

3. The apparatus of claim 1 wherein said feed stream computer means comprises an overall system computer and said yield and cost computer means and said new feed stream computer means comprises a local computer located at said one unit and remote from said system computer, said apparatus including means for transmitting said new feed stream signals from said local computer to said system computer.

4. The apparatus of claim 1 wherein said feed stream computer means, said new feed stream computer means and said yield and cost computer means comprises an overall system computer, said apparatus comprising means for transmitting said new yield and cost signals from said system computer to said unit control means.

5. The apparatus of claim 1 wherein said feed stream computer means includes means for computing the initial feed stream flow rates in accordance with the linear program model having an overall system objective function of maximizing the profitability

Z = c.sub.1 x.sub.1 + . . . . c.sub.n x.sub.n

subject to

p.sub.1 x.sub.1 + . . . . p.sub.n x.sub.n = Q

x.sub.i .gtoreq. 0 for all i

where p.sub.j = ##EQU36## c.sub.j is the marginal profit coefficient of the jth unit, x.sub.i is the flow rate of the ith stream, p.sub.j is the yield column of the jth unit and a function of operating conditions, Q is the demand constraint column, a.sub.ij is the yield coefficeint of the ith feed stream producing the jth product stream unit, and b.sub.i is the demand coefficient of the ith product stream.

6. The apparatus of claim 5 wherein said unit yield and cost computer means includes means for automatically computing a local objective function for said one unit, said local objective function defining changes in product yield and product cost or profit for increasing the profitability of the overall system.

7. The apparatus of claim 6 wherein said means for automatically computing the local objective function computes the local objective function

f.sub.m = (.DELTA.c.sub.m - .pi..degree..DELTA. p.sub.m) > 0

where .DELTA.c.sub.m is the change in the marginal profit coefficient for the mth unit, .DELTA.p.sub.m is the change in the mth yield column for the mth unit and .pi..degree. are the simplex multipliers [c.sub.1 ...c.sub. n ].beta..degree. at the previous operating conditions where .beta..degree. is the inverse of the basis matrix in [p.sub.1...p.sub. n ].

8. The apparatus of claim 7 wherein said other feed stream computer means includes means for computing said new feed stream flow rates ##EQU37##

9. An apparatus for controlling the operation of a processing system including a plurality of individual processing units supplied by a plurality of feed streams of materials being transformed by the units into a plurality of product streams flowing from the units, a plurality of local control means associated respectively with the units, each of the local control means including condition control means responsive to condition control signals generated by the local control means for controlling the operating of each unit, and central control means including flow rate control means associated with each of the feed streams and responsive to flow rate control signals generated by the central control means, the product yield and the product cost or profit of each unit being dependent upon the operating conditions of each unit, said apparatus comprising:

means for automatically computing the initial feed stream flow rates for a given product yield and a given product cost at each of the units so as to maximize the profitability of the overall system, said central control means generating initial flow rate control signals in response to the computed initial feed stream flow rates and applying the initial flow rate control signals to the flow rate control means to establish the initial feed stream flow rates;

means for substantially continuously and automatically computing a new product yield and new product cost or profit for at least one of the units corresponding to an increase in the profitability of the overall system, said local control means associated with said one unit substantially continuously generating new condition control signals corresponding to the new computed product yield and the new product cost or profit and applying the new condition control signals to the condition control means associated with said one unit to establish new operating conditions; and

means for substantially continuously and automatically computing new feed stream flow rates for all of the feed streams in response to the new product yield and the new product cost or profit for said one unit so as to maximize the profitability of the overall system, said central control means substantially continuously generating new flow control signals representing new feed stream flow rates and applying the new flow rate signals to the flow rate control means of the feed streams to establish the new feed stream flow rates.

10. The apparatus of claim 9 wherein

said means for automatically computing the new product yield and the new product cost or profit are located at the local control means associated with the unit; said apparatus further comprising means for generating signals representing the new product yield and the new product cost or profit at the local control means; and

means for transmitting the signals representing the new product yield and the new product cost or profit to said central control means.

11. The apparatus of claim 9 further comprising:

means for substantially continuously and automatically computing local objective functions for each of the units at the local control means, said local objective functions defining changes in product yield and product cost or profit for increasing the profitability of the overall system;

said means for substantially continuously and automatically computing the new product yield and the new product cost or profit being located at the local control means and being responsive to and consistent with the local objective function of said one unit.

12. The apparatus of claim 11 wherein

said means for automatically computing feed stream flow rates operates so as to satisfy and be consistent with a linear program model having an overall objective function of maximizing the profitability

Z = c.sub.1 x.sub.1 + .... c.sub.n x.sub.n

subject to

p.sub.1 x.sub.1 + .... p.sub.n x.sub.n = Q

x.sub.i .gtoreq. 0 for all i

where p.sub.j = ##EQU38## c.sub.j is the marginal profit coefficient of the jth unit, x.sub.i is the flow rate of the ith stream, p.sub.j is the yield column of the jth unit and a function of operating conditions, Q is the demand constraint column, a.sub.ij is the yield coefficient of the ith stream for the jth unit, and b.sub.i is the demand coefficient of the ith stream;

said means for substantially continuously and automatically computing the local objective functions operates so as to satisfy and be consistent with the equation

f.sub.m = (.DELTA.c.sub.m - .pi..degree. .DELTA. p.sub.m) > 0

where .DELTA.c.sub.m is the change in the marginal profit coefficient for the mth unit, .DELTA.p.sub.m is the change in the mth yield column for the mth unit, .pi..degree. are the simplex multipliers at the previous operating conditions; ; and

said means for substantially continuously and automatically computing the new feed stream flow rates operates so as to satisfy and be consistent with the equation ##EQU39## where s.sub.m is the mth sensitivity coefficient in the sensitivity coefficients s = .beta..degree. p.sub.m and .beta..degree. is the inverse of the base matrix .beta. formed of the columns p.sub.j.

13. The apparatus of claim 9 further comprising:

means for substantially continuously and automatically computing local objective functions for each of the units at the central control means, said local objective functions defining changes in product yield and product cost or profit for increasing the profitability of the overall system; and

said means for substantially continuously and automatically computing the new product yield and the new product cost or profit being located at the central control means and being responsive to and consistent with the signals representing the local objective function of said one unit

14. Apparatus for controlling a processing system including a plurality of individual processing units supplied by a plurality of feed streams of materials being transformed by the units into a plurality of product streams flowing from the units, said system being characterized by a linear program model having an overall objective function of maximizing the profitability

Z = c.sub.1 x.sub.1 + .... c.sub.n x.sub.n

subject to

p.sub.1 x.sub.1 + .... p.sub.n x.sub.n = Q

where c.sub.j is the marginal profit coefficient of the coefficient jth unit, p.sub.j is a column vector of yield coefficients for the jth unit and the function of operating conditions of said jth unit, Q is a column vector of demand constraints on the product streams and x.sub.i is the flow rate of the ith feed stream, said system comprising:

unit control means for controlling the operating conditions of said individual units;

first computing means coupled to said unit control means for computing yield coefficients and profit or cost coefficients for each jth unit at a given set of operating conditions;

second computing means for computing the flow rate x.sub.i for each said ith stream so as to maximize the profitability Z for the computed yield coefficients and profit or cost coefficients for each said jth unit at said given set of operating conditions;

feed stream control means coupled to said second computing means for controlling the flow rate of materials in said feed streams; and

third computing means coupled to said first computing means for computing a new flow rate x.sub.i.sup.1 for each said ith feed streams so as to maximize the profitability Z for a change in operating conditions at the mth one of said units, said change in operating conditions being represented by a substituted vector column p.sub.m.sup.1 for said mth unit and a substituted cost coefficient c.sub.m.sup.1, said third computing means comprising:

first storage means for storing signals representing each flow rate x.sub.i .degree. at a given set of operating conditions;

second storage means for storing signals representing the change in yield coefficients .DELTA.p.sub.m and a change in cost coefficients .DELTA.c.sub.m at said mth unit;

third storage means for storing signals representing the inverse basis matrix .beta..degree. of the basis matrix .beta..degree. of the vector columns p.sub.j .degree. at the given set of operating conditions;

means coupled to said second storage means and said third storage means for generating signals representing s = .beta..degree..DELTA. p.sub.m ; and

means coupled to said first storage means and said means for generating said s signals for generating signals representing ##EQU40## where s.sub.m is the sensitivity coefficient corresponding to said mth unit and x.sub.m .degree. is the feed stream flow to said mth unit;

said third computing means coupled to said feed stream control means so as to control said feed stream flow rate in response to said signals representing ##EQU41##

15. The apparatus of claim 14 further comprising fourth computing means coupled to said first computing means for determining the change in operating conditions for said mth unit so as to satisfy the local objective function f.sub.m = c.sub.m - .pi..degree. p.sub.m thereby maximizing the overall profitability of the system.

16. Apparatus for controlling a processing system including a plurality of individual processing units supplied by a plurality of feed streams of materials being transformed by the units into a plurality of product streams flowing from the units, said system being characterized by a linear program model having an overall objective function of maximizing profitability

Z = c.sub.1 x.sub.1 + .... c.sub.n x.sub.n

subject to

p.sub.1 x.sub.1 + .... p.sub.n x.sub.n = Q

where c.sub.j is the marginal profit coefficient of the jth unit, p.sub.j is a convector of yield coefficients for the jth unit and a function of operating conditions of said jth unit, Q is column vector of demand constraints on the product streams and x.sub.i is the flow rate of the ith feed stream, said apparatus comprising:

means for controlling the operating conditions at the local unit so as to correspond to column vectors p.sub.j for a given set of operating conditions at each unit;

means for controlling the feed stream flow rates x.sub.i to the local units so as to satisfy the overall objective function Z; and

means for computing the new feed stream flow rates x.sub.i.sup.1 to replace old feed stream flow rates x.degree..sub.i after a change in operating conditions at the mth unit resulting in the substitution of a new yield column vector p.sub.m.sup.1 for an old yield column vector p.sub.m.sup.1 and a new cost coefficient c.sub.m.sup.1 for an old cost coefficient c.degree..sub.m comprising:

first storage means for storing signals representing each flow rate x.degree..sub.i at a given set of operating conditions;

second storage means for storing signals representing the change in yield coefficients .DELTA.p.sub.m and a change in cost coefficients .DELTA.c.sub.m at said mth unit;

third storage means for storing signals representing the inverse basis matrix .beta..degree. of the basis matrix .beta..degree. of the vector columns p.degree..sub.j at said set of operating conditions;

means coupled to said second storage means and said third storage means for generating signals representing x = .beta..degree..DELTA.p.sub.m ; and

means coupled to said first storage means and said means for generating said s.sub.i signal for generating signals representing ##EQU42## where s.sub.m is the sensitivity coefficient corresponding to said mth unit and x.sub.m is the feed stream flow to said mth unit;

said third computing means coupled to said feed stream control means so as to control said feed stream flow rate in response to said signals representing ##EQU43##

17. The apparatus of claim 13 wherein

said means for automatically computing feed stream flow rates operates so as to satisfy and be consistent with a linear program model having an overall objective function of maximizing the profitability

Z = c.sub.1 x.sub.1 + .... c.sub. x.sub.n

subject to

p.sub.1 x.sub.1 + .... p.sub.n x.sub.n = Q

x.sub.i .gtoreq. 0 for all i

where p.sub.j = ##EQU44## c.sub.j is the marginal profit coefficient of the jth unit, x.sub.i is the flow rate of the ith stream, p.sub.j is the yield column of the jth unit and a function of operating conditions, Q is the demand constraint column, a.sub.ij is the yield coefficient of the ith stream for the jth unit, and b.sub.i is the demand coefficient of the ith stream;

said means for substantially continuously and automatically computing the local objective functions operates so as to satisfy and be consistent with the equation

f.sub.m = (.DELTA. c.sub.m - .pi..degree. .DELTA.p.sub.m) > 0

where .DELTA.c.sub.m is the change in the marginal profit coefficient for the mth unit, .DELTA.p.sub.m is the change in the mth yield column for the mth unit, .pi..degree. are the simplex multipliers at the previous operating conditions; and

said means for substantially continuously and automatically computing the new feed stream flow rates operates so as to satisfy and be consistent with the equation ##EQU45## where s.sub.m is the mth sensitivity coefficient in the sensitivity coefficients s = .beta..degree. p.sub.m and .beta..degree. is the inverse of the base matrix .beta. formed of the columns p.sub.j.

18. The apparatus of claim 1 wherein said unit yield and cost computer means includes means for automatically computing a local objective function for said one unit, said local objective function defining changes in product yield and product cost or profit for increasing the profitability of the overall system.