Optimal Simulation-Based Operations Planning ... - ACS Publications

2364

Ind. Eng. Chem. Res. 1999, 38, 2364-2374

Optimal Simulation-Based Operations Planning with Quantitative Safety Constraints Steven P. Asprey,* Rafael Batres, Tetsuo Fuchino, and Yuji Naka Process Systems Engineering Research Laboratory of Resources Utilization, Tokyo Institute of Technology, Midori-ku, Nagatsuta, Yokohama 226, Japan

A new two-layer optimization scheme is presented for operating procedure synthesis in the presence of quantitative safety path contraints. Layer 1 optimizes the overall operations time by varying the “operator time constant”sthe time that an operator is applied to the system. Layer 2 performs the path-constrained optimization by minimizing the sum-of-squared errors between the current system state vector and a prespecified goal state vector by applying the system operators. New system states are calculated using dynamic simulation of the equations describing the system. A detailed example involving mixture control of potentially flammable chemicals is described, illustrating all aspects of the proposed methodology. Introduction Operating procedure synthesis (OPS) can be viewed as searching for a set of sequenced primitive operations that transform a system from an initial state to a prespecified goal state through a series of intermediate states or path. These primitive operations must be carried out in such a way that no violations are made of any relevant operational safety constraints along the path due to either process equipment or process chemistry (i.e., the plan is “physically feasible”1). Such operating procedures are then used for plant start-up, normal operations, or shutdown. In light of this, an operational design methodology has been developed,2 of which operating procedures are an integral part. Classical work on OPS draws on developments in the artificial intelligence (AI) field. Either linear planning3,4 or nonlinear planning1,5,6 approaches have been applied to chemical engineering problems. Given the fact that AI search mechanisms are well-known to suffer from combinatorial explosion in larger-sized problems, the nonlinear planning approach is conceptually more appealing, as it attempts to mitigate the combinatorial aspect of the problem. However, these past approaches are fundamentally limited in their application because they only handle qualitative constraints on the path between system states (e.g., “chemicals A and B cannot be mixed in the unit together”). In addition, as pointed out by Li et al.,7 operating procedure synthesis in the later 1980s focussed on plan feasibility without due consideration to optimizing time or cost. Only recently have researchers been focusing on incorporating such quantative constraints as “chemicals A and B should be mixed in such a way so as to avoid the flammability region of concentrations”.7-9 Li et al.7 introduced a nonlinear programming (NLP) approach to search for an optimal process operating path. However, in this attempt, the time of application of system operators (or the operator time constant as defined in this paper) is not explicitly accounted for. In * Author to whom correspondence should be addressed. E-mail: [email protected]. Tel.: +81(45)924-5271. Fax: +81(45)924-5270.

addition, the time horizon was not properly discretized to allow the NLP routine to guide the process along a possibly time-optimal and feasible path. Furthermore, solving the problem in the “continuous” domain may produce unrealistic system operator actions. For example, at some point during the optimization, the NLP routine may require one valve to be open 0.05% for 0.3 ssan action that is not physically realizable. Mathematical representation of quantitative safety constraints was not discussed by Li et al.7 In more recent work, Galań and Barton9 approach the OPS problem using a hybrid discrete/continuous framework, as an optimal control problem. However, as they are quick to point out, such an approach cannot guarantee a globally optimal solution, having to settle for plan feasibility only. Galań and Barton9 present an example application of their methodology using a mixture control problem very similar to the one presented here, involing a tank change over operation from oxygen to methane by using nitrogen as a diluent. In this work, a new approach is taken to explicitly account for time optimality, while maintaining realizable operator actions. For solving the problem, we discretize the time dimension, however, not in an overall a priori fashionsthe actual number of time steps are not known a priori. By discretizing the problem, we can represent the system using evolving “states”, each represented by a state vector. As operators are applied to the system, new state vectors describing the system evolving through time are calculated using dynamic simulation. A nested optimization scheme is proposed, consisting of two distinct layers: (i) layer 1sa layer that optimizes the overall operations time by adjusting the “operator time constant”sthe time duration for which an operator is applied to the system; (ii) layer 2sa layer in which a path-constrained optimization minimizes the difference between the current state vector and the goal state vector by using the system operators. In this fashion, the two objective functions are optimized in a sequential, but iterative, manner. The details of this algorithm are outlined in the paper. A simple, continuous, one-dimensional case in layer 1 is investigated, in which we use the same operator time constants for each system operator, that because

10.1021/ie980708a CCC: $18.00 © 1999 American Chemical Society Published on Web 05/11/1999

Ind. Eng. Chem. Res., Vol. 38, No. 6, 1999 2365

Figure 1. Mixture control problem showing a pressure vessel, three feed lines with valves, and an exit line with a valve.

of the sequential nature of the overall optimization technique are kept constant throughout the optimization in layer 2. Optimization techniques for layer 1 include simply bounded global NLP algorithms (operator time constant g a physically realizable valve time). For this continuous simulated annealing is discussed. Optimization techniques for layer 2 include AI search algorithmssincluding brute-force methods such as depthfirst or breadth-first searches and more elegant methods such as A*, iterative-deepening A*, and depth-first branch-and-bound.10,11 The A* technique is used here and will be discussed in detail. As an alternative, the entire problem could be cast into a mixed integer nonlinear programming (MINLP) problem and solved by using such global methods as outer approximation12 or nonlinear branch-and-bound.13 Implementation of the MINLP could be done by use of a “superstructure” problem representation. However, the superstructure approach requires a priori knowledge of the problem “depth”,14 allowing for the proper number of valve manipulations to attain the goal state. To complicate matters, the problem depth changes with changes in the operator time constant. Further refinements of the problem include specifying nonideal valve responses, including more integer levels of valve openings and including constraint dependence on environment variables such as temperature. Implementation of these refinements will be pointed out, but are left as a future exercise. A detailed example is described, illustrating all aspects of the proposed methods. Problem Definition A good example of OPS in the presence of quantitative safety constraints is the synthesis of a sequence of valve operations to control the concentrations of chemical species, while avoiding a flammability limit during a changeover or start-up/shut-down process. Such an example can be found in Galań and Barton9 for a N2/ O2/CH4 mixture changeover, and also in CCPS15 for an air/steam/propylene mixture during reactor shutdown for acrylic acid. In either case, the mixture control setup can be depicted as in Figure 1. Here, we use the CCPS example, for the case of a process start-up. Feed lines for air, steam, and propylene are present, each with its own control valve (assumed to behave ideally). The initial conditions for start-up are atmospheric pressure, with the vessel filled with air. For illustrative purposes, the quantitative safety path constraint and a fictitious optimization solution scenario are shown in Figure 2. The percent by volume of air within the vessel is given by 100-%steam-%propylene. The describing equations for dynamic simulation can be written in terms of partial pressures of the individual components as

Figure 2. Planning in the presence of quantitative safety constraints. Shown here is the flammability path constraint and an illustrative optimization solution path.

dpj RT ) (δ N - δENE) dt V i i

(1)

∑j pj ) PT

(2)

δi, δE ∈ {0,1}

(3)

where pj represents the partial pressure of component j, PT is the total system pressure, R is the universal gas constant, T is the temperature of the tank, Ni is the molar flow rate of the ith feed line (NE for the tank exit line), and δi is the binary integer representing the state of valve i. In this model, we have assumed isothermal conditions, ideal gas mixtures, and operation of the control valves in a bang-bang fashion. The vessel pressure is controlled at a prespecified setpoint using δE by a simple on/off controller, while feed-line pressures are assumed to be high enough such that no backflow through the system occurs. A slightly more detailed model for the N2/O2/CH4 mixture changeover system, derived in terms of mole fractions,

yj )

pj PT

(4)

can be found in Galań and Barton,9 where nonideal valve responses (opening and closing speeds) are also accounted for. Table 1 shows the system values used in this study. In this work, we adopt the discretization approach by using binary valve operators. In this framework, the OPS problem is most easily represented in a state space1,3 or problem space. The vector of component partial pressures represents a state. A state-space operator acting on the system represents a manipulation of a valve. In a state- or problem-space graph, the states of the space are represented by nodes of the graph and the operators by (directed) edges between nodes. The sequence of operations for establishing desired process behavior is then represented as a path between the initial and goal states. An illustration of this is shown in Figure 3. Representation of Quantitative Safety ConstraintssThe Least-Squares Spline Quantitative safety constraints in chemical processing systems are formed on the basis of experimental data,

2366 Ind. Eng. Chem. Res., Vol. 38, No. 6, 1999

the following form: g

s(x) )

∑ ciMi,k+1(x)

(7)

i)-k

Figure 3. Planning problem or state-space graph. Shown is a state vector of three component partial pressures. The depth level of each state in the state-space is shown in square brackets.

in which Mi,k+1(x) denotes the B-spline with knots λi, λi+1, ..., λi+k+1. The spline function s(x) becomes a single polynomial on [a, b] if the discontinuities of its kth derivative at the interior knots λq all vanish. Given the data values fq at the points xq, where q ) 1, 2, ..., m1, xq < xq+1, the goal is to determine a spline approximation s(x), while trying to find a compromise between the following two objectives: (i) The prescribed values fq should be fitted “closely enough”. (ii) The approximating spline should be “smooth enough”, in the sense that the discontinuities in its kth derivative are as small as possible. Now, in order to be able to formulate this criterion mathematically, both some measure of smoothness and some measure of closeness of fit need to be proposed. For the latter, simply take the sum-of-squared residuals: m1

δ(cj) )

Table 1. Mixture Control System Values system variable

value

T V Ni

500 K 50.0 m3 0.5 mol/s

(8)

As a suitable smoothing norm, η, that is, a measure of the lack of smoothness, Dierckx16 proposes g

collected under small-scale, well-controlled environments. Mathematical representation of these constraints to a satisfactory degree of accuracy can be difficult, especially in the case of noisy experimental data. Galań and Barton9 express the flammability limit or path constraint using a simplified logical proposition. However, it should be pointed out that such a representation tends to be too simplified. From the problem description depicted in Figure 2, it is obvious that an optimal OPS solution will occur when the mixture within the holding tank is kept close to the path constraint. Thus, the more accurately the boundary is represented, the closer to true optimality the solution will be. In this study, we have chosen to represent the quantitative safety constraints by using a least-squares B-spline representation of eq 23 below. The following derivations are based on those presented in Dierckx.16 Consider the strictly increasing sequence of real numbers:

a ) λ0 < λ1 < ... < λg < λg+1 ) b

g

∑ (fq - i)-k ∑ ciMi,k+1(xq))2 q)1

(5)

Then the function s(x) is called a spline of degree k on [a, b], with knots λi, where i ) 1, 2, ..., g, if the following conditions are satisfied: (i) in each interval [λi, λi+1], where i ) 0, 1, ..., g, s(x) is given by some polynomial of degree k or less. (ii) s(x) and its derivatives of orders 1, 2, ..., k - 1 are continuous everywhere in [a, b]. If we introduce the additional knots,

λ-k ) λ-k+1 ) ... ) λ-1 ) a, b ) λg+2 ) ... ) λg+k ) λg+k+1 (6) then every such spline has a unique representation of

η(cj) )

g

∑ ∑ ciai,q)2 (

(9)

q)1 i)-k

where

ai,q )

{

0 (-)k+1k!(λi+k+1 - λi) ′ Πi,k (λq)

if i < q - k - 1 or i > q if q - k - 1 e i e q

}

(10) and

Πi,k(t) ) (t - λi)(t - λi+1)‚‚‚(t - λi+k+1)

Then the approximation criterion can be formulated as

minimize η(cj) subject to the constraint δ(cj) e S where S is a given, nonnegative constant which will control the extent of smoothing and therefore is called the smoothing factor. Using a Lagrange method, the problem of finding the coefficients of a spline sp(x), defined for positive values of p, can be cast as the leastsquares solution of the system: g

∑ ciMi,k+1(x) ) fq

q ) 1, 2, ..., m1

i)-k

(11) 1

g

∑ ciai,q ) 0

xpi)-k

q ) 1, 2, ..., g

when p is given the value of the positive root of F(p) ) S with

Ind. Eng. Chem. Res., Vol. 38, No. 6, 1999 2367 m1

F(p) )

∑

(fq - sp(xq))2

(12)

q)1

Once a set of knots is found such that

F(∞) e S < F(0)

(13)

there exists a single spline sp(x) with these knots for which F(p) ) S. The value of p can be determined by an iterative method. The choice of knots tries to take into account the specific behavior of the function underlying the data. The chosen set of knots must satisfy condition 13. Thus, first determine the least-squares polynomial s0(x). If F0(∞) e S, this polynomial is a solution to the problem. However, usually F0(∞) > S, in which case we determine successive least-squares splines sgj(x) with increasing number of knots gj, until condition 13 is satisfied. More on the details of a particular knot selection and placement can be found in Dierckx’s paper.16 Least-squares B-splines retain the underlying, nonlinear trend of the experimental data, without succumbing to the experimental noise. Subsequently, these spline representations have continuous derivatives, making them ideal for use within an optimization framework. As outlined by Dierckx,16 all spline-based methods require the specification of the degree of polynomial, as well as the smoothing factor to control the trade-off between closeness of fit and smoothness of fit. If nothing is known about the variance of the data to be approximated, the smoothing factor must be determined by trial and error. Adaptive degree polynomial filters have been devised, using various statistical measures of fit (i.e., the F test) as feedback in the determination of the polynomial degree17 and smoothing factor. An implementation of the smoothing spline method discussed above can be found in the FORTRAN package DIERCKX on the Internet at NETLIB (www.netlib.org/dierckx). Documentation on both implementation and more algorithm specifics can be found in Dierckx.16 An additive safety margin to account for model inaccuracies could easily be incorporated into the final spline representation. Proposed Optimization Scheme In this work, a new approach is taken to explicitly optimize the total operations time, while maintaining realizable operator actions for obtaining a prespecified goal state of the system. For simplicity, binary operators of close valve (0) and open valve (1) are used. In the state-space representation, as operators are applied to the system, new state vectors describing the system evolving through time are calculated using dynamic simulation, with states including the three-component partial pressures. The optimization scheme and solution methods are discussed in detail here. A nested optimization scheme is proposed, consisting of two distinct layers: (1) Layer 1 optimizes the overall operations time by adjusting the “operator time constant”sthe time duration for which an operator is applied to the system. Objective function values are supplied by the final solution resulting from layer 2. (2) Layer 2 is where a constrained optimization minimizes the difference between the current state vector and the goal state vector by using the system

operators. Nonlinear quantitative safety constraints as shown in Figure 2 are applied in this layer. In this fashion, the two objective functions are optimized in a sequential, but iterative, manner. From layer 1, a single value of the operator time constant is passed to layer 2 (and is held constant while the solution is determined in layer 2). Within layer 2, one complete discrete optimization is performed to bring the system to a prespecified goal state. In order to accomplish this, an undetermined number of system operators must be applied to the system, which in turn, results in a particular value for the overall time. Once the overall time is determined by the completion of the optimization in layer 2, it is passed back to layer 1 as the value for the objective function. An illustration of these concepts is depicted in Figure 4, while the algorithm is as follows: (i) Choose an initial value for ∆t in layer 1. (ii) Pass ∆t to layer 2 and attain the goal state by applying a sequence of operator actions to the system determined by the A* search technique. (iii) Once the goal state has been reached, compute the overall operations time (tf) from the sequence of operator actions. (iv) Pass the overall operations time (tf) to layer 1 for the objective function value. (v) Change ∆t via the simulated annealing method and return to step ii. (vi) Iterate steps ii-v to convergence. More flexibility in layer 1 could be added by allowing options for increasing/decreasing operator time constants as the goal state is approached in layer 2. This could be accomplished by using parametric polynomials that are functions of the components of the state vector (pj or PT). As a result, the layer 1 optimization would become multidimensionalsthe number of dimensions depending on the number of parameters in the polynomials. This is left as a future exercise. Layer 1sContinuous Simulated Annealing. The optimization problem in Layer 1 can be stated formally as follows:

min tf

(14)

subject to ∆t >

(15)

∆t

where tf is the total operational time to take the system from the initial state to the goal state, ∆t is the operator time constant, and is a predetermined realizable lower limit for the system operators and is a function of the performance capabilities of the control valves used in practice. In order to ensure a globally optimal solution to the OPS problem, continuous simulated annealing was chosen as the optimization method for layer 1. Several methods have been published on the use of simulated annealing for the minimization of multimodal, continuous objective functions.19-21 The implementations of Cardoso et al.19 and Press and Teukolsky20 both combine the discrete combinatorial simulated annealing method with the continuous-variable simplex optimizer of Nelder and Mead,22 while Corana et al.21 combined simulated annealing with a random search mechanism that operates along each coordinate axis. Given a terminating criterion ζ, a number of successive temperature reductions to test for termination Nζ, a test for step variation Ns, a step-varying criterion c, a test for temperature reduction NT, and a temperature reduc-


Figure 4. Two-layer optimization scheme.

tion coefficient rT, the algorithm presented by Corana et al.21 can be briefly described as follows: (i) Choose a starting point xi and evaluate the objective function fi ) f(xi). (ii) Generate a random point x* along one coordinate direction h as

x* ) xi + rvmheh

(16)

where eh is the vector of the hth coordinate direction, r is a random number generated in the range [-1, 1], and vmh is the component of the step vector vm along the same direction. (iii) If the hth coordinate of x* lies outside the definition domain of f (i.e., if xh* < ah or xh* > bh), then return to step ii. (iv) Compute f* ) f(x*); if (f* e fi) then accept the new point else (f* > fi) then accept or reject the point with acceptance probability p (the Metropolis move):

p ) exp

( ) fi - f* Tk

(17)

(v) Repeat for all coordinate directions in xi by going back to step ii. (vi) If the test for step variation fails (number of step variation iterations < NS), go back to step ii; else update the step vector vm according to the number of accepted points nu.

(

vmu ) vmu 1 + cu

)

nu/NS - 0.6 0.4

if nu > 0.6NS (18)

vmu )

vmu 0.4 - nu/NS 1 + cu 0.4 vmu ) vmu

if nu < 0.4NS

otherwise

(19)

(20)

(vii) If the test for temperature reduction fails (number of iterations at temperature Tk < NT), go back to step ii; else reduce the temperature to

Tk+1 ) rT × Tk

(21)

/ | e ζ, (viii) If the test for convergence |f /k - f k-u where u ) 1, ..., Nζ fails, go back to step ii; else stop.

Simulated annealing is a generalization of a Monte Carlo method for examining the equations of state and frozen states of n-body systems.18 The concept is based on the manner in which liquids freeze or metals recrystallize in the process of annealing. In an annealing process, a melt, initially at high temperature (T) and disordered, is slowly cooled so that the system at any time is approximately in thermodynamic equilibrium. As cooling proceeds, the system becomes more ordered and approaches a “frozen” ground state at T ) 0. Hence, the process can be thought of as an adiabatic approach to the lowest energy state. If the initial temperature of the system is too low or cooling is done insufficiently, slowly the system may become quenched, forming defects or freezing out in metastable states (i.e., trapped in a local minimum energy state). Starting at some high temperature T, a sequence of points is generated until a thermodynamic equilibrium is approached. In this


case, NT (the number of iterations before a temperature reduction) is used to give good sampling statistics for the current temperature. After equilibration, the temperature is then decreased and the entire process repeated until a frozen state is achieved at T ) 0. An implementation of the continuous simulated annealing optimization method by Goffe et al.23 based on the algorithm of Corana et al.21 can be obtained on the Internet at NETLIB (www.netlib.org/opt/simann.f). The success of the algorithm rests on finding an appropriate starting value for the temperature T and giving the routine a correct measure of the size of region of interest (vm). Comments on selecting a value of T are given in Goffe et al.23 More on the implementation of this algorithm for the specific case of mixture control will be presented in a subsequent section. Layer 2sA* Search. The optimization problem in layer 2 can be stated as follows:

min|pj0 - pjgoal|

(22)

subject to φ(p1(n), p3(n)) e 0

(23)

δi(n)

δ(n) ∈ {0,1} where pj and δi(n) are defined as before, and φ( ) represents the flammability constraint. The path constraint introduces nonconvexity into the objective function; thus a globally optimal AI search mechanism was chosen. Search is a universal problem-solving mechanism in artificial intelligence (AI). Although many brute-force search algorithms exist for finding a path from the initial state to the goal state, their time complexities grow exponentially with problem size. This is termed combinatorial explosion, and as a result, the size of problems that can be solved with these techniques is very limited.10 In order to solve larger and more practical problems, domain-specific knowledge must be added to improve search efficiencysin other words a heuristic evluation function is added to the search mechanism. For a fixed goal state, a heuristic evaluation is a function of a node, h(n), that estimates the distance from node n to the given goal state. The A* algorithm is a best-first search in which the cost associated with a node is f(n) ) g(n) + h(n), where g(n) is the cost of the path from the initial state to node n (and is given by -∑j(pj - pj0)2, and is negative so as to maximize the summation), while h(n) is the heuristic estimate of the cost of a path from node n to the goal (and is given by ∑j(pj - pjGOAL)2). This gives the A* search a “hill-climbing” ability. Thus, f(n) estimates the lowest total cost of any solution path going through node n. At each point a node with lowest f value is chosen for expansion. If two expanded nodes have equal f values, the node with a lower h value is chosen. The algorithm maintains a closed list of those nodes that have already been expanded and an open list of those nodes that have been generated but not yet expanded. The algorithm begins with just the initial state in the open list. At each cycle, a node on the open list with the minimum f(n) value is expanded, generating all of its children, and is placed on the closed list. The heuristic function is applied to the children, and they are placed on the open list in order of their heuristic values. The algorithm continues until a goal state is reached. The greatest benefit of the A* search is the fact that it finds a globally

Figure 5. Optimization solution showing the time-optimal trajectory around the safety path constraint.

optimal solution path to a goal if the heuristic function is admissible, meaning that it never overestimates the actual cost.11 To overcome the known, poor memory requirements of the A* search, a recursive minimin search can be used.10 The algorithm searches forward from the current state to a fixed depth. At the search horizon, the A* evaluation function f(n) ) g(n) + h(n) is applied to the frontier nodes. A single move is then made to the neighbor of the current state with the minimum value. Minimin’s memory requirement is much smaller than the original A* algorithm, is much easier to implement, and runs much faster than A* since it does not incur the overhead of managing the open and closed lists. In this study, a further adaptation of the A* search was implemented by using the augmented cost function f(n) ) w(n)g(n) + h(n) + P(n), where w(n) is an exponential decay function dependent upon the position of the current node n in the problem space and P(n) is a barrier function (Heavyside function) that is active when the flammability path constraint is active and is assigned a large, finite number. The function w(n) serves to relax the hill-climbing capability of the search as the solution progresses toward the goal state. The exponential decay term is problem-specific and was best when set to

w(n) ) exp(-n∆ty1GOALy2GOALy3GOAL)

(18)

where yjGOAL is the goal state value of the mole fraction of the jth component in the mixture. Results and Discussion Implementation of the two-layer optimization scheme was accomplished using Visual Basic and FORTRAN90. A graphical user interface (GUI) was built, in which the user is presented with a graph of the flammability path constraint using a B-spline representation, and the current solution path attempted by the optimization scheme. A screen capture of the GUI is presented in Figure 5, showing the final solution to the OPS problem in the context of the problem definition given in Figure 2. In order to maintain an interactive environment, the continuous simulated annealing calculations are performed in Visual Basic. The A* search mechanism and


Figure 6. Optimization surface for layer 2 showing the sum-ofsquared errors surface, the safety path constraint represented as a barrier function, and the true minimum located at (5.0, 11.0, 84.0).

code for dynamic simulation was implemented in a FORTRAN-90 Windows dynamic link library (DLL) to provide ample calculation speed. Two buttons exist on the main program form: (i) one for a complete optimization of the problem and (ii) one for a layer 2 optimization solution only, using the operator time constant manually entered by the user. A global optimizer was required in layer 2, as the presence of the flammability path constraint inherently introduces nonconvexities in the objective function. For example, inspection of Figure 5 shows that while starting at (p1,p2,p3) ) (0.0,0.0,0.0) with a goal state of (5.0,11.0,84.0), a local optimum in the sum-of-squared errors between the current state and the goal state will occur at (5.0,3.8,91.2), at the flammability constraint. This is clearly illustrated in the three-dimensional representation of the optimization surface in Figure 6, showing the true minimum at (5.0,11.0,84.0), while also exhibiting the local minimum at the safety path constraint. As shown, the safety path constraint is represented as a barrier function. A global optimization

algorithm with a hill-climbing capability has the potential to overcome this problem; however, nonglobal optimizers would not advance from this obviously suboptimal solution. With the augmented cost function described in the previous section, the A* search mechanism was able to find the global optimum for all realizable operator time constants. Depending on the operator time constant used, the A* algorithm took on average 20-30 s to reach a solution. Preliminary manual intervention using only the layer 2 optimization to provide values of tf showed that because of the discretization of the problem, the objective function for layer 1 turns out to be highly nonconvex. Figure 7 shows the objective function values for the optimization in layer 1 (where each value is given by a complete optimization performed in layer 2). As can be seen, many local optima exist; thus, obtaining a globally optimum solution required the use of a global optimization algorithm for continuous variables, such as the continuous simulated annealing methods described earlier. The implementation of the algorithm obtained from NETLIB has several parameters that must be set by the user a priori. These include the initial temperature for annealing (50.0), the temperature reduction factor (0.85), the number of iterations before temperature reduction (50), the error tolerance for convergence (1 × 10-5), the maximum number of iterations (2000), the upper and lower bounds on each independent variable (2.0 and 45.0, respectively), and a step length vector that must be initially set to encompass the region of interest of the independent variables (30.0). By using the values of the parameters shown in brackets, and performing a complete two-layer optimization, a solution was reached in 346 iterations in layer 1, with a final operator time constant (∆t) of 37.19 s and a final overall operations time (tf) of 1450 s. The total computation time was 2.9 h on a Pentium 266 MHz machine. Although simulated annealing has proven to be a very robust optimization method in the face of many local optima, improvement in solution time is most likely to come from use of other methods in layer 1 that require fewer function evaluations, as opposed to using a random

Figure 7. Optimization scenario for layer 1 showing the objective function values for the optimization problem in layer 1. As can be seen, the objective function is highly nonconvex, with a global optimum at 37.19 s for the operator time constant.


Figure 8. Optimal valve action time profiles showing the valve operations to give the time-optimal solution of Figure 5.

search mechanism that is known for requiring a high number of function evaluations. A graphical representation of the valve action time profiles for the time-optimal solution is shown in Figure 8. Inspection of Figure 8 shows that valve 3 oscillates in the time range 800-1200 s, possibly indicating that an optimal valve position would be 0.5 (or half-open). To investigate this further, we added a third valve level,

such that δi(n) ∈ {0, 0.5, 1}. Solution of such a problem increases the number of valve pattern positions (and thus the number of expanded nodes) to 33, as opposed to 23 in the binary case, for a single node in the statespace representation. As a result, the A* solution took more than twice as long to find a solution. Thus, for the sake of comparison, we did not carry out a layer 1 optimization, but instead selected the A* solution with


Figure 9. Optimal valve action time profiles with three valve levels showing the valve operations to give a comparable solution to that of Figure 5.

∆t ) 40 s, resulting in the valve action time profiles shown in Figure 9. As expected, the result of adding the third valve position improved the valve time profile for valve 3, not requiring as many valve manipulations. However, the

valve 2 time profile now shows some oscillatory behavior in the same time range. Again, further addition of integer valve levels may remedy the situation, with an asymptotic approach to the continuous case. In situations such as this, the hybrid approach of Galań and


Barton9 may prove to be beneficial if a globally optimal solution is not required and one can settle for plan feasibility. In a general case, where m valve levels or positions are used, there would be m3 child nodes for each node in the state-space. In practice, a finite upper limit would be reached for m in a computer-controlled environment. For cases where m is large, the present method would be able to find a solution; however, the computation time would become quite large. As an alternative, other elaborations of the A* technique, such as the iterative-deepening A* search algorithm, may tend to be more efficient and could be explored. Conclusions A new nested, two-layer optimization scheme has been proposed and tested for the purposes of OPS in the presence of quantitative safety path constraints. When the discretized formulation for such a problem was used, it was found that the objective functions in both layers were highly nonconvex. In order to ensure a globally optimal solution, continuous simulated annealing was used for layer 1, while a modified A* AI search mechanism was used for layer 2. The methods were demonstrated successfully on a problem in which the making of a mixture of potentially flammable chemicals was controlled in the presence of a flammability path constraint. The mathematical representation of the flammability path constraint was accomplished by using least-squares B-splines, while a barrier function method was used to ensure no constraint violations. Future refinements of the problem include specifying nonideal valve responses, including more integer levels of valve openings and constraint dependence on environment variables such as temperature. Future work in our laboratory will focus on the implementation and integration of these techniques in a plant wide situation, where these techniques can be used to aid global planners for safe plant start-up, steady-state operation, shutdown, and abnormal situation management.

NT ) test for temperature reduction for simulated annealing Nζ ) test for the number of successive temperature reductions for termination of the simulated annealing algorithm p ) positive root of F(p); Metropolis acceptance probability for simulated annealing P(n) ) barrier function (Heavyside function) pj ) partial pressure of the jth component PT ) total prssure R ) universal gas constant r ) random number for simulated annealing rT ) temperature reduction coefficient for simulated annealing S ) smoothing factor for spline determination s(x) ) spline function representation T ) temperature (both system temperature and simulated annealing temperature) tf ) final time or overall time of operations V ) vessel volume vmh ) hth coordinate direction of vm w(n) ) exponential decay function xq ) qth independent variable value for a spline approximation yj ) mole fraction of the jth component Vector Variables c ) step-varying criterion for simulated annealing eh ) vector of the hth coordinate direction x ) independent variable vector for simulated annealing vm ) step vector for simulated annealing Greek Variables ∆t ) operator time constant δi ) binary integer for the state of valve i ) lower limit for ∆t φ( ) ) functional representation of the flammability constraint λi ) ith knot of a spline δ( ) ) sum-of-squares function for spline determination η( ) ) smoothing norm for spline determination ζ ) terminating criterion for simulated annealing

Acknowledgment S. P. Asprey would like to acknowledge the Japan Society for Promotion of Science (JSPS) for funding. Nomenclature a, b ) lower and upper bounds for the range of independent variables for a spline representation of data and independent variables for simulated annealing ci ) ith weighting coefficient for spline representation f ) objective function value for simulated annealing fq ) qth data point for determining a spline approximation g(n) ) heuristic estimate of cost of the path from the initial node to node n h(n) ) heuristic estimate of cost of the path from node n to the goal node k ) degree of a spline m ) number of valve positions m1 ) number of data points for spline representation Mi,k+1 ) B-spline function n ) node number in the state-space representation nu ) number of points accepted by the simulated annealing algorithm Ni ) molar flow rate of the ith feed line NS ) test for step variation for simulated annealing

Literature Cited (1) Lakshmanan, R.; Stephanopoulos, G. Synthesis of Operating Procedures for Complete Chemical PlantssI. Hierarchical, Structured Modeling for Nonlinear Planning. Comput. Chem. Eng. 1988, 12, 985. (2) Naka, Y.; Lu, M. L.; Takiyama, H. Operational Design for Start-up of Chemical Processes. Comput. Chem. Eng. 1997, 21, 997. (3) Fusillo, R. H.; Powers, G. J. A Synthesis Method for Chemical Plant Operating Procedures. Comput. Chem. Eng. 1987, 11, 369. (4) Fusillo, R. H.; Powers, G. J. Operating Procedure Synthesis Using Local Models and Distributed Goals. Comput. Chem. Eng. 1988, 12, 1023. (5) Lakshmanan, R.; Stephanopoulos, G. Synthesis of Operating Procedures for Complete Chemical PlantssII. A Nonlinear Planning Methodology. Comput. Chem. Eng. 1988, 12, 1003. (6) Lakshmanan, R.; Stephanopoulos, G. Synthesis of Operating Procedures for Complete Chemical PlantssIII. Planning in the Presence of Qualitative Mixing Constraints. Comput. Chem. Eng. 1990, 14, 301. (7) Li, H. S.; Lu, M. L.; Naka, Y. A Two-Tier Methodology for Synthesis of Operating Procedures. Comput. Chem. Eng. 1997, 21S, S899. (8) Lee, J. J.; Norris, W. D., II; Fishwick, P. A. An ObjectOriented Multi-Model Design for Integrating Simulation and Planning Tasks. J. Sys. Eng. 1993, 3, 220.

2374 Ind. Eng. Chem. Res., Vol. 38, No. 6, 1999 (9) Galań, S.; Barton, P. I. Dynamic Optimization Formulations for Operating Procedure Synthesis. Presented at the AIChE Annual Meeting, Los Angeles, CA, Nov 1997. (10) Korf, R. E. Artificial Intelligence Search Algorithms. In The CRC Handbook of Algorithms and Theory of Computation; Atallah, M. J., Ed.; CRC Press: Boca Raton, FL, 1998. (11) McAllester, D. Graph Search and STRIPS Planning. Lecture Notes on Artificial Intelligence, available at http://www.research.att.com/∼dmac/notes/, 1994. (12) Kocis, G. R.; Grossman, I. E. Global Optimization of NonConvex MINLP Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27, 1407. (13) Ryoo, H. S.; Sahinidis, N. V. Gobal Optimization of NonConvex NLPs and MINLPs with Applications in Process Design. Comput. Chem. Eng. 1995, 19, 551. (14) Avraam, M. P.; Shah, N.; Pantelides, C. C. Modelling and Optimization of General Hybrid Systems in the Continuous Time Domain. Comput. Chem. Eng. 1998, 22S, S221. (15) Center for Chemical Process Safety (CCPS). Guidelines for Safe Process Operations and Maintenance; AIChE: New York, 1995. (16) Dierckx, P. A Fast Algorithm for Smoothing Data on a Rectangular Grid While using Spline Functions. SIAM J. Numer. Anal. 1982, 19, 1286. (17) Barak, P. Smoothing and Differentiation by an AdaptiveDegree Polynomial Filter. Anal. Chem. 1995, 67, 2758.

(18) Metropolis, N.; Rosenbluth, A.; Rosenbluth, M.; Teller, A.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087. (19) Cardoso, M. F.; Salcedo, R. L.; Feyo De Azevedo, S. The Simplex-Simulated Annealing Approach to Continuous Nonlinear Optimization. Comput. Chem. Eng. 1996 20, 1065. (20) Press, W. H.; Teukolsky, S. A. Simulated Annealing Optimization Over Continuous Spaces. Comput. Phys. 1991, 5, 426. (21) Corana, A.; Marchesi, M.; Martini, C.; Ridella, S. Minimizing Multimodal Functions of Continuous Variables With the Simulated Annealing Algorithm. ACM Trans. Math. Software 1987, 13, 263. (22) Nelder, J. A.; Mead, R. A Simplex Method for Function Minimization. Comput. J. 1964, 7, 308. (23) Goffe, W. L.; Ferrier, G. D.; Rogers, J. Global Optimization of Statistical Functions with Simulated Annealing. J. Econ. 1994, 60, 65.

Received for review November 9, 1998 Revised manuscript received March 18, 1999 Accepted March 19, 1999 IE980708A

Optimal Simulation-Based Operations Planning ... - ACS Publications

Recommend Documents