Feedback Control of Crystal Size Distribution for Cooling Batch Crystallization Using Deep Learning-Based Image Analysis

Abstract

The shape of the crystal size distribution directly determines the quality of crystal products. It is often assumed that distributional properties of crystal size conform to the Gaussian distribution or the log normal distribution. The mean and variance or relative crystal number are widely adopted to describe the crystal size distribution and taken as the control objectives. Therefore, the resulting control methods have difficulties in controlling the crystal size distribution with a general shape. In this article, a novel feedback control system of crystal size distribution based on image analysis is designed for the effective control of crystal size distribution with a general shape. First, a deep learning network-based image analysis method is adopted and implemented to extract the crystal size distribution. Second, the crystal size distribution is approximated by a radial basis function neural network. Consequently, a feedback controller is designed and the tracking control of the target crystal size distribution is finally realized. The results of crystallization experiments demonstrate the effectiveness of the proposed method.

1. Introduction

Cooling batch crystallization is widely applied in the production of raw material of nonferrous metals, drugs, and fine chemicals. Because the shape of the crystal size distribution directly determines the quality of crystal products, the description and control of the crystal size distribution (CSD) shape has attracted the interest of many researchers [1,2]. In a typical industrial application, feedback control of easily measured variables, such as temperature, is quite common, but feedback control of CSD is rarely done [3]. The reason stems from the online measurement of the full distribution and the distributed nature of the crystallization process.

The research of mathematical models is an effective method to describe the crystal growth dynamics. Hulburt et al. [3] and Randolph et al. [4] developed the population balance (PB) models to describe the temporal and spatial evolution of crystal particles. Grosso et al. [5,6,7] proposed the Fokker–Planck equation (FPE) and considered the natural fluctuations where the crystal size was selected as a random variable. These models provided effective tools to understand the dynamic and underlying phenomena of the cooling batch crystallization process. Semino and Ray [8] firstly addressed controllability analysis problem of crystal size distribution on the basis of the population balance model. However, the solution of these models can become complicated, time consuming, and inaccurate due to the inherent complexity. Hence, few experimental results have been reported in the literature before, in which the full mathematical models are adopted to implement feedback control. To deal with this problem, Chiu and Christofides [9] adopted the method of weighted residuals and an approximate inertial manifold to construct a low-order ordinary differential equation from the PB model. As a result, the overwhelming majority of real-time controllers adopted the reduced order moment model because of its decreased computational expense and high accuracy [10,11,12]. In the past years, the mean and variance or relative crystal number have been widely adopted to describe the crystal size distribution and taken as the control objectives [13]. Therefore, the published control methods cannot easily effectively control the distribution shape.

In addition, distribution shape control methods can also be found in many industrial systems other than crystallization process. Flores-Cerillo and MacGregor [14] controlled the full CSD of a styrene emulsion polymerization reactor adopting a PLS model to predict the weight-averaged CSD and Gaussian distribution to parameterize the distribution. Wang [15] and Guo [16] proposed stochastic distribution control algorithms to ensure the output probability density functions (PDFs) to track a desired PDF shape. Stochastic distribution control algorithms have been applied in many industrial processes such as the refining process [17] and foam size in the process of copper roughing [18]. In brief, the experimental implementation of a feedback control system that uses the crystal size distribution is still an open problem in the literature.

As mentioned above, another important issue for feedback control of CSD is online measurement of the full distribution. Process analytical technology (PAT) provides in situ information about the crystallization process and has been widely used in the real-time monitoring and control of the crystal size distribution [19,20]. FBRM is an effective PAT tool and enables the efficient design of feedback control systems for the crystallization process [21,22,23]. Szilágyi et al. [24] demonstrated an experimental implementation result of a nonlinear model predictive control system that used the chord length distribution (CLD) signal as feedback information. However, the laser light scattering based measured signal of FBRM is easily be influenced by crystal properties or operating conditions (e.g., mixing) [25]. Moreover, the chord length distribution (CLD) of crystals is statistically related to the CSD. For the control of CSD, it is necessary to restore the CSD from CLD. Hence, the time-consuming calculations [26,27] cannot be avoided and the geometry of the crystals must be known.

In recent years, image processing tools have been widely applied for crystallization processes [25]. Image processing algorithms, e.g., edge detection, region-based segmentation, and clustering segmentation have been developed to some extent [28,29,30,31,32]. These techniques can be used for the monitoring and control of the crystal size distribution. Borsos et al. [25] proposed an image analysis based real-time direct nucleation control method. Ghadipasha et al. [33,34] presented an image processing technique based online controller to reach a desired crystal mean size and standard deviation. However, the traditional image processing methods have difficulties in detecting crystals that are touching and overlapping. Moreover, another disadvantage of the traditional image analysis methods is that there are many parameters needed to be tuned according to the imaging conditions. Therefore, researchers has adopted the deep learning to extract the underlying patterns in the images [35,36]. Gao et al. [37] adopted Mask Regional Convolutional Neural Networks (R-CNNs) to measure the individual crystals. Manee et al. [38] proposed a modified deep learning neural network based on the RetinaNet, which mitigated the crystal detection problem in high-density solute. The deep learning-based image analysis method is referred to as an end-to-end method without manually tuning the threshold parameters, which makes the deep learning-based approach suitable for online feedback control of CSD for crystallization processes.

In this contribution, a novel image analysis-based CSD feedback control method is proposed for the cooling batch crystallization process. First, a deep learning network-based crystal image analysis method is adopted and implemented to extract crystal size distribution. Deep learning-based image analysis methods are becoming increasingly attractive due to the adaptability of this technique for feedback control application. Second, the crystal size distribution is approximated by a radial basis function neural network. Consequently, a feedback control system of crystal size distribution is designed and the tracking control of the target crystal size distribution is finally realized.

2. Image Analysis Based Feedback Control of CSD

2.1. Experimental Setup

The feedback control system of crystal size distribution with the experimental setup is presented in Figure 1. The crystallization experiment is carried out in a 2L jacketed glass crystallizer. The mixing in the crystallizer is realized by a variable speed stirrer. The temperature in the crystallizer is measured by a PT100 thermometer (produced by Lauda company, Lauda Koenigshofen, Germany). The temperature control device adopts the proline enhanced thermostatic controller of the Lauda company, Lauda Koenigshofen, Germany. The camera probe (produced by PharmaVision Nanosonic Technology Ltd., Qingdao, China.) is inserted into the solute, using an optical lens and a USB3.0 Vision Camera. The maximum frame rate of the camera is up to 30 fps and the minimum exposure time is up to 44 μs. The captured image has a resolution of 2000 × 1536 pixels. The imaging system is connected to an industrial computer. The crystals in the solute can be captured in real time through the matching image analysis software. At the same time, the online measurement software of crystal size distribution can read the pictures in real time to segment, extract, and compute the statistics of the crystals in the images. The data of the crystal size distribution curve are finally generated. The CSD control system is composed of an industrial computer. The distribution control software was designed and developed. The proposed distribution control method computes the next temperature profile as the target temperature and sends it to the Lauda thermostat system, which implements it through a chiller and heater.

2.2. Measurement of Crystal Size Distribution

2.2.1. Deep Learning Based Image Analysis Method

Mask R-CNN is a two-stage multi-task deep learning neural network [38,39], which can be used for target detection, target classification, semantic segmentation, instance segmentation, and other different tasks. The overall architecture of Mask R-CNN is presented in Figure 2. Mask R-CNN combines the ideas of Feature Pyramid Networks (FPN) and Fully Convolutional Network (FCN) with an FCN branch added to Fasters-RCNN, which can generate the corresponding mask. First of all, Resnet101 and FPN form the backbone network to extract feature maps from the images, and then the Region Proposal Network (RPN) generates proposals of crystal object locations. The second stage is the head of the adopted Mask R-CNN, which is composed of three branches. The box branch is used to generate the locations of crystal particles. The class branch is adopted to classify the crystal particles. The mask branch is applied to provide pixel-level predictions of crystal particles using a full convolutional network. Finally, the outputs of the three branches are integrated together to generate a pixel-level map, in which each crystal is represented by a different color. In the training process, we use the official Mask R-CNN project on the open-source framework TensorFlow and Keras (http://github.com/matterport/Mask_RCNN, accessed on 1 April 2019).

2.2.2. Kernel Density Estimation of Experimental Histogram

The number of crystals is a key issue to measure the crystal size distribution. Although a large number of crystals can increase the measurement accuracy of the crystal size distribution, it will take more time for image analysis. To test the effect of the crystal number on the results, two performance indexes [40,41] are introduced, which are reduced confidence interval and the coefficient of variation. A large number of images, comprising sufficient crystals, is considered for image analysis to calculate the crystal size distribution. The histogram can be drawn after the information of crystal size is obtained. Because the histogram can be regarded as a superposition of a series of step functions and the probability density is not continuous and differentiable, it is necessary to smooth the calculated histogram. In this work, kernel density estimation method is adopted to obtain the smooth crystal size distribution curve.

2.3. Feedback Control Algorithm of Crystal Size Distribution

As mentioned above, traditional stochastic control methods only consider the output mean and variance of crystal population. The object of the feedback control of the crystal size distribution is to control the shape of the output probability density functions (PDFs) [17,42]. For the crystallization dynamic system, the output crystal size is denoted as a continuous random variable y ∈ [a, ζ]. Let u(t) be the control input that controls the distribution of crystal size y. The probability of the output crystal size y can be characterized by its PDF γ(y, u(k)):

$$P(a<y<\zeta ,u\left(k\right))={{\displaystyle \int }}_a^{\zeta }\gamma \left(y,u\left(k\right)\right)\mathrm{d}y$$

(1)

In this contribution, it is assumed that γ(y, u(k)) is continuous and bounded with the interval [a, b], which is assumed to be known. According to the approximation principle of RBF-NN, the square root of the output PDF is approximated using the following RBF neural network:

$$\sqrt{\gamma \left(y,u\left(k\right)\right)}={\displaystyle \sum }_{i=1}^n{\omega }_i\left(u\left(k\right)\right)R_i\left(y\right)+\epsilon \left(y\right),i=1,2,\dots ,n$$

(2)

where γ(y, u(t)) represents the crystal size distribution in the crystallization process, y is the characteristic size variable of the crystal, Ri(y) represents the ith basis function, ${\omega }_i\left(u\left(t\right)\right)$ is the weight corresponding to the ith basis function, and $\epsilon$(y) represents the approximation error. The Gaussian basis function is expressed as follows.

$$R_i\left(y\right)=\exp \left(-\frac{{(x-{\mu }_i)}^2}{{\delta }_i{}^2}\right),i=1,2,\dots ,n$$

(3)

where μ represents the center value of the Gaussian basis function and δ represents the width of the Gaussian basis function.

By choosing the appropriate Gaussian basis function, the approximation error can be reduced infinitely. So, in the case of ignoring the approximation error, Equation (2) is changed to the following equation:

$$\sqrt{\gamma \left(y,u\left(t\right)\right)}={\displaystyle \sum }_{i=1}^n{\omega }_i\left(u\left(t\right)\right)R_i\left(y\right),i=1,2,\dots ,n$$

(4)

The target crystal size distribution$g\left(y\right)$ can be estimated by the same Gaussian basis functions as:

$$\sqrt{g\left(y\right)}={\displaystyle \sum }_{i=1}^n{\omega }_{gi}\left(t\right)R_i\left(y\right),i=1,2,\dots ,n$$

(5)

where ${\omega }_{gi}$ is the target weight vector.

Tracking error e(y,t) can be defined as follows:

$$e\left(y,t\right)=\sqrt{g\left(y\right)}-\sqrt{\gamma \left(y,u\left(t\right)\right)}=C\left(y\right)w\left(t\right)$$

(6)

where $C\left(y\right)=\left[R_1\left(y\right),R_2\left(y\right), \dots ,R_n\left(y\right)\right]$, $w\left(t\right)={\left[{\omega }_{g1}-{\omega }_1,{\omega }_{g2}-{\omega }_2,\cdots ,{\omega }_{gn}-{\omega }_n\right]}^T$.

The error w(t) is a function of the Gaussian basis function output y and time t. According to the continuity theory of the function, if and only if w(t) tends to 0, e(y,t) tends to 0. Therefore, the real-time crystal size distribution can be realized by controlling the real-time weight to track the target weight. Now we can approximate the crystal size distribution curve of the crystallization process by a set of selected Gaussian basis functions, and transform the control problem of the crystal size distribution into the control problem of the corresponding weights of a set of selected Gaussian basis functions.

The feedback control strategy applied is a PI algorithm by manipulating the crystallizer temperature, which is presented mathematically as follows:

$$u=K_{P}w+K_I{{\displaystyle \int }}_0^{t}wd\tau$$

(7)

where u denotes input variables, which provide the target temperature for the thermostatic control system. $K_P$ and $K_{I }$illustrate the proportional gain matrix and integral time constant matrix.

3. Results and Discussion

To exhibit the effectiveness of the presented control system, a cooling crystallization experiment with the growth of alum crystals from water was conducted. The alum crystals are produced by Aladdin company, Shanghai, China. The solubility of alum crystals is 10.8 g/100 g H₂O at 20 °C. A quantity of 1500 g distilled water and 165 g high-purity alum were added into the crystallizer and the rotating speed of the agitator was set to 150 rpm/min to evenly mix the solution.

In the training stage of the deep learning network, linear the cooling control method at the cooling rate of 0.1 °C/min is first implemented to gather enough images. Images captured from different stages can improve the accuracy and robustness of the deep learning network. In this work, a total of 250 pictures are selected from the total crystallization process as the training set, and another 70 pictures are selected as the evaluation set. Through the common data expansion method such as symmetry and rotation of the acquired images, the number of images in the training set is expanded from 250 to 1000 and number of images in the evaluation set is expanded from 70 to 280. Finally, all of the images are resized from 2000 × 1536 to 800 × 615 pixels to meet the network requirements of input. The VGG Image Annotator (VIA) image annotation tool is used for annotation of the crystals. In each image, the target areas are surround by enclosed areas created by the polygon tool of VIA along the boundary of the target area. Among them, regions marked for crystals shall not intersect with each other except the boundary. To deal with the problem of crystal overlapping and touching in the images, the modal labeling method [43] is adopted in this paper. The hidden edges due to crystal overlapping and touching can be reliably estimated by human labelers. The training process was performed on a hardware platform of R7-4800H CPU, 32 GB DDR4 RAM, 1 TB SSD, and RTX 2060 GPU. The number of training epochs was selected as 30. The network was trained for 3.7 h. The performance was evaluated by the mask average precision (AP) from the COCO evaluation metrics. The AP and AP50 were 24.17 and 59.72, respectively. The training loss curve is shown in Figure 3. Further improvement can be achieved by enhancing the training dataset.

The results of Mask R-CNN are presented in Figure 4. The raw images are padded and resized to 1024-pixel-wide squares by Mask R-CNN. After processing, the mask branch produces binary images with the same spatial size as the input images, as presented in Figure 4b. The crystals are separated from the background of the binary images. The binary images are combined with the outputs of the class branch and the box branch. The results are presented in Figure 4c. The pixel-level maps are obtained with each crystal object represented by a different color. Then, the pixel-level map is used to extract and compute the related parameters of crystals. Multidimensional information can be extracted from the crystallization process through the image analysis method. Crystals can be described in area, length, width, etc. The relative number of crystals can also be counted. Without loss of generality, the size of a crystal object is characterized by the width. This information is then used to represent the crystal size distribution by means of statistics and kernel density estimation.

In the inference stage, a test was performed on a set of 4000 crystals. Figure 5 presents the variations in the reduced confidence interval and the coefficient of variation. A minimum of 1000 crystals should be analyzed to ensure a reduced confidence interval smaller than 5%. As shown in Figure 5, the effect of the change in the slide around 2500 crystals can be seen on the coefficients of variation. Therefore, the number of crystals is recommended to exceed 2500. Figure 6 shows an illustrative example of the result of kernel density estimation method, i.e., the smoothed crystal size distribution is compared with the histogram of crystal particles.

As mentioned above, the radial basis function neural network is adopted to approximate the actual CSD. In the following experiment, we selected three Gaussian basis functions to validate the approximate effects on the actual CSD. The central value vector corresponding to the Gaussian basis function is: [335.4; 563; 816.6] and the width vector is: [22,782; 26,505; 27,000]. The target distribution and the estimated distribution of crystals are presented in the Figure 7. The obtained target weight vector is: [0.5133; 0.3019; 0.1413]. Hence, the control of the distribution of crystal size is transformed into the control of the target weight vector representing the crystal size distribution.

Before starting the online control of crystal size distribution, the open-loop control method is used to make the system close to the asymptotic region of target shape of crystal size distribution. Firstly, the temperature in the crystallizer is rapidly raised to 35 °C by the thermostatic system and maintained for 30 min until the crystal is completely dissolved. Then, the temperature in the crystallizer is reduced to 20 °C at the cooling rate of 0.1 °C/min and the temperature is kept for 10 min to obtain the stable saturated solution. The temperature in the crystallizer is further reduced to 19 °C at the cooling rate of 0.1 °C/min. At the same time, 7 g of evenly ground seed is added into the crystallizer and maintained at 19 °C for 10 min to promote nucleation and to speed up the growth of alum crystals. Finally, the temperature in the crystallizer is further reduced at a cooling rate of 0.1 °C/min. Images are started to be captured and crystal size distribution of alum crystals is measured from the captured images. Then, the proposed controller is put into use. The sampling time is selected as 10 min to ensure that the image processing steps can be completed. Figure 8 illustrates the size distribution evolution of the crystals during the whole batch. Figure 9 describes the changes in real-time weights. Figure 10 describes the set temperature and real-time temperature variation in the crystallizer. At beginning of the cooling crystallization experiment, the error between target curve and actual curve is clear, as can be seen from Figure 8 and Figure 9. The initial temperature is about 12.4 °C and then it is gradually reduced to enhance crystal growth under the action of the proposed control method. After 60 min, when the crystal size distribution curve exceeds its target crystal size distribution, the temperate increases and crystals begin to dissolve such that the sizes of crystals are reduced. Although there are two overshoots in the crystal size distribution curve, the system output is kept to the desired crystal size distribution curve after about 160 min. It can be seen that the CSD and the related weights gradually approach their targets. It can also be seen from Figure 10 that the distribution controller frequently regulates the temperature of the solution to ensure the crystal size distribution curve tracking the target crystal size distribution curve. Hence, the regulation of CSD is implemented by the circulation of growth and dissolution due to temperature control. The crystallization process reaches a steady state. In addition, the whole evolution trend of the CSD and the related weights is similar, as can be observed from the results of Figure 8 and Figure 9. These results are a further proof of the basic idea that the control problem of the crystal size distribution can be transformed into the control problem of the corresponding weights. It is worth noting that the proposed PI controller has constant parameters; therefore, it cannot perfectly deal with the complex nonlinear effect of crystallization process. It can be seen from Figure 8 and Figure 9 that there is an obvious overshoot in the CSD and the weights’ trajectory, and, finally, there is a small deviation from the target in the steady state. Nonetheless, physical experiments show that this method can effectively track the target crystal size distribution and achieve good experimental results.

4. Conclusions

In the current work, a feedback control system of crystal size distribution based on image analysis is proposed. The image analysis technology using deep learning and kernel density estimation is adopted to realize the online measurement of CSD. The weights of the CSD computed at different times are estimated using RBF basis functions, and the control problem of the CSD is transformed into the control problem of weights representing the crystal size distribution function. A feedback controller is designed to track the target CSD shape. Overall, the proposed deep learning-based image analysis provides a new direction for real-time characterization of CSD. Furthermore, the results obtained corroborated that the proposed control method represents a powerful tool for effective control of CSD with a general shape.