Pareto-Optimal Fronts for Simple Crystallization Systems Using

Article pubs.acs.org/IECR

Cite This: Ind. Eng. Chem. Res. 2019, 58, 14239−14251

Pareto-Optimal Fronts for Simple Crystallization Systems Using Pontryagin’s Minimum Principle Yu-Ti Tseng, Hao-Jen Pan, and Jeffrey D. Ward*

Downloaded via UNIV OF CAMBRIDGE on August 20, 2019 at 10:25:00 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Department of Chemical Engineering, National Taiwan University, Taipei 106-07, Taiwan ABSTRACT: Pontryagin’s minimum principle is applied to determine Pareto-optimal fronts for multiobjective optimization of seeded batch crystallization processes. For cases with two objectives, when both objectives are based on properties of the higher moments or lower moments of the crystal size distribution (for example, weight mean size and weight mean coefficient of variation), there is relatively little difference among the trajectories along the optimal front. By contrast, when one objective is based on the lower moments (e.g., minimizing the number of nuclei) and another objective is based on the higher moments (e.g., minimizing the nucleated mass), the difference in the shape of the trajectories along the front is more pronounced. In this case, a constant growth rate trajectory may represent a reasonable compromise between the competing objectives.

1. INTRODUCTION Optimization of batch crystallization processes has received considerable attention in the literature. Besides seed properties, the most important component of the operating recipe is the temperature or growth rate trajectory over the course of the batch. Many researchers have applied optimization routines to determine the best trajectory in some sense. Methods that have been applied can be divided into two categories: those based on control vector iteration1−11 and those based on orbital flatness12,13 and optimal control theory (Pontryagin’s minimum principle).14−22 General reviews on control and optimization of batch crystallization processes have also been published.24,25 Many different objective functions have been considered for batch crystallization. Since there are many different objective functions that might be of interest, batch crystallization is wellsuited to analysis using multiobjective optimization methods such as the construction of Pareto-optimal fronts. Several researchers have considered multiobjective optimization for batch crystallization using control vector iteration methods,26−28 and Bajcinca et al.20 constructed Pareto fronts using optimal control theory for the case of the batch time and the nucleated mass as objective functions. However, many other combinations of objective functions may be of interest to researchers, and to our knowledge there are no other reports in the open literature of researchers developing Pareto-optimal fronts for batch crystallization by using optimal control theory. Therefore, in this work, optimal control theory is applied to determine Pareto-optimal fronts for several combinations of objective functions for seeded batch crystallization.

used later. Readers are referred to refs 12−22 for additional information. The process model used in this work is quite simple: it accounts only for ordinary growth. Aggregation, breakage, size-dependent growth, imperfect mixing, and other phenomena are not considered. 2.1. Crystallization Model. The crystallization process model used in this work is based on the method of moments. The ith moment (mi/m3·s or mi/kg·s) of the crystal size distribution is given by μi =

∞

Lif (L) dL ,

i = 0, 1, 2, ...

(1)

where L is the characteristic length of particles and f(L) is the number crystal size distribution. The rate of change of the moments of the crystal size distribution is given by

dμ0 dt dμi dt

=B

= iGμi − 1 ,

(2)

i = 1, 2, ...

(3)

where G is the crystal growth rate (m/s) and B is the nucleation rate (no./m3·s or no./kg·s). The relative supersaturation is given by S=

2. THEORY This section presents a brief summary of the batch crystallization process model used in this work and a brief review of the formulation of the optimal control problem that will be © 2019 American Chemical Society

∫0

C − Csat Csat

Received: Revised: Accepted: Published: 14239

(4)

May 1, 2019 July 6, 2019 July 8, 2019 July 8, 2019 DOI: 10.1021/acs.iecr.9b02394 Ind. Eng. Chem. Res. 2019, 58, 14239−14251

Article

Industrial & Engineering Chemistry Research where C and Csat, which is a function of temperature (i.e., Csat = Csat(T)), are the solute concentration and saturation solute concentration in units of kg/kg solvent, respectively. The crystal nucleation and growth rates are given by the following empirical expressions:

Table 1. Parameters for Potassium Nitrate−Water and Pentaerythritol−Water Systems31,32 name mass of solvent density of crystals volumetric shape factor nucleation parameter nucleation parameter nucleation parameterc growth parameter growth parameter total batch time maximum growth rate constant for seed size distribution width of seed size distribution initial seed mean size initial third moment of seed crystals final third moment of seed crystals

G = kg S g

(5)

msolv ρc kv

B = kbS bμ3 j

(6)

kb

The relationship between concentration and time is given by the material balance: C(t ) = C0 − ρc k v(μ3 (t ) − μ3,s,ini )

b j

(7)

kg g tf Gmax

where kv is the shape factor that relates the volume of the crystals to the third moment, ρc is the density of crystals, and μ3,s,ini is the third moment of the seed crystals. The seed-grown (subscript “s”) and nucleus grown (subscript “n”) crystals can be tracked separately so that μi = μi,s + μi,n: dμ0,s dt dμ0,n dt

= 0, = B,

dμi ,s dt dμi ,n dt

= iGμi − 1,s

i = 1, 2, ...

a w

(8) L0

= iGμi − 1,n

i = 1, 2, ...

μ3,s,ini

(9)

The initial conditions of eq 8 are based on the initial seed moments (μi,s,ini), while the initial conditions of eq 9 are zero since there are no nucleated crystals at the start of the batch. In this work, we consider the number size distribution of seed f(L)seed as a parabolic form: f (L)seed = a(L − L0(1 − w))(L − L0(1 + w))

μ3,s,f,c a

units

1.65 2100 1

1 1400 1

kg kg/m3 dimensionless

3.47 × 104a

2.47 × 109b

1.78

3.8

dimensionless

1

0

dimensionless

6.97 1.32 360 5.71 × 10−3

3.44 1.9 420 5.71 × 10−3

mm/min dimensionless min mm/min

−1.10 × 1011

−5.72 × 108

1/mm3

0.2

0.2

dimensionless

0.1

0.2

mm

1.20 × 103

4.00 × 102

mm3

1.08 × 105

2.07 × 104

mm3

In no./mm3·min. bIn no./kg·min. cThe third moment term.

name zero moment of nuclei third moment of nuclei number mean size weight mean size number mean coefficient of variance weight mean coefficient of variance

abbrev

expression

units

mu0n

μ0,n,f

dimensionless

mu3n

μ3,n,f

mm3

NMS WMS

−μ1,f/μ0,f −μ4,f/μ3,f

mm mm

NMCOV

μ2,f μ0,f /μ1,f 2 − 1

dimensionless

WMCOV

μ5,f μ3,f /μ4,f 2 − 1

dimensionless

where α ∈ [0,1] and ϕ1 and ϕ2 are two objectives listed in Table 2 (subscript “f” indicates value at the end of a batch). A point on the Pareto-optimal front can be derived by specifying α before solving the optimization problem. Further discussion about whether the inequality constraint on batch time is active or not is shown in section 2.3. 2.3. Derivation of Pontryagin’s Minimum Principle Solution. The development in this section is based on the work of Hofmann, Raisch, and co-workers.12−21 Readers can also refer to ref 29 for mathematical details of Pontryagin’s minimum principle and examples of applying the theory. In order to apply Pontryagin’s minimum principle in this work, a transformed time τ is defined as

min ϕ = αϕ1(μ0,n (tf )...μ5,n (tf )) G(t )

(11)

s.t.

d τ = G (t ) d t ,

tf ≤ tf,c

τ(t = 0) = 0

(12)

The physical meaning of τ is the change in the characteristic length of the seed-grown particles from the beginning of the batch. By introducing the transformed time, eq 8 can be converted into the following ordinary differential equations:

μ3,s (tf ) = μ3,s,f,c G(t ) ∈ (0, Gmax ] ,

pentaerythritol

Table 2. Definition of Objective Functions Considered in This Work

(10)

where a is a constant depending on the seed loading (mass of seed added in per unit mass of volume of slurry), L0 is the initial seed mean size, and w is the width of the seed size distribution. Two crystallization systems (potassium nitrate and pentaerythritol with water as solvent) are considered in this work, and values of all parameters in these equations are given in Table 1. 2.2. General Formulation of Multiobjective Optimization Problem. In this work, we consider a cooling batch crystallization process with objective functions based on the values of moments at the end of batch. The optimization problems are subjected to a constraint on the final seed-grown volume (proportional to the third moment of seed-grown nuclei μ3,s,f,c), batch time tf,c, and interval of growth rate (0, Gmax]. Every single objective function considered in this work is listed in Table 2. The problem is a dynamic optimization problem with growth rate as input variable. The linear scalarization method is used to convert multiobjective function problems into single objective ones. The optimization problem can be stated as follows:

+ (1 − α)ϕ2(μ0,n (tf )...μ5,n (tf ))

potassium nitrate

variable

∀ τ ∈ [0, tf ] 14240

DOI: 10.1021/acs.iecr.9b02394 Ind. Eng. Chem. Res. 2019, 58, 14239−14251

Article

Industrial & Engineering Chemistry Research dμ0,s

dμ1,s

= 0,

dτ dμ3,s dτ

dτ

= μ0,s ,

dμ2,s dτ

B (μ (τ ), u(τ )) = kbkg −b / gu(τ )(g − b)/ g μ3,s (τ ) j G 3,s

= 2μ1,s ,

= 3μ2,s

Since there are no nuclei at the beginning of the batch, the initial values of x1−x6 are zero. The state x7 is also zero initially. For the value of x7 at the end of batch, if the batch time inequality constraint is active, x7(τf) = tf,c. The boundary conditions for eqs 19−21 can therefore be stated as

(13)

Combining eqs 12 and 13, values of seed moments (μ0,s to μ3,s) can be expressed as a function of τ:

x(0) = 0,

μ0,s (τ ) = μ0,s,ini

min ϕ = αϕ1(x1(τf ), ..., x6(τf ))

μ3,s (τ ) = μ0,s,ini τ 3 + 3μ1,s,ini τ 2 + 3μ2,s,ini τ + μ3,s,ini

u(τ )

+ (1 − α)ϕ2(x1(τf ), ..., x6(τf ))

The seed-grown volume constraint stated in eq 11 can thus be restated as μ3,s(τf) = μ3,s,f,c; τf is the value of τ at the end of batch, which can be solved explicitly. The control input u is defined as the inverse of the growth rate: 1 G (τ )

x(0) = 0,

d x = 5ψ1x 2 + 4ψ2x3 + 3ψ3x4 dτ B + 2ψ4x5 + ψ5x6 + ψ6 (μ3,s (τ ), u(τ )) + ψ7u G

H (x , u , Ψ , τ ) = ψ T

The states of the system are defined as follows:

x5 = μ1,n ,

x6 = μ0,n ,

d ∂ ψ=− H dτ ∂x

(17)

The key simplification is to neglect the effect of the nucleated mass on the nucleation rate, i.e., to assume that B(u , μ3 ) ≈ B(u , μ3,s )

dx 2 = 4x3 , dτ

d d d ψ = 0, ψ = −5ψ1 , ψ = −4ψ2 , dτ 1 dτ 2 dτ 3 d d d ψ4 = −3ψ3 , ψ5 = −2ψ4 , ψ = −ψ5 , dτ dτ dτ 6 d ψ =0 dτ 7

(18)

dx6 B = (μ3,s (τ ), u(τ )) dτ G dx 7 = u(τ ) dτ

dx 3 = 3x4 , dτ

(26)

Combining eqs 25 and 26 yields

Hofmann and Raisch18 showed that if the objective is to minimize μ3,n of the potassium nitrate−water system, the error caused by the simplification in solving the optimal control problem is less than 3% even when the volume ratio of nucleated nuclei to seed-grown at the end of batch (μ3,n(tf)/μ3,s,f,c) is 33%. Furthermore, the nucleation rate of the pentaerythritol− water system is independent of magma density (equals ρkvμ3), implying that the simplification does not have any effect on the result. Based on this assumption and the moment process model presented previously, rates of change of the states with respect to τ can be written as dx1 = 5x 2 , dτ dx5 = x6 dτ

(25)

where ψ is called the costate (or the adjoint state). Differential equations for solving costates are

x4 = μ2,n ,

x7 = t

∀ τ ∈ [0, τf ]

The Hamiltonian is defined as

(16)

x3 = μ3,n ,

x 7(τf ) = tf,c

u(τ ) ∈ U = [umin , ∞)

The real time t in terms of τ is tracked by the following equation: dt 1 = = u(τ ) dτ G (τ )

(24)

s.t.

(15)

x 2 = μ4,n ,

(23)

where x = [x1, x2, ..., x7] . Based on eqs 12−23, the optimization problem can be restated as

(14)

μ2,s (τ ) = μ0,s,ini τ 2 + 2μ1,s,ini τ + μ2,s,ini

x1 = μ5,n ,

x 7(τf ) = tf,c T

μ1,s (τ ) = μ0,s,ini τ + μ1,s,ini

u(τ ) =

(22)

(27)

In order to derive boundary conditions for eq 27, define the terminal constraint functions ε as eq 28. ε is a vector function whose components are constraint on the final value of state variables. In this work, the only component in ε is ε1 = x7(τf) − tf,c = 0. ε(x(τf ), τf ) = 0

(28)

Variational analysis, which investigates increments in objective with respect to states, costates, control input and time, shows that the necessary terminal condition for the optimal solution is

dx4 = 2x5 , dτ

ij ∂ϕT yz ∂ε T jj + ν − ψ zzz jj ∂x z ∂x k {

(19)

T

δ x(τf )

ij ∂ϕT yz ∂ε T + jjj + ν + H zzz j ∂τ z ∂τ k {

(20)

(21)

The expression of B/G in terms of the control input u can be derived from the process model:

=0 14241

τf

T

δτf τf

(29) DOI: 10.1021/acs.iecr.9b02394 Ind. Eng. Chem. Res. 2019, 58, 14239−14251

Article

Industrial & Engineering Chemistry Research

Figure 1. Results of prespecified trajectories of potassium nitrate system: (a) temperature trajectory; (b) growth rate trajectory; (c) NSDn,f; (d) VSDn,f.

l g /b o yz o o 1 ijj ψ6 b − g j o z o j z if u 0 ≥ umin kbμ3,s (τ )zz o jj k g ψ u* = o m g o k 7 { o o o o o o if u 0 < umin n umin

where vector ν is a Lagrange multiplier associated with the terminal constraints, δτf is the variation of the final transformed time, and δx(τf) is the variation of the final state. Notice that τf is fixed in this work, indicating that the second term in eq 29 is zero. A necessary condition for an optimal input trajectory is that the optimal control input u* must satisfy

After solving the optimal control problem, the temperature trajectory corresponding to optimal control input can be derived by combining the definition of supersaturation (eq 4), material balance (eq 7), and solubility function Csat(T). The final number size distribution of the nuclei (NSDn,f) and final volume size distribution of nucleated nuclei (VSDn,f) can be also derived by B B plotting G (τ ) versus (τf − τ) and G (τ )(τf − τ )3 versus (τf − τ), respectively.

u*(τ ) = argmin H(x , u , ψ , τ ) u∈U

B = argmin ψ6 (μ3,s (τ ), u(τ )) + ψ7u G u∈U = argmin ψ6kbkg −b / gu(τ )(g − b)/ g μ3,s (τ ) j + ψ7u u∈U

(30)

Suppose the constraint x7(τf) = tf,c is relaxed, then it can be deduced that ψ7(τ) = 0 for τ ∈ [0,τf] by eqs 27 and 29. Since b > g > 0 (nucleation rate is more sensitive to supersaturation than growth rate) in the case studies, if there exists a τ+ ∈ [0,τf] such that ψ6(τ+) > 0, eq 30 attains its minimum when u → ∞, leading to infinite value on x7(τf). This suggests that the batch time constraint is active whenever τ+ exists. Furthermore, we can conclude that ψ7(τ) = ν based on eqs 28 and 29. A candidate for the unconstrained minimizer u0 can be found by setting the partial derivative of H with respect to u to be zero:

3. RESULTS 3.1. Commonly Considered Temperature Trajectory. To compare results of optimized trajectories and commonly considered ones, four kinds of prespecified trajectories are discussed: linear cooling (Linear), cubic cooling (Cubic), constant growth rate (ConstG), and Mullin−Nývlt (M−N).30 All trajectories are subjected to the same final seed-grown volume and batch time constraints. Figures 1 and 2 show the NSDn,f, VSDn,f, temperature trajectory, and growth rate trajectory results for the potassium nitrate and pentaerythritol systems. The linear trajectory leads to an early growth result (growth rate is larger in the beginning of batch), while the cubic and M−N trajectories lead to a lategrowth result (growth rate is larger in the end of batch). For early growth trajectories, the number of nuclei formed will be relatively small while the total third moment will be relatively large and vice versa. This result is consistent with previous work.5,6,22,23 3.2. Problem 1. Third Moment of Nucleated Nuclei (mu3n) vs Zero Moment of Nucleated Nuclei (mu0n). The two objectives in problem 1 are to minimize the third moment and the zero moment of nucleated crystals at the end

b−g ∂ H = −ψ6 kbkg −b / g(u 0)−b / g μ3,s j (τ ) + ψ7 = 0 0 g ∂u u (31)

Solving eq 31 yields

yz 1 ij ψ b − g kbμ3,s j (τ )zzzz u = jjjj 6 kg k ψ7 g {

(33)

g /b

0

(32)

0

If u < umin, the constraint on the admissible range of control input is active. The expression of optimal control input u* is 14242


Article


Figure 2. Results of prespecified trajectories of pentaerythritol system: (a) temperature trajectory; (b) growth rate trajectory; (c) NSDn,f; (d) VSDn,f.

of the process with fixed final seed crystal mass and the fixed batch time. Problem 1 can be stated as min ϕ = αx3(τf ) + (1 − α)x6(τf )

(34)

u(τ )

s.t. x(0) = 0,

x 7(τf ) = tf,c

u(τ ) ∈ U = [umin , ∞)

∀ τ ∈ [0, τf ]

Combining eqs 27, 29, and 34, the following expressions for the costates can be obtained: ψ1(τ ) = 0,

ψ2(τ ) = 0,

ψ4(τ ) = −3α(τ − τf ),

ψ3(τ ) = α , ψ5(τ ) = 3α(τ − τf )2 ,

ψ6(τ ) = −α(τ − τf )3 + (1 − α),

ψ7(τ ) = ν

Figure 3. Pareto-optimal front for problem 1, potassium nitrate system. (35)

Notice that ψ6(0) = ατf + (1 − α) > 0. Taking τ = 0, the batch time constraint is guaranteed to be active. Substituting these expressions for ψ6(τ) and ψ7(τ) from eq 35 into eq 32, then the corresponding states x(τ) can be obtained by integrating eqs 19−21 and the value of the Lagrange multiplier ν can be determined by x7(τf) = tf,c. Figures 3 and 4 show the Pareto-optimal front for both crystallization systems of problem 1. Results for the four commonly considered trajectories described in section 3.1 subject to the same values of μ3,s,f,c and tf,c are also shown in Figures 3 and 4. The four points all lie above the Paretooptimal front as expected since they are suboptimal except for the constant growth rate trajectory for the pentaerythritol system, which has the same objective function values as the curve minimizing mu0n. Since the nucleation rate of the pentaerythritol is independent of the value of the third moment in the model, the trajectory that minimizes the number of nuclei is one of constant growth. The Pareto front is convex, and there are points along the Pareto front that are near the minimum value for each 4

+

Figure 4. Pareto-optimal front for problem 1, pentaerythritol system.

objective function separately. In order to identify a point that might be considered to represent the “best” trade-off between the two objective functions, we define a parameter θ as follows: 14243


Article


“middle” on the Pareto front in Figures 3 and 4 is the point on the front for which the value of θ is minimal. Values of θ can also be calculated for points not on the Pareto front, and the value of θ for four prespecified temperature trajectories and the middle point for each optimization problem are shown in Tables 3 and 4. Figures 5 and 6 show the growth rate and temperature trajectories as well as the final nucleated particle size distribution for individual single-objective optimizations and for the middle point on the Pareto-optimal front for both systems. For the potassium nitrate system, the trajectory that minimizes mu0n is an early growth trajectory while the trajectory that minimizes mu3n is a late-growth trajectory. This is consistent with previous findings22,23 and reflects that fact that if batch time and production rate (seed-grown volume) are constrained and the nucleation rate increases with increasing magma density, a higher supersaturation at the beginning of the batch results in fewer nuclei formed but a greater nucleated mass. This is because nuclei formed earlier have an opportunity to grow to a larger size. By contrast, when the supersaturation is greater near the end of batch, the number of nuclei formed is greater but the nucleated mass is smaller due to the small volume of these nuclei. The growth rate trajectory corresponding to the middle point in Figure 5 is relatively flat, suggesting that a constant growth rate policy may represent a reasonable tradeoff between objective functions in this case. The point corresponding to a constant growth rate is not far from the Pareto front and the middle point in Figure 3, also. For the pentaerythritol system, the model predicts that the nucleation rate does not depend on magma density, implying that neither early- nor late-growth trajectories are preferred when minimizing the number of nuclei. Thus, the constant growth trajectory minimizes mu0n. The middle point thought to represent the best trade-off between objective functions corresponds to a trajectory where the growth rate increases gradually during the batch.

Table 3. Values of Percentage Sum Function θ for Potassium Nitrate−Water System problem ϕ1 ϕ2

1

2

mu3n mu0n

−NMS mu3n

3

4

WMCOV −WMS θ value

NMCOV mu0n

trajectory

problem 1

problem 2

problem 3

problem 4

middle Linear Cubic ConstG M−N

17.23 41.78 33.20 19.49 25.51

17.04 43.03 26.94 18.16 21.74

0.04 11.36 8.22 4.24 6.08

5.98 7.40 37.83 30.83 34.67

Table 4. Values of Percentage Sum Function θ for Pentaerythritol−Water System problem ϕ1 ϕ2

1

2

mu3n mu0n

−NMS mu3n

3

4

WMCOV −WMS θ value

NMCOV mu0n

trajectory

problem 1

problem 2

problem 3

problem 4

middle Linear Cubic ConstG M−N

27.17 478.33 214.16 52.91 171.77

36.70 371.34 167.14 65.69 138.64

0.38 61.86 44.74 10.00 38.69

52.50 126.48 117.24 68.34 108.56

ij ϕ − ϕ ϕ2 − ϕ2,opt j 1 1,opt + θ = 100jjjj j ϕ1,opt ϕ2,opt k

yz zz zz zz {

(36)

where ϕ1,opt and ϕ2,opt refer to optimal values under optimization of the single objective function. The point marked

Figure 5. Results for individual single-objective optimizations and middle point on the Pareto-optimal front for problem 1, potassium nitrate system: (a) temperature trajectory; (b) growth rate trajectory; (c) NSDn,f; (d) VSDn,f. 14244


Article


Figure 6. Results for individual single-objective optimizations and middle point on the Pareto-optimal front for problem 1, pentaerythritol system: (a) temperature trajectory; (b) growth rate trajectory; (c) NSDn,f; (d) VSDn,f.

3.3. Problem 2. Number Mean Size (NMS) vs the Third Moment of Nucleated Nuclei (mu3n). Problem 2 can be stated as ji μ1,s,f + x5(τf ) zyz zz min ϕ = αx3(τf ) + (1 − α)jjjj− j μ u(τ ) + x6(τf ) zz 0,s,f k {

(37)

s.t. x(0) = 0,

x 7(τf ) = tf,c

u(τ ) ∈ U = [umin , ∞)

∀ τ ∈ [0, τf ]

Expressions for the costates can be derived from eqs 27, 29, and 37: ψ1(τ ) = 0,

ψ2(τ ) = 0,

Figure 7. Pareto-optimal front for problem 2, potassium nitrate system.

ψ3(τ ) = α ,

ψ4(τ ) = −3α(τ − τf ), ψ5(τ ) = 3α(τ − τf )2 + (1 − α)k1 ,

(38)

ψ6(τ ) = −α(τ − τf )3 + (1 − α)[−k1(τ − τf ) + k 2],

ψ7(τ ) = ν

where the constant k1 and k2 are μ0,f 1 k1 = − , k2 = μ0,f μ1,f 2

(39)

Since k2 > 0, we can take τ = τf so that the batch time constraint is guaranteed to be active. Figures 7 and 8 show the Pareto-optimal front for problem 2. Notice that the y-axis is the negative of the value of NMS because NMS should be maximized. Again the Pareto front is a convex and monotonically decreasing function and all four prespecified temperature trajectories lie above and to the right of the front because they are suboptimal. +


Figures 9 and 10 show the growth rate and temperature trajectories as well as the final nucleated particle size distribution for individual single-objective optimizations and 14245


Article


Figure 9. Results for individual single-objective optimizations and middle point on the Pareto-optimal front for problem 2, potassium nitrate system: (a) temperature trajectory; (b) growth rate trajectory; (c) NSDn,f; (d) VSDn,f.


for the middle point on the Pareto-optimal front for both systems for problem 2. The linear cooling trajectory results in a relatively large value of the NMS because it is an early growth trajectory and NMS is based on the lower moments at the end of batch (μ0,f and μ1,f). The constant growth rate trajectory has the smallest θ value among the four prespecified trajectories again since maximizing the NMS requires an early growth trajectory while minimizing mu3n requires a late-growth trajectory. Furthermore, trajectories corresponding to the middle points for both systems show relatively small change in growth

rate during the batch. The result again suggests that the constant growth rate trajectory represents a reasonable tradeoff between early growth and late-growth objective functions. Figure 10a shows that the temperature trajectory corresponding to optimizing the NMS is not monotonically decreasing but rather increases abruptly at about 25 min. This is because the optimal growth rate trajectory changes very rapidly at this time. This suggests that it may not be practical or even possible to implement the optimal trajectory for this single objective function. 14246


Article

Industrial & Engineering Chemistry Research min ϕ = α u(τ )

(μ3,s,f + x3(τf )) (μ5,s,f + x1(τf )) (μ4,s,f + x 2(τf ))2

ij μ + x 2(τf ) yz 4,s,f zz + (1 − α)jjjj− z j μ + x3(τf ) zz k 3,s,f {

−1

(40)

s.t. x(0) = 0,

x 7(τf ) = tf,c

u(τ ) ∈ U = [umin , ∞)

∀ τ ∈ [0, τf ]

Expressions for the costate trajectories for problem 3 can be derived from eqs 27, 29, and 40: Figure 11. Pareto-optimal front for problem 3, potassium nitrate system.

ψ1(τ ) = αk1 ψ2(τ ) = α( −5k1(τ − τf ) + k 2) + (1 − α)k4 ψ3(τ ) = α[10k1(τ − τf )2 − 4k 2(τ − τf ) + k 3] + (1 − α)[−4k4(τ − τf ) + k5] ψ4(τ ) = α[−10k1(τ − τf )3 + 6k 2(τ − τf )2 − 3k 3(τ − τf )] + (1 − α)[6k4(τ − τf )2 − 3k5(τ − τf )] ψ5(τ ) = α[5k1(τ − τf )4 − 4k 2(τ − τf )3 + 3k 3(τ − τf )2 ] + (1 − α)[−4k4(τ − τf )3 + 3k5(τ − τf )2 ]

(41)

ψ6(τ ) = α[−k1(τ − τf )5 + k 2(τ − τf )4 − k 3 (τ − τf )3 ] + (1 − α)[k4(τ − τf )4 − k5(τ − τf )3 ]


3.4. Problem 3. Weight Mean Coefficient of Variation (WMCOV) vs Weight Mean Size (WMS). Problem 3 can be stated as

ψ7(τ ) = ν



Article



where the constants k1−k5 are yz 1 jij μ3,f μ5,f z − 1zzz k = jjj 2 zz 2 jj μ4,f k { μ3,f k1 = k 2 μ4,f

−1/2

k 2 = −2k k3 = k

μ4,f 3

μ5,f

(42)

μ4,f 2

k4 = − k5 =

μ3,f μ5,f

Figure 15. Pareto-optimal front for problem 4, potassium nitrate system.

1 μ3,f

μ4,f μ3,f

(see Tables 3 and 4) for the middle point on the Pareto-optimal front. The reason is that both objectives are defined based on the higher moments (from μ3,f to μ5,f). In this case, the middle trajectory is also very similar to the single objective cases. 3.5. Problem 4. Number Mean Coefficient of Variance (NMCOV) vs Zero Moment of Nucleated Nuclei (mu0n). Problem 4 is defined as

2

To prove the existence of τ+, the expression for ψ6(τ) can be rewritten as ψ6(τ ) = α(τf − τ )3 [k1(τ − τf )2 − k 2(τ − τf )1 + k 3] + (1 − α)(τf − τ )3 [−k4(τ − τf ) + k5]

(43)

min ϕ = α u(τ )

Since (τf − τ)3 > 0 for τ ∈ [0, τf) and k3, k5 > 0, it can be deduced that ψ6 > 0 when τ − τf is sufficiently small, which indicates that τ+ exists and the batch time constraint is guaranteed to be active. Figures 11−14 show the Pareto-optimal front, trajectories, and final nucleated particle size distributions for problem 3. For this problem, the trade-off between the two objectives is not significant, which can be deduced from the similarity of growth rate trajectories, final nucleated particle size distributions (see Figures 13 and 14), and relatively small value of θ

(μ0,s,f + x6(τf )) × (μ2,s,f + x4(τf )) (μ1,s,f + x5(τf ))2

+ (1 − α)(x6(τf ))

−1 (44)

s.t. x(0) = 0,

x 7(τf ) = tf,c

u(τ ) ∈ U = [umin , ∞)

∀ τ ∈ [0, τf ]

Based on eqs 27, 29, and 44, the corresponding costate equations are 14248


Article


Since k3 > 0, τ+ can be chosen as τf and the batch time constraint is guaranteed to be active. Figures 15−18 show the Pareto-optimal front, trajectories, and crystal size distribution of nucleated crystals for simultaneous minimization of number mean coefficient of variance (NMCOV) and the zero moment of nuclei (mu0n). For the potassium nitrate system, since both objective functions are early growth ones, the linear trajectory, which is an early growth trajectory, leads to the smallest value of θ. As for the pentaerythritol system, since the trajectory that minimizes mu0n in the pentaerythritol system is the constant growth rate trajectory, the linear trajectory leads to relatively large values of mu0n and θ, as shown in Figure 16. Figure 17 shows that maximum growth rate constraint is active at the beginning of the batch for the optimal trajectories for the single objective functions and the middle trajectory for the potassium nitrate system. The time duration that the maximum growth rate constraint is active for the middle trajectory is between the values for the two single objective function cases, which is also reasonable since the middle trajectory represents the trade-off between two single objectives. 3.6. Multiobjective Optimization Involving Total Batch Time. Batch time is also a crucial objective that should be taken into consideration when optimizing a batch process. In this work, the batch time is represented as the value of a state variable at the final transformed time (see eq 11). To treat the batch time as an objective, the final state constraint is changed and the second objective function is set to be one of those listed in Table 2. Two Pareto-optimal fronts were constructed using this method with mu3n and WMS as the second objective function. The results are shown in Figures 19 and 20. To normalize the value of mu3n for comparison of two crystallization systems, the y-axes in Figures 19a and 20a are mu3n divided by μ3,s,f,c. For the potassium nitrate system, the value of mu3n is less than 10% of μ3,s,f,c when the batch time is longer than 100 min. For the pentaerythritol system, the value of mu3n increases


ψ1(τ ) = 0,

ψ2(τ ) = 0,

ψ4(τ ) = αk1,

ψ3(τ ) = 0,

ψ5(τ ) = − 2αk1(τ − τf ) + k 2 ,

ψ6(τ ) = α(k1(τ − τf )2 − k 2(τ − τf ) + k 3) + (1 − α),

(45)

ψ7(τ ) = ν

where the constants k1−k3 are yz 1 ijj μ0,f μ2,f zz − k = jjj 1 zz zz 2 jj μ1,f 2 k { μ0,f k1 = k 2 μ1,f

−1/2

k 2 = − 2k k3 = k

μ0,f μ2,f

(46)

μ1,f 3

μ2,f μ0,f 2



Article



Figure 19. Pareto-optimal front involving batch time, potassium nitrate system: (a) mu3n versus batch time; (b) WMS versus batch time.

Figure 20. Pareto-optimal front involving batch time, pentaerythritol system: (a) mu3n versus batch time; (b) WMS versus batch time.

observed: The objective function based on the lower moments favored an early growth trajectory and the objective function based on the higher moments favored a late-growth trajectory. Furthermore, the Pareto fronts were quite convex so that there was a segment along the optimal front that seems to represent a good trade-off between the objective functions; i.e., the value of both objective functions was reasonably close to the singleobjective optimal value. Growth rate trajectories for points in these segments were nearly constant, suggesting that a nearly constant trajectory may represent a reasonable trade-off between early growth and late-growth objectives. When both objectives were optimized by the same kind of trajectories (for example, WMCOV and WMS), the trade-off was much less significant and the range of values spanned by the Pareto front was much smaller, indicating, for example, that

more rapidly as the batch time decreases and even exceeds the seed-grown mass for batch times less than about 100 min. Furthermore, the WMS of the pentaerythritol system is also more sensitive to batch time. The result suggests that the sensitivity of objective functions based on moments to batch time is strongly case-dependent.

4. CONCLUSIONS Pareto-optimal fronts for several combinations of objectives widely considered in seeded batch crystallization were constructed by application of Pontryagin’s minimum principle. For two-objective cases when one objective function was optimized by early growth trajectory and the other one was optimized by late-growth trajectory, a clear trade-off was 14250


Article


(14) Jones, A. G. Optimal operation of a batch cooling crystallizer. Chem. Eng. Sci. 1974, 29, 1075−1087. (15) Ajinkya, M. B.; Ray, W. H. On the optimal operation of crystallization processes. Chem. Eng. Commun. 1974, 1, 181−186. (16) Morari, M. Some comments on the optimal operation of batch crystallizers. Chem. Eng. Commun. 1980, 4, 167−171. (17) Corriou, J. P.; Rohani, S. A new look at optimal control of a batch crystallizer. AIChE J. 2008, 54, 3188−3206. (18) Hofmann, S., Raisch, J. Application of optimal control theory to a batch crystallizer using orbital flatness. In 16th Nordic Process Control Workshop, Lund, Sweden; Lund Institute of Technology: 2010; pp 25−27. (19) Bajcinca, N. Analytic solutions to optimal control problems in crystal growth processes. J. Process Control 2013, 23, 224−241. (20) Bajcinca, N.; Hofmann, S.; Raisch, J.; Sundmacher, K. Robust and optimal control scenarios for batch crystallization processes. In Proceedings of ESCAPE-20, Ischia, Naples, Italy; Elsevier: 2010; pp 1605−1610. (21) Bajcinca, N.; Hofmann, S. Optimal control for batch crystallization with size-dependent growth kinetics. In Proceedings of the 2011 American Control Conference, San Francisco, CA; IEEE: 2011; pp 2558−2565. (22) Tseng, Y. T.; Ward, J. D. Comparison of objective functions for batch crystallization using a simple process model and Pontryagin’s minimum principle. Comput. Chem. Eng. 2017, 99, 271−279. (23) Hsu, C. W.; Ward, J. D. The Best Objective Function for Seeded Batch Crystallization. AIChE J. 2013, 59 (2), 390−398. (24) Nagy, Z. K.; Braatz, R. D. Advances and New Directions in Crystallization Control. Annu. Rev. Chem. Biomol. Eng. 2012, 3 (3), 55−75. (25) Nagy, Z. K.; Fevotte, G.; Kramer, H.; Simon, L. L. Recent advances in the monitoring, modelling and control of crystallization systems. Chem. Eng. Res. Des. 2013, 91, 1903−1922. (26) Sarkar, D.; Rohani, S.; Jutan, A. Multi-objective optimization of seeded batch crystallization processes. Chem. Eng. Sci. 2006, 61 (16), 5282−5295. (27) Trifkovic, M.; Sheikhzadeh, M.; Rohani, S. Kinetics estimation and single and multi-objective optimization of a seeded, anti-solvent, isothermal batch crystallizer. Ind. Eng. Chem. Res. 2008, 47 (5), 1586− 1595. (28) Acevedo, D.; Tandy, Y.; Nagy, Z. K. Multiobjective Optimization of an Unseeded Batch Cooling Crystallizer for Shape and Size Manipulation. Ind. Eng. Chem. Res. 2015, 54 (7), 2156− 2166. (29) Lewis, F. L.; Vrabie, D. L.; Syrmos, V. L. Optimal Control, 3rd ed.; Wiley: Hoboken, NJ, 2012. (30) Mullin, J. W.; Nývlt, J. Programmed cooling of batch crystallizers. Chem. Eng. Sci. 1971, 26 (3), 369−377. (31) Miller, S. M.; Rawlings, J. B. Model identification and control strategies for batch cooling crystallizers. AIChE J. 1994, 40, 1312− 1327. (32) Bernardo, A.; Giulietti, M. Modeling of crystal growth and nucleation rates for pentaerythritol batch crystallization. Chem. Eng. Res. Des. 2010, 88 (10), 1356−1364.

the sacrifice in the value of the WMS required to minimize the WMCOV was relatively small. If minimizing the batch time is also included as an objective, analysis shows that the Pareto front for a two-objective problem including minimizing the batch time as one objective can be constructed by solving the single optimization problem for different final values of a state variable. Result shows that if the other objective function is WMS or mu3n, its value always improves with increasing batch time, which is in accordance with the conclusion that the batch time inequality constraint is always active. However, the extent of improvement depends strongly on the crystallization system. These results are based on very simple crystallization models that account only for nucleation and ordinary growth. Nevertheless, they capture the essential effect of the objective function on the optimization problem and are useful for understanding the inherent trade-offs that may be encountered when considering multiobjective optimization of batch crystallization processes.

■

AUTHOR INFORMATION

Corresponding Author

*Fax: +886-2-2369-1314. E-mail: jeff[email protected]. ORCID

Jeffrey D. Ward: 0000-0003-0727-7689 Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Kraft, D. On Converting Optimal Control Problems into Nonlinear Programming Problems. In Computational Mathematical Programming; Schittkowski, K., Ed.; Springer: Berlin, 1985; pp 261− 280. (2) Schlegel, M.; Stockmann, K.; Binder, T.; Marquardt, W. Dynamic optimization using adaptive control vector parameterization. Comput. Chem. Eng. 2005, 29, 1731−1751. (3) Worlitschek, J.; Mazzotti, M. Model-Based Optimization of Particle Size Distribution in Batch-Cooling Crystallization of Paracetamol. Cryst. Growth Des. 2004, 4 (5), 891−903. (4) Paengjuntuek, W.; Kittisupakorn, P.; Arpornwichanop, A. Batchto-batch Optimization of Batch Crystallization Processes. Chin. J. Chem. Eng. 2008, 16 (1), 26−29. (5) Chung, S. H.; Ma, D. L.; Braatz, R. D. Optimal seeding in batch crystallization. Can. J. Chem. Eng. 1999, 77 (3), 590−596. (6) Ma, D. L.; Tafti, D. K.; Braatz, R. D. Optimal control and simulation of multidimensional crystallization processes. Comput. Chem. Eng. 2002, 26 (7−8), 1103−1116. (7) Ma, D. L.; Braatz, R. D. Robust identification and control of batch processes. Comput. Chem. Eng. 2003, 27 (8−9), 1175−1184. (8) Nagy, Z. K.; Fujiwara, M.; Braatz, R. D. Modelling and control of combined cooling and antisolvent crystallization processes. J. Process Control 2008, 18 (9), 856−864. (9) Nagy, Z. K.; Braatz, R. D. Robust nonlinear model predictive control of batch processes. AIChE J. 2003, 49 (7), 1776−1786. (10) Lang, Y.-d.; Cervantes, A. M.; Biegler, L. T. Dynamic Optimization of a Batch Cooling Crystallization Process. Ind. Eng. Chem. Res. 1999, 38 (4), 1469−1477. (11) Hu, Q.; Rohani, S.; Jutan, A. Modelling and optimization of seeded batch crystallizers. Comput. Chem. Eng. 2005, 29 (4), 911− 918. (12) Vollmer, U.; Raisch, J. Control of batch cooling crystallization processes based on orbital flatness. Int. J. Control 2003, 76, 1635− 1643. (13) Vollmer, U.; Raisch, J. Control of batch crystallization - A system inversion approach. Chem. Eng. Process. 2006, 45, 874−885. 14251


Pareto-Optimal Fronts for Simple Crystallization Systems Using

Recommend Documents