Multiobjective Optimization of an Unseeded Batch Cooling Crystallizer

Article pubs.acs.org/IECR

Multiobjective Optimization of an Unseeded Batch Cooling Crystallizer for Shape and Size Manipulation David Acevedo, Yanssen Tandy, and Zoltan K. Nagy* School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907, United States ABSTRACT: The optimization of unseeded batch cooling crystallization systems is studied in a multiobjective framework. The length mean size and target aspect ratio (AR) of final crystals were considered as objectives, which give rise to a set of optimal solutions known as Pareto-optimal solutions. Two model compounds were considered: paracetamol and potassium dyhidrogen phosphate (KDP), nucleation and growth dominated systems, respectively. The optimization of KDP showed a dependence between the optimal profile and the target objectives weights. The optimal profile for paracetamol did not vary from a cubic cooling profile due to the low sensitivity of the AR to temperature changes. Furthermore, the longer batch time has a greater impact on the optimal AR that could be achieved due to the variation in supersaturation dependence. Moreover, it was shown that there is a relation between the classification of crystallization kinetics based on growth and nucleation mechanisms and the set of optimal cooling profiles for shape and size optimization.

1. INTRODUCTION Batch crystallization is an important unit operation in the pharmaceutical industry because of the high purity that can be obtained through the liquid−solid separation. Over 90% of the active pharmaceutical ingredients (API) are crystals of small organic molecules.1 The product quality and downstream process operations can be strongly influenced by the crystallization performance. Variations in size distribution, purity, and morphology of the crystals can affect downstream processes such as filtration, milling, and eventually, formulation. Extensive work has been done in the development of control approaches using online characterization techniques.2−4 Various control techniques have been developed in which the supersaturation and temperature profiles were optimized and controlled for a seeded batch crystallization process.5,6 Monitoring, control, and optimization of crystallization processes to achieve a uniform size distribution has been the main focus due to the technology limitations to monitor other properties such as crystal shape. During the past decade considerable progress has been achieved in the development of fundamental models to describe the shape evolution of crystals.7,8 A feedback control strategy was implemented to a cooling crystallization process to achieve a desired shape through the manipulation of additives.7 Furthermore, Yang et al. showed that the aspect ratio (AR) of crystals can be controlled by manipulating the supersaturation profile by variations in batch time and seed levels.9 Therefore, to ensure the efficiency of the process it is necessary to find the optimal operation in order to achieve the desired shape and size distribution. However, crystallization processes are complex in nature and are known to be highly nonlinear which impose greater difficulties in their modeling and control.8 Optimization of a crystallization process is used as a tool to achieve products with desired properties in a more efficient manner. Traditionally, the optimization of a batch cooling crystallizer is performed with respect to the cooling profile and seeding characteristics, which keep supersaturation at an optimum level through the process. The crystal size distribution © XXXX American Chemical Society

(CSD) has been commonly the objective considered for various types of crystallization processes due to its strong influence on the quality of the product.10−13 However, recent works have focused on the optimization of crystallization processes considering shape based objective functions. The optimal temperature profile was obtained for the cooling crystallization of L-glutamic acid (L-GA) using a morphological population balance model (PBM).14 It was shown that the supersaturation level and duration of the generation of short, medium and long particles are significantly different. An optimization algorithm was applied to a protein crystallization process, whereby the growth of individual faces was optimized with the aim of obtaining desired crystal size and shape distributions.15 In the context of multiobjective optimization, a trade-off between the performance objectives exists, in which a desirable change in one objective function results in an unfavorable change in another objective function. Therefore, the presence of multiple objectives usually gives rise to a family of nondominated solutions known as Pareto-optimal solutions because there may not exist one solution that is better with respect to other objectives.16 Extensive amounts of work have been performed for various crystallization systems in which the applicability of different algorithms has been shown considering objective functions related to size.16−19 A multiobjective optimization framework was applied to a multisegmented, multiaddition plug-flow crystallizer (MSMA-PFC) to determine the antisolvent flow rates into each segment that simultaneously maximize the average crystal size, and minimize the coefficient of variation.20 Our previous work developed a classification method by simultaneous consideration of the effect of temperature profiles on nucleation and growth rates of two different characteristic crystal dimensions.21 The various systems were classified with Received: September 25, 2014 Accepted: February 9, 2015

A

DOI: 10.1021/acs.iecr.5b00173 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research respect to the magnitudes of the nucleation rate and growth rates differences between the two characteristic dimensions, into high nucleation rate − high growth rate difference, high nucleation rate − low growth rate difference, low nucleation rate − high growth rate difference, and low nucleation rate − low growth rate difference systems. It was shown that generally in crystallization, an equilibrium shape is achieved, hence the variation in the aspect ratio is limited by only manipulating the supersaturation through cooling. The objective of this article is to evaluate the variation in size and shape for various types of systems from these classes through the development of Paretofronts.

crystals. Furthermore, a constant volumetric shape factor, kv, is assumed. The mass balance equation is given by25

2. MULTIDIMENSIONAL PB MODEL AND OPTIMIZATION PROBLEM 2.1. Model Formulation. The important characteristics of particulates products such as size and shape obtained through different types of processes can be represented through a population balance model (PBM).22,23 For a well-mixed batch crystallization processes in which it is assume that agglomeration and breakage of crystals are negligible, the population balance (PB) equation is given by

where GL and GW are the size independent growth rates along the length and width of the crystals, B being the primary nucleation rate, and S being the relative supersaturation. The kinetic parameters kn, kb, gn, and b are often sensitive to process conditions. Supersaturation is the driving force for growth and nucleation, which is given by the solute concentration, C, and the equilibrium concentration, Cs. The temperature dependence of the equilibrium concentration can be expressed using a polynomial expression fit to experimental data, given by

∂ n(t , X) + ∇X ·[Gn(t , X)] = Bδ(X − X 0) ∂t

(1)

I.C.: n(0, X) = n0(X)

(2)

∞

∫0 ∫0

∞

n(t , X)LiW jdLdW

dt dμij dt

=B

(7)

n = L, W

(8)

C − Cs(T ) Cs(T )

(9)

Cs = aT 4 + bT 3 + cT 2 + dT + e

(10)

2.2. Growth Kinetics Estimation Problem. The twodimensional growth kinetics are necessary in order to describe the shape evolution of a crystal. The estimation of the twodimensional growth kinetics parameters have shown to be a challenging problem due to the limitation of online monitoring techniques. The recent development of online methods to monitor the evolution of the shape of crystals such as the Process Vision and Measurement (PVM, Mettler-Toledo) system have increased the possibilities to control shape. The PVM has been utilized on a wide range of applications in crystallization studies such as for the determination of growth kinetics, polymorphic form, and characterization of metastable zone width.26,27 Through the analysis of the online images obtained it is possible to monitor the growth of different characteristic lengths. Therefore, it is possible to estimate the dynamic variation of the AR throughout the crystallization process, which coupled with the solute concentration can be used as model variables in order to obtain the growth kinetic parameters. The determination of the growth kinetics can be formulated as a nonlinear optimization problem.28,29 Let θ be the vector of kinetic parameters from eq 7 θ T = [kLk WgL g W ]

(11)

Then the estimation problem is given by l

(3)

2

min ∑ ∑ (yij − yij*)2 θ

j=1 i=1

(12)

where yij and yij* are the model predicted and measured variables for the variable i and sampling interval j. The variables considered in the optimization problem were AR and solute concentration. The growth kinetic parameters were obtained for various initial conditions that showed similar results. The confidence intervals were estimated because the parameter estimates are stochastic variables with corresponding probability distributions.30 The intervals were obtained through the projection of uncertainty via linear error propagation.31 2.3. Multiobjective Optimization Problem. The final crystal properties in an unseeded, well-mixed batch crystal-

(4)

= iG Lμ(i − 1)j + jG W μi(j − 1) + BL0i W0j ,

i , j = 1, 2..., n

Gn = knS gn

S=

The 2D moment equations are given by

dμ00

(6)

B = k bS b

where n(t,X) is the density distribution, X = [X1, X2, ..., XN] is the vector containing the various characteristic lengths, G = [G1, G2, ..., GN] is the vector of growth rates, B is the nucleation rate, X0 is the size of the nuclei, δ is the Dirac δ-function acting at X = X0, and n0(X) is the initial (seed) distribution. Previous work demonstrated the use of a 1-D PB model to describe the size and shape of crystals in a seeded batch crystallization of an organic, industrial photochemical; however, the habit dynamics only applied to the seed crystals due to the assumptions of the model.24 This work focuses on size and shape of crystals in an unseeded batch crystallization process thus a 2-D PB model is implemented. Furthermore, a classification framework introduced in previous work is studied and it involves two characteristic growth dimensions which can be considered as the length and width of the crystals, L and W, respectively.21 The method of moments (MOM) is a common technique used to solve the PBM by transforming it into a set of ordinary differential equations (ODEs), which are solved numerically. The moments corresponding to a 2D CSD are known as the cross-moment equation, which are given by25 μij =

dc = −2ρc k vG L(μ11 − μ20 ) − ρc k vG W μ20 dt

(5)

A solute mass balance equation in terms of moments completes the multidimensional PBM for an unseeded batch crystallizer. The amount of solute leaving the solution depends on crystal growth with the assumption of negligible size of nucleation B


Article

Industrial & Engineering Chemistry Research

various optimization steps in order to evaluate and construct the Pareto-optimal front.

lization are determined by the supersaturation trajectory, which is considered to be dependent on temperature, T(t); the supersaturation is generated by reducing the temperature throughout the crystallization process. Therefore, the crystallizer temperature can be used to derive an optimal temperature profile, which can lead to the desired shape and size distribution. The multidimensional PB model was applied to derive the optimal temperature profile for a defined objective function, which is subjected to several constraints. The temperature range and ramp are considered as constraints to ensure that the temperature profiles stay within the operating conditions. Furthermore, the product yield and batch time are also constraint related to the productivity and efficiency of the crystallization process. Therefore, the optimization problem formulated with equality and nonequality constraints is given by

3. MATERIALS AND METHODS 3.1. Materials. Acetaminophen (paracetamol, SigmaAldrich, MO), potassium dihydrogen phosphate (KH2PO4, KDP, Macron Fine Chemicals, PA), and acetylsalicylic acid (ASA or aspirin, Sigma-Aldrich, MO), were used as model systems. Potash alum, (KAl(SO4)2·12H2O,), potassium nitrate (KNO3), and L-glutamic acid (L-GA) were analyzed for the simulation case studies. The parameters for eq 10 to estimate solubility information of the various systems analyzed through this study were obtained from the literature.14,21,32,33 3.2. Experimental Setup. A 500 mL lab scale glass jacketed vessel was employed for all experiments conducted in this work. The agitations speed was set to be 350 rpm using an overhead stirrer with a three blade retreat curve impeller. The temperature was controlled with a PT100 thermocouple using a Julabo F25-ME refrigerated and heating circulator controlled via the Crystallization Process Informatics Systems (CryPRINS) interface. A G400 FBRM Lasentec probe (MettlerToledo Autochem Ltd.) was positioned in the crystallizer for all experiments performed in this work in order to monitor the chord length distribution (CLD). The solution concentration was monitored during the growth kinetics estimation experiments using a Carl Zeiss MCS621 UV/vis spectrometer. The FBRM, UV/vis and temperature data were recorded every 10 s. A V819 PVM Lasentec probe (Mettler-Toledo Autochem Ltd.) was used for the growth kinetics estimation experiments in order to monitor the shape evolution of the crystals studied. Pictures were obtained every second during the various growth kinetics estimation experiments performed. 3.3. Methods. The batch cooling experiments were performed by preparing saturated solutions for the various systems studied by adding 300 mL of ethanol to the amount of crystals needed to reach the equilibrium concentration. The equilibrium concentration was determined from the solubility information for each system at a fixed saturation temperature of 40 °C. The experiments performed can be divided depending on their aim: growth kinetic parameter estimation, and optimization experiments. 3.3.1. Parameter Estimation Experimental Method. The growth kinetic parameter estimation experiments were performed for paracetamol and ASA. A systematic experimental procedure was designed that uses a carefully designed set of experiments that trigger different crystallization mechanisms, such as growth and primary nucleation, allowing the sequential estimation of the various kinetic parameters. This is an efficient procedure that enables the robust estimation of highly correlated kinetic parameters with minimum number of experiments. In this work, primary nucleation and growth were only considered as crystallization mechanisms. First, seeds of each compound were prepared in order to obtain undamaged and nonagglomerated crystals and suppress nucleation during the growth kinetic estimation experiments. The solution was equilibrated at 50 °C and held constant for 60 min at 350 rpm to ensure complete dissolution of crystals. The system was linearly cooled to the initial temperature and held constant for other 60 min. The seed preparation and growth kinetic estimation experiments were performed following a cubic cooling from an initial temperature of 40 to 20 °C. The cubic cooling profile implemented is given by

min J(T ) T (t )

subject to: c1(t ) = Tmin ≤ T(t ) ≤ Tmax. c 2(t ) =

dT ≤0 dt

c3(t ) = C(tfinal) − 0.4C0 ≤ 0 c4(t ) = tfinal = 90min

(13)

where Tmax. is the initial temperature, Tmin is the final temperature, C0 is the initial concentration, and tfinal is the batch time. Temperature is the manipulated variable and it was parametrized as a piecewise linear function with respect to time. The time was discretized with respect to the total batch time as given by

t = 0:

tfinal − t 0 : tfinal N

(14)

where N is a parameter that determines the number of discretization. A sequential procedure was applied in which the number of discretization was increased consecutively. The optimal temperature profile solution was used as initial condition for the next optimization in which the number of discretization is increased to find a smoother optimal solution. The objective function, J, considered in this optimization problem is a vector consisting of two objective functions that are the size and shape descriptors. In this work, the average crystal size for only one of the characteristics dimensions and AR were considered as the objective functions. The length number mean crystal size is given by μ Ln = 10 μ00 (15) whereas the AR is determined as follows μ AR = 10 μ01

(16)

The optimization problem formulated in eq 13 is a constrained nonlinear programming problem (NLP) due to the nonlinear PBM, and algebraic equality and nonequality constraints. Therefore, the successive quadratic programming (SQP) method implemented in the ‘fgoalattain’ function in the MATLAB optimization toolbox was utilized to calculate the operating profiles. The weighting vector was varied throughout C


Article

Industrial & Engineering Chemistry Research T (t ) =

(Tmin − Tmax.) tfinal 3

× t 3 + Tmax.

(17)

The seeds were added to the saturated solution 30 min after the system reached the initial temperature. The seed loading used was 12% (w/w) of the initial mass utilized to prepare the saturated solution. It was observed through various experiments performed at various seed loadings (5−12%) that this amount of seeds was enough to ensure negligible secondary nucleation. The wet seeds and crystals were filtered and washed with acetone to inhibit agglomeration and growth of crystals. Samples were taken after the seeds and crystals were dried in order to obtain offline image data for paracetamol and aspirin. ImageJ was utilized to obtain the dynamic variation of the AR from the PVM images. The AR and mean size of the seeds crystals were obtained using Adobe Photoshop CS 5.1 software; a significant amount of seed crystal images (>200) were analyzed in order to achieve a representative sample of seed crystals for the analysis. The images were analyzed using the customized ruler tool, based on the pixel scale from the microscope. The length and width of the crystals were measured and the AR is determined by dividing the length by the width mean size. The moments of the distribution were obtained using the mean size values from the offline images. 3.3.2. Experimental Method for Validation of the Optimization Results. To validate the multiobjective optimization method, separate unseeded batch cooling experiments involving KDP were conducted. The crystals were initially dissolved at 50 °C and cooled to the desired initial temperature of 40 °C. The temperature profile followed for each set of experiments was obtained from the optimization results, while maintaining the initial and final temperature at 40 and 20 °C, respectively. The crystals were dried and offline image analysis was used to determine the AR following the similar image analysis method explained in section 3.3.1.The results obtained from these experiment were then compared to the calculated data obtained using the multiobjective optimization.

4. RESULTS AND DISCUSSION 4.1. Evaluation of the Crystallization Classification Framework for Active Pharmaceutical Ingredients. The two-dimensional growth and nucleation kinetic parameters were estimated for paracetamol and aspirin in order to evaluate the application of the classification framework in the case of these model pharmaceutical compounds based on their crystallization kinetics. A classification scheme was developed in previous work considering mainly inorganic compounds as model systems.21 The aim of this initial work was the validation of the applicability of the systematic classification framework to active pharmaceutical ingredients. The crystallization kinetics for both systems was determined using online data obtained from PVM images. Figure 1 shows the PVM images obtained during the seeded batch cooling crystallization of paracetamol used to determine the growth kinetics for two crystal characteristic dimensions. The shape of the seeds shown in Figure 1a changed between the first 60 min of experimentation since it reached a prism-like habit (Figure 1b,c). The aspect ratio (AR) through the cooling crystallization was determined using the images obtained online as a measure of the growth variation of the two characteristic crystal dimensions (length, and width). The moments, AR, and mean size of the initial seeds were determined in order to have a better initialization for the least-squares optimization problem. The optimization

Figure 1. PVM images obtained for the determination of paracetamol two-dimensional growth kinetics at time: (a) 0, (b) 60, and (c) 120 min.

problem was solved utilizing the SQP method named “fmincon”, which calls the simulation model iteratively in MATLAB. Therefore, the MOM was used since it reduces the simulation time considerably compared to other methods. The measured solution concentration and AR for a series of time steps were used for the ordinary least-squares calculation throughout the optimization. The same method was utilized to determine the two-dimensional growth kinetics parameters of aspirin. The resulting growth kinetics parameters with their confidence intervals are given in Table 1. The estimated two-dimensional growth kinetics were used for the determination of the classification of the two systems with respect to the size independent growth differences and primary nucleation sensitivity for a set of 30 different temperature profiles. Previous work done using KDP, LGA, potash alum, and KNO3 as model systems demonstrated the effect and interaction of the kinetics in the achievable shape and size through a proposed classification framework;21 this provides the tools to investigate the practical controllability of the attainable shape and size in a unseeded crystallization D


Article


Table 1. Two-Dimensional Growth Kinetic Parameters with 95% Confidence Intervals Determined for Paracetamol (PCM) and Aspirin (ASA) from Cooling Crystallization Experiments from 40 to 20 °Ca

a

systems

kL × 10−2

kW × 10−2

gL

gW

C0 (g/g)

kv

ρc (g/cm3)

PCM ASA

1.80 ± 0.08 6.91 ± 0.05

2.06 ± 0.09 0.87 ± 0.04

1.19 ± 0.14 1.70 ± 0.51

1.00 ± 0.19 1.36 ± 0.45

0.27 0.39

0.24 0.67

1.29 1.40

kL and kW are in cm/min. Parameters utilized for simulation and optimization of paracetamol and aspirin.

Figure 2. (a) Growth rates difference and primary nucleation following a natural, cubic, and linear temperature profiles and (b) length number mean size and final aspect ratio variation obtained following 30 different temperature profiles for aspirin (ASA), paracetamol (PCM), KNO3, L-GA, and KDP. For the systems, PCM and ASA kinetic parameter sensitivity analysis results are also included to illustrate the variations of the properties due to uncertainties in the kinetic parameters.

compared to compounds in the same region (e.g., potash). Therefore, it is clear that changing the temperature trajectories will not have a significant effect on the achievable AR for a highly nucleation dominated systems, which shows similar supersaturation dependence for different facet growth rates.21 A significant variation in the growth rate difference was obtained for ASA compared to paracetamol, as shown in Figure 2a. This trend is expected due to the difference in magnitude between the growth kinetic parameters for aspirin compared to paracetamol. The system can be classified as a system with low nucleation rate and high growth rate difference. A significant variation in the AR and larger mean size should be expected for ASA, similarly as observed for systems that have similar behavior, such as KDP and L-GA. Figure 2b demonstrates that ASA does show a higher sensitivity in the AR by temperature variations. Furthermore, different ranges of variation of the growth differences was observed for the various growth kinetic combinations due to the larger range in the confidence intervals of the growth orders. However, it can be seen that the AR variability is significant since the nucleation rate dependence on supersaturation is low; the supersaturation will be generated and consumed mostly due to the growth of the two characteristic dimensions. It was shown experimentally in previous studies that in the case of systems that belong to this class, supersaturation profile during the cooling crystallization processes can have a significant impact on the final shape of crystals.9,21 The supersaturation dependence of each system has a direct impact on the region in which the system can be classified due to the strong interplay between the growth and nucleation mechanisms and their competitive dependence on supersaturation. This classification scheme can be used as an initial step to develop control approaches to tailor the unseeded cooling crystallization process to achieve a desirable size and shape. The results also demonstrate the applicability of the classification scheme to analyze the attainable size and shape of active pharmaceutical ingredients, such as paracetamol and aspirin. Each class of system shows a different sensitivity of the AR and mean size to the temperature changes. Therefore, an optimization problem can be formulated in order to evaluate

system with respect to variations in temperature profiles. A set of 30 temperature profiles was considered for the base case study in order to simulate a wide operational range. Profiles between natural and programmed (estimated with a cubic profile) were simulated by formulating a piecewise linear function. The parameters utilized in this simulation study are shown in Table 1. A set of 10 different growth kinetic parameter combinations were simulated for paracetamol and aspirin in order to evaluate the sensitivity of the model and classification framework. The aim of the uncertainty analysis was to investigate the variations of attainable size and shape regions in realistic industrial conditions when the kinetics at the industrial setup may change compared to the laboratory measurements used to classify the compound, due to effects caused by scale up, mixing or impurities in the solvent or raw material. The parameters were chosen within the confidence intervals. Figure 2a shows the nucleation and growth rate difference variation for paracetamol, aspirin, with respect to the various systems studied in the previous work. The trend obtained for paracetamol shows a significantly higher sensitivity to the nucleation rate and lower variation in growth rate differences via temperature variations compared to the other systems. The low sensitivity of the growth rate difference demonstrates that the temperature profile has a low effect on the achievable AR. Paracetamol shows a lower variability in the growth rate difference and nucleation range compared to the other systems, hence it can be classified as a compound with high nucleation rate and low growth rate difference (see Figure 2a). The final crystals obtained for this system should not achieve a large mean size and be of low AR. Figure 2b shows the variation of the final AR and mean size obtained for 30 different cooling trajectories. The final AR obtained did not show significant variation due to changes in temperature profiles with respect to the same growth kinetic combination; the variations in the kinetics only affects the magnitude of the AR, which can be related to the high supersaturation effect on the nucleation rate Furthermore, it can be seen that large mean size was observed for paracetamol since the nucleation rate range obtained is much lower E


Article


final crystals. Region II can be related to a transition zone between regions I and II.16 A distinct temperature profile can be associated with each point of the Pareto-optimal front. Three cooling profiles representing different regions of the Pareto-optimal front in Figure 3 are shown in Figure 4a. All

the optimal temperature profile that can be obtained for different system classes. 4.2. Multiobjective Optimization. For the multiobjective optimization approach, paracetamol and KDP crystallization were chosen from the classification scheme (Figure 2) because the two systems are classified in two different regions, which showed strongly different interplay between the crystallization kinetics (refer to section 4.1). Furthermore, KDP and paracetamol were chosen in order to simulate inorganic and organic molecules. The crystallization kinetics for KDP shown in Table 2 was obtained from the literature.21,25 Choosing these Table 2. Parameters Used in KDP Optimization Studies parameters

values

kL gL kw gw kb b

0.60 1.74 0.07 1.48 4.49 × 106 2.04

units cm/min cm/min #/cm3 min

parameters

values

units

kv ρc Tmin Tmax. tbatch tstep

0.67 2.34 40 20 90 4.5

g/cm3 °C °C min min

two systems also enables to provide a comparison of the multiobjective optimization results for an inorganic and an organic compound. The temperature effect on the mean size and aspect ratio are significantly different for the two systems. Therefore, the problem of achieving a target shape, while growing the nucleated crystals can be formulated and solved as a multiobjective optimization problem. 4.2.1. KDP Optimization Case Study. The objective of the optimization problem was to maximize the final mean size that can be achieved while maintaining a target shape (AR = 1). The optimization problem was formulated initially for KDP since it shows the highest sensitivity in the AR due to temperature profile changes. The parameters used in the optimization study of KDP are shown in Table 2. Figure 3 shows the Pareto-

Figure 4. (a) Temperature, and (b) relative supersaturation profiles obtained from optimization of an unseeded batch cooling crystallization at various objective function weights (regions): I (solid line), II (dash line), and III (dots).

three temperature profiles can be characterized by progressively increasing cooling rates. However, the temperature trajectory obtained for region I is closer to linear temperature profile compared to the profiles obtained for the other regions. Various works have shown that a linear temperature profile is obtained when the AR is minimized.15,34 Moreover, it was found that a linear cooling profile caused crystals change from needle-like to plate-like.35 Region I also corresponds to the solution of the single objective optimization problem where the final AR is minimized. The supersaturation profile obtained from the solution of this optimization problem (region I) is shown in Figure 4b. A significant increase in the supersaturation level occurs during the initial part of the process, followed by a steady decrease in the supersaturation throughout the rest of the batch. For unseeded crystallization, it is typical that the optimal supersaturation profile shows an initial peak, which is required for the generation of initial “seed” nuclei after which rapidly drops to lower levels to avoid nucleation and maximize growth of the this is the objective.36 A similar trend was obtained by Yang et.al. in which the effects of different supersaturation profiles, batch times, and seed loadings were studied. It was observed that fast increase in the supersaturation level during the initial part of the cooling batch crystallization process resulted in crystals with the lowest AR.9 The temperature profile selected from region III, as shown in Figure 4a, describes the solution that satisfies the maximization of the number-based mean size. The cooling profile follows a

Figure 3. Pareto optimal-front obtained for the multiobjective optimization problem of maximizing the number mean length and minimizing the aspect ratio to a target shape for an unseeded batch crystallization process of KDP.

optimal front obtained for the multiobjective optimization problem using KDP as model system. The trade-off between the two objectives considered in the optimization problem can be clearly observed from the figure. The final shape that can be obtained can be improved (AR = 1) only at the expense of decreasing the mean size of the final crystals. The Paretooptimal front can be divided in three distinct regions that show different sensitivity to the manipulated variable. Region I shows low sensitivity in the variation of the AR when the mean size is increases. Region III is characterized by a rapid increase in the AR for any attempt to increase the number mean size of the F


Article


Figure 5. Dynamic evolution of aspect ratio and number base mean length obtained from the optimization of an unseeded batch cooling crystallization following the temperature profiles obtained from the three regions of the Pareto-optimal front: I (solid line), II (dash line), and III (dots).

level throughout the batch can be considered almost constant after primary nucleation occurs. Therefore, the number base length mean size shows a uniform increase throughout the entire batch. These observations demonstrate the direct impact of the supersaturation dynamics in the evolution of the AR and number based mean length during the cooling batch crystallization of a highly nucleating system such as KDP. The optimal temperature profiles obtained for the various points of the three regions considered in Figure 4a were implemented in various unseeded batch crystallization experiments, as explained in section 3.3. A tetragonal prism combined with a tetragonal bipyramid shape was obtained for the final crystals of KDP for each scenario, as shown in Figure 6. The increase in the length characteristic dimension is evident in the case of the final KDP crystals obtained following the optimal temperature profile from region III compared to the one from region I (Figure 6a,c). The mean AR and mean length obtained for each scenario using image analysis is shown in Table 3. The lowest AR was achieved for the set of experiments in which the optimal temperature profile implemented was obtained from region I in the Pareto-optimal front shown in Figure 3. The largest crystals were obtained for the cooling crystallization experiments using a temperature profile from region III. Smaller crystals were obtained using region I scenario and vice versa. Table 3 also shows the expected trade-off between the two properties and the effect of choosing a temperature profile from a certain operating region. As expected, for the experiment using a temperature profile from region II, the AR and mean size obtained were between the values of corresponding to the experiments with temperature profiles from regions I and III. 4.2.2. Effect of Batch Time on KDP Pareto Front. Batch time can affect not only the final properties of crystals but also the optimal trajectory of the process; this was demonstrated in previous work for seeded cooling batch systems considering growth rate dispersion only in which the batch time was formulated as objective function and size as constraint.37 However, in our work the effect of batch time in the final properties, AR and length number mean size was studied through the evolution of the Pareto-optimal front for two different batch times (90 and 180 min). The optimal temperature profiles for the unseeded batch cooling crystallization of KDP were determined; the optimal profiles obtained for the 180 min batch scenario follow a similar shape as the one obtained for 90 min (refer to Figure 4(a)). The Pareto-optimal front obtained for the two scenarios are shown in Figure 7. The magnitude of the achievable AR decreased considerably

cubic cooling trend, in which an initial slow temperature decrease is observed followed by a region of very slow decrease. The initial temperature decrease is to generate the necessary supersaturation level for primary nucleation to occur as it can be seen from Figure 4b. Subsequently, this is followed by a region of constant supersaturation in which the small amount of nucleated crystals will growth constantly. The final linear decrease in the temperature after 70 min occurs due to the final temperature constraint. The process is maintained at a low supersaturation level in order to avoid significant nucleation and maximize the growth of crystals. Previous work demonstrated a similar cooling profile tendency for unseeded cooling crystallization processes, in which the objective was to maximize the average crystal size.11 The fast cooling during the final part of the cooling crystallization process can be associated with the final temperature constraint considered in the optimization problem, which forces the temperature profile to decreases to the desired final value within the fixed batch time. However, the solution of the optimization problem considered by region III yields to larger mean size at the expense of higher AR. The temperature profile selected from region II yields a compromised value of AR and number base mean size since this region is characterized by a transition zone between regions I and III. Therefore, the optimal supersaturation trajectory obtained for region II shows low peaks compared to the other profiles throughout the crystallization. It can be seen that the time at which the supersaturation peak occurs changes with different objective functions.16 The supersaturation peak has an effect in the dynamic behavior of the mean size and AR. Figure 5 shows the dynamic evolution of the AR and number base length mean size obtained for the conditions considered for each region. The dynamic evolution of the AR for the three scenarios shows a sudden increase in the value that occurs approximately during the same time as the supersaturation peaks shown in Figure 4b. For region I, nucleation occurs initially due to the sudden increase in the supersaturation level and it continuously decreases throughout the batch. The crystals formed show high initial AR values which continuously decreases subsequently throughout the rest of the batch. A similar tendency can be observed for the dynamic evolution of the number mean length shown in Figure 5b. The sudden increase in the supersaturation level for the scenario in which the number base length mean size is only maximized (region III) leads to a sudden increase in both AR and mean length, due to the fast cooling. The AR trend observed for region II shows the expected relatively stable behavior because the supersaturation G


Article


Figure 7. Pareto-optimal fronts obtained for the multiobjective optimization problem of minimizing aspect ratio and maximizing length number mean size of KDP crystals for total batch times of 90 and 180 min.

be no significant variation in the final number mean length achieved. This behavior can be related to the supersaturation trajectory and level achieved during the batch. The increase of batch time allows the system to be controlled at lower supersaturation levels due to the longer growth period (after the supersaturation peak) for the minimization of AR. 4.2.3. Paracetamol Optimization Case Study. In the second case study, paracetamol was used as a model compound for the multiobjective optimization framework because it shows a nucleation dominated behavior, as shown in section 4.1. The kinetic parameters used in the optimization study are shown in Table 1. The initial and final temperature, batch time, and time step were set to be the ones utilized in the KDP optimization case study (Table 2). The objective of the optimization problem was to maximize the final mean size that can be achieved while maintaining a target shape (AR = 1). Figure 8

Figure 6. Microscopy images of the final KDP crystals obtained from an unseeded batch cooling crystallization from 40 to 20 °C following the optimal temperature profiles obtained from region: (a) I, (b) II, and (c) III. Figure 8. Pareto optimal-front obtained for the multiobjective optimization problem of maximizing the number mean length and minimizing the deviation of the aspect ratio from a target shape for an unseeded batch crystallization process of paracetamol.

Table 3. Aspect Ratio and Mean Length Obtained for the Unseeded Batch Cooling Crystallization of KDP Following the Optimal Cooling Profiles Obtained at Regions I, II, and III from 40 to 20 °C variable

AR

Ln (μm)

I II III

1.36 ± 0.13 1.85 ± 0.27 1.99 ± 0.29

374.43 ± 117 407.90 ± 169 498.65 157

shows the Pareto-optimal front obtained for the multiobjective optimization of paracetamol. It can be seen that the AR does not vary significantly; the difference between the single objectives scenarios is only of 0.01. The optimal profiles obtained using various objectives did not vary from a cubic cooling profile significantly. These results indicate that the final crystal shape obtained from a cooling crystallization of paracetamol will not vary significantly. Paracetamol is a nucleation dominated system, as described in section 4.1. The multiobjective framework shows that for a nucleation dominated system the final crystal shape practically cannot be manipulated by supersaturation (temperature profile) only. Previous work demonstrated that this type of systems show low sensitivity of the final achievable shape to variations in the temperature profile.21 For system of high growth rate variation different final crystal shapes could be

between the two scenarios. However, the number mean length obtained does not show significant increase, the range of variation of the Pareto-optimal front also being similar. The Pareto-optimal fronts obtained for the two scenarios show the compromise between achieving larger mean size and lower aspect ratio. It can be seen that for an optimal profile in which the objective is to achieve higher mean size the effect of increasing batch time is only to achieve larger mean size; the final AR does not change significantly. Similarly, if the objective is to achieve lower AR by increasing the batch time, there will H


Article


crystallization system does not belong to class I, it is recommended to study other alternatives to modify and control the size and shape of crystals such as seed crystal size distribution and amount and/or additives. Nonetheless, if the system shows high variability in the growth kinetics difference and low primary nucleation variation, significant variability in the final shape and size should be expected, hence supersaturation can be used to manipulate these properties with a compromise. Further analysis can be performed through a multiobjective optimization framework that can be used as open loop control study for later development of the temperature control approach. The framework presented provides a clear and fast method for analysis of unseeded batch cooling crystallization systems that can be used as screening of the controllability of the final size and shape due to temperature profile variations and for the estimation of the operating temperature profile that stirs the crystallization process toward desired final product properties. While the multiobjective optimization-based analysis approach is presented here for size and shape control, similar approach could be applied for the simultaneous consideration of other crystalline product property indicators, such as coefficient of variation of size and shape. Additionally, although the model compounds used in this work and the classification approach are simple organic and inorganic molecules, the approach and systematic framework proposed for the analysis of attainable regions of property control based on the classification approach, can be used for industrially more relevant complex molecules if the approximate nucleation rate and 2D growth rates are determined.

obtained by changing supersaturation trajectory through temperature profiles, as observed in section 4.2.1 for KDP. This was observed in Figure 2 in which the achievable controllability of the final size and shape for this class of system was demonstrated. The methodology presented throughout this work demonstrates that there is a clear correlation between the achievable final shape and size in a multiobjective optimization framework and the kinetic classification observed in a multidimensional cooling crystallization system. The multiobjective optimization framework studied demonstrated the expected behavior for the two classes of crystallization systems studied. Therefore, the results obtained throughout this work validates the relation between the kinetic classification and expected controllability for size and shape for inorganic and organic molecules, as described in section 4.1 and previous work.21 Additionally, the approximate characteristic form of the temperature trajectory that would lead to minimum aspect ratio versus maximum crystal size can be estimated based on the kinetic class a particular compound belongs to. These results with the classification approach can be used for the rapid initial design of crystallization processes as outlined in Figure 9. The

5. CONCLUSIONS The kinetic parameters for the crystallization of paracetamol and aspirin were identified in order to evaluate the application of a systematic classification framework to active pharmaceutical ingredients for the analysis of attainable size and shape. A multiobjective optimization framework for size and shape optimization was implemented for two model compounds: a growth dominated inorganic compound (KDP) and a nucleation dominated pharmaceutical ingredient (paracetamol). The optimization results showed the correlation between the crystallization kinetics and achievable final shape and size. Higher variability in the final shape is expected for growth difference dominated systems such as KDP. The experimental investigations confirmed the conclusions of the model-based optimization results. The proposed method can be applied in order to evaluate the design space of the achievable shape and size for unseeded batch cooling crystallization systems.

■

AUTHOR INFORMATION

Corresponding Author

Figure 9. Framework proposed to study the controllability of size and shape for unseeded batch cooling systems.

*Z. K. Nagy. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

classification framework requires that the crystallization kinetics, primary nucleation and two-dimensional growth kinetics are known for the system of interest. This could be a limitation of the method because the estimation of the twodimensional growth kinetics is time and labor consuming. However, after estimation of the crystallization kinetics, it is possible to study not only the controllability of size and shape but also the sensitivity of primary nucleation and growth kinetics with respect to the temperature profile. If a

ACKNOWLEDGMENTS This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant no. DGE-1333468. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. The support from MettlerI


Article


(10) Mostafa-Nowee, S.; Abbas, A.; Romagnoli, J. A. Optimization in seeded cooling crystallization: A parameter estimation and dynamic optimization study. Chem. Eng. Process. 2007, 46, 1096−1106. (11) Choong, K. L.; Smith, R. Optimization of batch cooling crystallization. Chem. Eng. Sci. 2004, 59, 313−327. (12) Ma, D. L.; Tafti, D. K.; Braatz, R. D. Optimal control and simulation of multidimensional crystallization processes. Comput. Chem. Eng. 2002, 26, 1103−1116. (13) Yang, Y.; Nagy, Z. K. Model-based systematic design and analysis approach for unseeded combined cooling and antisolvent crystallization (CCAC) systems. Cryst. Growth Des. 2014, 14, 687− 698. (14) Ma, C. Y.; Wang, X. Z. Closed-loop control of crystal shape in cooling crystallization of L-glutamic acid. J. Process Control 2012, 22, 72−81. (15) Liu, J. J.; Hu, Y. D.; Wang, X. Z. Optimization and control of crystal shape and size in protein crystallization process. Comput. Chem. Eng. 2013, 57, 133−140. (16) Sarkar, D.; Rohani, S.; Jutan, A. Multi-objective optimization of seeded batch crystallization processes. Chem. Eng. Sci. 2006, 61 (16), 5282−5295. (17) Sarkar, D.; Modak, J. M. Pareto-optimal solutions for multiobjective optimization of fed-batch bioreactors using nondominated sorting genetic algorithm. Chem. Eng. Sci. 2005, 60 (2), 481−492. (18) Sarkar, D.; Rohani, S.; Jutan, A. Multi-objective Optimization of Semibatch Reactive Crystallization Processes. AIChE J. 2007, 53 (5), 1164−1177. (19) 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. (20) Ridder, B. J.; Majumder, A.; Nagy, Z. K. Population balance model based multi-objective optimization of a multi-segment multiaddition (MSMA) continuous plug flow antisolvent crystallizer. Ind. Eng. Chem. Res. 2014, 53, 4387−4397. (21) Acevedo, D.; Nagy, Z. K. Systematic classification of unseeded batch crystallization systems for achievable shape and size analysis. J. Cryst. Growth 2014, 394, 97−105. (22) Hulbert, H. M.; Katz, S. Some problems in particle technology. A satistical mechanical formulation. Chem. Eng. Sci. 1964, 19, 555− 574. (23) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering, 1st ed.; Academic Press: Waltham, MA, 2000; p 355. (24) Matthews, H. B.; Rawlings, J. B. Batch crystallization of a photochemical: Modeling, control, and filtration. AIChE J. 1998, 44 (5), 1119−1127. (25) Majumder, A.; Nagy, Z. K. Prediction and control of crystal shape distribution in the presence of crystal growth modifiers. Chem. Eng. Sci. 2013, 101, 593−602. (26) Barrett, M.; McNamara, M.; Hao, H.; Barrett, P.; Glennon, B. Supersaturation tracking for the development, optimization and control of crystallization processes. Chem. Eng. Res. Des. 2010, 88 (8), 1108−1119. (27) Barrett, P.; Glennon, B. Characterizing the metastable zone width and solubility curve using Lasentec FBRM and PVM. Chem. Eng. Res. Des. 2002, 80 (7), 799−805. (28) Rawlings, J. B.; Miller, S. M.; Witkowski, W. R. Model identification and control of solution crystallization processes: A review. Ind. Eng. Chem. Res. 1993, 32 (7), 1275−1296. (29) Kempkes, M.; Eggers, J.; Mazzotti, M. Measurement of particle size and shape by FBRM and in situ microscopy. Chem. Eng. Sci. 2008, 63 (19), 4656−4675. (30) Gunawan, R.; Ma, D. L.; Fujiwara, M.; Braatz, R. D. Identification of kinetic parameters in multidimensional crystallization processes. Int. J. Mod. Phys. 2002, 16 (1), 367−374. (31) Seber, G. A. F.; Wild, C. J. Nonlinear Regression; John Wiley & Sons: New York, 1989.

Toledo Inc. for making accessible the use of the PVM is acknowledged.

■

NOMENCLATURE AR = aspect ratio b = nucleation kinetic order [dimensionless] B = birth terms in the PBE [cm4 s−1] C = solute concentration [g cm−3] C0 = initial concentration [g cm−3] Cs = solubility [g cm−3] gn = growth kinetic order of n dimension [dimensionless] Gn = crystal growth rate of n dimension [cm s−1] J = objective function vector kb = nucleation kinetic parameter [cm4 s−1] kn = growth kinetic parameter of n dimension [cm s−1] kv = volumetric shape factor [dimensionless] L = length dimension crystal size [cm] Ln = length number mean crystal size [μm] n = number density distribution [cm−4] n0 = initial number density distribution [cm−4] S = relative supsaturation [dimensionless] t = time [s] tstep = time step [min] tfinal = batch time [min] T = temperature [K] Tmin = final temperature [°C] Tmax. = initial temperature [°C] W = width dimension crystal size [cm] X = characteristic dimensions X0 = initial characteristic dimensions ρc = crystal density [g cm−3] θ = vector of kinetic parameters μij = i,j moment of number density distribution [g cmi−3]

■

REFERENCES

(1) Alvarez, A. J.; Myerson, A. S. Continuous plug flow crystallization of pharmaceutical compounds. Cryst. Growth Des. 2010, 10.5, 2219− 2228. (2) Braatz, R. D. Advanced control of crystallization processes. Ann. Rev. Control 2002, 25 (1), 87−99. (3) Fujiwara, M.; Nagy, Z. K.; Chew, J. W.; Braatz, R. D. Firstprinciples and direct design approaches for the control of pharmaceutical crystallization. J. Process Control 2005, 15 (5), 493− 504. (4) Saleemi, A. N.; Rielly, C. D.; Nagy, Z. K. Monitoring of the combined cooling and antisolvent crystallisation of mixtures of aminobenzoic acid isomers using ATR-UV/vis spectroscopy and FBRM. Chem. Eng. Sci. 2012, 77, 122−129. (5) Nagy, Z. K.; Aamir, E. Systematic design of supersaturation controlled crystallization processes for shaping the crystal size distribution using an analytical estimator. Chem. Eng. Sci. 2012, 84, 656−670. (6) Aamir, E.; Rielly, C. D.; Nagy, Z. K. Experimental evaluation of the targeted direct design of temperature trajectories for growthdominated crystallization processes using an analytical crystal size distribution estimator. Ind. Eng. Chem. Res. 2012, 51, 16677−16687. (7) Wan, J.; Wang, X. Z.; Ma, C. Y. Particle shape manipulation and optimization in cooling crystallization involving multiple crystal morphological forms. AIChE J. 2009, 55 (8), 2049−2061. (8) Patience, D. B.; Rawlings, J. B. Particle-shape monitoring and control in crystallization processes. AIChE J. 2001, 47, 2125−2130. (9) Yang, G.; Kubota, N.; Sha, Z.; Louhi-Kultanen, M.; Wang, J. Crystal shape control by manipulating supersaturation in batch cooling crystallization. Cryst. Growth Des. 2006, 6 (12), 2799−2903. J


Article

Industrial & Engineering Chemistry Research (32) Seidell, A. D. Handbook of Solubility Data for Pharmaceuticals; Van Nostrand Company: New York, 1919. (33) Mitchell, N. A.; Frawley, P. J.; Ó ’Ciardhá, C. T. Nucleation kinetics of paracetamol-ethanols solutions from induction time experiments using Lasentec FBRM. J. Cryst. Growth 2011, 321, 91−99. (34) Wan, J.; Wang, X. Z.; Ma, C. Y. Particle shape manipulation and optimization in cooling crystallization involving multiple crystal morphological forms. AIChE J. 2009, 55 (8), 2049−2061. (35) Liu, J. J.; Ma, C. Y.; Hu, Y. D.; Wang, X. Z. Modelling proteing crystallisation using morphological population balance models. Chem. Eng. Res. Des. 2010, 88, 437−446. (36) Nagy, Z. K.; Chew, J. W.; Fujiwara, M.; Braatz, R. D. Comparative performance of concentration and temperature controlled batch crystallization. J. Process Control 2008, 18 (3), 399−407. (37) Patience, D. B.; Dell’Orco, P. C.; Rawlings, J. B. Optimal operation of a seeded pharmaceutical crystallization with growthdependent dispersion. Org. Process Res. Dev. 2004, 8, 609−615.

K


Multiobjective Optimization of an Unseeded Batch Cooling Crystallizer

Recommend Documents