ARTICLE IN PRESS

Available online at www.sciencedirect.com

ScienceDirect
Resource-Efficient Technologies ■■ (2017) ■■–■■
www.elsevier.com/locate/reffit

Research paper
Modelling and control of crystallization process
S.K. Jha 1, S. Karthika 2, T.K. Radhakrishnan 2,*
1
Department of Instrumentation and Control Engineering, National Institute of Technology, Tiruchirappalli 620 015, India
2
Department of Chemical Engineering, National Institute of Technology, Tiruchirappalli 620 015, India
Received 30 June 2016; received in revised form 6 January 2017; accepted 9 January 2017
Available online

Abstract
Batch crystallizers are predominantly used in chemical industries like pharmaceuticals, food industries and specialty chemicals. The nonlinear
nature of the batch process leads to difficulties when the objective is to obtain a uniform Crystal Size Distribution (CSD). In this study, a linear
PI controller is designed using classical controller tuning methods for controlling the crystallizer outlet temperature by manipulating the inlet
jacket temperature; however, the response is not satisfactory. A simple PID controller cannot guarantee a satisfactory response that is why an
optimal controller is designed to keep the concentration and temperature in a range that suits our needs. Any typical process operation has
constraints on states, inputs and outputs. So, a nonlinear process needs to be operated satisfying the constraints. Hence, a nonlinear controller like
Generic Model Controller (GMC) which is similar in structure to the PI controller is implemented. It minimizes the derivative of the squared error,
thus improving the output response of the process. Minimization of crystal size variation is considered as an objective function in this study. Model
predictive control is also designed that uses advanced optimization algorithm to minimize the error while linearizing the process. Constraints are
fed into the MPC toolbox in MATLAB and Prediction, Control horizons and Performance weights are tuned using Sridhar and Cooper Method.
Performances of all the three controllers (PID, GMC and MPC) are compared and it is found that MPC is the most superior one in terms of settling
time and percentage overshoot.
© 2017 Tomsk Polytechnic University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license
(http://creativecommons.org/licenses/by-nc-nd/4.0/).
Keywords: Batch crystallization; Crystal size distribution; Crystal size variation; Optimal control; Constraints; Model predictive control

1. Introduction concentration has been possible and the advanced control con-
cepts related to model free controller design and robust control
Crystallization is one of the oldest and widely used unit
of crystal shape and size have been developed [1].
operations in chemical industry. It is commonly used in phar-
maceuticals, fine chemicals and food industries. It is usually the
first step in any separation process and produces high-purity 1.1. Operations of batch crystallizer
products. The batch crystallization mode is widely adopted in A pictorial representation of batch cooling crystallizer is
industries for the production highly pure products in small shown in Fig. 1, in which the crystals are formed when the
quantities. The methodologies of batch crystallization model- solubility of a solute in a solution decreases as a result of
ling are well established and the control strategies have under- cooling the solution. Decrease in the solubility can also be
gone many refinements over the past decades. Due to achieved by evaporating the solvent from the solution or by
advancement of sensor technology, measurement of solution adding another solvent which precipitates the solute. The crys-
tallization process consists of two major events, nucleation and
crystal growth which continue to occur simultaneously, and
International Conference on Separation Technologies in Chemical, their rate is driven by the existing supersaturation in the solu-
Biochemical, Petroleum and Environmental Engineering (TECHNOSCAPE tion. Supersaturation can be altered by changing the operating
2016). conditions [2]. Depending upon the conditions, either nucle-
* Corresponding author: Department of Chemical Engineering, National
Institute of Technology, Tiruchirappalli 620 015, India. Tel.: 0431 2503104; fax:
ation or growth may be predominant over the other, and as a
91-431-250133. result, crystals with different sizes and shapes are obtained. As
E-mail address: [email protected] (T.K. Radhakrishnan). the time progresses, supersaturation level decreases and a stage
http://dx.doi.org/10.1016/j.reffit.2017.01.002
2405-6537/© 2017 Tomsk Polytechnic University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license
(http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer review under responsibility of Tomsk Polytechnic University.

Please cite this article in press as: S.K. Jha, S. Karthika, T.K. Radhakrishnan, Modelling and control of crystallization process, Resource-Efficient Technologies (2017), doi:
10.1016/j.reffit.2017.01.002
 ARTICLE IN PRESS
2 S.K. Jha et al. / Resource-Efficient Technologies ■■ (2017) ■■–■■

modelling of batch cooling crystallizer, and method of

moments is discussed. In section 3, the design principles of
conventional PID controller, Generic Model Controller (GMC)
and the Model Predictive Controller (MPC) are presented and
the performance of these controllers is discussed in the subse-
quent section 4.

2. Population balance modelling of batch cooling

crystallizer
Modelling of the batch crystallizer is done by first-principles
method by solving Population Balance Equation (PBE), mass
and energy balance equations simultaneously [6–8]. PB equa-
tion describes the behaviour of particulate system like crystal-
Fig. 1. Jacketed batch cooling crystallizer. lizer. It consists of partial differential equation which describes
the evolution of CSD as a function of time and the crystal size.
One dimensional PBE for a batch crystallization system is
given as
of equilibrium is reached where the crystallization process
reaches a steady state. After the batch time, the slurry is dis- ∂n ( L, t ) ∂Gn ( L, t )
+ = Bo (1)
charged from the crystallizer and crystals are separated by ∂t ∂L
filtration and dried. The size of the product formed is very
important since it helps in increasing the efficacy of down- where n(L, t) is the number density function (1/m4) which
stream processes like filtration, washing etc. Hence, the prime describes the number of crystals per crystal length and volume
objective of crystallizer is to obtain the crystals with narrower of slurry, t is the time (s), L is the characteristic size of crystals
size distribution. (m), G is the crystal growth rate (m/s) and Bo is the nucleation
The operating region in the crystallizer can be identified rate of crystals of infinitesimal size (1/m3s). The initial and
between a series of clear points (solubility) and cloud points boundary conditions for solving the PBE are given as
(metastable limit) at various temperatures. For every solute– B0
solvent combination, this can be obtained by experimentation n ( L, 0 ) = 0, n (0, t ) = (2)
G
and the region bounded by these curves is termed as the meta-
stable zone. In this zone, spontaneous nucleation does not occur If growth rate is assumed to be independent of crystal length
but it can be controlled. Hence, in cooling crystallization the and initial birth rate is considered to be negligible, the simpli-
solute concentration profile is maintained well within the meta- fied form of PBE can be written as
stable zone and close to the solubility curve in order to promote
∂n ( L, t ) ∂n ( L, t )
crystal growth and avoid spontaneous nucleation. +G =0 (3)
The classical approach to control is based on first principles ∂t ∂L
models in which a model is constructed from material, energy The solute mass balance equation provides the change of
and population balances ignoring agglomeration and breakage concentration inside the crystallizer as the time progresses. The
phenomenon. The growth and nucleation kinetics need to be mass balance is written in the form of ordinary differential
accurately determined and the model should give the robust and equation
accurate description of process. The direct design approach to
dC ∞
control is based on following a desired trajectory within the = −3ρc K v ∫ nGL2 dL (4)
operating region. The desired trajectory of the crystallizer falls dt 0

in the metastable zone which is bounded by solubility and where C is the solute concentration, ρc is the density of crystals
metastable curves. Here, a concentration trajectory as a func- and Kv is the volumetric shape factor. The energy balances for
tion of temperature is followed with simple feedback control the suspension and the cooling jacket are given below
system which is very simple to adopt and does not require the
accurate prediction of nucleation and growth kinetics [3–5]. dT − UA 3ρ K ΔH ∞
The objective of this study is to analyse and compare the
= (T − Ti ) − c v ∫0 nGL2dL (5)
dt MC p Cp
performance of the conventional and the advanced controllers
which track the desired cooling profile. The temperature control dTi Fi UA
strategy is the most widely adopted technique because of its = (Tisp − Ti ) + (T − Ti ) (6)
dt Vi ` i Vi C pi
simplicity since it requires only temperature of the crystallizer
to be controlled and is used in the present study. Potassium where T is the temperature of the crystal suspension at any time
nitrate–water system is taken as a model system and the inlet t, Ti is the time dependent temperature of fluid in the jacket
jacket temperature is taken as the manipulated variable. The (water), Fi, Vi, ρi and Cpi are respectively the flow rate, volume,
paper is organized as follows. In section 2, population balance density and specific heat of the inlet jacket fluid respectively,

Please cite this article in press as: S.K. Jha, S. Karthika, T.K. Radhakrishnan, Modelling and control of crystallization process, Resource-Efficient Technologies (2017), doi:
10.1016/j.reffit.2017.01.002
 ARTICLE IN PRESS
S.K. Jha et al. / Resource-Efficient Technologies ■■ (2017) ■■–■■ 3

Table 1
Parameters used for simulation of crystallization process.
Parameters Definition Value Units
b Nucleation rate exponent 1.78
Kb Nucleation rate constant 4.64*1011 1/m3.s
g Growth rate exponent 1.32
Kg Growth rate constant 1.1612*10−4 m/s
C0 Initial concentration of solution 18 kg KNO3/kg water
ρc Density of crystals 2110 kg/m3
ρi Density of water 1000 kg/m3
Kv Volume shape factor 1
M Mass of slurry 30 kg
Cp Specific heat capacity of solution 95.06 J/kg.K
Cpi Specific heat of water 4184 J/kg.K
Vi Volume shape factor 0.905 m3
Fi Flow rate of water 0.001 kg/s
U.A. Areas times overall heat transfer coefficient 0.125 J/K.s
μ0 (0) Zeroth moments (total number of crystals) 10−6
μ1 (0) First moment (total length of crystals) 10−6
μ2 (0) Second moment (total area of crystals) 10−6
μ3 (0) Third moment (total volume of crystals) 10−6
n (L, t) number density function which describes the number of crystals per crystal (1/m4)
length and volume of slurry
S Supersaturation which defines the driving potential for the nucleation and growth
to occur

Tisp is the set point temperature provided, U and A are the the model equations, which have the following physical
overall heat transfer coefficient and area available for heat significance; μ0 represent total crystal number per unit volume,
transfer respectively, ρ is the density of the slurry and ΔH is the μ1 represent total crystal length per unit volume, μ2 represent
heat of crystallization. The growth rate G and birth rate B rate total crystal surface area per unit volume, μ3 represent total
at any time t in the crystallizer are obtained from empirical crystal volume per unit volume. The rate equations for different
equations moments can be written as
G = K g .Sg (7) dμ 0
= B (t) (10)
dt
B = K b .Sb (8)
dμ i
where Kg and g are growth rate parameters and Kb and b are = i ∗ G ( t ) ∗ μ in−1 (11)
dt
birth rate parameters and S is the supersaturation which is
defined as the driving potential for nucleation and growth to All the above mentioned ordinary differential equations,
occur. In terms of concentration, S = (C − Cs)/Cs, where C is the mass and energy balance equations are solved in SIMULINK
solution concentration and Cs is the concentration at saturated environment of MATLAB software using ODE45 solver as well
conditions. as forward Euler method. The parameters are taken from Miller
[9], for carrying out the simulation and is shown in Table 1.
2.1. Method of moments Metastable zone width is assumed to be constant having a value
of ΔT = 2.2 K for temperature range used. The saturation con-
PBE as described in equation (1) is in the form of a partial centration (Cs) for KNO3–water system is given by the equation
differential equation. The numerical schemes for solving the
PBE are very complex and time consuming. In order to reduce Cs = 0.1286 + 0.00588.T2 + 0.0001721.T3 (12)
the model complexity, it is advisable to convert this partial
differential equation into a set of ordinary differential equa- 3. Design of controllers
tions. The most widely used method to reduce the complexity is 3.1. PID controller design
to use the method of moments. In principle, the method of
moments defines the ith moment in terms of the population Crystallization process can be controlled by two ways; either
density function by by controlling the crystallizer temperature or by controlling the
crystallizer concentration. Inlet jacket temperature is used as
∞
μ i = ∫ n ( L, t ) Li ( t ) dL (9) the manipulated variable. As the jacket temperature is changed,
0
concentration and temperature of the crystallizer is changed.
where μi is the ith moment and n(L, t) is the population density Since, concentration control requires some expensive sensors,
function. Usually, first four moments are needed for solving temperature control is widely adopted in literatures. A basic

Please cite this article in press as: S.K. Jha, S. Karthika, T.K. Radhakrishnan, Modelling and control of crystallization process, Resource-Efficient Technologies (2017), doi:
10.1016/j.reffit.2017.01.002
 ARTICLE IN PRESS
4 S.K. Jha et al. / Resource-Efficient Technologies ■■ (2017) ■■–■■

PID controller is designed for temperature control strategy. I = ∫ ( y d − y ) dt

2
(15)
Output temperature of the crystallizer is fed back and difference
between the set point and measured output is calculated. This where y and y d denote the derivative of the process output and
error signal is then fed to the controller and an actuating signal the desired trajectory. To minimize the performance index, let
is then generated from the controller so that the output reaches dI/dt equals zero, then
the set point. Zeigler–Nichols tuning settings are used for
tuning the PID parameters. ( y d − y ) = 0 (16)

3.2. Generic model control (GMC) Now, y can be written as,

A PI controller is also designed for this case. But the dy dy dx

= ∗
response was not satisfactory. It took a lot of time for the system dt dx dt
to reach steady state value. Controller efficiency was not dy dh ( x )
improving even on changing the parameter values. So there was Now, = and,
dx dx
a need of designing a different controller which is as simple in
dy dh dh
structure as a PI controller but which will give a good output = ∗ f (x) + ∗ g (x) ∗ u (t) (17)
response. Hence, a Generic Model Controller for temperature dt dx dx
control of exothermic batch crystallizer is designed in which
Combining the above equations and putting in Equation
rate of change of temperature is used as a manipulated variable.
(16), we will get
3.2.1. Generic model controller algorithm formulation
GMC has several advantages that make it better from other
control algorithms. The process model appears directly in the
K1 ( y d − y ) + K 2 ∫ ( y d − y ) dt − { dh
dx
dh
∗ f (x) + ∗ g (x) ∗ u (t) = 0
dx }
(18)
control algorithm and it does not need to be linearized before
use, allowing for the inherent nonlinearity of the operation to be From the above equation, u (t) can be calculated in terms of
taken into account. It can be used for discrete as well as con- K1 and K2.
tinuous process plant model. By design, GMC provides feed-
back control of the rate of change of the process output. This 3.2.2. Tuning of the parameters
suggests that the rate of temperature change (output) can be Several methods have been proposed to find the values of K1
used directly as a manipulated variable. The relationship and K2. The common methods used are
between feed forward and feedback control is explicitly stated 1. Finding the ultimate gain and ultimate period of the
in the GMC algorithm. Finally and importantly, the GMC output response and multiplying by a scalar.
framework permits us to develop a control algorithm that can be 2. By calculating the damping coefficient (ƺ) and time con-
used for both heat-up and temperature maintenance and there- stant (τ) of the desired trajectory. K1 and K2 can be
fore eliminates the need for a switching criterion between dif- written as:
ferent algorithms; which makes it very robust in nature.
The Generic Model Control (GMC) algorithm derived for
the control of batch reactor (Lee and Sullivan, 1988) is used for 3. By trial and error method.
control of crystallizer in this study [10]. GMC requires a
dynamic model of the process written in standard state variable 3.3. Model predictive control (MPC)
form. The controller is formulated by solving the dynamic
process model for the derivative of the controlled variable, y, In chemical industries, there are certain inequality con-
and letting it equal what is, in effect, a proportional integral straints on input and output variables that need to be considered
term operating on the difference between the current value of y for optimal operation. So, in order to follow those constraints,
and its desired value, yd. The desired trajectory can be assumed some advanced controllers need to be designed that follow
as, advanced optimization algorithm to minimize the error. Model
predictive controller predict the future output if an accurate
y d i = K1 ( y d − y ) + K 2 ∫ ( y d − y ) dt (13) dynamic model of process is available and then depending on
predicted and measured outputs, appropriate change in mea-
where K1 and K2 are tunable parameters. Consider a non-linear sured input can be calculated.
system Basic structure of MPC is given in Fig. 2. Current values of
the output variables are calculated using a process model
x = f ( x ) + g ( x ) u ( t ) (14)
described in section 2 and then difference between predicted
y = h (x) and actual outputs are used as a feedback signal to a prediction
block. The predicted outputs are used in controlled calculation
where f(x), g(x) and h(x) are vector fields. The scalar field and set point calculation after considering constraints on the
describing the output function is simply the state itself. The input and output variables. MPC configuration is analogous to
performance index for GMC is given by [11], both internal model control configuration and smith predictor

Please cite this article in press as: S.K. Jha, S. Karthika, T.K. Radhakrishnan, Modelling and control of crystallization process, Resource-Efficient Technologies (2017), doi:
10.1016/j.reffit.2017.01.002
 ARTICLE IN PRESS
S.K. Jha et al. / Resource-Efficient Technologies ■■ (2017) ■■–■■ 5

Fig. 2. Block diagram of MPC.

configuration because model and process are parallel acted and MPC toolbox model advisor is used to calculate the perfor-
difference act as feedback control signal. But coordination of mance weights. MPC toolbox of the SIMULINK environment
the control and set point calculations makes MPC superior than is used to carry out the simulation. Input and output constraints
others. Moreover, MPC is largely used in MIMO control prob- are provided and weights and gain are tuned to find out the
lems than IMC or Smith predictor. optimal performance.
The main objective of the MPC is to determine the set of
control moves so that predicted model output reaches the 4. Results
desired set point. At any sampling time, MPC decides a set of M
A classical PID controller to control the temperature of
input values containing a present input u(k) and future (M-1)
the crystallizer has been proposed in which temperature of the
inputs. Set of inputs are calculated in such a way that a set of P
jacket fluid is taken as the manipulated variable to control the
predicted output should reach the optimum set point. Objective
crystallizer temperature. The parameter values for the PID con-
function is optimized to find out control moves. Number of
trollers were obtained minimizing the integral square error
control moves M is called the control horizon while number of
(ISE). The response of PID controller for the set point tracking
prediction P is referred to as the prediction horizon.
is shown in Fig. 3. Ultimate period and gain are calculated by
The least square objective function for a prediction horizon
increasing the value of Kc until a sustained oscillation curve is
P and a control horizon M is written as follows,
obtained. Zeigler–Nichols tuning settings are given by Kc = 0.6
P 2 M −1
∅ = Q.∑ i =1 ( rk +1 − yˆ k + i ) + R.∑ i = 0 Δu 2 K + i (19) Ku, τI = PU/2,τD = PU/8. For this case, it is found out to be
Kc = 1.8, τI = 0.6, τD = 1.35. The controller is tracking the set
where ŷ denotes the predicted model outputs, ‘r’ is the set point changes for different types of step changes provided. But,
point, Δu is the change in manipulated input from one sample the overshoot appears to be higher and it takes a minimum of
time to the next, ‘Q’ and ‘R’ is a weight for the change in the 53 min for the response to settle at steady value. The response
output and manipulated input respectively. The optimization of GMC for set point tracking is shown in Fig. 4. It provides a
problem deals with a minimization of the objective function by better response when compared with PID controller. The over-
manipulating control moves. shoot has reduced and the response settles after 49 min. Both
controllers were compared (Fig. 5, Table 2) in terms of percent-
3.3.1. Tuning strategy review age overshoot, settling times and it was found that GMC
response was quite better than that of the PID. By designing a
1. Approximate the process dynamics of all controller MPC controller with prediction horizon as 51 and control
output–process variable pairs with FOPDT models horizon as 11 with sampling times of 0.5 s for 2000 s, con-
(Shridhar and Cooper, 1997) [12]. trolled state variables of batch crystallization process are
obtained and are given in Fig. 6. There is a drastic decrease in
y (s ) K.e − θs
=
u (s ) τs + 1

2. Select the sample time as close as possible to T = Min Table 2

(0.1τ, 0.5θ) Comparison between GMC, PID and MPC.
3. Compute the prediction horizon (P) and control horizon Controller Percent Settling ISE
overshoot time (mins)
(C) using
PID 18% 53 594.55
5τ θ τ θ GMC 10% 49 284.39
P = + +1 c = + +1 MPC 1% 16.67 50.73
T T T T

Please cite this article in press as: S.K. Jha, S. Karthika, T.K. Radhakrishnan, Modelling and control of crystallization process, Resource-Efficient Technologies (2017), doi:
10.1016/j.reffit.2017.01.002
 ARTICLE IN PRESS
6 S.K. Jha et al. / Resource-Efficient Technologies ■■ (2017) ■■–■■

Fig. 3. PID response for set point tracking.

the percentage overshoot and the response settles sooner when response is satisfactory. Temperature set point tracking is
compared with other controllers. The MPC response is quite achieved for every case by keeping the thermodynamics of the
better than other controllers in tracking the set point. system in mind. Simulation is done for quite a number of
operating points. The MPC toolbox of MATLAB simulation
block is used for implementing the MPC controller. Several
5. Conclusions
constraints on input and output variables are provided. Perfor-
The primary objective of this study is to evaluate the oppor- mance of MPC is then compared with that of PID and Generic
tunities for model based control to improve product quality Model controller (GMC) and it is found that MPC gives a better
and process productivity of industrial crystallizers. All the con- response based on percent overshoot, settling time and integral
trollers have been successfully implemented and the output squared error.

Fig. 4. GMC response for set point tracking.

Please cite this article in press as: S.K. Jha, S. Karthika, T.K. Radhakrishnan, Modelling and control of crystallization process, Resource-Efficient Technologies (2017), doi:
10.1016/j.reffit.2017.01.002
 ARTICLE IN PRESS
S.K. Jha et al. / Resource-Efficient Technologies ■■ (2017) ■■–■■ 7

Fig. 5. Comparison of PID and GMC responses for temperature set point tracking.

Fig. 6. MPC response for set point tracking.

References [6] J.W. Mullin, Crystallization, vol. 7, third ed., Butterworth-Heinemann,

London, 1993. 354.
[1] Z.K. Nagy, R.D. Braatz, Advances and new directions in crystallization [7] A.D. Randolph, M.A. Larson, Theory of Particulate Processes, Academic
control, Ann Rev Chem Biomol Eng 3 (2012) 55–75. Press, New York, 1998.
[2] A.G. Jones, J.W. Mullin, Programmed cooling crystallization of [8] D. Ramakrishna, Population Balances-Theory and Applications to
potassium sulfate solutions, Chem. Eng. Sci. 29 (1974) 105–117. Particulate Systems in Chemical Engineering, Academic Press, San
[3] J.B. Rawlings, S.M. Miller, W.R. Witkowski, Model identification and Diego, CA, 2000.
control of solution crystallization process: a review, Ind. Eng. Chem. Res. [9] S.M. Miller, Modeling and quality control strategies for batch cooling
32 (1993) 1275–1296. crystallizers, Ph.D. thesis, University of Texas at Austin 1993.
[4] Z.K. Nagy, J.W. Chew, M. Fujiwara, R.D. Braatz, Advances in the [10] P.I. Lee, G.R. Sullivan, Generic model control, Comp. Chem.Eng. 12 (6)
modeling and control of batch crystallizers, in: Proc. of the IFAC Symp. (1988) 573–580.
on Advanced Control of Chemical Processes, Elsevier Scientific, Oxford, [11] B.J. Cott, S. Macchietto, Temperature control of exothermic batch reactors
UK, 2004, pp. 83–90. using generic model control, Ind. Eng. Chem. Res. 28 (1989) 1177–1184.
[5] M. Fujiwara, Z.K. Nagy, J.W. Chew, R.D. Braatz, First-principles and [12] R. Shridhar, D.J. Cooper, A tuning strategy for unconstrained SISO model
direct design approaches for the control of pharmaceutical crystallization, predictive control, Ind. Eng. Chem. Res. 36 (3) (1997) 729–746.
J Proc Contr 15 (5) (2005) 493–504.

Please cite this article in press as: S.K. Jha, S. Karthika, T.K. Radhakrishnan, Modelling and control of crystallization process, Resource-Efficient Technologies (2017), doi:
10.1016/j.reffit.2017.01.002