Inversion method of particle size distribution of milk fat based on improved MPGA

Milk fat’s particle size and distribution not only affect product quality, but also have great impacts on food safety in the economy and society. Based on total light scattering method, this paper has studied the inversion method of particle size distribution under dependent mode condition by combining multi-population genetic algorithm (MPGA) with Tikhonov smooth function. It has minimized the influence from light-absorb medium to improve the inversion accuracy. The approach introduces Tikhonov smooth function and apparent optical parameters to build an objective fitness function and weaken the ill condition of the particle size inversion equation. It also introduces multi-population genetic algorithm to solve the premature convergence of genetic algorithms. The results show that the relative error of the milk fat simulation solution with a nominal diameter is -3.52%, which meets the national standard of ±8% and better than the relative error of -5.01% of the standard genetic algorithm. Thus, the improved MPGA can reconstruct particle size distribution, with a good reliability and stability.

Introduction

Particle size and distribution not only affect product performance and quality, but also influence the economy and society in the aspects of environmental protection and human health (Cai et al., 2010). In recent years, the demands for online monitoring of particulate matter have increased, especially in environmental testing and food safety. More scholars from global leading institutions focus on particle size inversion methods.

According to Wang Li’s finding, the method of using regularization to build objective function and establish an improved mode searching model showed a good inversion performance. She also proposed using generalized eikonal approximation to effectively calculate the coefficient of light extinction for measuring the particle size (Wang, 2013; Wang and Sun, 2013). Based on angular scattering method, Mroczka et al. measured particle size distribution using regularization least squares algorithm with constraints. Although they obtained a good inversion result, the team failed to reconstruct the bimodal distribution function (Mroczka and Szczuczyński, 2013). Clementi et al. reversed polystyrene latex particles with different shape, size distribution width and diameter range. They used multi-angle dynamic light scattering and capered Tikhonov and Bayesian methods to solve the ill non-linear problem. The paper concluded that Bayesian method produced better results than Tikhonov regularization (Igushi and Yoshida, 2011). Compared with classical methods, intelligent optimization algorithm showed better whole search ability in measuring particle size distribution and can efficiently overcome the effects from the noise. However, one algorithm cannot have all the advantages of other algorithms. It can only manage to reduce one or another disadvantage, improving the iteration speed or calculation accuracy. The detection of particulate matter across the world mainly focuses on inorganic molecules distributed in the air, rather than organic ones dissolved in the liquid media (He et al., 2018).

The following issues arise when particle matter in milk is detected: during the detection of organic molecules in a liquid medium, the medium may absorb light and affect the outcomes; in the total light scattering method, the first kind of Fredholm integral equation is a typical ill-posed problem; when calculating particle size, the algorithm exhibits delayed convergence, inaccurate inversion and poor noise resistance (Maguire et al., 2018).

This research suggests an enhanced MPGA based on the total light scattering approach to address these issues. It simplifies the objective function using Tikhonov smooth function (Wang and Sun, 2013) and reduces ill conditions of the inversion equation. Based on the extinction spectra calculated by known specific parameters, an inversion algorithm is built to generate a more precise particle size distribution. Apparent optical parameters (Fu and Sun, 2001; Yang et al., 2002; Fu and Sun, 2006; Ding et al., 2020; Wrana et al., 2020) are introduced based on total light scattering to solve the extinction coefficient and create an objective function. This effectively overcomes the influence of the light-absorbing medium. The study also achieves the particle size inversion of milk fat in medium of liquid absorption, addresses the typical ill-posed problem in the total light scattering method, and provides a reliable method for determining particle size of milk fat (Lienert et al., 2001; Wang et al., 2021; Zhou et al., 2021; Cai et al., 2022).

Materials and methods

Experimental materials

The size of the milk fat particles in this study, which satisfied the unimodal R-R distribution function, was 0.1 $\mu m$ ∼10 $\mu m$. Ethanolamine, of which the refractive index was $m=1.454+1.59\times {10}^{-5}i$ and could be classified as liquid-absorbing medium, was selected as the medium to dissolve milk fat particles. Anhydrous milk fat was dissolved in the ethanolamine solvent at a concentration of 6% to create milk fat simulated liquid. The liquid was then thoroughly churned using a JJ-1 precision timed mixer while being maintained at a constant 30 °C temperature. The prepared simulated liquid was then measured with a Malvern particle size analyzer. According to the analyzer, the size of the measured particle system was 5.075 $\mu m$ (Ding et al., 2020).

Design of optical experimental system

Figure 1 shows the schematic diagram of the optical experimental system. As can be seen from Figure 1, the whole system included photoelectric sensor, light source, data receiving unit and sample pool. Lasers of 0.3838, 0.4598, 0.5163, 0.6269 and 0.77 $\mu m$ were selected as the light source (Ding et al., 2020). To improve accuracy, the study used quartz cuvette with low light consumption as the sample pool, the thickness of which was 12.33 mm. The light intensity of each wavelength was collected by the optical power detection device (optical power meter) (Ding et al., 2020).

FIGURE 1

FIGURE 1. Schematic diagram of particle optical experimental system.

Measurement method

First, the ethanolamine medium was added into the sample pool. The measuring area was irradiated with a laser at a preset wavelength. Then, the scattered light intensity was collected by a photoelectric sensor positioned at a 90° angle. The collected signal was input to a computer for processing and was set as the incident light intensity. Then the simulated liquid containing milk fat particles was added to the sample pool in place of the original liquid. Using the same method above, the signal value was recorded as the scattered light intensity. To avoid errors, the ratio of the first result and the second result was processed and calculated by the computer as the extinction value (Wang, 2013). The excitation experiment was repeated for 20 times on the simulated liquid of milk fat particles with lasers of different wavelengths, and the average of the results was taken as the final light intensity (Cai et al., 2010; Wang, 2013; Wang and Sun, 2013; Ding et al., 2020).

Basic concept and improvement of total light scattering method

The measurement of the particle size followed the Lambert-Beer law.

The monochromatic parallel lights’ incident intensity was $I_0$ and their wavelength was $\text{λ}$. When the lights were incident on the particle system, they was scattered and absorbed, and the transmitted intensity $I$ would attenuate. The change can be expressed as (Cai et al., 2010):

$-dI=I\tau dL(1)$

Where: $\tau$ was the turbidity of the particle medium and $L$ was the thickness of the medium to be measured. Assume that $\tau$ in the particle system was independent of the optical path L. In other words, the spatial distribution in the medium was disordered and uniform. Eq. 1 can be integrated along $L$:

$-{\displaystyle \underset{I_0}{\overset{I}{\int }}\frac{1}{I}dI}={\displaystyle \underset{0}{\overset{L}{\int }}\tau dL}(2)$

The relationship between the incident light intensity and the transmitted light intensity was changed into:

$ln\frac{I\left(\lambda \right)}{I_0\left(\lambda \right)}=-\tau \left(\lambda \right)L(3)$

Where: $I\left(\lambda \right)/I_0\left(\lambda \right)$ was the extinction value with a wavelength of $\lambda$, $\tau \left(\lambda \right)$ was the turbidity with a wavelength of $\lambda$.

When the to-be-measured particles satisfy the condition of irrelevant single scattering and were all spherical particles with the size of D. In the single scattering system, the $N_D$ single particle’s extinction value was:

$ln\frac{I\left(\lambda \right)}{I_0\left(\lambda \right)}=-\frac{\pi }{4}L{N}_D{D}^2{Q}_{ext}\left(\lambda ,m,D\right)(4)$

Where: $Q_{ext}\left(\lambda ,m,D\right)$ was the extinction coefficient, in which $m$ was the complex refractive index of the relative medium. $Q_{ext}$ can be deduced based on Mie scattering theory: (Cai et al., 2010)

$Q_{ext}=\frac{2}{{\alpha }^2}{\displaystyle \sum _{n=1}^{\infty }\left(2n+1\right)}Re\left(a_n+b_n\right)(5)$

Where: $\text{α}=\pi D/\lambda$ was the size parameter of zero-dimensional particles; $a_n$ and $b_n$ were the Mie coefficients related to Bessel and Hankel functions.

$a_n=\frac{{\psi }_n^{\prime }\left(m\alpha \right){\psi }_n\left(\alpha \right)-m{\psi }_n\left(m\alpha \right){\psi }_n^{\prime }\left(\alpha \right)}{{\psi }_n^{\prime }\left(m\alpha \right){\xi }_n\left(\alpha \right)-m{\psi }_n\left(m\alpha \right){\xi }_n^{\prime }\left(\alpha \right)}(6)$

$b_n=\frac{m{\psi }_n^{\prime }\left(m\alpha \right){\psi }_n\left(\alpha \right)-{\psi }_n\left(m\alpha \right){\psi }_n^{\prime }\left(\alpha \right)}{m{\psi }_n^{\prime }\left(m\alpha \right){\xi }_n\left(\alpha \right)-{\psi }_n\left(m\alpha \right){\xi }_n^{\prime }\left(\alpha \right)}(7)$

${\psi }_n\left(\alpha \right)=\sqrt{\frac{\pi \alpha }{2}}{J}_{n+1/2}\left(\alpha \right)(8)$

${\xi }_n\left(\alpha \right)=\sqrt{\frac{\pi \alpha }{2}}{H}_{n+1/2}^2\left(\alpha \right)(9)$

In Eqs 8, 9, $J_{n+1/2}$ and $H_{n+1/2}$ were the Bessel function of half-integer order and the Hankel function of first kind, respectively. $m=m_r+m_{i}i$ was the complex refractive index of the particles relative to their surrounding medium. The sum of infinite series was considered when we calculated the extinction coefficients using Mie scattering theory. The sum upper limit $N_{stop}$ was usually derived by the empirical formula of Wiscombe (Wang, 2013; Wang and Sun, 2013):

$N_{stop}=\{\begin{array}{c}\alpha +4{\alpha }^{1/3}+1\mathrm{,}0.02\le \alpha \le 8\\ \alpha +4.05{\alpha }^{1/3}+2\mathrm{,}8\le \alpha \le 4200\\ \alpha +4{\alpha }^{1/3}+2\mathrm{,}4200\le \alpha \le 20000\end{array}(10)$

The calculation formula of the extinction coefficient $Q_{ext}$ in Mie scattering theory was:

$Q_{ext}=\frac{2}{{\alpha }^2}{\displaystyle \sum _{n=1}^{N_{stop}}\left(2n+1\right)}Re\left(a_n+b_n\right)(11)$

In measurement, the particle system was a dispersion with multiple scattered matters rather than a monodisperse system. Thus, when we had the wavelengths, the extinction value of spherical particles was calculated with the following equation:

$ln\frac{I\left(\lambda \right)}{I_0\left(\lambda \right)}=-\frac{\pi }{4}L{\displaystyle \underset{D_{min}}{\overset{D_{max}}{\int }}{Q}_{ext}}\left(\lambda ,m,D\right)N\left(D\right)D^{2}dD(12)$

Where: $D_{max}$ and $D_{min}$ were the maximum and minimum values that meet the particle size distribution in the system; $N\left(D\right)$ was the number of particles in $\left[D+dD\right]$. Volume distribution was used for calculation, because quantity distribution was relatively complex:

$ln\frac{I\left(\lambda \right)}{I_0\left(\lambda \right)}=-\frac{3}{2}L{N}_D{\displaystyle \underset{D_{min}}{\overset{D_{max}}{\int }}\frac{Q_{ext}\left(\lambda ,m,D\right)}{D}}f\left(D\right)dD(13)$

Where: $N_D$ was the total number of to-be-measured particles; $f\left(D\right)$ was the volume frequency distribution function of the particle system. This formula was the Fredholm integral equation of the first kind. As it is difficult to calculate directly, the equation can be further discretized, as shown in Eq. 14:

$ln\frac{I\left(\lambda \right)}{I_0\left(\lambda \right)}=-\frac{3}{2}L{N}_D{\displaystyle \sum _{j=1}^S{c}_j}\frac{Q_{ext}\left(\lambda ,m,\widetilde{D_j}\right)}{\widetilde{D_j}}f\left(D_j\right)(14)$

Where: $S$ denoted the grading number of particle size in $\left[D_{min},D_{max}\right]$; $c_j$ was the numerical integration coefficient; $f_j\left({\tilde{D}}_j\right)=\underset{j}{\overset{j+1}{{\displaystyle \int }}}f\left(D\right)d\left(D\right)$; and ${\tilde{D}}_j$ represented the midpoint of subinterval $\left[D_j,D_{j+1}\right]$. There were S types of particle size distribution to be solved, and they were measured using multiple wavelengths. As a result, the linear equation set was obtained (Wang, 2013; Wang and Sun, 2013):

$E=Af(15)$

Where: $E=\left[ln\left(I_1/I_{10}\right);\cdots ln\left(I_U/I_{U0}\right)\right]$; U was the wavelength; and $A=\left[A_{ij}\right]$ indicated the $U\times S$ weight matrix, with all the elements expressed as: $A_{ij}=-3C_{j}L{N}_D{Q}_{ext}\left({\lambda }_i,m,{\tilde{D}}_j\right)/\left(2{\tilde{D}}_J\right),i=1,\cdots ,U,j=1,\cdots ,S$; $f=\left[f_1\left({\tilde{D}}_1\right),\cdots ,f_S\left({\tilde{D}}_S\right)\right]$.

Inversion of particle size distribution by total light scattering method used algorithms to inverse particle size, with given extinction values of incident light wavelengths. Based on whether the distribution function was preset in data processing, the inversion can be divided into dependent algorithm and independent algorithm. After optical measurement of the simulated liquid of milk fat particle, the extinction value was obtained. Then, total light scattering method was used to inverse particle size. (The experiment and simulation results were obtained by dependent algorithm).

Dependent algorithm, also known as function restriction method, required presetting a distribution function during the calculation. Its concept was simple. The measured extinction values of multiple incident light wavelengths need to be substituted into the preset distribution function to solve the characteristic parameters of the function. This method was widely used for its simple steps. The specific steps in this paper were as follows:

The distribution function of the to-be-measured particle system was preset as follows:

$f\left(D\right)=f\left(D,D_1,D_2\right)(16)$

Where: $D_1$ and $D_2$ were the characteristic parameters of the distribution function $f\left(D\right)$.

In inversion, Eq. 16 was substituted into Eq. 14 by referring to extinction values:

$ln\frac{I{\left(\lambda \right)}_i}{I_0{\left(\lambda \right)}_i}=-\frac{3}{2}L{N}_D{\displaystyle \sum _{j=1}^S{c}_j}\frac{Q_{ext}\left({\lambda }_i,m,{\tilde{D}}_j\right)}{{\tilde{D}}_j}{f}_j\left({\tilde{D}}_j,D_1,D_2\right)(17)$

Where: $i$ indicated the measured wavelength ${\lambda }_i$. The objective function was:

$\text{OBE}={\left|\frac{ln{\left(I/I_0\right)}_i}{ln{\left(I/I_0\right)}_k}-\frac{-\frac{3}{2}L{N}_D{\displaystyle \sum _{j=1}^S{c}_j}\frac{Q_{ext}\left({\lambda }_i,m,{\tilde{D}}_j\right)}{{\tilde{D}}_j}{f}_j\left({\tilde{D}}_j,D_1,D_2\right)}{-\frac{3}{2}L{N}_D{\displaystyle \sum _{j=1}^S{c}_j}\frac{Q_{ext}\left({\lambda }_k,m,{\tilde{D}}_j\right)}{{\tilde{D}}_j}{f}_j\left({\tilde{D}}_j,D_1,D_2\right)}\right|}^2\mathrm{,}i\ne k(18)$

Thus, the problem of inversing to-be-measured particle size distribution transformed into solving parameters $D_1$ and $D_2$ corresponding to the minimum objective function. By substituting parameters $D_1$ and $D_2$ into Eq. 19, the Sauter average particle size was obtained and compared with the actual particle size (Wang, 2013; Wang and Sun, 2013).

$D_{32}=\frac{{\displaystyle \underset{D_{min}}{\overset{D_{max}}{\int }}{D}^3}f\left(D\right)dD}{{\displaystyle \underset{D_{min}}{\overset{D_{max}}{\int }}{D}^2}f\left(D\right)dD}(19)$

Where: $f\left(D\right)$ was the distribution function of the particle size; $D$ denoted the single particle size; $D_{max}$ and $D_{min}$ represented the maximum and minimum particle sizes in the particle system.

Based on the method of total light scattering, Tikhonov regularization was introduced into the objective function. The equation was:

$min=\Vert E-Af{\Vert }^2+\alpha \Vert f{\Vert }^2(20)$

Where: $\alpha$ denoted the regularization parameter which adjusted the relative weight between the residual $E-A{f}^2$ and the regularized $f^2$. The Tikhonov smooth objective function was:

$\phi \left(f\right)=min\left\{\begin{array}{c}{\displaystyle \sum }_{i=1}^s{\left\{ln{\left(\frac{I}{I_0}\right)}_{{\lambda }_i}-\left[-\frac{3}{2}L{N}_D{\displaystyle \sum _{j=1}^N{c}_j}\frac{Q_{ext}\left({\lambda }_i,m,D_j\right)}{D_j}f\left(D_j\right)\right]\right\}}^2\\ +\alpha {\displaystyle \sum _{j=1}^N{\left[f\left(D_j\right)\right]}^2}\end{array}\right\}(21)$

Where: when calculating the extinction value, the apparent optical scattering principle was used, so light absorption character of the liquid medium should be considered. The regularization parameter was solved by the L curve method (Yang, 2016; Zhou et al., 20212021), in which the maximum curvature point was the most appropriate regularization parameter. Let $\rho =lgb-A{x}_2^{\alpha }$, $\theta =lg{x}_2^{\alpha }$, then the curvature can be defined as (Yang, 2016; Wang et al., 2019; Zhou et al., 20212021):

$K\left(\alpha \right)=\frac{{\rho }^{\prime }{\theta }^{″}-{\rho }^{″}{\theta }^{\prime }}{{\left[{\left({\rho }^{\prime }\right)}^2+{\left({\theta }^{\prime }\right)}^2\right]}^{\frac{3}{2}}}(22)$

Where: The $\text{α}$, which corresponded to the maximum value of $K\left(\alpha \right)$, was the desired regularization parameter. Thus, the problem of inversing particle size distribution was transformed to finding the minimum value of the objective function.

After the objective function was constructed by Tikhonov (Yang, 2016; Wang et al., 2019; Zhou et al., 20212021) regularization, MPGA (Guo et al., 2020; Shi et al., 2021; Yang et al., 2022) was introduced to find the minimum value.

Solving the objective function minimum value with MPGA

1) Taking the unimodal R-R distribution function as an example, the first step was to determine $\overline{D}$ and K, the characteristic parameter values of the R-R distribution function. Then, n initial populations were randomly generated. Each population contained m individuals. Each individual, which was coded in binary to form a gene coding of a chromosome, represented different set of $\overline{ D}$ and K, which were the characteristic parameter values of a R-R distribution function.

2) The objective fitness function was determined. The function in this paper was $\text{y}=E-A{f}^2$ and $\text{y}=E-A{f}^2+\alpha f^2$.

3) The fitness of each individual in each community was calculated using the substituted individuals, and the findings were rated within each population from high to poor.

4) Excellent individuals were moved from each subpopulation to adjacent subpopulations to replace low-fitness individuals.

5) Selection operation was carried out. To scale the fitness value of individuals in each sub-population, Roulette Wheel Selection was adopted. The chance of elimination would increase as the fitness rating decreased. Then the individual was assessed to see if it met the optimization criteria. If it was qualified, the results would be decoded and exported to obtain the optimal solution and the iteration stopped; Otherwise, the next step would be taken.

6) The fitness function values were compared to determine the regenerated individuals of the population. Excellent individuals were picked and reproduced manually and added to the elite population to be saved. It was also necessary to determine whether they were repeated. If so, the repeated one would be eliminated and the selection operation would start.

7) Crossover and mutation were carried out based on the randomly generated probabilities of each subpopulation. This was to obtain next generation of populations and individuals.

8) A new generation of populations were generated by crossover and mutation, and returned to step (3).

9) Iteration ceased and the results were exported when the algorithm reached the predetermined accuracy. Otherwise, crossover and mutation would continue in step (5) to produce new individuals and populations. The least preserver algebra in the elite population was determined until the preset iterative cutoff value was met.

Figure 2 shows the MPGA structure diagram. Each population was an independent evolution of standard genetic algorithm (SGA), in which single-point crossover and site mutation were used for global search, and roulette for selection (Bouvier et al., 2019).

FIGURE 2

FIGURE 2. MPGA structure diagram.

Data analysis

Excel 2016 was used to process the experiment’s preliminary statistical data. Data and graphics were analyzed and processed with Matlab 2016a, including the calculation of apparent extinction coefficient (Fu and Sun, 2001; Yang et al., 2002; Fu and Sun, 2006; Ding et al., 2020; Wrana et al., 2020), regularization parameter and Mie scattering coefficient. Matlab 2016a was used to inverse different algorithms for calculating particle size distribution. Sauter average particle size $D_{32}$ (Han et al., 2021; Romanov et al., 2021; Tolosana-Moranchel et al., 2021) was applied to assess the accuracy of the measured samples.

Results and analysis

Model establishment and analysis of the simulation experiment

First, numerous simulation experiments were used to verify the accuracy and reliability of the algorithm. The commonly used R-R function was adopted as the particle size distribution function for the inversion with a dependent algorithm (Xu and Nieto-Vesperinas, 2019; Liu et al., 2021; Li et al., 2022). In other words, the characteristic parameters in the known equation set was determined with optimal solution. The volume frequency distribution function under unimodal and bimodal conditions was as follows (Cai et al., 2010):

$f_{R\text{R}-\text{s}}\left(D\right)=\frac{K}{\overline{D}}\cdot {\left(\frac{D}{\overline{D}}\right)}^{K-1}\cdot e^{\left[-{\left(\frac{D}{\overline{D}}\right)}^K\right]}(23)$

$f_{RR-b}\left(D\right)=n\left[\frac{K_1}{{\overline{D}}_1}\cdot {\left(\frac{D}{{\overline{D}}_1}\right)}^{K_1-1}\cdot e^{\left[-{\left(\frac{D}{{\overline{D}}_1}\right)}^{K_1}\right]}\right]+\left(1-n\right)\left[\frac{K_2}{{\overline{D}}_2}\cdot {\left(\frac{D}{{\overline{D}}_2}\right)}^{K_2-1}\cdot e^{\left[-{\left(\frac{D}{{\overline{D}}_2}\right)}^{K_2}\right]}\right](24)$

Where: $\text{RR}-\text{s}$ and $\text{RR}-\text{b}$ represented unimodal and bimodal RR distributions, respectively. Moreover, $\overline{D}$, $\text{k}$, ${\overline{D}}_1$, $k_1$, ${\overline{D}}_2$, $k_2$ were the characteristic parameters, $0\le \text{n}\le 1$.

The complex refractive index was set at $m=1.33+0.01i$. The simulation extinction values of incident wavelengths used 0.3838, 0.4598, 0.5163, 0.6269 and 0.77 $\text{ μm}$. In SGA, the population size, probability of selection, and mutation probability was set at 50, 0.7, and 0.05, respectively. The iteration times of unimodal and bimodal distribution functions were set at 200 and 300, respectively. In MPGA, the initial population number was set at 10, population size 50, selection probability [0.7, 0.9], and mutation probability [0.001.0.05]. The iteration times were the same as in the SGA. The approximate range of distribution parameters following RR distribution was $\text{R}-\text{R}:1<\text{k}<9;1<\overline{\text{D}}<9$, which was inferred from the particle size range of $0.1\sim 10\text{ μm}$. The inversion error was defined as follows (Doicu et al., 2019; Shcherbakov et al., 2019; Zhang et al., 2019):

$\text{δ}=\frac{\sqrt{{\displaystyle \sum _{j=1}^S{\left[f\left(D_j\right)-F\left(D_j\right)\right]}^2}}}{\sqrt{{\displaystyle \sum _{j=1}^S{\left[f\left(D_j\right)\right]}^2}}}(25)$

First, the research verified the effect of regularization and non-regularization objective functions on the reconstruction of particle size distribution. It compared the results with 5% noise or without noise with MPGA, respectively. Then, the least squares error formula and the LS error with Tikhonov smooth function constraint were used as the objective fitness functions. As a result, the following regularization settings were specified for the latter:

The regularization parameters without noise and with 5% noise were $9.329\times {10}^{-6}$ and $6.142\times {10}^{-4}$ under unimodal distribution function, respectively. As for bimodal distribution function, the parameters were $5.231\times {10}^{-5}$ and 0.0028, respectively.

Figure 3 depicts the inversion results for the unimodal R-R distribution, where Figure 3A under the condition of without noise and Figure 3B under the condition of with 5% noise. When no noise is added, the peak position and height of the volume frequency distribution function curves generated by MPGA inversion based on regularization and non-regularization objective functions are very close to the set values, particularly when the peak position of the regularization basically matches the preset curve, despite the variation of the peak height. However, the identification of particle size distribution is unaffected, although the curve distribution of non-regularization moves smoothly. Despite the fact that the peak position is deviated, resulting in judgment errors in the particle size distribution, the total effect is within the error range. When 5% noise is added, the regularization curve modifications are not apparent, which is consistent with the preset distribution curve. Although the peak height varied slightly, it is still within an acceptable range, demonstrating that the inversion algorithm with regularization is strongly noise resistant. The function curve with non-regularization has a great peak position and height offset as the noise increases, while the reconstruction effect is low, showing that the noise resistance of the inversion algorithm without regularization is inadequate.

FIGURE 3

FIGURE 3. Inversion results following unimodal R-R distribution without noise and with 5% noise.

Figure 4 illustrates the inversion results following for the bimodal R-R distribution, where Figure 4A under the condition of without noise and Figure 4B under the condition of with 5% noise. As can be observed, the volume frequency distribution function curves produced by MPGA inversion based on regularization objective function are highly consistent with one another in the absence of noise. The error is acceptable despite slightly higher peaks. The peak positions are basically identical, which does not affect the particle size distribution identification. However, the overall curve deviation of non-regularization is quite large, and the height deviation of the first peak and the position deviation of the second peak are obvious. Additionally, the inversion impact was reduced due to the complexity of the bimodal distribution function, demonstrating that an inversion algorithm with regularization outperforms the one with non-regularization in complex function inversion. Even when 5% noise is added, the overall effect of the regularization curve is still better than that of non-regularization. As the noise rises, the peak heights of the inversion curve for both objective functions grow, indicating that the overall inversion impact is proportional to the bimodal distribution function’s complexity. However, the curve produced by the inversion algorithm with regularization is closer to the preset distribution curve, suggesting that the inversion algorithm with regularization is still more noise-resistant than that with non-regularization even under complex bimodal distribution function.

FIGURE 4

FIGURE 4. Inversion results following bimodal R-R distribution without noise and with 5% noise.

Next, the research then calculated the inversion of the particle size distribution with MPGA and SGA with effects compared to verify the reliability of improved genetic algorithm based MPGA. Figure 5 shows the changing rules of the fitness function value as it varies with iteration times under the unimodal distribution function for SGA and MPGA, where the fitness function is regularization objective function. Figure 5 shows that iteration results using SGA stabilize after 65 iterations, despite the fact that the findings are not unique. Comparatively, MPGA iteration outcomes stabilized after 25 rounds, yielding consistent final findings.

FIGURE 5

FIGURE 5. Change rule of the fitness function value varying with iteration times with SGA and MPGA.

Figures 6, 7 show the volume frequency distribution function curve calculated based on the inversion of particle size distribution after unimodal and bimodal R-R functions with MPGA and SGA, where Figures 6A, 7A under the condition of without noise and Figures 6B, 7B under the condition of with 5% noise. Without noise, the inversion results of both algorithms for unimodal and bimodal functions agree with the predetermined values, and the peak position and height could be rebuilt effectively. Among them, the early peak heights of the bimodal distribution function deviates to a certain degree, however, the peak height deviates further with SGA due to the immature convergence of SGA in the complex objective function situation. The immature convergence of SGA exacerbated by the addition of 5 percent noise, and both the unimodal and bimodal functions exhibit substantial deviations. The inversion effect of MPGA is slightly superior to that of SGA. Although the peak location and height vary, they are still acceptable and the overall effect is superior to that of SGA, indicating that the anti-noise ability of MPGA is outperforms the capability of SGA.

FIGURE 6

FIGURE 6. Inversion results following unimodal distribution with different algorithms.

FIGURE 7

FIGURE 7. Inversion results following bimodal distribution with different algorithms.

Experimental results and analysis

Considering the light absorption of liquid medium, milk fat particles were employed as the research object in actual measurements to further validate the paper’s algorithm. Table 1 displays the extinction levels measured at various wavelengths determined by the mean value of 20 experimental data groups.

TABLE 1

TABLE 1. Extinction values measured at different wavelengths.

Table 2 Displays the inversion data obtained from five inversion tests with SGA and MPGA. According to Table 1, the k values of both algorithms are relatively substantial, indicating that the particle size distribution of the simulation solution of milk fat in the experiment is relatively concentrated and consistent with narrow band distribution. The results of the inversion indicate that the final average values of the two algorithms are quite similar. In contrast, in many trials, the final iteration’s SGA findings are variable, whereas MPGA results are consistent. This demonstrates that MPGA is more robust and effective at optimizing systems and addressing complicated issues.

TABLE 2

TABLE 2. Comparison of inversion results between SGA and MPGA.

Figure 8 illustrates the particle size results obtained with SGA and MPGA. In Figure 8, the Sauter average particle sizes $D_{32}$ are 4.8206 $\mu m$ and 4.8962 $\mu m$, respectively, while the relative errors are -5.01% and -3.52%, respectively, meeting the national standard of ±8%.

FIGURE 8

FIGURE 8. Comparison of volume frequency distribution maps obtained from inversion with SGA and MPGA.

Discussion

Comparative analysis of the inversion results of different objective functions

The least square error is frequently used as the objective function in particle size inversion calculations. In their book (Cai et al., 2010), Cai et al. elaborated on the solution of objective functions and emphasized that the least squares error is used as the objective function for inversion calculations. In practical calculation, however, the extinction coefficient determines the linear equation of the extinction spectrum, primarily because the linear equation resulting from the highly ill-conditioned coefficient matrix is unstable as a result of the glitches generated during the initial calculation process. If we directly solve the least-squares solution as an objective fitness function, the calculation results may not be unified with a reduced level of precision. Therefore, the introduction of Tikhonov regularization into the creation of objective functions can effectively solve the problem of ill solutions. The main goal of Tikhonov regularization is to modify the relative weight between residual and regularization terms and to reduce the ill-condition of the objective function by the introduction of smooth constraints, hence optimizing the stability of the function to be measured. The approach changes the problem into a solution for least-squares problems with penalty factors. By introducing regularization parameters, the literature (Yang, 2016; Wang et al., 2019; Zhou et al., 20212021) efficiently addresses the problem of coefficient matrix pathology encountered during the development of objective functions.

Table 3 compares the characteristic parameters and inversion errors (using MPGA for inversion) for unimodal and bimodal distribution functions calculated using the inversion algorithm with regularization, non-regularization, and 5% noise. As shown in Table 3, the inversion errors of the objective function built by regularization and non-regularization under the unimodal R-R distribution function are 0.0402 and 0.0974, respectively, without noise, and 0.0799 and 0.1450, respectively, with 5% noise. Non-regularization under the bimodal R-R distribution function yields the following results: the inversion errors are 0.0539 and 0.1085 without noise, and 0.0917 and 0.1473, respectively, with 5% noise.

TABLE 3

TABLE 3. Comparison of the characteristic parameters and the inversion errors under unimodal and bimodal distribution functions.

The simulation experiment compares the inversion errors under two different objective functions. By comparing the regularization and non-regularization objective functions, it was discovered that the inversion error with the regularization objective function was significantly smaller than that with the non-regularization objective function, regardless of whether the distribution functions were unimodal or bimodal. The regularization inversion error with 5% noise was even lower than the non-regularization inversion error without noise, indicating that the regularization had a positive effect on alleviating the ill condition of the discrete complex coefficient matrix, enhancing the accuracy of particle size distribution inversion results, and effectively reducing the complexity of the objective function and enhancing the stability. In addition, research data (Wang, 2013; Zhang, 2015; Zhao et al., 2015) indicate that the inversion of the regularization goal function yields excellent results.

Comparative analysis of different inversion algorithms

This work compares the inversion outcomes of SGA and MPGA by calculating the inversion of the particle size of milk fat. We initially examined the iteration times. MPGA times were significantly less than those of SGA because MPGA utilised multiple populations simultaneously in the initial population setting, and the populations were both independent and interconnected. In addition, the optimal individuals from each population were preserved in the elite population for repeat selection and to avoid loss, thereby drastically lowering the number of iterations required to locate the optimal individuals. However, SGA only had one population per iteration for optimal individual selection, which was insufficient to assure loss and repeated selection. Through MPGA, the literature (Yang et al., 2022) significantly enhances the algorithm’s speed. In addition, MPGA’s final iteration findings were consistent, indicating that MPGA’s target parameter search was more accurate and exhaustive. It was mostly owing to the fact that the control parameters in optimal individual selection were distinct. In SGA, the crossover and mutation probabilities were fixed. The crossover and mutation operations determined the global and local search capabilities of the algorithm, while the probability setting influenced the algorithm’s overall search capability. The optimization results for various selections differed considerably. MPGA made up for SGA’s deficiency by randomly generating probability parameters within the specified range to distinguish the crossover and mutation probabilities for each subpopulation. In the meantime, populations evolved synchronously, covering both global and local searches of the algorithm. The problem of inconsistent experimental results is solved by MPGA in the literature (Guo et al., 2020; Shi et al., 2021).

Table 4 compares the characteristic characteristics and inversion errors (using regularization goal function for inversion) of particle size distribution functions produced by MPGA and SGA inversion without noise and with 5% noise. According to Table 3, the inversion results of MPGA and SGA under unimodal R-R distribution function are as follows: the inversion errors are 0.0402 and 0.0563, respectively, without noise; 0.0799 and 0.1158, respectively, with 5% noise; and under bimodal R-R distribution function, the inversion errors are 0.0539 and 0.0665, respectively, without noise; and 0.0917 and 0.1384, respectively, with 5% noise.

TABLE 4

TABLE 4. Comparison of characteristic parameters and inversion errors of particle size distribution under unimodal and bimodal distribution functions.

According to the analysis, MPGA has better inversion effects than SGA, indicating that MPGA has stronger anti-noise capabilities. The inversion error comparison suggests that MPGA reconstruction is more precise than SGA reconstruction. Errors in SGA reconstruction could be decreased further, but doing so would require additional iterations and increase computation costs. SGA requires more search into random points, but MPGA effectively eliminates this drawback. Based on the outcome of this single operation, MPGA was able to better match the specified particle size distribution with high reconstruction accuracy under the same conditions, resolving the premature convergence of SGA. MPGA has faster iteration speed and higher inversion accuracy than SGA, hence avoiding the issue of SGA’s premature convergence.

Conclusion

1) The combination of MPGA and Tikhonov smooth function supports particle size inversion calculations with the total light scattering method in an efficient manner. Due to their high dependability and stability, both functions reconstruct the particle size distribution in a liquid absorption medium.

2) Integrating Tikhonov regularization and apparent optical parameters into the construction of the objective function can effectively resolve the high ill condition of the objective function caused by the first type of Fredholm integral equation with the total light scattering method, as well as reduce the medium’s light absorption.

3) Compared to conventional genetic algorithms, the approach suggested in this study achieves a detection error of -3.52% with a greater detection accuracy, thereby providing a dependable method for milk molecular identification. The algorithm is ubiquitous and applicable to other domains as well.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

GD and ZZ: critical review of the manuscript; YW and PJ: statistical data analyses and writing of the manuscript; supervision of the statistical data analyses; All authors contributed to the article and approved the submitted version.

Funding

This study was funded by Heilongjiang Natural Science Joint Foundation of China (LH 2021C068) and Heilongjiang Bayi Agricultural University’s post-doctoral funding.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Bouvier, M., Yon, J., Lefevre, G., and Grisch, F. (2019). A novel approach for in-situ soot size distribution measurement based on spectrally resolved light scattering. J. Quantitative Spectrosc. Radiat. Transf. 225, 58–68. doi:10.1016/j.jqsrt.2018.12.018

CrossRef Full Text | Google Scholar

Cai, T. T., Yu, D. M., Liu, H. N., and Gao, F. K. (2022). Computational analysis of variational inequalities using mean extra-gradient approach. Mathematics 10 (13), 2318. doi:10.3390/math10132318

CrossRef Full Text | Google Scholar

Cai, X. S., Su, M. X., and Shen, J. Q. (2010). Particle granularity measurement technology and application. Beijing: Chemical Industry Press.

Google Scholar

Ding, G., Duan, H., Wang, X., and Zhou, Z. (2020). Analysis of extinction spectrum apparent property and extraction of specific wavelength of milk fat particles in liquid absorbing medium. Int. J. Perform. Eng. 16 (3), 430–437. doi:10.23940/ijpe.20.03.p12.430437

CrossRef Full Text | Google Scholar

Doicu, A., Vasilyeva, A. A., Efremenko, D. S., Wirth, C. L., and Wriedt, T. (2019). A light scattering model for total internal reflection microscopy of geometrically anisotropic particles. J. Mod. Opt. 66 (10), 1139–1151. doi:10.1080/09500340.2019.1605005

CrossRef Full Text | Google Scholar

Fu, Q., and Sun, W. (2006). Apparent optical properties of spherical particles in absorbing medium. J. Quant. Spectrosc. Radiat. Transf. 100 (1-3), 137–142. doi:10.1016/j.jqsrt.2005.11.031

CrossRef Full Text | Google Scholar

Fu, Q., and Sun, W. (2001). Mie theory for light scattering by a spherical particle in an absorbing medium. Appl. Opt. 40 (9), 1354–1361. doi:10.1364/ao.40.001354

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, C., Yang, Z., Li, S., and Lou, J. (2020). Predicting the water-conducting fracture zone (WCFZ) height using an MPGA-SVR approach. Sustainability 12 (5), 1809. doi:10.3390/su12051809

CrossRef Full Text | Google Scholar

Han, J., Liu, X., Jiang, M., Wang, Z., and Xu, M. (2021). A novel light scattering method with size analysis and correction for on-line measurement of particulate matter concentration. J. Hazard. Mater. 401, 123721. doi:10.1016/j.jhazmat.2020.123721

PubMed Abstract | CrossRef Full Text | Google Scholar

He, Z., Mao, J., and Han, X. (2018). Non-parametric estimation of particle size distribution from spectral extinction data with PCA approach. Powder Technol. 325, 510–518. doi:10.1016/j.powtec.2017.11.050

CrossRef Full Text | Google Scholar

Igushi, T., and Yoshida, H. (2011). Investigation of low-angle laser light scattering patterns using the modified Twomey iterative method for particle sizing. Rev. Sci. Instrum. 82 (1), 015111. doi:10.1063/1.3520136

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, P., Li, Y., Gao, R., Xu, C., and Shang, Y. (2022). New exploration on bifurcation in fractional-order genetic regulatory networks incorporating both type delays. Eur. Phys. J. Plus 137 (5), 598. doi:10.1140/epjp/s13360-022-02726-3

CrossRef Full Text | Google Scholar

Lienert, B. R., Porter, J. N., and Sharma, S. K. (2001). Repetitive genetic inversion of optical extinction data. Appl. Opt. 40 (21), 3476–3482. doi:10.1364/ao.40.003476

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Mao, S., Song, S., Huang, L., Belfiore, L. A., and Tang, J. (2021). Towards applicable photoacoustic micro-fluidic pumps: Tunable excitation wavelength and improved stability by fabrication of Ag-Au alloying nanoparticles. J. Alloys Compd., 161091. doi:10.1016/j.jallcom.2021.161091

CrossRef Full Text | Google Scholar

Maguire, C. M., Rösslein, M., Wick, P., and Prina-Mello, A. (2018). Characterisation of particles in solution-a perspective on light scattering and comparative technologies. Sci. Technol. Adv. Mater. 19 (1), 732–745. doi:10.1080/14686996.2018.1517587

PubMed Abstract | CrossRef Full Text | Google Scholar

Mroczka, J., and Szczuczyński, D. (2013). Improved technique of retrieving particle size distribution from angular scattering measurements. J. Quantitative Spectrosc. Radiat. Transf. 129, 12948–12959. doi:10.1016/j.jqsrt.2013.05.030

CrossRef Full Text | Google Scholar

Romanov, A. V., Yurkin, M. a. J. L., and Reviews, P. (2021). Single‐particle characterization by elastic light scattering. Laser & Photonics Rev. 15 (2), 2000368. doi:10.1002/lpor.202000368

CrossRef Full Text | Google Scholar

Shcherbakov, M. R., Lemasters, R., Fan, Z., Song, J., Lian, T., Harutyunyan, H., et al. (2019). Time-variant metasurfaces enable tunable spectral bands of negative extinction. Optica 6 (11), 1441–1442. doi:10.1364/optica.6.001441

CrossRef Full Text | Google Scholar

Shi, B., Chang, X., Liu, Z., Pang, Y., Xu, Y., Mao, L., et al. (2021). Intelligent identification of sequence stratigraphy constrained by multipopulation genetic algorithm and dynamic time warping technique: A case study of lower cretaceous qingshuihe formation in hinterland of junggar basin (NW China). Basin Res. 33 (4), 2517–2544. doi:10.1111/bre.12567

CrossRef Full Text | Google Scholar

Tolosana-Moranchel, A., Pecharromán, C., Faraldos, M., and Bahamonde, A. (2021). Strong effect of light scattering by distribution of TiO2 particle aggregates on photocatalytic efficiency in aqueous suspensions. Chem. Eng. J. 403, 126186. doi:10.1016/j.cej.2020.126186

CrossRef Full Text | Google Scholar

Wang, L. (2013). Research on particle size distribution reconstruction algorithm based on spectral extinction method. Harbin, China: Harbin Institute of Technology.

Google Scholar

Wang, L., and Sun, X. G. (2013). Research on pattern search method for inversion of particle size distribution in spectral extinction technique. Spectrosc. Spectr. Analysis 33 (03), 618–622.

PubMed Abstract | Google Scholar

Wang, P., Wang, L., Leung, H., and Zhang, G. (2021). Super-resolution mapping based on spatial-spectral correlation for spectral imagery. IEEE Trans. Geosci. Remote Sens. 59 (3), 2256–2268. doi:10.1109/tgrs.2020.3004353

CrossRef Full Text | Google Scholar

Wang, Y. F., Liu, H. M., and Fu, Z. W. (2019). Low-light image enhancement via the absorption light scattering model. IEEE Trans. Image Process. 28 (11), 5679–5690. doi:10.1109/tip.2019.2922106

PubMed Abstract | CrossRef Full Text | Google Scholar

Wrana, F., Von Savigny, C., Zalach, J., and Thomason, L. W. (2020). Retrieval of stratospheric aerosol size distribution parameters using SAGE-III/ISS extinction measurements at three wavelengths. Atmos. Meas. Tech. Discuss., 1–20.

Google Scholar

Xu, X., and Nieto-Vesperinas, M. (2019). Azimuthal imaginary poynting momentum density. Phys. Rev. Lett. 123 (23), 233902. doi:10.1103/physrevlett.123.233902

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, P., Gao, B. C., Wiscombe, W. J., Mishchenko, M. I., Platnick, S. E., Huang, H. L., et al. (2002). Inherent and apparent scattering properties of coated or uncoated spheres embedded in an absorbing host medium. Appl. Opt. 41 (15), 2740–2759. doi:10.1364/ao.41.002740

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Yang, R., Tan, S., Yu, X., and Fang, L. (2022). MPGA-based-ECMS for energy optimization of a hybrid electric city bus with dual planetary gear. Proc. Institution Mech. Eng. Part D J. Automob. Eng. 236 (8), 1889–1909. doi:10.1177/09544070211041074

CrossRef Full Text | Google Scholar

Yang, Y. (2016). Study on detection technology of PM2.5 particle size distribution. Hangzhou, China: China Jiliang University.

Google Scholar

Zhang, W., Shen, J., Thomas, J. C., Mu, T., Xu, Y., Xiu, W., et al. (2019). Particle size distribution recovery in dynamic light scattering by optimized multi-parameter regularization based on the singular value distribution. Powder Technol. 353, 320–329. doi:10.1016/j.powtec.2019.05.040

CrossRef Full Text | Google Scholar

Zhang, Y. M. (2015). Research on the optimization method of inversion algorithm in polydisperse nanoparticle size based on regularization. Jinan, China: Qilu University of Technology.

Google Scholar

Zhao, Y. J., Gao, C. B., et al. (2015). Improving computational method of mine dust distribution by eikonal approximation. Coal Mine Saf. 46 (08), 178–180.

Google Scholar

Zhou, G., Li, C., Zhang, D., Liu, D., Zhou, X., and Zhan, J. (2021). Overview of underwater transmission characteristics of oceanic LiDAR. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 14, 8144–8159. doi:10.1109/jstars.2021.3100395

CrossRef Full Text | Google Scholar

Zhou, G. Q., Li, W. H., Zhou, X., Tan, Y. Z., Lin, G. C., Li, X. X., et al. (2021). An innovative echo detection system with STM32 gated and PMT adjustable gain for airborne LiDAR. Int. J. Remote Sens. 42 (24), 9187–9211. doi:10.1080/01431161.2021.1975844

CrossRef Full Text | Google Scholar