Hybrid Dynamic Optimization Using Artificial Chemical Process

8400

Ind. Eng. Chem. Res. 2006, 45, 8400-8412

Hybrid Dynamic Optimization Using Artificial Chemical Process: Extended LARES-PR Roberto Irizarry* DuPont Electronic Technologies, 14 T.W. Alexander DriVe, P.O. Box 13999, Research Triangle Park, North Carolina 27709-3999

This paper presents a methodology for the dynamic optimization of hybrid processes. The proposed methodology considers the most general class of hybrid dynamic systems consisting of a finite set of dynamic subsystems, in which sequence and duration is not established a priori. The objective is to determine the optimal sequence of dynamic subsystems, the transition times from one subsystem to another, and the optimal profile of the control laws. LARES-PR algorithm for the solution of dynamic optimization problems (Irizarry, R. Chem. Eng. Sci. 2005, 60, 5663-5681) is extended to solve this type of problem. The efficiency and robustness of the method is demonstrated with a set of challenging benchmark problems. The proposed methodology is then applied to a multicomponent nanoparticle alloy synthesis. This prototype model is a very challenging one in which a mechanism switch not changes not only the model structure but also the solution methodology for each mode. For instance, in one phase, the particles grown by coalescence are modeled using discretized population-balance models. When a mechanism switch is induced, a new type of particle is introduced into the system. The subsequent growth of the binary system is modeled using the Monte Carlo method. The proposed methodology can handle this situation effectively. In general, the method is very robust, is simple to implement, and can be used to solve a wide range of hybrid optimization instances with different types of complexities. 1. Introduction Most industrial processes are operated in the batch mode, combining periods of continuous operations with discrete changes (events). These events may be triggered when a logical condition is satisfied or may be imposed by design. An example of the first case is when a valve is closed in response to a temperature change in the reactor to avoid a runaway reaction. In the second case, the transition times and dynamic subsystems are determined as part of the optimal solution. For example, a process is designed with valves opening and closing optimally to maximize yield. In addition to chemical processes, many biological processes consist of hybrid dynamics with switching mechanisms as well (Ahtchi-Ali and Pedersenb1). The dynamic optimization of these hybrid systems is a very difficult task because of the discontinuities embedded in the continuous dynamics. Furthermore, most chemical and biological systems exhibit highly nonlinear dynamics described by complex models that need to be solved numerically, resulting in a large set of ordinary differential equations or even an equationless simulation like Monte Carlo methods. The literature on hybrid dynamic optimization is limited because of the mentioned difficulty, especially when deterministic optimization methods (Floudas2) are used to solve the optimization problem. Progress on deterministic methods is limited to small systems and problem types. For example, in the study by Manon et al.,3 a sequential controller that generates this hybrid sequence, which consists of discrete controllable events impulsed at precisely precalculated time instants, was developed using extended Pontryagin’s maximum principle. This strategy is limited to very small systems. Deterministic solutions of global optimization problems with hybrid systems emdedded have recently been proposed (Barton and Lee4 and Lee and Barton5). Different from previous methodologies, these algorithms guarantee the location of the global solution within some * Phone: (191) 248-5246. E-mail: [email protected].

user-supplied optimality tolerance. In the study by Barton and Lee,4 a deterministic branch-and-bound method was developed for the dynamic optimization of linear hybrid systems with fixed transition times and a fixed sequence of modes. The work was extended by Lee and Barton5 to determine the optimal mode sequence while keeping the transition times fixed. A second algorithm was developed to determine the optimal transition times while keeping the mode sequence fixed. The systems studied with these algorithms are linear models of small to medium size. In the study by Barton et al.,6 a gas system was solved successfully using the ICRS stochastic optimization algorithm. In this case, the state variable triggers the mechanism switch. This algorithm cannot be used for the case in which the switching from one mechanism to another is part of the decision variables because the ICRS algorithm only considers real variables as decision variables. In this work, LARES-PR (Irizarry7) is modified to consider the most general type of hybrid optimization problems. That is, problems with a set of control laws, transition times, and mode sequences to be optimally determined. LARES-PR is a general-purpose optimization algorithm to solve dynamic optimization problems.The method uses the LARES algorithm to solve the optimization problem. LARES is a global stochastic optimization algorithm based on a new paradigm called “artificial chemical process”.8 In general, stochastic optimization methods (LARES, genetic algorithms, evolutionary algorithms, simulated annealing, etc.) are direct-search methods in which the only information used during the search is the value of the objective function to be minimized. They have the capability of effectively escaping from local minima and eventually reaching a near-global optimum. Many of these algorithms are based on heuristics with convergence properties guaranteed only in a probabilistic sense. In practice, they have proven to be very robust in a wide range of applications. The objective of this work is to show that LARES-PR and the extended LARES-PR are effective tools to solve a wide

10.1021/ie060463z CCC: $33.50 © 2006 American Chemical Society Published on Web 11/08/2006

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8401

Figure 1. Schematic of LARES emulating a chemical process consisting of four unit operations: L ) load tank, AR ) activation-reaction reactor, E ) etranscio unit, and S ) separation unit.

range of hybrid dynamic optimization problems. After solving a set of benchmark problems, the methodology is applied to the hybrid optimization of nanoparticle synthesis. This is a novel application of the proposed methodology in which mechanism switches are induced by design to further control the particle size and distribution of a binary particulate system. The design considers two phases. In the first phase synthesis, A-type particles are grown by coalescence, then B-type particles are added to make the final product, an AB-alloy particulate system with a desired particle size distribution and composition. During the hybrid dynamic optimization, the optimal conditions for phase 1 and phase 2 together with the optimal transition time are determined as solutions to the optimization problem. In Section 2, the artificial chemical process paradigm is briefly reviewed. The algorithm is described in detail in Appendix A (see ref 8 for a more detailed description). In Section 3, the LARES-PR methodology for dynamic optimization introduced by Irizarry7 is reviewed. The hybrid dynamic optimization problem is stated and the nomenclature used is introduced in Section 4. In Section 5, the extended LARESPR is presented for the solution to hybrid dynamic optimization problems in its most general form. In Section 6, the performance of the current algorithm is studied with a set of case studies. The robustness of the algorithm is shown along with the importance of considering the problem in its most general form to find better designs. The production of a multicomponent particulate process is considered in Section 7. This challenging problem, solved for the first time, shows the power of the current methodology in which, in addition to the discrete event embedded in the continuous dynamics, part of the dynamics is described by a large set of ordinary differential equations (ODEs) while other parts are solved by an equation-free Monte Carlo method. 2. Concept of Artificial Chemical Process and LARES Algorithm The optimization algorithm used to solve the optimization problem is based on a new paradigm called “artificial chemical process” (see refs 7 and 8). This paradigm considers a chemical process like the one shown in Figure 1, in which a finite set of molecules are distributed between four sets (called L, AR, E, and S) representing “unit operations”. Each molecule can be in one of the possible finite states allowed (the set of possible states can be different for each molecule). Similar to a real chemical

process, a sequence of steps consisting of transfer-reaction operations is executed. As the “chemical process” is operated, the state of some of the molecules may change. The molecules in the artificial chemical plant represent an encoding of the actual decision variables, which can be of any type (real, integer, logical, combinatorial, etc.). For example, if a given decision variable is a real variable, it can be encoded into a binary string. For each bit in the binary string, a twostate molecule is assigned. The state of the molecule sets the bit value (one state is associated with one, and the other is associated with zero). When the molecule in the artificial plant changes its state, the corresponding bit changes its value accordingly. The same procedure is used for all other decision variables that are real. As another example, for a logical variable, a single molecule is needed; the molecule can have only two possible states, one that will correspond to the value TRUE and the other that will correspond to the value FALSE. For other type of variables, molecules with multiple states (>2) may be needed (i.e., combinatorial variables). This artificial chemical plant is a sequential algorithm with one function evaluation at a time (process step or iteration). The type of transfer-reaction operation to be executed next depends on the previous value of the objective function. See the schematic representation of the artificial chemical process shown in Figure 1. At some iteration during the algorithm, a set of molecules can be transferred from the load tank to the activation reactor in which the molecules are in a new state induced by a “chemical reaction”. In another type of transferreaction operation, some of the reacted molecules can be sent to the extraction unit where the molecules are deactivated back into the previous state entering the reactor. These extracted molecules could be sent to the separation unit or recycled back into the activation reactor where the molecules are reactivated to a new state. After a “good batch” is accomplished (better objective function is found), the activation reactor can be emptied into the separator unit; then a new “batch” is started. While this “plant” is in operation, the actual optimization problem is converging to a near-global optimal value. The actual algorithm is presented in Appendix A and described in the Flowchart in Figure 2. Although this is a new algorithm, it has proven to be very effective in terms of robustness and speed in a wide range of problems (Irizarry7-9). It has the capability of escaping from local minima and converges to a near-global optimal solution. The encoding capabilities allow for the solution of a large class of problems in which the decision variables can be virtually of any type (real, integer, nonnumerical, etc.). 3. Dynamic Optimization Using LARES-PR The dynamic optimization problem, DOP, can be formulated as

Find u(t) over t ∈ [to, tf] such that minF(x(tf))

(1a)

dx ) Ψ(x(t),u(t),t) dt

(1b)

u(t)

subject to:

8402

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006

Figure 2. LARES algorithm.

x(to) ) xo

(1c)

h(x(t),u(t)) ) 0

(1d)

c(x(t),u(t)) e 0

(1e)

uL e u(t) e uU

(1f)

where F is the performance index, u is the vector of control variables (control law), and x is the vector of state variables. The set of constraints consists of the dynamic model (eq 1b), the initial conditions of the state variables (eq 1c), the equality constraints (eq 1d), the inequality constraints (eq 1e), and the bounds on the control variables (eq 1f). A new algorithm for dynamic optimization was developed by Irizarry.7 The algorithm proved to be very powerful in solving dynamic optimization problems with different degrees of model complexity, size, and structure. The main steps of the algorithm are reviewed here, and then the modifications needed to consider hybrid dynamic optimization are presented in the next section. This procedure decodes the LARES decision variables (molecules) into a flexible representation of the control profile based on three key elements: (a) the use of finite elements, (b) variable-length elements/subdomains, and (c) the use of piece-

wise multifunctions. Using this representation of the control law, smooth regions, drastic changes in functionality, singularities, and discontinuities can be found as part of the solution of the optimization problem with a reduced number of decision variables. Section 3.1 describes the representation of the control law. A novel encoding of this representation of the control law is presented in Section 3.2. 3.1. Approximation of the Control Variable. In this algorithm, the control law is represented by continuous subsystems, each of them defined over variable-length nonoverlapping time intervals. For each subdomain, a function is either represented by a finite element trial function or by a userspecified function resulting in a hybrid formulation. 3.1.1. Variable-Length Finite Element Representation. Finite elements form part of the base functions that span the subspace searched by LARES. It is not used for discretization as in the simultaneous optimization and solution method. It is also a building block for the hybrid method discussed in the next section. In this formulation, the time-space is discretized with variablelength finite elements, Ei, described by the nodes τi and τi+1 where τi e τi+1. The element, Ei, is described with Me nodes. The first two nodes are (τi, τi+1) and are determined as part of


each element Ei, which identifies the functional form used in the element from a list of possible functions,

{

u ) f(t; Ri,Ri+1) if δi ) 1 u(t) ) u ) h(t; Ri,Ri+1) if δi ) 2 u ) g(t; Ri,Ri+1) if δi ) 3

Figure 3. Internal and external node numbering for the finite element representation.

Figure 4. Hybrid finite element/user-specified profile representation.

the solution of the DOP. The positions of any additional internal nodes are fixed with respect to (τi, τi+1). The value of the control law at the node points is υi, i ) 1, N. The internal element node values are υIi,k, where k is the counter of the internal element number. See Figure 3 for external and internal node numbering. The parametrization of the control using the finite elements trial functions and the min/max operators is defined as Me

u(t) ) max(umin, min{υiφ1 + υi+1φ2 +

I υi,k φk(l(t)), umax}) ∑ k)3

(2)

υIi,k

are unknown constant coefficients to be where υi and determined by the algorithm, l(t) is a simple linear transformation between the element time interval and the interval [-1, 1], and φi(ξ) is the finite element trial function for node i, defined in the interval ξ ∈ [-1,1]. If higher-order trial functions are used, each internal node will add more of an undetermined internal parameter to the set of unknown coefficients. 3.1.2. Variable-Length Combinatorial Representation. This formulation starts with a finite-element discretization using the highest order of elements to be considered in the analysis. Then a combinatorial variable is defined for each element, used as a switching mechanism. The combinatorial variable can replace the element with a lower-order element or with a new continuous subfunction defined over the same interval as that of the eliminated element. Figure 4 illustrates the hybrid formulation. In this example, a quadratic finite-element discretization is modified by the combinatorial variable into a hybrid formulation. The first two elements are replaced with lower-order elements (eliminating the middle node as a variable), the third element is replaced with a user-specified function, and element 3 is replaced with a second user-specified function. The last element remains as a quadratic element. It is also important when there are drastic changes in the functional form in different regions of time-space. This approach will allow a high-order approximation with a reduced number of undetermined coefficients. Notice that the word hybrid formulation is used in a different context: a control function composed of different type of trial functions. The hybrid dynamic optimization problem consists of discrete events embedded in the continuous dynamics of the problem. The multiple function capability is included by adding a combinatorial variable, δi (used as the switching variable), for

(3)

where τi e t e τi+1; f, h, and g are possible functional forms; and Ri is a vector of unknown parameters associated with element Ei. Ri ) {υi, υIi,k, other adjustable parameters for userdefined function}. For example, if there is a bang-bang optimal control, the constant functions h ) umin and g ) umax can be added to the set. In this case, when δi ) 2, if t is in element Ei, the value of u is constant, u(t) ) umin. Figure 4 shows schematically a control profile built with the hybrid approach. Notice that the same parameters can be used as adjustable parameters for different functions (multipurpose). 3.2. Coding the Control Function. Given the representation of the control law in its most general form in eq 3, the next step is to encode the unknown coefficients and combinatorial variables into molecule variables. The strategy for each type of variable is discussed next. 3.2.1. Unknown Node Values and Combinatorial Parameters. To encode each real variable, υi, i ) 1, N, K molecules (ηij, j ) 1, K) are assigned with two possible states, ηij ∈ {0,1}, for the binary encoding of υi in the range of [umin, umax]. The combinatorial variables for each finite element, δi, i ) 1, N 1, is an integer number whose range equals the number of possible functions to be used. There is one molecule ωi assigned for the combinatorial variable δi. The decoding of this variable is trivial: δi ) ωi (the number of states equals the number of possible functions with each state associated with a functional form). 3.2.2. Node Location Parameters. The most difficult variable to encode into LARES are the node location variables because of the constraint τi e τi+1 for i ) 2, N - 1. If standard binary encoding of these parameters in the range τi ∈ [0, tf] is utilized, this constraint will be violated frequently during LARES iterations. To solve this problem, two-step encoding was introduced by Irizarry.7 First, instead of encoding the nodelocation variables τi, computational variables si, i ) 2, N - 1, are used as decision variables defined in a fixed interval [0, 1]. These variables are then represented by R binary molecules to encode the variable in the range [0, 1], (ψij, j ) 1, R). In the second step, a transformation is defined that maps the computational variables into the node-location variables:

τi ) Ti(s2, s2, ..., sN-1)

(4)

This transformation is defined in such a way that, for any value of the vector s, the node order constraint is satisfied automatically. 3.2.2.1. Moving Partition (MP) Transformation. This transformation involves dynamic partitioning of time-space into disjoint spaces with a small gap between the spaces. Each element of the partition is assigned to each variable τi; in this way, each search vector will satisfy the constraint automatically. The boundaries used by Banga et al.10 are used to generate the partition

τbi

b b b τL,b i ≡ τi - β(τi - τi-1)/2

(5)

b ≡ τbi + β(τi+1 - τbi )/2 τU,b i

(6)

where is the best value found so far in the current iteration and the parameter β is used to control the gap between the intervals. A value of 0.65 is used in this work. This partitioning

8404


was used in the study by Banga et al.,10 in which the controlled random search method was customized to restrict sampling of the time variables. Here, this partitioning is used for the same purpose but in a different way. It is used to define a transformation of the computational variables by a mapping called moving U,b partition (MP), Ti: [0, 1] f [τL,b i ,τi ] defined as follows: U,b τi ) τL,b - τL,b i + si(τi i )

(7)

To use this transformation, the initial point needs to be feasible. 3.2.3. Inner Node Variables. When the quadratic or higherorder trial functions are part of the set of possible functions to approximate u, the variables υIi,k, k ) 3, Me, i ) 1, N - 1, also form part of the decision variables to be determined. In this work, special care is also taken for this variable in order to control the level of curvature allowed for the nonlinear trial function within the element. Two-step encoding is also used for this variable as follows: A computational variable is added to the set of decision variables, zi,k, k ) 3, Me, i ) 1, N - 1, defined in the interval [zL, zU], where zL and zU control the degree of curvature allowed. In this work, these parameters were set to 0 and 1, respectively. These variables are represented in the LARES artificial chemical plant with molecules (ζikj, j ) 1, K) with two possible states, ζikj ∈ {0,1}. The variables υIi,k are generated from the computational variables using the following transformation: I ) υi + zi,k(υi+1 - υi) υi,k

(8)

4. Hybrid Dynamic Optimization Problem In the description of the hybrid dynamic optimization problem, the nomenclature introduced by Barton and Lee4 is utilized. Let ti (i ) 1, ..., Ne + 1) be an increasing finite series of times where the system makes an instantaneous transition (changes from one continuous system to a new one like a mechanism change, for example). The transition can be model switches or state jumps. In the interval Ii ) [ti, ti+1], called the epoch, the system evolves continuously. It is assumed that there are Ne epochs. The continuous subsystems describing the dynamics in each epoch are called modes. The hybrid system is determined by the sequence of epochs, Tτ ){Ii, i ) 1, Ne} and the sequence of modes Tµ ) {mode for epoch 1, mode for epoch 2, ..., mode for epoch Ne}. The instantaneous transition from one mode to another mode is called an eVent, and the time, ti, when the instantaneous transition occurs is called the transition time. We consider three possible variations of the hybrid sequence and identify them as problem types for easy reference in the subsequent discussion. (i) Type 1: The optimal mode sequence, Tµ, is determined as part of the solution of the hybrid dynamic optimization problem while keeping the transition times fixed. (ii) Type 2: The algorithm determines the optimal transition times, Tτ, while keeping the mode sequence fixed. (iii) Type 3: The most general problem is solved where the transitions times and the sequence of modes (Tµ, Tτ) are part of the decision variables. (iv) Types 4 and 5: In these cases, the mode sequence and the transition times are fixed (type 4) or triggered by time or state variables (type 5). In both cases, discontinuities are embedded in the continuous dynamics but are not independently manipulated by the optimization algorithm. When the transition from current mode to a new mode is triggered by the state variables (types 4 and 5), the LARES-

Figure 5. Generalized profile/epoch representation.

PR algorithm can be utilized without any modifications, because the model is a black box for the algorithm. On the other hand, when the mode sequence or the transition times are part of the control strategy (Types 1-3), the LARES-PR methodology can be easily extended to consider these types of problems as described in the next section. 5. Extension of LARES-PR to Include Hybrid Dynamics The strategy is to have two independent discretizations for the time domain. One consists of the hybrid finite element profile of the original LARES-PR methodology, and the second discretization consists of discrete event times and mode specification for each epoch as shown in Figure 5. Notice that, for a problem with Mc control laws to be determined, the hybrid system becomes Mc + 1 independent discretizations. The first Mc corresponds to the hybrid finite element profiles, and the last one corresponds to the mode sequence/transition time discretization. Also notice that the word hybrid is used in two different contexts in this work: hybrid finite element profile is the generalized profile representation for a control law introduced by Irizarry,7 and hybrid dynamic optimization refers to the type of optimization problem to be solved (types 1-5). To encode the second profile, a similar methodology used for the control profile is utilized. Let us assume the most general case where Tτ and Tµ are part of the optimization problem. The mode for each epoch, mi, i ) 1, Ne, is an integer number (pointer) whose range equals the number of possible continuous subsystems. There is one molecule σi assigned for each epoch whose value represents the mode (continuous subsystem) for the epoch. The decoding of this variable is trivial: mi ) σi (the number of states equals the number of possible subsystems with each state associated with a functional form). For Tτ, the same two-step encoding in Section 3.2.2 is used. First, computational variables are defined for each discrete event time ri, i ) 2, Ne are used as decision variables defined in a fixed interval [0, 1]. These variables are then represented by R binary molecules to encode the variable in the range [0, 1], (θij, j ) 1, R). In the second step, a transformation is defined that maps the computational variables into the node-location variables: ti ) tL,b + i L,b ri(tU,b t ), where the upper and lower moving bounds are i i calculated using similar formulas as eqs 5 and 6. 5.1. Expanded LARES-PR Procedure. The extended LARES-PR algorithm is shown schematically in Figure 6 and described in detail in the following steps: 1. For the current LARES iteration, a new trial molecular state is generated in the artificial chemical process, ξht ) (ηtij, ψtij, θtij, ζtikj, ωti, σti). 2. Decode to generate the real parameters:

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8405 NC

Fp(x,u) ) F(x,u) +

pj max(0,cj(x,u)) ∑ j)1

(9)

where pj is the penalty parameter (a large positive number). The modified LARES-PR algorithm is solved using the expanded objective function in eq 9 instead of the problem objective function, F(x,u). As the algorithm searches for a nearglobal minimum, the terms in the summation will tend to zero, thus satisfying the constraints. 6. Benchmark Solutions

Figure 6. Extended LARES-PR algorithm.

Figure 7. Discrete event integration. t (a) Node values: υti ) decode(ηi1 , ..., ηtiK), i ) 1, N. (b) Computational variables used for node location: sti ) t , ..., ψtiR), i ) 2, N - 1. decode(ψi1 (c) Computational variables used for node location: rti ) t , ..., θtiR), i ) 2, Ne. decode(θi1 (d) Computational variables used for internal element interpolation: ztik ) decode(ζtik1, ..., ζtikP), i ) 1, N - 1. 3. Apply MP to generate the node location, τti, i ) 2, N - 1, and the transition times, tti, i ) 2, Ne. 4. Apply the transformation for the internal-node variables (if needed): υIik, i ) 1, N - 1, k ) 3, Me. 5. Set the combinatorial parameters: δti ) ωti, i ) 1, N - 1. 6. Set the mode for each epoch: mti ) σti, i ) 1, Ne. 7. Build the trial control law using the parameters in steps 2-5: ut(t) ) u(t; υt,τt,δt,υI,t). 8. Integrate the dynamic model using the trial control, ut(t), transition times, tti, i ) 2, Ne, and the mode for each epoch, mti, i ) 1, Ne. 9. Evaluate the performance index. Return the performance index to LARES. If the LARES stopping criterion is not satisfied, LARES will generate a new trial vector; return to step 1. The other modification of the algorithm is that the integration of the dynamic model is performed in event-driven intervals, as shown in Figure 7. When the terminal time is part of the solution of the OCP, the position of the node N becomes part of the decision variables. For this type of problem, a large value of time T is selected, for the MP, where τbN e T. When inequality constraints are added to the optimization problem (i.e., cj(x,u) < 0, j ) 1, NC), they are implemented using the penalty method. The approach is to simply substitute the objective function of the original problem with the following expanded objective function,

Reported solutions of hybrid dynamic optimization problems are limited. In most cases, the systems consist of small-size problems with linear dynamical systems. In this section, we consider the problems solved in the works of Barton and Lee4 and Lee and Barton.5 Here, results are presented for three example problems ranging from very small to average size. This is an interesting set of problems because the problem size can be parametrized to study the performance of the algorithm vs problem size. Furthermore, the global optimal solution is available. The analysis of examples 2 and 3 are further extended to consider a type-3 problem. All calculations were performed on an Intel pentium M 1.6 Ghz machine with 504 MB RAM. The termination criterion is the user-specified maximum number of iterations. The time to reach the global optimum is studied to understand the convergence properties of this algorithm and to suggest the maximum number of iterations for problems of similar size and structure. 6.1. Example 1. This example was solved by Barton and Lee4 with their deterministic method for the type-5 problem (fixed-time and fixed-mode sequence). This example consists of the following hybrid-dynamic optimization:

min F(p) ) p

∫03 - x(p,t)2 dt

(10a)

where p ∈ [-4,4] and x(p,t) is the solution of the following hybrid system

Mode 1:

dx ) -2tx + p dt

(10b)

dx x + p ) dt t + 10

(10c)

Mode 2:

The problem was solved by a determistic method with x(p,0) ) 1 and fixed Tµ ) {1, 2, 1}. The fixed transition times are t1 ) 1 and t2 ) 2. The global solution found was p ) 4 with F ) -15.360. The same problem was solved using LARES-PR, finding the same solution. The simulation was performed using four constant elements for the control law. Then the problem was relaxed to a type-2 problem, and the constant control function was relaxed to a piecewise linear finite element with 13 linear elements to represent the control law. The global solution was improved to F ) -29.9 with new transition times t1 ) 0.64 and t2 ) 2.997, and the same constant profile (p ) 4) was found. The solution was 0.2% from the final solution in 6 010 iterations, and the best value was reached in 56 000 iterations. This particular model was solved using the RungeKutta method with constant time steps of 0.001 and 0.000 05 to confirm a discretization-independent solution. 6.2. Example 2. This example was introduced by Lee and Barton5 to test their deterministic algorithm developed for linear hybrid systems of type 1 or 2. The problem consists of a plugflow reactor operating at steady state. The reactor is loaded with

8406


Table 1. Solution of Catalyst Problem Using LARES-PR (Hybrid Problem Type 1) Ne

F* (exact)

F* (LARES-PR)

N* × 1000

5 8 12 17 20 22

296.6 295.0 300.5 306.4 308.6 309.5

296.6 295.0 300.5 306.4 308.6 309.5

0.03 0.16 1.0 2.0 1.8 2.9

three possible choices of solid catalyst that can be loaded into sections of the reactor. Initially, the reactor is loaded with pure A. The reaction scheme is

A9 8 I; A 9 8 W1; I 9 8P; I 9 8W2 k k k k 1

2

3

(11)

4

The mass balances are

dxA ) -(k1 + k2)xA dl

(12)

dxW1 ) k2xA dl

(13)

dxI ) k1xA - (k3 + k4)xI dl

(14)

dxW2

) k4xI

(15)

dxP ) k3xI dl

(16)

dl

with initial conditions xA ) 1000; xW1 ) xW2 ) xI ) xP ) 0; and reactor length ) 1. In these expressions, x is the molar concentration, l is the position in the reactor, and kj is the rate constant of reaction j. The rate constant in each section of the reactor depends on the type of catalyst:

Catalyst 1: k1 ) 2.098, k2 ) 1.317, k3 ) 0.021, k4 ) 0.033 Catalyst 2: k1 ) 29.53, k2 ) 110.2, k3 ) 0.295, k4 ) 0.079 Catalyst 3: k1 ) 182.6, k2 ) 2325, k3 ) 1.826, k4 ) 0.143 The objective is to find the optimal distribution of catalyst to maximize profit, which is the value of the product minus the waste treatment of byproducts:

min - F ) 0.01xW1(1) + 0.1xW2(1) - xP(1) Tµ

(17)

The reactor is divided into Ne equal segments (epochs). The type of catalyst in each segment is determined as part of the optimization problem. To compare the results with the analytical ones, the same problem is solved using the extended LARESPR using a fixed grid and only the mode changes (type 1). This is coded very easily by assigning one molecule to each epoch with three states (one for each possible catalyst). For this problem, the maximum number of function evaluations was set to 5 000. As shown in Table 1, in all cases the exact solution was found. The best-so-far performance as a function of function evaluations is shown in Figure 8. The number of function

Figure 8. Performance for example 2 for Ne ) (a) 5, (b) 12, (c) 17, and (d) 22. Table 2. Solution of Catalyst Problem Using LARES-PR (Hybrid Problem Type 3) Ne

F* (LARES-PR)

N* × 1000

3 4 5 8 9 10

314.2 296.9 314.2 314.2 314.2 314.2

8.1 0.16 11.6 1.0 4.1 2.7

evaluations required to reach the global optimum for the largest problem was 3 000 iterations. For a midsize problem with Ne ) 8, the number of iterations was as low as 160 to find the global optimum. Notice that in most conditions the solution is very close to the global optimum after 1 000 function evaluations. The dynamic models were integrated using LSODA using the discrete event integration scheme in Figure 7. These results show the robustness of the proposed algorithm in finding a near-global optimal solution for the different types of hybrid optimization problems. Although the convergence to a near-global optimum for a stochastic optimization algorithm is guaranteed only in a probabilistic sense, in this particular problem, the exact solution was found in all numerical experiments. For cases in which the size of the catalyst is limited to a fixed length, the solutions in Table 1 are the global optimum. If in practice there are no constraints in the length of catalyst that can be used, then solving the problem as a type-1 problem could result in a suboptimal design because of the artificial constraint imposed by fixing the transition times, as shown in the next set of results. The same problem is relaxed to consider the possibility of a variable-length section (type 3). This type problem was not considered by Lee and Barton,5 because their algorithm considers a type 1 or 2 type problem. As shown in Table 2, LARESPR was able to find a global optimum of 314.2, a better value than all the results in Table 1 for a fixed transition times. These results demonstrate the importance of having the flexibility of solving the hybrid problem in its most general form to find a true optimal design. Figure 9 shows the distribution of catalyst as a function of different design types. Design #1 correspond to the true global optimum design found considering a type-3 problem with Ne ) 5. Designs 2-5 correspond to solutions of type-2 problems with Ne ) 5, 12, 20, and 22, respectively. As Ne is increased, it approaches the true optimal design at the expense of increasing the problem size. Even at Ne ) 22, a suboptimal design is obtained (309.5) when compared with the optimal design of a type-3 problem for Ne ) 5 (314.2). The performance of the relaxed problem is shown in Figure 10. It is interesting to note that, for a medium problem size, the solution converged faster to a near-optimum value in ∼1 0004 000 iterations, while for smaller problems, Ne ) 5 it required


Figure 11. Computation time for example 2. Figure 9. Optimal catalyst distribution for different designs. Design #: (1) distribution obtained by solving type 3 problem with 5 epochs (best design), (2) type 1 problem with five epochs, (3) type 1 problem with 12 epochs, (4) type 1 problem with 20 epochs, and (5) type 1 problem with 22 epoch.

{

N˙ A ) 4NB - 2NA N˙ B ) 2NA - 5NB N˙ P ) NB

Mode 3:

(18d)

Constraints: Ne

δ2σ FB e 2Ne ∑ i)1

(18e)

NP(1) g 0.6

(18f)

log(NA(1) + NB(1) + 1) e 1.2

(18g)

1 e FA, FB e 5

(18h)

i

Figure 10. Performance for the catalyst problem (type-3 problem).

more iterations (∼10 000). For an even smaller value Ne ) 4, a suboptimal value was found; it required 170 000 iterations to find the global optimum. This indicates that, if a too-small number of elements are used, the search is made more difficult by constraining the trial vectors during the search. Numerical experience with LARES-PR and the extended LARES-PR shows that the best results are found when the elements are of the order 5-15. In practice, it is a good idea to solve the same problem using different levels of discretization. Figure 11 shows an estimate of CPU time to reach the optimum solution as a function of the number of epochs. The increase in computational cost in going from a type-1 to a type-3 problem is not significant. In 4 s, the global optimum was found for a type-3 problem in the same order as when the problem is solved with 22 epochs. 6.3. Example 3. The second example used by Lee and Barton5 consists of a well-mixed reactor. The physical description is given by Barton and not repeated here. The hybrid optimization problem consists of

max

Tµ,FA(t),FB(t)

S)

Np(1) NA(1) + NB(1) + 1

(18a)

ST:

Mode 1:

Mode 2:

{ {

Models:

N˙ A ) 4NB - 2NA + FA N˙ B ) 2NA - 5NB N˙ P ) NB

(18b)

N˙ A ) 4NB - 2NA N˙ B ) 2NA - 5NB + FB N˙ P ) NB

(18c)

Initial conditions: NA(0) ) NB(0) ) NC(0) ) 0

(18i)

where Ni is the concentration of reactant i, FA is the flow rate of A, FB is the flow rate of B, and δ2σi has a value of 1 for epochs with a “molecule” state of 2 (σi ) 2) and 0 otherwise. Similar to example 2, the type-1 problem considered by Lee and Barton5 is considered first; then the same problem is relaxed to consider the generalized problem to show the flexibility and robustness of LARES-PR. As shown in Table 3, in all cases the exact solution was found. The relaxed problem (type 3) is shown in Table 4. In this case, by considering the relaxed problem in addition to having a better cost function, the design was simplified (an added bonus); the optimal profile consisted simply of adding FB ) 5 in a time period t ) [0, 0.59]. For the fixed case, all optimal solutions have a lower cost function and also require the use of FA (an added design cost and process complexity). Figure 12 shows the performance of the type-3 problem vs the type-1 problem. The relaxed problem converged faster to a better value than the problem with fixed transition times. Similar to a previous case, using a too-small number of epochs (Ne < 4) is not recommended for the type-3 problem. Using a larger number of epochs, the natural epoch sequence will emerge faster. Figure 13 shows the CPU time needed to reach the optimum value. It is interesting to see that, when solving a type-3 problem, in addition to finding a better solution, the algorithm approaches the global optimum in 1/10 the number of iterations required for a type-1 problem. In general, these examples demonstrated that the extended LARES-PR is a robust method to find the near-global optimum of different types of hybrid problems. In particular, it is very powerful in considering the most difficult one (type 3) with no increase in computational cost, since the maximum number of iterations can safely be set the same used for type-1 problems, as suggested by the results in Figure 13.

8408


Table 3. Solution of Isothermal Reactor Problem Using LARES-PR (Hybrid Problem Type 1) Ne

S* (exact)

S* (LARES-PR)

N* × 1000

5 6 7 8 9 10 11 12

0.238 0.236 0.240 0.240 0.241 0.243 0.242 0.243

0.238 0.236 0.240 0.240 0.241 0.243 0.240 0.243

27.5 66.7 44.1 115.2 172.0 151.2 37.5 45.1

Table 4. Solution of Isothermal Problem Using LARES-PR (Hybrid Problem Type 3) Ne

S* (LARES-PR)

N* × 1000

5 7 10 12

0.259 0.259 0.259 0.259

1.0 2.1 6.3 4.8

7. Application to Population Balance Models: Determination of Optimal Nucleation Profile and Mechanism Switches In this section, the hybrid-dynamic optimization methodology is applied to the synthesis of binary nanoparticles. The synthesis is made in two phases. Initially, the reactor is charged with nuclei with a volume Vo of particles of type A. More nuclei of A particles can be added to the system in a fed-batch mode. Then, at a time Tswitch, the B type particles are added to the reactor in a very fast mode (instantaneous addition). After Tswitch, the binary system is allowed to grow by coalescence. At the end of the process, particles with total volume V are composed of a partial volume of type A, VA, and a partial volume of type B, VB, where V ) VA + VB. It is further assumed that the average particle size can be monitored and the residence time, RT, is extended until the desired average particle size is achieved. A schematic of the process is shown in Figure 14. A multiobjective problem is stated as a single-objective optimization problem. First, to grow the particles to a desired

size, the system is allowed to evolve until the desired size is reached. Thus, this objective is included as an equality constraint impacting the process RT. The main objective is to grow particles with a desired volume composition of B/A and a desired average particle size with a minimal RT. The hybrid optimal control problem is formulated as follows: where the

min w1|E[VB/V] - Xset| + w2RT + p

J(t),Tswitch

(19a)

ST PSD of the binary particles is given by the following hybrid system:

Mode 1: ∂f(V,t) ) J(t)δ(V - Vo) + ∂t 1 V q(V - V′,V′)f(V - V′,t)f(V′,t) dV′ 2 0

∫

∫V∞q(V,V′)f(V,t)f(V′,t) dV′

(19b)

(This model is solved by the discretization method.) Mode 2: -

∂n(V,VA,t) V min (V,VA) ) 0 0 β(V - ,,VA - γ,γ) × ∂t n(V - ,VA - γ,t)n(,γ,t) dγ d -

∫∫

∫0∞ ∫0β(V,,VA,γ)n(V,VA,t)n(,γ,t) dγ d

(19c)

(Instead of solving this model, the system dynamics is simulated using the Monte Carlo method.) Equality constraint: ravg(RT) ) rset

(19d)

Combinatorial constraints: J(t) )

{

Jset t ∈ Ii 0 otherwise

Count(I) e Nmax

(19e) (19f)

Inequality constraint: 0 e Tswitch e Ts_max

Figure 12. Performance for example 3 (type 1 and type 3).

Figure 13. Computation time for example 3.

(19g)

where, in the objective function, E[x] is the expected value of x, Xset is the target value for the average volume fraction of B, w1 and w2 are weighing factors, and p is a penalty parameter. The constraint in eq 19b is the PBE describing the dynamic of the coalescing particle population of A particles when the system is in phase 1, where f is the particle density function, J is the nuclei addition rate of A particles, and q is the aggregation rate. V and V′ are particle volumes, and t is the time. The first term on the right-hand side is the nucleation term and is used as the control parameter. The profiles of the nuclei addition and switching time are optimally determined by solving eq 19. The constraint in eq 19c is a 2D population balance equation describing the dynamic of the coalescing binary particle system when the system is in phase 2, n(V, VA) is the particle density of the A-B binary particle system, and β is the 2D aggregation rate. The equality constraint in eq 19d determines the duration of the synthesis, where ravg(RT) is the average particle radius at time RT and rset is the target average particle radius.


Figure 14. Binary particulate process.

The constraints in eqs 19e and 19f are operational constraints added to the system, where Ii ) [ti0, tif] is the ith time interval for adding nuclei to the system at a fixed rate Jset, and Count(I) is the number of addition intervals limited by Nmax. The constraint in eq 19g limits the optimal switching time to the interval [0, Ts_max]. 7.1. Coalescence Mechanisms. When the system is in phase 1, the Brownian collision kernels with a DVLO correction is utilized,

qB(V,V′) )

8kT 1/3 (V + V′1/3)(V-1/3 + V′-1/3) 3µ

(20)

where k is Boltzmann’s constant, T is the reactor temperature, and µ is the viscosity. The Brownian collision kernel is corrected with the so-called stability factor, using standard DVLO theory to account for the attractive and repulsive forces between particles (Russel et al.11 and Schwarzer and Peukert12),

qB

q(V,V′) )

∫0

2

∞

exp(VT(x,V,V′)/kT) (x + 2)2

(21) dx

where x is the shortest distance between particle surfaces scaled by the average radius between the two particles. The interaction energy VT is the sum of all forces between particles. In particular, the effects of van der Waals attractive forces and electrostatic repulsive forces are considered (VT ) Vvdw + Velec). The van der Waals potential is given by

Vvdw ) -

[

AH 2a1a2 2a1a2 + + 6 R2 - (a - a )2 R2 - (a + a )2 1 2 1 2 ln

]

R2 - (a1 + a2)2 R2 - (a1 - a2)2

(22)

where ai is the radius of particle i, i ) 1, 2, and R is the distance between two particle centers; all other parameters and values are defined in Table 5. There are several formulations for the electrostatic interaction. To simplify the problem, the following expression is used,

Velec ) 4πψ02

a1a2 ln[1 + exp(-kx)] a1 + a2

(23)

7.2. Solution of the PBE Equations when the System is in Mode 1. The PBE is written in dimensionless forms using the following scales: the time is scaled with 1/(qofo), the particle density function is scaled with fo, the nuclei addition rate is scaled with qo fo fo, and the particle volume is scaled with Vo, where qo ) 8kT/3µ, fo is the initial particle density in the reactor,

Table 5. DVLO and Process Parameters variable

description

value

AH T m e k yo ao aAo aBo fo

Hamaker constant temperature viscosity dielectric constant of medium inverse Debye length surface potential nuclei particle radius initial A particle radius in phase I initial B particle radius in phase II initial particle concentration in phase I

1 × 10-19 J 300 K 0.001 Pa‚s 2.17 × 10-10 C/(V m) 2 × 108 (1/m) 50 mV 2 nm 2 nm 20 nm 1017 par/m3

and Vo is the nuclei volume. The PBE is discretized by the method of Kumar and Ramkrishna13 with 50 volume nodes (V1 ) Vo; Vi ) 1.3Vi-1, i ) 2, 10; Vi ) 1.5Vi-1, i ) 11, 50). 7.3. Solution of the PBE Equations when the System is in Mode 2. Although mode 2 could be solved directly using finite elements methods, the dynamic of this coalescing of a binary system is modeled by the Monte Carlo (MC) method using an event-driven simulation (Garcia et al.14 and Kruis et al.15). First, a simulation box Vs is selected. The number of particles in the simulation box is then determined from the particle concentration. The MC simulation consists of the following steps, executed iteratively until a final time is reached. In the first step, the time for the next collision, tcol, is determined (Garcia et al.14). The process time is moved to t ) t + tcol. The two colliding particles (i, j) are then chosen using the method described by Kruis et al.15 to coalesce and create the new particle k, with volume Vk ) Vi + Vj. The composition of the new particle is given from the composition of the mother particles (VAk ) VAi + VAj and VBk ) VBi + VBj). In this case, to simplify the MC simulation, a simple constant collision kernel is considered (β ) qo). All variables are made dimensionless using the same scales used in phase 1. 7.4. Switch from Mode 1 to Mode 2. At time Tswitch, the B particles are added rapidly (instant addition) to the reactor, and from that point on, the binary systems start to evolve. The transition from the numerical solution of the discretized population-balance equation to the Monte Carlo simulation of the binary systems consists of defining the initial ensemble of particles and volume of the simulation box. First, the number of simulation A-particles is set, NA. The initial simulation volume is determined from Vs ) NA/µo, where µo is the total particle concentration of A particles at the end of phase I (zero moment of f(VA, Tswitch)). The initial size for each of the NA A-particles in the simulation box is determined by sampling from the particle-size distribution obtained in phase 1, f(VA, Tswitch). It is assumed that B-particles are concentrated and that the reactor volume change after the addition of B-particles is negligible. The addition is made to have a one-to-one ratio of A and B particles in the reactor. With these assumptions, NB )

8410


NA particles with volume VBo are added to the simulation box, keeping the simulation box volume constant. The constraint in eqs19e and 19f forces the shape of the optimal profile to a set of constant addition intervals. The number of intervals and the duration of those intervals are determined as part of the solution of the optimization problem. The maximum number of possible addition intervals is Nmax. These constraints convert this optimization problem into a complex combinatorial one. LARES-PR can implement this constraint automatically by just using the combinatorial profile representation with two possible functions for δi ) {1, 2} ) {u(t) ) 0, u(t) ) Jset}. This is a very complex hybrid optimization algorithm of type 2 that cannot be solved by deterministic methods because of the following complications: (1) The dynamics in mode 2 is solved by Monte Carlo (equationless model). (2) The transition from phase 1 to phase 2 requires generating the initial population of particles using stochastic sampling (again, an equationless process). (3) The time horizon is not known a priori; it is determined by constraint eq 19d. 7.5. Results. Consider the following design with process conditions summarized in Table 5. In Phase 1, the initial conditions consist of a monodisperse system with A-particles with a unit dimensionless volume. Nuclei particles can be added while the process is in phase 1. The maximum dimensionless time allowed for phase one is Ts_max ) 6000. The dimensionless addition rate is constant Jset ) 0.005, but the addition intervals are decision variables. It is allowed to have a maximum of three addition intervals (Nmax ) 3). The Brownian aggregation kernel is utilized including DVLO effects with parameters shown in Table 5. The initial population of B-particles is monodispersed with a dimensionless volume of 1 000. The solution determines the optimal intervals for nuclei addition, the switching time from phase 1 to phase 2, and the residence time. All results are also presented in terms of volume fraction density function (Φ(V,t) ) V f(V,t)) vs particle radius. The target is to grow an ensemble average dimensionless radius of rset ) 60 (120 nm) binary particles at the end of phase 2 with an ensemble volume composition of VB/V ) 0.4(Xset ) 0.4). This problem was solved using the extended LARESPR set for type 2 with two epochs. The control profile consisted of three combinatorial segments. The termination criterion was set to 10 000 function evaluations. In the objective function, w1 ) 1, w2 ) 1/20 000, and p ) 0. Also, the penalty strategy was described by Irizarry,7 for combinatorial optimization. To avoid redundant search (element size adjustment of a constant profile), when all intervals have the same value (δi ) 1 ∀ i or δi ) 0 ∀ i), the penalty parameter is set to p ) 1000; otherwise, the parameter is 0 (after confirming that the two extreme solutions are suboptimal). The optimal solution has one addition interval (Count(I)* ) / 1, I/1[0, 4 949], Tswitch ) 5 977, the optimal residence time RT* ) 17 817. The solution to the optimization problem is the following process synthesis: (1) Add nuclei at constant rate of Jset ) 0.005 of a period of time [0, 4 949]. Stop the addition and let phase 1 continue until the switching time is reached, Tswitch ) 5 977. The B-particles are added very fast to start Phase 2, which continues for a period [5 977, 17 817]. This schedule produces the optimal conditions of an ensemble average radius of 120 nm and volume fraction of B-particles of 0.4, the exact specifications. The optimal residence time was 17 817, while some suboptimal solutions consisted of residence times as large as 150 000. Figure 15 shows the volume density function of A-particles after phase 1 and the binary particle after phase 2.

Figure 15. Volume fraction distribution vs dimensionless particle radius: (a) optimal profile of phase 1 and (b) optimal profile of constrained phase 2.

The objective function can be modified as a function of the specifications needed, for example, a low size variability, etc. The results presented in this section show the power of the current methodology to solve very complex hybrid dynamic optimization problems. Because designs are highly dependent on product specifications, this can be a valuable tool to study optimal conditions for a multipurpose facility. 8. Conclusions The LARES-PR algorithm was extended and applied to the solution of hybrid-dynamic optimization problems in its most general form. For problems of types 4 and 5 (as defined in this paper), the original LARES-PR can be utilized without any modifications. The model is a black box, and the discontinuities are captured during integration of the model. For problems of types 1-3, the modifications introduced in this paper allow results in a robust methodology to solve these types of problems. This can be used to solve hybrid-dynamic optimization problems in their most general form. This feature eliminates artificial constraints added when the solution is restricted to certain types of functionalities of the control law or a specified transition time for the switches. The structure of this algorithm can be exploited to find near-global optimal solutions with simplified discrete event sequences and solution features. Numerical solutions of benchmark problems demonstrate the robustness and effectiveness of the extended LARES-PR to solve a wide range of problems under the same framework. The hybrid optimization problem of a production of an alloy particle system was considered for the first time. This was a complex problem because, in addition to a switching in dynamics, the solution methodology switches during the integration. Notice that this type of problem cannot be solved by the deterministic method given the equation-free nature of the Monte Carlo simulation. These results demonstrate the usefulness of LARES-PR in solving optimal control problems resulting in the production of nanoparticles. As the complexity of the system increases, so does the solution methodology, ranging from a large set of ODEs to equation-free Monte Carlo type simulations. Since LARES-PR solves the model many times, an area of further research is to combine a detailed model with a reduced-order model to increase the speed for very large models. Appendix A. LARES Algorithm A flowchart of the algorithm is in Figure 1. Let L, AR, E, and S be four disjunctive sets whose elements are the molecule variables. Let F be the objective function to be minimized, expressed in terms of the vector of molecule variables, F ) F(x). Let the value of the best solution found so far be xg )


(xg1, ..., xgV). The algorithm consists of redistributing the molecule variables among the sets to generate new trial vectors using the following rule: whenever a molecule variable xj becomes a member of the set AR, a new value different from the best value found is assigned to it. Let xaj * xgj be the new value, which will not change until the variable is moved out of the set AR. For each iteration, the trial vector xt ) (xt1, ..., xtV) is then constructed using the following formula:

xtj )

{

xgj if xj ∉ AR xaj if xj ∈ AR

(24)

With these definitions, the algorithm is described next. Initialization: 1. The algorithm starts by initializing xg randomly and placing all variables in L (set AR ) E ) S ) Ø, L ) {x1, ..., xV}). Outer loop: Perturbation to form AR 2. Select a random number, |Trx| (|Trx| e |L|), from a uniform PDF,

|Trx| ) min({int}(F‚[V × co] + 1), |L|)

RR )

(25)

where F is uniformly distributed in (0, 1) and co is an adjustable parameter used to select the average fraction of elements to be selected from V, which is the total number of molecules representing the possible solutions to the problem. 3. Select |Trx| elements randomly from L to form the subset Trx, (Trx ⊂ L). 4. Transfer the subset to AR: L ) L\Trx; AR ) AR ∪ Trx. 5. Select a random new value for each molecule variable in Trx: xaj * xgj , ∀j ∈ Trx. 6. The new trial vector, xt, is generated using eq 24. 7. If F(xt) < F(xg), the trial vector is accepted as a new best solution found, xg ) xt. In this case, all of the elements in AR are sent to the set S: S ) S ∪ AR; AR ) Ø. If the algorithm termination criterion is achieved, exit the algorithm and return xg as the solution to the optimization problem. 8. If a better solution was found in step 7, skip the inner loop and perform another outer loop iteration: Go to step 17. Otherwise, continue with step 9. 9.Initialize the parameter RP: RP ) F(xt). This parameter is used and modified in the “goodness” test in the inner loop. Also set |AR|0 ) |AR|, which is the initial number of molecules in AR before starting the next inner loop. Inner loop: Iterative improvement of AR 10. Select a random number, |E|, from a prescribed PDF, |E| e |AR|. The formula used is

|E| ) min[int(F‚[|AR|0 × ci] + 1), |AR|]

in E will prefer to stay in their ground state (xj ) xgj ) to generate better solutions. In this case, the elements in E are transferred to S: S ) S ∪ E; E ) Ø and the metric RP is updated, RP ) F(xt). When F(xt) > RP, the hypothesis is that there is a high probability that most elements in E will induce a better solution if they are in a different state from their ground state. In this case, a new activated state is generated for all elements in E (xj ) xaj * xgj , ∀j ∈ E), and all of the elements in E are transferred back to AR (AR ) AR ∪ E; E ) Ø). Check conditions to exit or continue the inner loop: 16. If one of the following conditions is satisfied, the algorithm will exit the inner loop (go to step 17): • If a better solution was found in step 14. • The number of elements in AR is less than or equal to one, |AR| e 1. For some particular problems, this restriction could be generalized to a generic threshold parameter (|AR| < c), but this case is not considered in this work. • The large recycle-ratio (RR) is defined as

(26)

where |AR|0 is defined in step 9 and ci is another adjustable parameter. 11. Select |E| elements randomly from AR to form the subset E. 12. Extract the subset E from AR: AR ) AR\E. 13. The new trial vector, xt, is generated using eq 24. 14. If F(xt) < F(xg), the trial vector is accepted as a new best solution found, xg ) xt. All of the elements in AR are sent to the set S: S ) S ∪ AR; AR ) Ø. If the algorithm termination criterion is achieved, exit the algorithm and return xg as the solution to the optimization problem. 15. Improvement criterion for AR: If F(xt) e RP, the hypothesis is that there is a high probability that most elements

rec g RRT |AR|0

(27)

where rec is the counter of the number of times that E is sent back to AR in the current inner loop and |AR|0 is defined in step 9 as the initial number of molecules in AR generated in the outer loop. The parameter RRT is adjustable. Otherwise, a new inner loop iteration is started: go to step 10. Note that most of the inner loops should be terminated naturally either by a better solution or by |AR| e 1. The RRT is set large enough to let most of the inner loops be terminated naturally by a better solution or |AR| e 1. The RRT should not be too large, because then many iterations will be allowed for an undesired inner loop, affecting the algorithm’s performance. Check the number of elements in L and AR, and the algorithm termination criterion. 17. Check the number of elements in the sets L. If the number of elements in L is below a prescribed value, LT, all of the elements in the set S are transferred: L ) L ∪ S; S ) Ø. If the number of elements in L is still too low (|L| e LT), then all of the elements in AR are transferred to L: L ) L ∪ AR; AR ) Ø. 18. Start an outer loop iteration, returning to step 2. The following parameters are used: RRT ) 1.0, co ) 0.3, ci ) 0.25, LT ) V/2. Literature Cited (1) Ahtchi-Ali, B.; Pedersenb, H. Dynamic Simulation of Operational and Biological Switches in Biochemical Reactors. Process Biochem. 1992, 27, 3. (2) Floudas C. A. Deterministic global optimization: Theory, methods and applications. NonconVex optimization and its applications; Kluwer Academic Publishers: Norwell, MA, 2000. (3) Manon, P.; Claire, V. R.; Gilles, G. Optimal control of hybrid dynamical systems: application in process engineering. Control Eng. Pract. 2002, 10, 133. (4) Barton, P. I.; Lee, C. K. Design of process operations using hybrid dynamic optimization. Comput. Chem. Eng. 2004, 28, 955. (5) Lee, C. K.; Barton, P. I. Optimal mode sequences for linear hybrid systems. IEEE Trans. Autom. Control 2005, submitted for publication. (6) Barton, P. I.; Banga, J. R.; Gala’n, S. Optimization of hybrid discrete: Continuous dynamic systems. Comput. Chem. Eng. 2000, 24, 2171. (7) Irizarry, R. A generalized framework for solving dynamic optimization problems using the artificial chemical process paradigm: Applications to particulate processes and discrete dynamic systems. Chem. Eng. Sci. 2005, 60, 5663.

8412


(8) Irizarry, R. LARES: An Artificial Chemical Process Approach for Optimization. EVol. Comp. J. 2004, 12 (4), 435. (9) Irizarry, R. Fuzzy classification with an artificial chemical process. Chem. Eng. Sci. 2005, 60, 399. (10) Banga, J. R.; Irizarry-Rivera, R.; Seider, W. D. Stochastic optimization for optimal and model-predictive control. Comput. Chem. Eng. 1998, 22, 603. (11) Russel, W. B.; Saville, D. A.; Schowalter, W. R. Colloidal dispersions; Cambridge University Press: Cambridge, U.K., 1989. (12) Schwarzer, H. C.; Peukert, W. Prediction of aggregation kinetics based on surface properties of nanoparticles. Chem. Eng. Sci. 2005, 60, 11.

(13) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization: I. A fixed point technique. Chem. Eng. Sci. 1996, 51, 1311. (14) Garcia, A. L.; van den Broek, C.; Aertsens, M.; Serneels, R. A Monte Carlo method of coagulation. Physica 1987, 143A, 535. (15) Kruis, F. E.; Maisels, A.; Fissan, H. Direct simulation Monte Carlo method for particle coagulation and aggregation. AIChE J. 2000, 46, 1735.

ReceiVed for reView April 12, 2006 ReVised manuscript receiVed October 4, 2006 Accepted October 4, 2006 IE060463Z

Hybrid Dynamic Optimization Using Artificial Chemical Process

Recommend Documents