Identification of a non-linear dynamic model of the bubble size distribution in a pilot flotation column

International Journal of Mineral Processing 145 (2015) 7–16

Contents lists available at ScienceDirect

International Journal of Mineral Processing journal homepage: www.elsevier.com/locate/ijminpro

Identification of a non-linear dynamic model of the bubble size distribution in a pilot flotation column A. Riquelme a, A. Desbiens a,⁎, R. del Villar b, M. Maldonado c a b c

Department of Electrical and Computer Engineering, LOOP (Laboratoire d'Observation et d'Optimisation des Procédés), Université Laval, Pavillon Pouliot, Québec City G1V 0A6, Canada Department of Mining, Metallurgical and Materials Engineering, LOOP (Laboratoire d'Observation et d'Optimisation des Procédés), Université Laval, Pavillon Pouliot, Québec City G1V 0A6, Canada Department of Metallurgical Engineering, University of Santiago, Chile

a r t i c l e

i n f o

Article history: Received 22 August 2014 Received in revised form 24 September 2015 Accepted 3 November 2015 Available online 6 November 2015 Keywords: Flotation column Wiener model Non-linear model Bubble size distribution Image processing Hough transform

a b s t r a c t Gas dispersion properties play an important role in flotation as they partially determine the metallurgical performance. The objective of this work is to obtain a dynamic model of the bubble size distribution in a twophase (air and water) pilot flotation column. The steps are: 1) to measure and count the bubbles from digital images taken by a camera, 2) to estimate a log–normal distribution of the bubble sizes, 3) to estimate a Wiener model whose outputs are the mean and standard deviation of the distribution while the inputs are the superficial gas velocity set-point and the superficial shearing water velocity set-point. The size and number of bubbles in each image are evaluated by a bubble detection technique based on circular Hough transform (CHT). This technique allows overcoming issues related to the detection of large single bubbles as well as clusters. Tests are carried out in the laboratory flotation column using different concentrations of frother and air flow rates. Results are compared with visual counting, as well as with a commonly used bubble detection method based on circular particle detection (CPD). The estimated number of bubbles is very similar to what is obtained with a visual count. Compared to CPD algorithm, CHT significantly improves D32 estimation (error of 3% instead of 18% with the former) with a comparable processing time. Then, the mean and standard deviation of a log–normal distribution are estimated by maximizing a likelihood function thereby leading to very good fits for the distributions. Finally, a series of model identification tests are conducted by manipulating the superficial gas velocity setpoint and the superficial shearing water velocity set-point (i.e. air and shearing water flow rate set-points through the sparger). A Wiener model is then estimated to predict the corresponding mean and standard deviation of the log–normal distribution. Validation tests confirm the quality of the non-linear model. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Column flotation is a mineral processing method used to separate minerals based on the difference in surface hydrophobicity. The process consists in injecting air bubbles into a cylindrical vessel (typically 12 to 14 m height) where ground ore slurry is introduced at approximately 2 m from the top. Air is dispersed into bubbles at the bottom of the column through a sparger manifold. Prior to its introduction into the flotation column, the slurry is usually conditioned with proper chemicals rendering the surface of valuable minerals hydrophobic, in other words inducing their attraction to air bubbles. Particles exhibiting hydrophobic surfaces thus attach to the rising air bubbles. These bubble-particle aggregates then rises to the top of the column forming a froth phase rich in the valuable minerals. The overflowing froth at the top of the cell (usually the valuable product) is called the concentrate. Hydrophilic particles, which do not attach to the rising bubbles,

⁎ Corresponding author. E-mail address: [email protected] (A. Desbiens).

http://dx.doi.org/10.1016/j.minpro.2015.11.003 0301-7516/© 2015 Elsevier B.V. All rights reserved.

settle down and exit the column at the bottom port, forming what is called the tailings. Gas dispersion properties play an important role in flotation as they partially determine the metallurgical performance. One of the most relevant variable is the bubble surface area flux (Sb), which is usually calculated as a function of the superficial gas velocity (JG) and the Sauter bubble mean diameter (D32) determined from the bubble size distribution (BSD). The Sauter bubble diameter is defined as the diameter of a set of identical size bubbles having the same total surface and volume of the original BSD. After some simplifications the following formula is obtained:

D32 ¼

i X

i X 3 2 di = di

n

n

ð1Þ

where di is the ith observed bubble diameter and n is the number of counted bubbles. Sb actually represents the rate of bubble surface available for particle collection. On this regard, it is worth noting that Yianatos and Contreras (2010) developed a model for estimating the

8

A. Riquelme et al. / International Journal of Mineral Processing 145 (2015) 7–16

Fig. 1. Sample image of bubbles.

carrying capacity (maximum bubble coverage by particles) of a given flotation machine as a function of the bubble Sauter mean diameter and Sb, among some other variables. However, a given D32 can be obtained with dissimilar shapes of bubble size distributions (Maldonado et al., 2008), thus leading to multiple possible flotation recoveries. This is in accordance with the work by Heiskanen (2000) suggesting that matching the bubble size distribution to the particle size distribution of the valuable mineral would increase its recovery. The same concept also applies to the Sb, as a same value can be obtained from different JG and D32 combinations. To control the entire BSD thus becomes essential to investigate the relationship between the flotation kinetics and hydrodynamic characteristics of the process. In order to design such a control strategy a proper measurement for BSD is required. Once this measurement is obtained, it would then be possible to estimate a dynamic model related to the manipulated variables available for control purposes. It is thus clear that a reliable technique for measuring bubble size diameters is necessary. Rodrigues and Rubio (2003) reviewed different methods presently used to estimate the BSD, which include drift flux

Fig. 2. Original image and results after each step of bubble detection processing.

A. Riquelme et al. / International Journal of Mineral Processing 145 (2015) 7–16

9

Fig. 3. Accumulator matrix: calculation of the matrix elements for three points (for a target radius). The elements take the value 0, 1 or 2.

analysis, electro-resistivity, optical, and image analysis. These methods are also categorized as intrusive or non-intrusive. Image analysis is probably nowadays the most widely used method for bubble size measurement, since it is the most straightforward technique. However, it should be emphasized that vision-based techniques require a high computing capacity. Recent developments in digital image processing provide new opportunities for developing more specific and efficient algorithms to detect and measure bubbles in dispersions. Automated processing methods generally exhibit three distinct stages: filtering, edge detection and solid particle detection. Particle detection aims at identifying circular elements (Han et al., 2002; Lin et al., 2007; Zhu et al., 2001), thus excluding bubble clusters and large bubbles which are often non-spherical. Other methods include an off-line manual stage to take into account clustered and larger bubbles (Vinnett et al., 2012), thus excluding possibilities of on-line process control applications. The Hough transform (HT) is frequently employed to detect edges. It can isolate features of a particular shape within an image (Duda and Hart, 1972). The calculation complexity of applying the HT is proportional to the number of parameters needed to describe the pattern (Duda and Hart, 1972). The circular Hough transform (CHT) is a variation of the general HT, where the target shape is a circle. This transform is recognized as a robust technique for curve detection, and it is widely used to detect overlapped circles (Daraei et al., 2011; Goulermas and Liatsis, 2012; Nostrati and Karimi, 2011). This method will be applied in this paper, with some image processing techniques, to properly estimate the BSD. To eventually control the BSD in flotation processes, it is required to correctly measure the BSD. Density estimation deals with the problem of estimating the probability density function (PDF) of a set of data samples (bubble sizes in this case). The different existing methods for density estimation can be classified into three groups: 1) parametric methods, 2) non-parametric methods, and 3) semi-parametric methods. Parametric methods are based on a defined and restricted function, leaving just a few parameters to be determined (depending on the objective function). The most common PDFs are Gaussian, exponential,

Fig. 4. A Wiener model consist of a linear time-invariant dynamic system with a unitary gain in series with a static non-linear element.

Fig. 5. Laboratory column and its instrument.

Rayleigh and Chi-square. The estimation problem consists of finding the PDF model parameters such that the probability that the data have been generated by the resulting PDF is maximized (maximum likelihood). Because the PDF is pre-defined, it may fail to adequately represent data in many cases.

Fig. 6. Frit-and sleeve sparger prototype (Kratch et al., 2008).

10

A. Riquelme et al. / International Journal of Mineral Processing 145 (2015) 7–16

Fig. 7. Examples of bubble detection results using CHT for small bubble clusters (left) and large bubble clusters (right).

Non-parametric methods do not assume any fixed function, which allows modelling any data set without restrictions. So far, histograms have been used for quick visualization of the BSD (Gomez and Finch, 2007). The critical parameter in this particular case is the bin width. When it is very small, the resulting density model is spiky. On the other hand, when the bin width is too large, the resulting model is too smooth. In both cases, an inappropriate bin choice could not represent the distribution correctly (Bishop, 2006). One important feature for this method is that the estimated density exhibits some discontinuities. This is an important disadvantage for control purposes. A semi-parametric method combines the advantages of both previous techniques and can comprise linear combinations of parametric functions. Maldonado et al. (2008) proposed, for the bubble size density estimation, the use of a Gaussian Mixture Models (GMM) whose parameters are estimated based on an Expectation Maximization (EM) algorithm (Titterington et al., 1985). (Maldonado et al. (2009) also suggested a control strategy formulated as an optimization problem which is based on GMM. However, the proposed cost function minimization algorithm does not guarantee a global minimum and thus the use of other methods should be investigated (i.e. free-derivative optimization algorithms such as genetic algorithms). According to previous works on flotation processes, bubble size distributions can be reasonably well fitted with log–normal distributions (Grau and Heiskanen, 2005; Majumder et al., 2006). In this paper, log– normal functions are computed by the EM method for BSD modelling purposes. Process dynamics can often be adequately modelled with linear time invariant system around an operating point. However, this strategy may

have limitations in describing non-linear processes. The computational complexity of modelling non-linear processes can be reduced by using some specific model structures. Particularly, some non-linear models are composed of linear dynamic subsystem and memoryless nonlinear functions such as Wiener systems (Park and Lee, 2004; Iqbal and Aziz, 2011; Gomez and Baeyens, 2005; Shiotani and Kobayashi, 2009), Hammerstein systems (Sung, 2002; Gomez and Baeyens, 2005), and their variations (Ogunfunmi, 2007). Other non-linear approaches can be found in various articles (Haber and Unbehauen, 1990; Nelles, 2001). A Wiener model is a convenient technique for representing any time-invariant system with fading memory (Boyd and Chua, 1985). Specifically it is a simple choice for a multi-input multi-output (MIMO) system where the non-linear element depends on more than one output variable of the linear element. Since it appears from experimental data that only the static gains of the BSD model are significantly non-linear, a Wiener model is a good structure for this system. This work presents a technique based on image processing to analyse pictures of bubbles taken from a laboratory flotation column operating with a two-phase (water–air) system. From the measurements, BSD is estimated with a parametric method (log–normal distribution). Then, a non-linear Wiener process is estimated to model the dynamics of BSD. Section 2 describes the algorithm to detect and measure the BSD using the Hough transform. Section 3 details the process for estimating a parametric representation for the BSD. In Section 4, the Wiener model identification steps are explained. Section 5 presents the laboratory setup. Finally, Section 6 shows the experimental results for bubble size measurement, BSD parametrization and Wiener model identification.

Fig. 8. Detection of cluster bubbles with a) circular particle detection, and b) circular Hough algorithm.

A. Riquelme et al. / International Journal of Mineral Processing 145 (2015) 7–16

2. Bubble detection algorithm This section presents the fundamentals of an image processing application for bubble detection based on CHT. More details about the bubble detection algorithm can be found in the paper authored by Riquelme et al. (2013). CHT offers better performances than circular particle detection (CPD) since it is able to detect not perfectly circular shapes. The images used are 800 × 600 grayscale pictures. Fig. 1 shows a raw sample image which should undergo a series of preconditioning steps before applying CHT: filtering, background subtraction, contrast enhancement and image inversion. Fig. 2a is a zoom of the sample image, with overlapped bubbles (clusters), to better see the results of each processing step. 2.1. Image preconditioning 2.1.1. Median filter The first stage of the process consists in applying a 10 × 10 median filter (Russ, 1998). This filter is a non-linear operation. It reduces salt and pepper noise while preserving the edges of bubbles. In the present case, the noise sources are light reflection on bubbles and in the background. Fig. 2b shows the result of applying a median filter. 2.1.2. Background subtraction This operation creates a new image by applying an erosion operator, followed by a dilatation using a zone of 5 × 5 pixel disc into the previously filtered image (Russ, 1998). This step removes snowflakes having a radius less than five pixels. The resulting image is later subtracted from the previously obtained image. This step allows extracting the gradient of the background produced by light diffusion and imperfections. Fig. 2c is the resulting image after elimination of the background. 2.1.3. Contrast enhancement The goal is to enhance bubble contour lines. The contrast is improved by an adaptive histogram equalization step (Russ, 1998). In addition, the process also highlights bright inner parts. This generates a circle in the centre of each bubble, which must be later eliminated. Fig. 2d is the resulting contrasted image. 2.1.4. Filling holes As explained above, a bright circle is produced in the centre of each bubble as if they behave like a mirror reflecting the light source. This circle is filled with a grey surrounded tone before inverting the greyscale image. A hole is defined as a set of background pixels which cannot be reached by filling in the background from the edge of the image after applying a flood-fill operation (Soille, 2003). Fig. 2e shows the resulting image after filling the holes inside the circles.

Fig. 9. Cumulative functions obtained from visual (manual) count, CHT and CPD.

11

Table 1 Comparison for each method. Method

Bubble detected

D32 (pixels)

D32 error (%)

Visual CHT CPD

537 525 377

33 32 27

– 3 18

2.2. Edge-map construction and Hough transform Circular Hough transform is a variation of the Hough transform aiming at detecting partial or complete circular shapes. CHT transforms the edged image into a set of accumulated values in a parameter space. The basic idea is described in (Kimme et al., 1975). Centred on each point of the edge, a circle with a target radius is drawn. Fig. 3 shows three of these points and the corresponding circles. At the coordinates where the perimeter of such a circle passes, the value of an accumulator matrix is incremented. In this example, the elements of the accumulator matrix would be 0 everywhere except on the three circles where it would take the value of 1, 2 (where 2 circles cross) or 3 (where 3 circles cross). In this example, as the edge corresponds to the same shape as the selected circle, a large value of the accumulator matrix would appear at the centre of the edge if several points around the edge were used. It would then be concluded that the edge corresponds to a circle having the radius of the selected circle. The accumulator matrix is calculated for circles with different radius. Local maxima are then searched in the accumulator array, indicating the presence of the centre of a circle for the given radius. The threshold for finding this maximum can be tuned to detect less perfect circles. A modified version of CHT, described in (Peng et al., 2007), is based on the gradient field of an image to compute the accumulator. The maximum intensity in the accumulation matrix represents the centre position of a sphere (Fig. 2f). The detection of peaks in the accumulation array is done in three steps: 1) a Laplacian of Gaussian (LoG) filter to fill the transition between local peaks and to make the peak separation easier, 2) a threshold to eliminate all the values which are not local maxima (Figs. 2g), and 3) a region-growing algorithm to detect the centroid of the isolated region, to determine the centre of each circle (Fig. 2h). Detected circles are accumulated until the process is completed. The final results are the position and the radius of each detected circle.

3. BSD parametrization Previous researchers have suggested that bubble size distributions in flotation processes can be reasonably well fit with log–normal

Fig. 10. Histogram of BSD and estimated log–normal distribution.

12

A. Riquelme et al. / International Journal of Mineral Processing 145 (2015) 7–16

Table 2 Log–normal estimation (10 ppm F150).

Table 3 Log–normal estimation (JG =1 cm/s).

JG (cm/s)

0.4

1

1.6

F150 (ppm)

10

20

30

Bubbles counted D32 from Eq. (3) (pixels) D32 from CHT data (pixels) D32 error (%) μ^ Std. Dev. of μ^ ^ σ ^ Std. Dev. of σ

4242 17.84 17.98 0.8 2.7334 0.003734 0.2432 0.002641

3356 21.42 21.5 0.4 2.9306 0.003991 0.2312 0.002823

2785 25.19 25.31 0.5 3.015 0.005514 0.2910 0.003901

Bubbles counted D32 from Eq. (3) (pixels) D32 from CHT data (pixels) D32 error (%) μ^ Std. Dev. of μ^ ^ σ ^ Std. Dev. of σ

3356 21.42 21.51 0.4 2.9306 0.003991 0.2312 0.002823

7966 15.55 16.3 4.6 2.617 0.002529 0.2258 0.001788

8700 13.84 14.27 3.0 2.556 0.001819 0.1697 0.001286

distributions (Grau and Heiskanen, 2005; Majumder et al., 2006). A log– normal probability density function is defined as: P ðxÞ ¼

2 2 1 pffiffiffiffiffiffi e−ð ln ½x−μ Þ =ð2σ Þ σ 2π x

ð2Þ

function, it does not modify the optimal solution. The log-likelihood is defined as: ln ðLL ðθ; x1 ; …; xn ÞÞ ¼ −

n X

n X ln X i þ ½ f ð ln ðxi ÞjθÞ:

i¼1

ð6Þ

i¼1

where σ is the standard deviation and μ is the mean. The Sauter mean diameter can be calculated from the log–normal distribution using directly the first and second moments of the function as (Vinnett et al., 2012):

The estimation of the expected log-likelihood of a single observation in the model is:

D32 ¼ eμþ5=2σ :

The standard approach to estimate the parameters of a defined PDF (θ), is the MLE, which requires maximization of the following loglikelihood function:

2

ð3Þ

The maximum likelihood estimation (MLE) method was used to estimate the log–normal distribution best fitting the measured BSD (Dempster et al., 1977). The MLE aims at finding the values of the log– normal distribution parameters (mean and variance) which makes the modelled data most likely (maximizing the likelihood function). Lets define a set of N independent observations xn coming from a distribution with an unknown PDF f, where θ is the vector of parameters (σ 2 and μ). The likelihood function is defined as: n

Lðθ; x1 ; …; xn Þ ¼ ∏ f ðxi jθÞ

"

#) ̂

fθmle g⊆ arg max l L ðθjx1 ; …; xn Þ θ∈Θ

:

ð8Þ

Taking derivatives with respect to σ 2 and μ, and solving the resulting system of first order conditions, the MLE for the parameters are: μ̂ ¼

" n X

# ln xi

ð9Þ

n

i¼1

and.

This estimation is the fundamental likelihood for a Gaussian distribution. To estimate a log–normal distribution, the following transform is made: n

ð7Þ

(

ð4Þ

i¼1

LL ðθ; x1 ; …; xn Þ ¼ ∏ ½ð1=xi Þf ð ln ðxi ÞjθÞ:

̂

l L ¼ ð1=nÞ ln LL :

ð5Þ

i¼1

" 2

σ̂ ¼

n X

ln ðxi −μ Þ2

# n:

ð10Þ

i¼1

4. Identification of a non-linear model

It is mathematically convenient to work with the logarithm of the likelihood function, called the log-likelihood. Since this is a monotone

Developing an accurate model is a fundamental step toward the design of a successful control strategy. Many types of non-linear models are available. A Wiener model is one of the simplest non-linear representation. It consists of a linear time-invariant (LTI) dynamic system

Fig. 11. Log–normal estimations for different flow rates at the same frother concentration.

Fig. 12. Log–normal estimated for different frother concentrations.

A. Riquelme et al. / International Journal of Mineral Processing 145 (2015) 7–16

13

Table 4 Estimated transfer function parameters.

Fig. 13. Wiener model to explain the dynamic behaviour.

with a unitary gain in series with a static non-linear (memory-less) element (Fig. 4). Usually the non-linear element is selected as an invertible function for eventual control design purposes. Many techniques have been suggested to identify Wiener models. Iqbal and Aziz (2011) compare different methods such as the blind approach (Bai, 2002), the frequency domain approach (Bai, 2003), the support vector approach (Totterman and Toivonen, 2002) and the maximum-likelihood approach (Hagenblad et al., 2008). There are three possible ways to identify a Wiener model (Cervantes et al., 2003): • Simultaneous: in this case, both elements (LTI system and non-linear gain) are estimated at the same time. • Linear–Non-linear (L–N): the parameters of the linear transfer function are first identified. An intermediate signal is then generated to calculate the non-linear static element. • Non-linear–Linear (N–L): the static non-linear parameters are first identified from the steady-state data. Then, the linear dynamic parameters are obtained from dynamic data. This last approach will be further described as it will be selected.

The steps for N–L identification are the following: 1. Collect a set of steady-state data. Assuming that the LTI system has a gain equal to 1, the static non-linearity can be identified from this set of data. This non-linear function must be invertible for the subsequent estimation of v(t) (Fig. 4). 2. Collect a set of dynamic data. Once the non-linearity has been identified and assuming that the non-linear function is invertible, it is possible to calculate v(t) (Fig. 4) by applying the inverse of the non-linearity. Then, the linear system can be identified using the signals u(t) and v(t). 5. Experimental set up A schematic representation of the experimental set-up is shown in Fig. 5. The laboratory flotation column is made of 5.08 cm diameter

Transfer function

τ (sec.)

θ (sec.)

MGμ(s) MGσ(s) MLμ(s) MLσ(s)

25 ± 2.6 15 ± 3.4 61 ± 3.6 43 ± 3.4

57 ± 1.9 53 ± 2.5 117 ± 2.6 113 ± 2

acrylic tubes for a total height of 7 m. Two peristaltic pumps inject the feed at the middle of the column and extract the tailings at the bottom flow. Tailings and concentrate streams are recirculated to a common reservoir containing the feed solution made of tap water and F150 (a frothing agent commonly used in flotation to stabilize the froth). Air bubbles are generated with a stainless-steel frit-and-sleeve sparger located at the bottom of the column. The frit-and-sleeve sparger (Fig. 6), consists of a porous cylinder surrounded by a sleeve forming a gap through which water is injected, shearing the air injected through the sparger into tiny bubbles (Kratch et al., 2008). This sparger allows the modification of the bubble size by manipulating the shear-water flow rate through the gap. This means that bubble size can be varied independently of the air flow rate. i.e. it provides an extra degree of freedom for controlling the bubble size distribution. The gas flow rate is controlled by a mass flow controller and the shearing water flow is regulated by manipulating the speed of a gear pump with a PI controller. Bubble images are taken with a modified version of the McGill bubble viewer (Gomez and Finch, 2007). The bubble viewer is connected to the column by a 1.27 cm (0.5 in. sampling tube at 2 m from the top. Bubbles are returned to the column after having been exposed in the viewing chamber and photographed with a computer controlled charge-coupled device (CCD) camera. The camera is scaled at approximately 20 pixels/mm. The viewing chamber is uniformly illuminated with a diffused white-light source positioned opposite to the camera and aligned to this latter to have uniform and clear bubble contours, avoiding shaded areas. A computer running with Matlab (R2010a) was used to take the pictures, analyse the images and estimate the log–normal PDF. 6. Experimental results 6.1. Bubble size measurement Tests were carried out with different aeration conditions to generate a wide range of bubble sizes and shapes. Fig. 7 shows two cases where most of the bubbles in the foreground are detected by applying CHT, irrespective of their size (including large ones), state of clustering, and

Fig. 14. Surface estimation results. Left: μ estimation, right: σ estimation.

14

A. Riquelme et al. / International Journal of Mineral Processing 145 (2015) 7–16

bubbles larger than 30 pixels in diameter. With visual counting, a total of 537 bubbles were identified. Table 1 summarizes the results for each method. Clearly, CHT is superior to the circular shape analysis algorithm, with a total number of detected bubbles comparable to the visual count. The average diameter (D32) was also calculated, its value being more accurate for CHT method, with a relative error of 3% with respect to the visual counting.

Fig. 15. Monovariable validation test (top graph: model output in blue and measured value in red).

completion (e.g. at the boarder of the frame). Bubbles in the background do not have a sharp enough boarder, thus they are not detected. Thus detection only applies for the first layer of bubbles in the viewing chamber. The processing time for each picture is around two seconds.

6.1.1. Comparison between different measuring methods Three different methods for detecting bubbles were used: CHT, CPD (with ImageJ software), and visual counting. Fig. 8 compares the results obtained with (a) CPD (black circles: detected, fuzzy circles: undetected) and (b) CHT (detected bubbles are those with blue circles and a red cross marking their centre). It quite evident that CHT is superior in detecting clustered bubbles. Detected shapes from the use CHT can only be circles. Therefore, elliptical bubbles are detected as circles fitting within the ellipses. To compare the overall performance of CHT, a set of 25 pictures was analysed using the three methods mentioned above, and their cumulative function was estimated. Fig. 9 shows the results. The cumulative function obtained with CHT is very similar to what was obtained with visual counting. The performance of the CPD algorithm decreases for

6.1.2. Parametrization of bubble size distribution Three tests were carried out to validate the performance of the BSD parametrization. The objective of the first test is to validate the choice of a log–normal function by comparing with the histogram of the resulting measurements. In the second series of tests, the air flow rate was varied while keeping the frother concentration constant. Log–normal parameters were estimated and D32 was computed for comparison purposes. In the last tests, the air flow rate was kept constant while the frother concentration was changed. D32 and log–normal parameters were computed and compared for validation. Test 1. A set of 50 pictures/experiment were taken in steady state conditions. Each experiment corresponds to different frother concentrations and superficial gas rates. The maximum likelihood estimation algorithm was used to estimate the log–normal function with the 50 pictures. Fig. 10 is a typical result. It super-imposes the calibrated function over the histogram for an experiment done with a solution containing 10 ppm of F150 and a superficial gas velocity of 0.4 cm/s. The histogram was computed with equal bins of 1 pixel width each. Test 2. Different air flow rates were injected into the column through the sparger. Frother concentration was kept constant at 10 ppm for all tests. The methodology consisted in collecting 50 bubble pictures in steady-state operating conditions. Table 2 summarizes the experimental results, showing the total bubble count, the estimated D32 (calculated with Eq. (3)), the reference D32 estimated with the histogram (raw data from CHT detection method), the percentage of error with respect to the reference D32, and the estimated log–normal function parameters μ and σ. It is quite evident that the error of the D32 estimation is very low for all three cases. Fig. 11 presents the estimated log–normal function for the three superficial air velocities at the same frother concentration. Test 3. In these experiments, three different frother concentrations were used (10, 20 and 30 ppm), while keeping constant the superficial air velocity ( JG = 1 cm/s). Table 3 shows the results. D32 is correctly estimated by the log–normal model. Fig. 12 presents the three estimated log–normal functions. Not surprisingly, increasing frother concentration reduces the mean diameter and the distribution gets narrower.

Fig. 16. Multivariable validation tests (two top graphs: model output in blue and measured value in red).

A. Riquelme et al. / International Journal of Mineral Processing 145 (2015) 7–16

6.2. Wiener model identification Fig. 13 is the proposed Wiener model to explain the dynamic behaviour of the BSD parameters. JGs (cm/s) is the superficial velocity setpoint of the air injected into through the sparger, and JLs (cm/s) is the superficial velocity set-point of the shearing water. In steady state, JG and JL are respectively equal to JGs and JLs. The outputs μ and σ are the parameters that characterize the log–normal function representing the bubble size distribution in the column. MGμ, MGσ, MLμ and MLσ are the transfer functions for each output and input combination (all first−θs

order with dead time and unitary static gain: 1e 1þτs ), while fμ(vGμ, vLμ) and fσ(vGσ, vLσ) are the non-linear elements for each output. In order to estimate the Wiener model, the N–L method is used. As described in Section 4, the identification process is separated in two stages: determination of the non-linear elements and estimation of the transfer functions. 6.2.1. BSD static non-linearities In order to obtain the non-linear model, a set of data must be collected under steady-state conditions. Several experimental tests were conducted in the laboratory flotation column. A constant concentration of 10 ppm of frother was used in all tests. JLs was varied in steps of 0.2 cm/s, between 0 and 1, and JGs was modified in steps of 0.2 cm/s between 0.5 and 1.5 cm/s. A total of 36 steady-state points were collected. Fig. 14 depicts the results for μ and for σ. After the experiments were carried out, the static models were estimated using a non-linear least squares algorithm. Various mathematical models were tested. For μ, the estimated function is:   f μ vGμ ; vLμ ¼ 2:815 þ 0:1366vGμ −0:1245vLμ

ð12Þ

with R2 = 98.5%. The resulting fitting is shown in Fig. 14. Applying the natural logarithm to Eq. (12) leads to: ln ð f σ ðvGσ ; vLσ Þ−0:1975Þ ¼ −2:688 þ 1:105vGσ −2:011vLσ

imperfect circles, and bubbles at the borders of the frame). CHT leads to results consistent with a visual count, for the Sauter mean diameter, within a 3% error. The processing time is comparable to that of the commonly used circular shape based method. It must be emphasized that CHT requires automated pre-conditioning stages. A pre-tuning of CHT is thus necessary to determine the threshold to find the local maxima of the cumulative function in order to detect less-circular shapes. Fitting the BSD by a log–normal distribution leads to good results as suggested by previous works. Given the characteristics of the system, a non-linear relationship is found between the manipulated variables (JGs and JLs) and the log–normal parameters (μ and σ). The BSD dynamics are modelled by a Wiener structure. The good agreement between the response of the process and the response of the identified model confirms this choice. Since the non-linear elements are invertible, a linear control method could then be easily designed. This will be the objective of the future work. It would then be possible to select the desired bubble size distribution that optimizes the concentrate recovery for a given particle size or particle size distribution. Another challenge is the industrial application of the bubble size measurement algorithm proposed in this research. A complete study with a three-phase system (bubbles, water and particles) needs to be done. The bubble viewer will have to be adapted for continuously work with solids such as implementing an automated system to wash the viewing chamber, as proposed by Roy (2013).

References

ð11Þ

with R2 = 93.3%. The estimated function for σ is: f σ ðvGσ ; vLσ Þ ¼ 0:1975 þ 0:068eð1:105vGσ −2:011vLσ Þ

15

ð13Þ

which is a linear function. 6.2.2. BSD linear dynamics Once the non-linearities were calculated, a series of tests for the identification of the transfer functions were carried out (as detailed in Section 4). One input was kept constant, while the other was varied in steps. In the first test, JGs was equal 0.5 cm/s while steps were applied to JLs. Two similar experiments were recorded with JGs = 0.9 cm/s and JGs = 1.3 cm/s. Three other data sets were collected by keeping constant JLs (0.0 cm/s, 0.4 cm/s and 0.9 cm/s) and varying JGs. The transfer function parameters (τ and θ) were then estimated using a maximum likelihood estimation algorithm. The results are given in Table 4. Fig. 15 shows the comparison between the measured value of μ and its value predicted by the model, using validation data (i.e data that were not used for identification). Another validation test is illustrated in Fig. 16 where both inputs are manipulated. These validation tests confirm that the proposed model is able to adequately predict the dynamic behaviour of the BSD parameters. A controller could certainly be designed based on this model. 7. Conclusions and future work A new application of CHT for detecting bubbles in flotation columns is developed in this work. This method, based on the circular Hough transform, allows the detection of all sorts of bubbles (clusters,

Bai, E.-W., 2002. A blind approach to the Hamerstein–Wiener model identification. Automatica 38, 967–979. Bai, E.-W., 2003. Frequency domain identification of Wiener models. Automatica 39, 1521–1530. Bishop, C., 2006. Pattern Recognition and Machine Learning. Springer-Verlag. Boyd, S., Chua, L., 1985. Fading memory and the problem of approximating nonlinear operators with Volterra series. IEEE Trans. Circuits Syst. 320 (11), 1150–1161. Cervantes, A., Agamennoni, O., Figueroa, J., 2003. A nonlinear model predictive control system based on Wiener piecewise linear models. J. Process Control 13, 655–666. http://dx.doi.org/10.1016/S0959-1524(02)00,121-X. Daraei, M., Rezaie, J., Norouzi, N., Khajooeizadeh, A., 2011. An innovative implementation of circular Hough transform using eigenvalues of covariance matrix for detecting circles. Int. Symp. Electronic. Mar. 1, 397–400. Dempster, A.P., Laird, N.M., Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. 390 (1), 1–38. Duda, R., Hart, E., 1972. Use of the Hough transformationto detect lines and curves in pictures. Commun. ACM 1, 11–15. http://dx.doi.org/10.1145/361,237.361242. Gomez, J., Baeyens, E., 2005. Subspace identification of multivariable Hammerstein and Wiener models. Eur. J. Control. 11, 1–10. Gomez, C.O., Finch, J.A., 2007. Gas dispersion measurements in flotation cells. Int. J. Miner. Process. 84, 51–58. http://dx.doi.org/10.1016/j.minpro.2007.03.009. Goulermas, J., Liatsis, P., 2012. Novel combinatorial probabilistic Hough transform technique for detection of underwater bubbles. Proceedings of Machine Vision Applications in Industrial Inspection 3029, pp. 147–156. http://dx.doi.org/10.1117/12. 271237. Grau, R., Heiskanen, K., 2005. Bubble size distribution in laboratory scale flotation cells. Miner. Eng. 180 (12), 1164–1172. http://dx.doi.org/10.1016/j.mineng.2005.06.011. Haber, R., Unbehauen, H., 1990. Structure identification of nonlinear dynamic system, a survey on input/output approaches. Automatica 260 (4), 651–677. Hagenblad, A., Ljung, L., Wills, A., 2008. Maximum likelihood identification of Wiener models. Automatica 44, 2697–2705. Han, M., Park, Y., Yu, T., 2002. Development of a new method of measuring bubble size. Water Sci. Technol. Water Supply 20 (2), 77–83. Heiskanen, K., 2000. On the relationship between flotation rate and bubble surface area flux. Miner. Eng. 130 (2), 141–149. http://dx.doi.org/10.1016/S0892-6875(99)00,160–0. Iqbal, I., Aziz, N., 2011. Comparison of various Wiener model identification approach in modeling nonlinear process. 3rd Conference on Data Mining and Optimization 1, pp. 134–140. http://dx.doi.org/10.1109/DMO.2011.5976517. Kimme, C., Ballard, D., Sklansky, J., 1975. Finding circles by an array of accumulators. Commun. ACM 180 (2), 120–122. Kratch, W., Gomez, C., Finch, J., 2008. Controlling bubble size using a frit-and-sleeve sparger. Miner. Eng. 210 (9), 660–663. Lin, B., Recke, B., Knudsen, J., Jorgensen, S., 2007. Bubble size estimation for flotation processes. Miner. Eng. 210 (7), 539–548. http://dx.doi.org/10.1016/j.mineng.2007.11.004. Majumder, S., Kundu, G., Mukherjee, D., 2006. Bubble size distribution and gas–liquid interfacial area in a modified downflow bubble column. Chem. Eng. J. 122, 1–10. http:// dx.doi.org/10.1016/j.cej.2006.04.007.

16

A. Riquelme et al. / International Journal of Mineral Processing 145 (2015) 7–16

Maldonado, M., Desbiens, A., del Villar, R., Girgin, E., Gomez, C., 2008. On-line estimation of bubble size distributions using Gaussian mixture models. Proceedings of 5th International Mineral Processing Seminar: Procemin 2008 1, pp. 389–398. Maldonado, M., Desbiens, A., Del Villar, R., 2009. Potential use of model predictive control for optimizing the column flotation process. Int. J. Miner. Process. 93, 26–33. Nelles, O., 2001. Nonlinear System Identification: From Classical Approaches to Neural Networks and Fuzzy Models. Springer editors. Nostrati, M., Karimi, R., 2011. Detection of circular shapes from impulse noisy images using median and laplacian filter and circular Hough transform. International Conference on Electrical Engineering, Computing Science and Automatic Control 1, pp. 1–5. http://dx.doi.org/10.1109/ICEEE.2011.6106710. Ogunfunmi, T., 2007. AdaptiveNnonlinear System Identification: The Volterra and Wiener Model Approaches. Springer editors. Park, H., Lee, J., 2004. Simple approaches for the identification of Wiener-type nonlinear process. 5th Asian Control Conference 1, pp. 93–103. Peng, T., Balijepalli, A., Gupta, S., LeBrun, T., 2007. Algorithms for on-line monitoring of micro spheres in an optical tweezers-based assembly cell. J. Comput. Inf. Sci. Eng. 70 (4), 330–339. http://dx.doi.org/10.1115/1.2795306. Riquelme, A., Bouchard, J., Desbiens, A., del Villar, R., 2013. Bubble detection in flotation columns based on circular Hough transform. World Min. Congr. 1, 2–10. Rodrigues, R., Rubio, J., 2003. New basis for measuring the size distribution of bubbles. Miner. Eng. 16, 757–765. http://dx.doi.org/10.1016/S0892-6875(03)00,181-X.

Roy, J., 2013. Investigation sur l'impact de la concentration de moussant dans un système solide-liquide-gaz et automatisation d'un capteur de dimension des bulles (Master's thesis) Université Laval. Russ, J., 1998. The Image Processing Handbook. 3rd ed. CRC press. Shiotani, Y., Kobayashi, 2009. Identification of multi-input multi-output Wiener-type nonlinear system. ICCAS-SICE International Joint Conference 1, pp. 5244–5249. Soille, P., 2003. Morphological Image Analysis: Principles and Applications. SpringerVelag, NY. Sung, S., 2002. System identification method for Hammerstein processes. Ind. Eng. Chem. Res. 410 (17), 4295–4302. Titterington, D.M., Smith, A.F., Makov, U.E., 1985. Statistical Analysis of Finite Mixture Distributions. John Wiley & Sons Inc., UK. Totterman, S., Toivonen, H., 2002. Support vector method for identification of wiener models. J. Process Control 19, 1174–1181. Vinnett, L., Contreras, F., Yianatos, J., 2012. Gas dispersion pattern in mechanical flotation cells. Miner. Eng. 26, 80–85. http://dx.doi.org/10.1016/j.mineng.2011.11.003. Yianatos, J., Contreras, F., 2010. On the carrying capacity limitation in large flotation cells. Can. Metall. Q. 40 (8), 345–352. http://dx.doi.org/10.1179/000,844,310,795,937,578. Zhu, Y., Wu, Y., Manasseh, R., 2001. Rapid measurement of bubble size in gas-liquid flows using a bubble detection technique. Australasian Fluid Mechanics Conference 1, pp. 541–544.