Global and Dynamic Optimization using the Artificial Chemical Process Paradigm and Fast Monte Carlo Methods for the Solution of Population Balance Models

Global and dynamic optimization of engineering problems usually involves complex physico-chemical models as constraints. These models are in general highly non-linear, resulting in multimodal optimization problems. The model may have discontinuous behavior and/or include a very large set of variables. As the complexity of the systems increases, equation-free modeling is becoming more common (Kevrekidis, Gear and Hummer, 2004). For example, in particle dynamics, population balance models are sometimes more effectively solved by the Monte Carlo method. Stochastic global optimization methods are very important algorithms for the solution of these types of problems. They have been successfully applied to solve challenging problems that cannot be solved using gradient based methods. Stochastic optimization methods have also been used in many algorithms, in which solution of optimization problems is part of the algorithm. Global stochastic optimization strategies have been utilized in learning phase of pattern recognition algorithms using fuzzy logic (Irizarry, 2005b) and neuro-fuzzy systems (Lin, 2008). These methods have been used for the optimization of complex engineering designs involving computational fluid mechanics such as aerodynamics applications (Duvigneau and Visonneau, 2004). Other applications include the determination of molecular structures, including protein structure prediction and protein-small molecule interactions among others (Sahinis, 2009). Batch scheduling problems are another type of problem were stochastic optimization can be very efficient (Liu et al., 2010). In particular, the solution of dynamic optimization problems is also of great industrial importance for process development and process optimization, since most processes are dynamic. In this type of problem an optimal profile function is sought (vs. an optimal value for a set of variables). For example, in a fed-batch fermenter, the feed-rate schedule is optimized to maximize production of antibiotics, vitamins, enzymes, and other products (Banga et al., 2003). Another example is the determination of optimal temperature profiles in crystallization processes to control crystal size distribution (Ma, Tafti and Braatz, 2002). Dynamic optimization is also of central importance to the application of process control

Global and Dynamic Optimization using the Artificial Chemical Process Paradigm and Fast Monte Carlo Methods for the Solution of Population Balance Models 389 called molecules, and their respective values are called states. As a motivating example consider a case were θ is a vector of real variables. The real decision variables can be encoded using a binary representation (similar to GA). Unlike GA, in this encoding each bit of the string is associated with a molecule variable with two possible values (0 and 1), which represents the value of a specific bit in the binary encoding (see Figure 1). In general, any type of decision variables (including real, integer, logical, and combinatorial) can be encoded into a set of molecules, making the algorithm very flexible. Given the decision vector represented as a set of molecules, the LARES algorithm operates on these molecules to create new trial states. At each iteration of the algorithm, a subset of molecules will change state (value) to generate a new trial vector. The artificial chemical plant concept is based on the fact that chemical reactors convert a low-quality material into a high-value product by a series of reactions, feedback loops, and separation steps. The following algorithm is based on an abstraction of those steps.

The LARES artificial chemical plant paradigm
The LARES algorithm is an iterative improvement methodology, which considers one solution at a time. Given the decision variables encoded into molecules, the algorithm generates a movement of some molecules between four compartments or sets (called L, AR, E, and S). The four compartments are shown in Figure 2 (panel 1). The algorithm starts with all molecules assigned to a set L with an initial state. At each iteration, a set of rules determines the event to be triggered next. Each event is a stochastic subprocess whereby a subset of molecules is moved from one compartment to another. When molecules reach one specific compartment (AR), their state is changed (similar to a reactor in a real chemical plant). The set of rules for the next event selection are based on the previous values of the objective function and the objective function of the best value found so far. Before describing the algorithm in detail, the different types of possible triggered events are described. Let g j x be the state of the molecule j for the best value found so far (or initial trial), and 1 ( ,..., ) gg g V xxx = the vector of molecular states for the best value found. Let F be the objective function to be minimized. Figure 2 panel 2 shows an example of the state of six molecules for the best value found so far. One possible event consists of a set of molecules being transferred from the Load tank (L set) to the Activation Reactor (AR set), in which the molecules change state to a new random state, not change while they are inside AR. An example of this type of event is shown in Figure 2 panel 2, where molecules 1 and 4 were moved from L to AR and their states changed from 0 to 1 and 1 to 0, respectively (compare panel 1 and panel 2). This event generates a new trial vector as shown in Figure 2 panel 3, the performance of which is evaluated. In another type of event in a different iteration, some reacted molecules can be sent to the Extraction unit (E set) where they are deactivated back into their previous state upon entering the reactor ( g t j j xxj E =∀ ∈ ). These extracted molecules could be sent to the Separation unit (S set) or recycled back into the Activation Reactor, where the molecules are reactivated to a new state g a j j xxj E ≠∀ ∈ . At each iteration, a trial vector consists of the activated molecules in AR and the deactivated molecules in the other three sets: xxj AR =∀ ∉ . If, after any event in the current iteration, a "good batch" is accomplished (i.e., a better objective function is found, () () tb FF θ θ ≤ , the activation reactor can be emptied into the separator unit, S. In this case all molecules conserve the new state ( g a j j xxj AR =∀ ∈ ), and a new "batch" is then started. The algorithm is described in detail in the next section.  xxx j E = ≠∀ ∈ ) and transfer all molecules in E to AR ( E =∅) 6.5 Exit the inner loop if AR is too small or if the ratio of the number of iterations relative to the initial size of AR exceeds a given parameter, RRT. Otherwise go to Step 6.1. 7. If the size of L is less than a parameter LT, transfer all S molecules to L. In Step 1, 10 FV c =×; in Step 6.1 2 0 i FA Rc = × . In both cases, if the number of molecules selected is larger than the set, all molecules in the set are selected. The parameters used for this algorithm are: RRT = 1.0, c o = 0.3, c i = 0.25, LT = V/2. The algorithm was shown to be fast and robust when tested with problems of different degrees of multi-modality, discontinuity and flatness. The molecular representation allows the solution of a large class of problems. This structure is general in purpose but also has the flexibility to add problem-specific features. For example, the "locality" of these operators allows the inclusion of bias in the sub-set formation (Steps 3 and 11) or the transformation rule (Step 5).

Algorithm performance
This algorithm has been tested and extensively utilized to solve many optimization problems. Its performance in some of the test problems is reviewed in this section. The multi-modal random problem generator of Spears (Spears, 1998) was utilized to test LARES over various degrees of modality for binary representation. The problem generator generates a set of P random V-bit strings representing the location of the P peaks in space.
To evaluate the performance of an arbitrary string, the nearest peak is located (in Hamming space). Then the fitness of the bit string c is calculated as the number of bits the string has in common with that nearest peak, divided by V. The optimum fitness for an individual is 1.0.
The objective function used in LARES was () 1 -() Fc f c = , while -F(c) was used for the fitness function in GA simulations. Table 1 shows the results for four study cases. For each set of parameters V and P, 20 random problems were generated in each case. Each algorithm was run on each problem generated. LARES found the global maximum in all cases ( (* ) 0 fc = ), while GA failed to find the global maximum for cases 3 and 4, and μGA failed to find the global maximum in three out of four cases. For the first case, LARES converged to a global optimum in 78 function evaluations on average, while GA converged in 900 function evaluations and μGA converged in less than 350 function evaluations. For the second case, LARES found the global maximum in 647 evaluations on average, while GA converged in approximately 3,700 evaluations and μGA failed to find the global maximum in 20,000 function evaluations (see Figure 3). For the third and fourth study cases, LARES was the only algorithm that converged to the global maximum in nearly 30,000 function evaluations. This behavior was explored systematically by De Jong et al. (1997). In their analysis, the authors found that for V = 20, the simple GA will converge in less than 5,000 function evaluations. For V = 100, many trials failed to find the global optimum after 20,000 evaluations.
The algorithm has also being tested with Boolean satisfiability problems (SAT), which refers to the task of finding a truth assignment that makes a Boolean expression true. The Boolean satisfiability problem generator of Mitchel et al. (1992) was used to test the performance of LARES in solving random problems with different levels of epistasis. The model assumes a conjunctive normal form of the Boolean expression with C clauses. All clauses are also assumed to consist of the same number of literals, L. The vector of variables V is represented as a binary string. A random problem is generated to create C random clauses. Each clause is generated by randomly selecting L variables, and then each variable is negated with probability 0.5. Once a random L-SAT problem is defined, the fitness function, f, is given by the fraction of clauses that are satisfied by the assignment. Note that the main goal of this section is to study LARES with different levels of epistasis. For practical solution of this type of problem, methods such as GSAT (Selman and Kautz, 1993) and WSAT (Gottlieb et al., 2002) Table 2 shows the solution of a series of L-SAT random problems using LARES and GA. Each test consists of an average of over 20 randomly generated problems. In all simulations, the length of clauses, L, had a fixed value of 3. The number of variables, V, was also fixed at a value of 100. The number of clauses was used as a parameter in the simulation, ranging from 200 to 2400. LARES was faster than GA in all cases, but in the last two cases GA found a slightly better solution while μGA found a slightly worse solution to the L-SAT problem (see Figure 9). These results indicate that LARES also behaves very well with problems involving different levels of epistasis. The LARES algorithm was also applied to a very challenging test bed used by many authors to test real function optimization algorithms. Binary encoding was used to represent real variables. The algorithm was compared with GA using the same test bed, starting with the same initial guesses, and performing the same number of iterations. Comparisons are also made with other methods specifically designed for real-function optimization reported in the literature. Although literature in this field is extensive, few studies involve methods that are efficient for real function optimization. Reported algorithms include Differential Evolution (DE) (Storn and Price, 1997), the Breeder Genetic Algorithm (BGA) (12 Mühlenbein and Schlierkamp-Voosen, 1993), Evolutionary Algorithm with Soft Genetic operators (EASY) (Voigt, 1995), the Line-up Competition Algorithm (LCA) (Yan and Ma, 2001), Continuous Ant Colony Optimization (CACO) (Mathur et al., 2000), Adaptive Simulated Annealing (ASA) (Ingber and Rosen, 1992), Very Fast Simulated Annealing (VFSA) (Ingber, 1993), Guided Evolutionary Simulated Annealing (GESA) (Yip and Pao, 1995) and Covariance Matrix Adaptation Evolution Strategy (CMA-ES) (Hansen, Muller and Koumoutsakos, 2003). In summary, LARES had better performance than GA in most instances, and in many cases the speed of LARES is comparable to that of methods specially designed to operate with real-value optimization problems.

Solution of dynamic optimization problems using LARES
As discussed in the introduction, dynamic optimization is a very important type of optimization problem in operation research and engineering, since many systems of interest are dynamic. In particular, most processes in the chemical industry are batch processes in which optimal reactant addition (and/or temperature) profiles determine product quality. Dynamic optimization is also used in model predictive control systems. These types of problems are more effectively solved using stochastic optimization methods that can escape from local minima and are not affected by singularities. In Banga, Irizarry-Rivera, and Seider (1998), an efficient and robust algorithm was developed to solve dynamic optimization problems using a flexible parametrization of the control law, consisting of a piecewise variable-length linear function. This method resulted in a big improvement over the more traditional piecewise constant approximations (Roubos et al., 1999;Luus, 2000). A generalized algorithm that uses LARES with a very flexible control law representation has also been considered (Irizarry, 2005a). This algorithm is reviewed in this section. Unlike standard optimization to determine the optimal value of a set of real variables, in dynamic optimization, we seek an optimal function u(t) or a set of functions u i (t), i= 1,M. The dynamic optimization problem for a single control law can be formulated as: www.intechopen.com

Global and Dynamic Optimization using the Artificial Chemical Process Paradigm and Fast Monte Carlo Methods for the Solution of Population Balance Models 395
Find u(t) over subject to: where F is the performance index, u is the control law, and x the vector of state variables. The set of constraints consists of the dynamic model (Eqn. 3b), the initial conditions of the state variables (Eqn. 3c), the equality constraints (Eqn. 3d), the inequality constraints (Eqn. 3e), and bounds on the control variables (Eqn. 3f). Different methods to solve this type of problem have been reviewed recently by Banga et al. (2003).

LARES-PR algorithm
The previously introduced algorithm (Irizarry, 2005;Irizarry, 2006) consists of interfacing the LARES algorithm with a generalized representation of the control law. This procedure decodes the LARES decision variables (molecules) into a flexible representation of the control law based on three key elements: (a) variable-length segments, (b) the use of finite element trial functions to represent the control function in each segment (Zienkiewics, 1977), and (c) switching between different representations to model each segment with different functions. Figure 5 shows an example in which the possible profile is represented by three segments of different lengths. In the first segment, the control function is modeled with a quadratic finite element. The second segment is modeled with a constant function (step function). The third segment is modeled with a linear finite element. In this representation, the segment sizes, the type of function representing each segment, and the adjustable parameters of the selected function for each segment are all decision variables of the optimization problem to be solved. This representation spans a large functional space in which smooth regions, drastic changes in functionality, singularities, and discontinuities of the control function can be found as part of the solution for the optimization problem with a reduced number of decision variables. The unknown control profile is encoded according to the following procedure. For the segments represented with finite elements, the node values of the finite element are part of the decision variables. Molecules are assigned to encode each of these variables using binary encoding. The function selection for each segment is then performed as follows. The data structure starts with all segments represented by finite-element function using the highest order of elements to be considered in the analysis. Then, a logical variable is defined for each segment, which is used as a switching mechanism. The logical variable can replace the

Variable-length segments
Node location is part of the problem Time u(t) element with a lower-order element or with user-defined functions over the same segment interval. Figure 6 illustrates the hybrid formulation. This example consists of quadratic finite elements with three nodes representing each element. The logical variable with state (1,1,2,3,0) δ = replaces the first two elements with lower-order linear elements (eliminating the middle node as a variable), the third and fourth elements are replaced with user-specified functions, and element 5 is represented with a quadratic element. The combinatorial variables for each finite element, i δ , is an integer number whose range equals the number of possible functions to be used. One molecule is selected for each combinatorial variable. The number of states for the molecule equals the number of possible functions that can represent a segment.

Master finite elements
). These variables cannot be encoded directly into LARES, and if standard binary encoding of these parameters in the range is utilized, this constraint will be violated frequently during LARES iterations. This problem is avoided by the following two-step procedure, called Moving Partition (MP) transformation (Irizarry, 2005a where the parameter β is used to control the gap between the disjoint segments. With this boundary for each partition node and the trial vector, t i s , the actual trial partition is calculated from the following MP mapping, ,, The MP is shown schematically in Figure  With this description of the control law, the LARES-PR algorithm can be described as follows: LARES-PR algorithm. The overall algorithm is discussed in Irizarry (2005). It consists of interfacing this profile representation with LARES. After a new trial molecular state from LARES, the procedure described in this section consists of (1) decoding, (2) applying MP transformation, (3) building the control law determined by element type, size, and parameters, (4) integrating the model, and (5) feedback to LARES regarding the performance of the control

Simulation results
LARES-PR performance has been studied with a set of benchmark problems with low sensitivity of the objective function, bang-bang behaviors, singular arc, and discontinuities in the optimal profile. In all cases, the algorithm has proven to be efficient and robust. Figure 8 shows the solution of four optimization problems used by several authors as benchmark problems. The Van der Pol oscillator problem has been studied by Vassiliadis (1993), Tanartkit and Biegler (1995), Banga, Irizarry-Rivera and Seider (1998), and Vassiliadis et al. (1999). This problem was solved using the generalized control function, where each element can be represented by either linear or quadratic Lagrange polynomials. The optimal profile is shown in Figure 8a, with an improved performance index over other methods using only four elements and a smoother profile compared to previous results reported in the literature.

www.intechopen.com
Stochastic Optimization -Seeing the Optimal for the Uncertain 398 u(t) The second case shown in Figure 8b is a plug-flow reactor with singular arc. In this problem, a plug-flow reactor packed with a mixture of catalysts is used to perform the reaction ABC ⇔⇔. The fraction of catalyst is adjusted throughout the reactor to maximize the product C. This problem was solved using the generalized control law approximation with three possible functions for the element E i : δ i = {1, 2, 3} = {u(t) = constant finite element, u(t) = u max , u(t) = u min }. The optimal control law is shown in Figure 8b. Figure 8c shows the optimal production of secreted protein in a fed-batch reactor. This problem consists of a bioreactor operated in fed-batch mode studied by Park and Ramirez (1988), Luus (1992), Banga, Irizarry-Rivera and Seider (1988), Vassiliadis et al. (1999), and Sarkar and Modak (2003). This problem shows very low performance index sensitivity of the control profile, often leading to computational difficulties particularly when gradient-based algorithms are used. The fed-batch reactor problem was also solved using the generalized control law approximation, with three possible functions for the element E i : δ i = {1, 2, 3} = {u(t) = quadratic finite element, u(t) = u max , u(t) = u min }with 10 elements. The control law gives a global optimum for this problem. Figure 8d shows the optimal profile for the bang-bang control problem (see Irizarry 2005a for a detailed description of the problem). To solve this problem, eight elements were used, and the function approximation of each element is either a linear interpolation function or the bang-bang constant functions δ i ={1, 2, 3}={u(t)=linear trial function, u(t)= u max , u(t)= u min }. LARES-PR found the correct bang-bang feature as part of the solution, that is, δ*= {4,3,4,3,4,3,4}. This is an important element of the proposed method, which can be used with general-purpose approximation functions or with problem-specific functionalities for which the proposed algorithm will identify problem features in addition to the solution.

Molecules in trial state
In most cases, the near optimum value (less than 0.5% of global optimum) was found in less than 1,000-10,000 function evaluations. The algorithm continues to refine the solution at a slower rate, resulting in very accurate solutions. In most cases, a very accurate solution can be found in 50,000-100,000 iterations. The results demonstrate that LARES-PR is robust and has fast convergence properties when compared with other stochastic optimization methods. LARES-PR has also being applied to large-scale optimal control problems with discrete-time dynamics and multiple control laws. The dynamic integrated climate-economy (DICE) model for global warming (Nordhaus, 1994) is a model of a very important problem, which posed several challenges in finding the optimal profile. This model consists of maximizing the discounted sum of per capita utilities consumption subject to the dynamics of emissions, economic impact, and economic cost of policies to control global warming. Moles, Banga and Keller (2004) made an extensive study of optimal policy with a modified version of this model using different global optimization algorithms (ICRS, LJ, DE, SRES, GLOBAL, GCSOLVE). As discussed in Moles, Banga and Keller (2004), the numerical solution of this multimodal NLP is very challenging, due to the non-convexities and discontinuous nature of the dynamics. As the time horizon is discrete, the dynamic optimization problem can be formulated as a standard NLP problem with the value of the control laws at each discrete time as a decision variable. Using this approach, the number of decision variables increases as the time horizon, N t , increases (number of decision variables = N t *N u ), resulting in a large-scale nonlinear optimization problem. Alternatively, the profile representation of LARES-PR can be utilized to represent the control law with a very small set of decision variables. Figure 9 shows the performance of LARES, DE, CRS, and LARES-PR in solving the original DICE model for a time horizon of 50 decades. As shown in Figure 9, LARES and LARES-PR are faster than DE and CRS in converging to a near global optimum. In particular, LARES-PR was much faster than all methods with a high-quality solution: The best value found for each algorithm was: LARES-PR effectively solved this problem, which consisted of finding the optimal functionality of two simultaneous control laws. To implement multiple control laws, a representation is defined for each unknown profile (PR 1 and PR 2 ).

Monte Carlo methods for population balance problems
A great majority of products are composed of finely divided solids or contain finely divided particles as part of their composition. One example is metallic microparticles and nanoparticles used in electronic ink compositions. Another example is the dispersion of pigments in paints. Most pharmaceutical products include a crystallization of organic powders. The particle size distribution, shape, and composition of these finely divided particles control the properties of the final product. Therefore, the understanding of these particles and how they form is of great importance (Irizarry, 2010a). The macroscopic modeling of particle formation consists of formulating population balance equations for the problem at hand. When optimization of these systems is pursued, the population balance equations appear as part of the constraints (Irizarry, 2005;Irizarry, 2006).
Population balance models are continuity equations of a particle population evolving by different mechanisms (such as aggregation, breakage, nucleation, and growth). The continuous population balance equation, PBE, is a deterministic integro-differential equation that describes the dynamics of a particle density function as a function of continuous particle properties (i.e. volume, particle radius, surface area, etc.). As the dimensionality of the PBE increases, the direct numerical solution of these equations becomes more difficult. For a multidimensional population balance equation, the Monte Carlo (MC) solution is an attractive alternative (and in many cases the only option). In these methods, the system evolution is modeled by a simple stochastic game, which is robust and easy to implement (Gillespie, 1975;Garcia et al., 1987). For systems close to the thermodynamic limit, both the MC solution and the direct numerical solution of the PBE converge to the same results. In many situations of practical interest, the MC solution may become very slow. Several optimization approaches have been developed to increase the MC simulation speed. The point ensemble Monte Carlo (PEMC) algorithm and the τ−PEMC algorithm developed in 2007 (Irizarry, 2007a;Irizarry, 2007b) are approximated MC methods that increase simulation speed by orders of magnitude when compared with existing MC methods. Population balance models can be formulated as discrete events Markov processes (also known as a jump Markov process). The standard exact simulation method (exact MC) for jump Markov processes consists of selecting the time for the next event and the type of event sequentially until a final time is reached (one trajectory). Many trajectories are calculated to generate the probability distribution function of the Markov process variables. This simulation method uses the propensity functions of each event, s E , defined as follows: R(E s ) dt ≡ the probability that the event E s occurs in the time interval (t, t + dt).
The time for the next event, τ, is sampled from an exponential distribution: where R Σ is the total propensity of all possible events ( ). The probability of the event, i E , occurring next (i.e. that particles will aggregate) is proportional to the event propensity, () i RE : In the inverse method, this distribution is sampled using a uniform random variable (0,1) r ∈ and then solving the following equation: is the indexed list of all possible events, and T is the total number of events. The solution of this equation, f E , is the next event to be executed at time t τ + .
Alternatively, the acceptance-rejection method can be used to sample the next event [5].
This simulation procedure has been used to develop the stochastic simulation algorithm, SSA, for chemical kinetics (Gillespie, 1976;Gillespie, 1977). In this method the firing of a chemical reaction represents a discrete chemical kinetic event of the Markov process. This algorithm is also known as kinetic Monte Carlo (KMC). The exact MC has also been used for the MC solution of population balance models. In many situations of practical interest, the MC solution may become very slow, especially when the number of particles in the simulation box is increased and the total number of events becomes very large or when the computational cost of calculating all the rates, R Σ is large. In these cases the calculation of Eqs. (1) and (3) becomes very computationally expensive, slowing down the generation of trajectories. To further accelerate the MC simulation of population balance models, a new approach was introduced in 2007 (Irizarry, 2007a;Irizarry, 2007b). These algorithms are based on the construction of a jump Markov process called PERP, which approximates the actual jump Markov model. These algorithms are shown to reduce CPU time by orders of magnitude without sacrificing simulation accuracy, when compared with optimized exact MC methods. Unlike other coarse graining (or lumping) strategies in which information and identity is lost, in these algorithms, the history of each particle is retained while a coarse view of the process is taken. These two algorithms are summarized in the next section.

Fast Monte Carlo algorithms
The PEMC and τ-PEMC algorithms are based on the simulation of an approximated jump Markov process called PERP (Irizarry, 2007a;Irizarry, 2007b). This approximated Markov process is based on three ideas. First, the total population is "discretized" into subpopulations of particles with sizes of specified intervals. Each subpopulation is viewed as a "chemical species" with the number of particles in the subpopulation representing the number of molecules of that species in the simulation volume. Second, the inter-particle interactions (i.e. aggregation, nucleation, breakage) are viewed as a set of special types of reaction, in which the reaction products are allocated stochastically to the existing species using probability functions that are mass conserving on average. Third, the original set of subpopulations is coupled with the system of "chemical species". The PERP Markov process is described next.

Point ensembles
Reacting pseudo-species. RPC example: In this approximated process, the discrete events consist of a set of "reaction channels" where the reaction part involves the pseudo-species to mimic the actual particle event (i.e. aggregation between two particles). Unlike standard reaction channels of chemical kinetics, the product component consists of several steps, some of them also stochastic processes.
In the exact Markov process, an event, E, is defined in terms of the possible interactions between the simulation particles. For example, an aggregation between two particles i and j is an event that creates a new particle in the simulation volume ( : new i j Ex x x =+) and eliminates the mother particles from the simulation volume. The reactant component of the RPC mimics the same event, but between the pseudo-species instead of actual particles. For example, the events representing the aggregation mechanism in the original Markov process are replaced by an "aggregation reaction" of pseudo-species ( : si j ES S + ) in the PERP Markov process. In the case of aggregation, the propensity function for these "reaction channels" is given by: The propensity functions for other mechanisms have also been described (Irizarry, 2007a). Notice that the propensity of each event in the PERP Markov process is given only in terms of N. For the aggregation case discussed here, the PERP event is summarized in the following algorithm: . The product probability parameter, s P , is calculated using the mass conservation equation for this landing interval: All these events in the product component of the RPC are simply called PERP events. RPCs can be defined for any mechanism (Irizarry, 2007a) The PEMC algorithm is the exact MC simulation of the approximated PERP Markov process. A detailed description of these steps has been published (Irizarry, 2007b). Let ,1 , . . , i Ei T = be the list of all possible RPCs describing the population balance model at hand. The PEMC method is summarized in the following algorithm: Algorithm M2: PEMC (one iteration). 1. Find the time to fire the next RPC using Eq. (7). Update the time ttτ = + . 2. Find the RPC to be fired next solving Eq. (9). 3. Fire the selected PERP event (Algorithm M1) These steps are repeated until the final time is reached. This algorithm is very fast because the set of RPCs is more compact than the set of all possible events between particles, making the calculations of Eq. (7) and (9) very fast.

τ-PEMC
The τ−PEMC algorithm is a τ-leap solution of the approximated PERP Markov process. The τ-leap method is an approximated stochastic simulation, were many events are fired at once over the time interval. Consider a coarser time interval, τ, such that many events occur in this interval, but small enough that the propensity functions will not change appreciably during τ. When this condition is satisfied, all reaction channels can be considered as independent events (Gillespie, 2001), and the number of firings for each reaction, E j , is a Poisson random variable with distribution ( ) The accuracy and speed of the method depends on the selection of time τ during the simulation (Gillespie and Petzold, 2003). Since the Poisson distribution is not bounded, it could generate negative values for the concentration.
Another improvement to the method is to replace the Poisson distribution with a binomial distribution (Tian and Burrage, 2004;Chatterjee, Vlachos, and Katsoulakis, 2005): In this case the number of firings for each reaction, k j , is sampled from a binomial random The PERP process can be simulated using the τ-leap method, as previously described (Irizarry, 2007b). A pseudo-code for this algorithm is described as follows: Algorithm M3: τ-PEMC (one iteration) 1. Select the time parameter, τ .
2. For each RPC, s E , take a sample s k from a binomial distribution (Eq. (12) 3. Fire each PERP event, s E , s k times (execute algorithm M1 s k consecutive times).
Continue steps 1-3 until the final time is reached.

Performance of the PEMC algorithm with complex kernels
The numerical accuracy of these algorithms has been studied with complex coagulation kernels of physical relevance. Numerical results are compared with the generalized approximation method (GA) developed by Piskunov and Golubev (2002) and Piskunov et al. (2002). The values for the second moment by the GA method are considered the most accurate existing numerical results in the literature. They are used here as benchmark values. The following kernels are considered: i. The Brownian kernel, ii. The coagulation kernel simulating the process of migration and coalescence of particles on a heated substrate, iii. The gravitational coagulation of particles in the Stokes regime, ( ) For the gravitational kernel, some moments diverge after a critical point that depends on initial conditions. The critical point for the initial conditions used here is in the range of 0.5-0.8 (Piskunov et al., 2002).

www.intechopen.com
Stochastic Optimization -Seeing the Optimal for the Uncertain

406
The initial conditions used in the simulation for q B and q + kernels are a mono-disperse solution with unit particle volume and unit total particle concentration. For the case of q C , the following initial conditions are utilized: Table 3, the CPU time and final second moments of the simulations for different methods are compared as a function of kernel type and number of particles. For the τ-PEMC, two representative coarse-graining factors (f = 100 and f = 1,000) were examined. All calculations were performed using Visual Fortran 6.6 on an Intel Pentium M 1.6 GHz machine with 504 MB RAM. Table 3 shows that τ-PEMC can be 14-44 times faster than PEMC. For kernel q B , the τ-PEMC method gives accurate results in all cases with a coarse-graining factor f = 100. For the case f = 1,000, the number of particles had to be increased to 50,000 to generate accurate second moments (see Table 3). The q + kernel is more computationally challenging than the Brownian kernel, resulting in a wider distribution. Numerical results presented previously (Irizarry, 2007a) show that PEMC also quickly converges to the GA values. Table 3 shows the solution convergence to the exact value as a function of increasing the number of particles. For the factor f = 1000, large deviations of the second moments were observed even when a large number of particles was utilized. Good results were obtained for f = 100.
For the q C kernel, some moments diverge at a critical value (gelling point). For the initial conditions used here, this critical value is believed to be between 0.5 and 0.8 (Piskunov et al., 2002). Unlike most numerical methods, the PEMC gives accurate results for all times, especially at time 0.6, where all discretization methods diverge (see Irizarry, 2007a For simulations near gelling points, the performance of the τ-PEMC deteriorates, and the PEMC should be used in this case. Table 3 also shows results using the acceptance rejection (AR) method of Garcia et. al. (1987) for comparison. For the coagulation case, this method is very attractive because it is simple to implement and avoids calculation of all interactions between particles, making the method very computationally efficient. As noticed by other authors, the CPU time and accuracy of this method depends drastically on the problem (kernel and initial conditions). As shown in Table 1, for some cases the AR method can be very efficient (q B ) while in other cases the CPU time is very high (q + ). Additionally, the AR method diverges for the q C kernel, also shown in Table 3.

Hybrid Monte Carlo algorithm for population balance models with stochastic and deterministic variables
Most MC implementations of population balance models have focused on the solution of the PBE to approximate macroscopic variables. Less attention has been focused on the solution of population balance models where some species are far from the thermodynamic limit (very dilute or finite) and other species can be considered deterministic (high concentration).
In this case the MC is more accurate than direct numerical solution, which ignores the inherent fluctuations of the system. This type of problem often results in a stochastic system that contains both stochastic and deterministic variables with multiple timescales for the different mechanisms. In this case, the direct MC simulation will be accurate but very ineffective in terms of CPU time. Furthermore, most of the computational time is spent sampling the fast events. This type of situation has been considered in the case of chemical kinetics and biological systems for which efficient hybrid algorithms have been developed to solve multi-scale problems (Salis and Kaznessis, 2005;Kaznessis, 2006;Haseltine and Rawlings, 2002;Haseltine and Rawlings, 2005). In these algorithms, the fast processes are approximated by continuous models, and the slow processes are solved by the exact MC method in a hybrid algorithm. Disparate scales in population balance models may arise because some species are concentrated while others are very dilute. For accurate simulation of the dilute species, a large number of simulation particles are needed. In this case, the exact MC methods become very slow. A recently introduced hybrid strategy (Irizarry, 2010b) is reviewed. In this strategy, the τ-PEMC is used for the parts of the system than can be considered large, and the PEMC is used for the stochastic events.

Hybrid algorithms in chemical kinetics
The hybrid strategy for chemical kinetic problems (Salis and Kaznessis, 2005;Kaznessis, 2006;Haseltine and Rawlings, 2002) is based on partitioning between fast and slow reactions. To split the system between slow and fast events, the following criteria for fast events are utilized: 1. The fast events occur many times in a small time interval.
2. The effect of these events on the number of particles and the propensity functions is small relative to the total propensity function and the total number of particles. The slow processes are simulated using SSA, while the fast processes are integrated using the Langevin equations. To make this hybrid simulation self-consistent, the coupling between both processes must be considered. If we look at the interval between the previous slow event at time o t and the time for the next slow event ( where r is a random number from the uniform distribution (0, 1). This equation replaces Eq. (7) for the time of the next slow event. In the interval [, ] oo tt τ + , the dynamics of the fast system can be integrated in a seamless way since by definition no slow events are present in this interval. Thus, Eq. (16) becomes a constraint in the hybrid strategy that needs to be monitored while integrating the fast process. As previously discussed (Salis and Kaznessis, 2005), one way to implement this constraint is to notice that the integral term is monotonically increases, and the second term is a negative term. Therefore, the time for the next slow event can by found by monitoring the zero crossing of the residual equation:

Hybrid algorithm in multi-scale MC simulation of particulate processes
The hybrid strategy described in Section 6.1 is utilized with the slow mechanisms simulated using the PEMC algorithm. Instead of approximating the fast mechanism with a continuous model, as in the case of chemical kinetics, the τ-PEMC algorithm is utilized to model the fast mechanisms. The τ-PEMC method allows for coarse simulation in time while maintaining individual particle properties, in contrast to continuous models such as Langevin equations in which particle integrity is lost. The hybrid algorithm consists of the τ-leap integration of the PERP Markov process for fast events while monitoring the residual Eq. (17) for the firing of the next slow event of the PERP Markov process. This process is shown schematically in Figure 11. The detailed description of the algorithm is given in Irizarry (2010b).

Simulation results
Consider a system with two type of particles, A and B. B particles are much smaller than A particles. A particles can grow by an aggregation mechanism. B particles are stable from aggregation with other B particles but can condense on the surface of growing A particles. Furthermore, it is assumed that B particles are very dilute compared to A particles. These conditions make A-B condensation events a stochastic process, while A-A aggregation events can be approximated as continuous events. Figure 12 shows five instances of the test problem for the case W = 100 (See Irizarry, 2010b for details). As shown in Figures 12a and 12b  algorithm can be utilized. As the rate of aggregation is reduced with time, the hybrid algorithm correctly switches to a PEMC simulation of the entire system. The hybrid algorithm can simulate stochastic variables (A-B) at speeds approaching τ-PEMC. The statistics of the condensation of B monomers for the case W = 100 is summarized in Table 4 for 1,000 simulations. The histogram of the parameter x 1 (which measures the concentration of B particles on A particles) for the hybrid-PEMC and PEMC solutions is shown in Figure 13. The hybrid and PEMC solutions are in excellent agreement. This result is remarkable considering that condensation events are very rare (~30 condensation events vs. ~60,000 aggregation events). The box plot for these parameters is shown in Figure 14. An analysis of variance was performed to compare the population generated by both methods.
For the case of x 1 , both simulations are statistically equivalent. The hybrid algorithm is 17 times faster than the PEMC algorithm per trajectory. This increase in speed is very important, since many trajectories (~10 3 ) are needed to generate good statistics for the slow process. If the model is used for optimization, many design instances are needed (~10 4 ), each consisting of many trajectories, resulting in a very large number of simulations (~10 7 ). In this case any increase in the simulation speed of a trajectory will have a tremendous impact on the total simulation speed.

PEMC
Hybrid-PEMC Average value of x 1 5.20 × 10 -4 5.13 × 10 -4 Standard deviation of x 1 8.9 × 10 -5 9.0 × 10 -5 Average number of slow aggregation events 58,623 58,629 Average number of condensation events 31.2 30.8 Table 4. A comparison of hybrid-PEMC and PEMC for the test problem with W = 100 (Average of 1,000 trajectories).  Hybrid and MC simulations produce statistically equivalent results, but the hybrid's increase in computational speed allows optimization problems involving these types of models to be solved very efficiently.

Conclusions and discussion
This chapter reviewed the artificial chemical process paradigm for global optimization. The LARES algorithm is very robust and efficient, converging to near-global optimal solutions when solving different classes of problems with different degrees of multi-modality, epistasis, flatness, and discontinuities. Future research will consider the use of the ACP paradigm in the development of new problem-specific algorithms. The algorithm was utilized to develop dynamic optimization strategies, LARES-PR and hybrid LARES-PR. The power of the algorithm lies in its utilization of a generalized approximation of the control function, composed of variable-length segments of finite elements of different orders or using specified functions. This generalized representation of the trial control law is possible due to the two-step encoding of the decision variables and the capability of LARES for multiple encoding. Multiple encoding allows the inclusion of different types of problemspecific finite elements (constant, linear, quadratic, etc.) and/or specialized functions to approximate the control law without any tailoring of the optimization algorithm. This approach is particularly effective for the solution of problems in which manipulated variables experience transition from smooth variations over time to discrete changes. Numerical experiments demonstrate that this algorithm is robust in finding global optimums for the different types of problems and definitions of the generalized control law introduced in this work. To accelerate optimization of systems that use MC simulations as part of their constraints, a new general-purpose MC algorithm for the dynamics of the particulate process was proposed, PEMC. The method has been shown to reduce CPU time without sacrificing simulation accuracy. While a coarse view of the process is taken, particle history is retained. The method was extended with the τ-PEMC method, proposed to further improve CPU time with negligible simulation errors. As with the original PEMC, internal coordinates can be handled effectively. The CPU times reported here show that accurate results can be achieved with simulation times less than a second using a low-end PC. These results demonstrate the feasibility of stochastic optimization using PEMC and τ-PEMC. A new hybrid strategy was developed to solve stochastic population balance models with multiple time scales. This self-consistent hybrid method combines the PEMC and τ-PEMC algorithms to accelerate simulation time while capturing the stochastic nature of the slow process. The simulation speed and accuracy of the hybrid strategy depends on the selection of the τ parameter, the criteria for the partition between slow and fast events, and the grid quality of the point ensembles.