Quasi-Random Sampling-Enhanced Metaheuristic Algorithms for CAMD-Based Solvent Selection in Octacosanol Extraction

1. Introduction

Natural products often occur as complex mixtures, making the isolation of target compounds an important objective because of their specific applications. In industry, extraction is one of the most widely used separation operations, and solvent selection is critical to its effectiveness. In an extraction system, the solvent should selectively attract the desired compound into the extract phase while leaving impurities in the raffinate phase. For a given system, solvents are generally selected using one of three approaches [1,2]:

(a)

Experimental screening of candidate solvents for the system,

(b)

Database screening based on physical, chemical, and related properties, and

(c)

Computer-aided molecular design (CAMD) to generate optimal solvents.

The first two approaches may not identify the best solvent and can lead to incomplete evaluation. Although experimental screening provides reliable results, it is often expensive and time-consuming. In contrast, CAMD uses thermodynamic models to design solvents tailored to target performance requirements.[1] In CAMD, molecules are generated by combining available functional groups to ensure that the resulting structures meet desired performance metrics. Because the number of possible group combinations and attachment patterns is extremely large, combinatorial optimization is central to this problem.

Papadimitriou and Steiglitz [3] define a combinatorial optimization problem as one that seeks an optimal solution, such as an integer, permutation, or graph structure, from a finite or countably infinite set. Common examples include the traveling salesman and scheduling problems.[4] CAMD can also be formulated as a combinatorial optimization problem. Metaheuristic algorithms, including ant colony optimization (ACO), evolutionary computation, genetic algorithms, and simulated annealing (SA), use random number generation during the search process, and their performance strongly depends on the quality of the sampling. In this work, we enhance SA and ACO by replacing the pseudo-random numbers used in conventional implementations with quasi-Monte Carlo samples generated through Hammersley sequence sampling. The resulting algorithms, efficient ant-colony optimization (EACO) and efficient simulated annealing (ESA), preserve k-dimensional sampling uniformity.

The optimization problem considered in this study focuses on extracting nutraceuticals from multicomponent sugarcane wax. Nutraceuticals are health-promoting compounds derived from food sources such as animal and plant waxes, as well as apple and grape seeds, and they provide medicinal benefits beyond basic nutrition.[5] As sugarcane is produced on a large scale worldwide, extracting nutraceuticals from sugarcane waste has become an important separation operation.[6] Sugarcane wax contains various long-chain fatty compounds, including alcohols, esters, aldehydes, acids, and saturated hydrocarbons. Among these, saturated alcohols known as policosanols are important nutraceuticals. Policosanols typically contain 18–34 carbon atoms and a hydroxyl group (-OH) on the alpha carbon. Octacosanol, a 28-carbon policosanol, is particularly promising because of its reported health benefits, including immune-supporting, anti-Parkinsonian, and anti-fatigue effects.[7,8] Therefore, this study uses CAMD to identify optimal solvents for extracting octacosanol from the multicomponent mixture present in sugarcane wax.

To provide context, the next section discusses the importance of octacosanol and reviews reported extraction studies for obtaining it from sugarcane wax. This is followed by a brief overview of CAMD, molecular construction, and the group contribution method. The formulation and development of the solvent selection optimization problem for extracting octacosanol from multicomponent sugarcane wax are then presented. Section 3 describes the algorithms used in this study, beginning with ACO and SA, followed by previously reported enhancements to these methods. It then introduces Hammersley sequence sampling (HSS), explains its advantages, and describes its incorporation into ACO and SA to develop EACO and ESA. The algorithmic frameworks and parameter values used in this work are also provided. Section 4 presents the solvents generated by the four solvers and compares the efficiency of EACO and ESA with that of ACO and SA, respectively. Finally, the study concludes by identifying the best solvents for extracting octacosanol from sugarcane wax and comparing the four solvers for this process-systems optimization problem.

2. The Solvent Selection Problem

Policosanols offer several health benefits, including lipid regulation, anti-atherosclerotic effects, and support in managing obesity, etc.[8,9] Among them, octacosanol is a particularly promising nutraceutical. Numerous studies have attempted to extract policosanols or octacosanol-rich fractions from sugarcane using methods such as Soxhlet extraction, supercritical CO₂ (99.99%), and subcritical dimethyl ether extraction.[10,11,12,13]

However, in most cases, solvent selection has been ad-hoc, guided largely by prior experience. This approach risks overlooking solvents with superior performance, especially, since experimentally screening a large solvent set is costly and time-consuming. A systematic, performance-based solvent selection strategy is therefore essential.

Solvent choice in extraction depends on properties such as the distribution coefficient, selectivity, recoverability, and environmental safety.[1,5] In this study, we focus on distribution coefficient and solvent selectivity.

• Distribution coefficient is the ratio of the desired solute’s concentration in the solvent phase to its concentration in the original mixture.
• Solvent selectivity measures how effectively a solvent attracts the desired solute relative to undesired components.

These metrics form the basis of the solvent selection optimization problem developed in the following sections.

2.1. Computer-Aided Molecular Design

The group contribution method estimates molecular properties—melting and boiling points, density, viscosity, and activity coefficients (γ)—from the functional groups that make up a molecule. As distribution coefficient and selectivity can be determined using activity coefficients, this method is central to our analysis.

CAMD reverses this logic: instead of predicting properties from a molecule, it constructs molecules from available groups to optimize performance metrics.[1] CAMD has been widely applied in designing solvents for separation processes, refrigerants, polymers, pharmaceuticals, and more.[1,2,14,15,16,17] Recent advances even integrate machine learning to generate molecules with drug-like characteristics.[18]

In CAMD, molecules are built by selecting:

• N₁: the total number of functional groups
• N₂⁽ⁱ⁾: the index of each group (i = 1,…,N₁)

Feasible molecular structures must satisfy chemical bonding rules — the octet rule for covalent systems.[5] Thus, N₁ and N₂⁽ⁱ⁾ serve as decision variables, while bonding constraints ensure chemically valid structures.

Figure 1 outlines the CAMD workflow for identifying optimal solvents for sugarcane wax extraction. The process begins by identifying mixture components and estimating their activity coefficients, followed by formulating and solving the solvent-performance optimization problem.

2.2. Group Contribution Method: UNIFAC

The group contribution method (GCM) predicts thermodynamic properties by decomposing molecules into functional groups and summing their contributions.[19] This approach is especially valuable when experimental data are unavailable, enabling property estimation across chemical engineering, molecular design, and process simulation.

Among the many properties GCM can predict, activity coefficients are particularly important for mixture behavior. They quantify deviations from ideality, with γ = 1 representing ideal solutions. For solvent design, the most widely used GCM for estimating activity coefficients is UNIFAC, which provides reliable predictions for non-ideal liquid mixturesFor molecule i present in a solution, its activity coefficient, ${\gamma }_i$, is defined as shown in Equation 1.[19]

$$U\left({\gamma }_i\right)=ln {\gamma }_i-\left(ln {\gamma }_i^C+ln {\gamma }_i^R\right)=0$$

(1)

Equation 1 is the standard UNIFAC equation to estimate the activity coefficient. The estimation of ${\gamma }_i$ is done in two parts namely, the combinatorial part (${\gamma }_i^C$) and the residual part (${\gamma }_i^R$). Such a distinction allows us to understand the effects of size and shape, and the group energy interactions. The combinatorial part takes care of the effects caused by size and shape while on the other hand, the residual part takes care of the energy interactions among the groups of molecules present in the solution.

Firstly, the combinatorial part is estimated as shown in Equation 2(a-f). [19]

$$ln{\gamma }_i^C=\ln \left({\displaystyle \frac{{\Phi }_i}{x_i^b}}\right)+5q_i\ln \left({\displaystyle \frac{{\theta }_i}{{\Phi }_i}}\right)+l_i-{\displaystyle \frac{{\Phi }_i}{x_i^b}}\sum _j{x}_j^b{l}_j$$

(2a)

where

${\Phi }_i={\displaystyle \frac{r_i{x}_i^b}{{\displaystyle {\sum }_j^{N_C}}{r}_j{x}_j^b}}$

(2b),

${\theta }_i={\displaystyle \frac{q_i{x}_i^b}{{\displaystyle {\sum }_j^{N_C}}{q}_j{x}_j^b}}$

(2c),

$l_i=5\left(r_i-q_i\right)-\left(r_i-1\right)$

(2d)

where ${\Phi }_i$ and ${\theta }_i$ are the molecular volume fraction and the molecular surface area fraction of molecule i, ranging from 1 to N_C molecules. This is followed by the estimation of van der Waals’ volume, $r_i$, and surface area, $q_i$, using Equation 2(e) and Equation 2(f), respectively.

$r_i={\displaystyle \sum _g^{N_g}}{\nu }_g^{\left(i\right)}{R}_g$

(2e),

$q_i={\displaystyle \sum _g^{N_g}}{\nu }_g^{\left(i\right)}{Q}_g$

(2f)

where $R_g$ and $Q_g$ are van der Waals’ group volume and group surface area values ranging from $i$ to $N_g$ groups, and ${\nu }_g^{\left(i\right)}$ is the total number of functional groups present in the molecule. The values of $R_g$ and $Q_g$ are obtained from the literature for this study.

While the estimation of the residual part is based on the Equation 3(a-e). [19]

$$ln{\gamma }_i^R=\sum _k{v}_k^i\left(\ln {\Gamma }_k-\ln {\Gamma }_k^i\right)$$

(3a)

where ${\Gamma }_k$ is defined as the residual activity coefficient of group k in a solution, shown in Equation 3(b). Similarly, ${\Gamma }_k^i$ is defined as the residual activity coefficient of group k in a solution containing only molecules of type i.

$$ln{\Gamma }_k=Q_k\left[1-ln\left(\sum _m{\Theta }_m{\Psi }_{mk}\right)-\sum _m{\displaystyle \frac{{\Theta }_m{\Psi }_{km}}{\sum _n{\Theta }_n{\Psi }_{nm}}}\right]$$

(3b)

where m and n represent all the functional groups in the molecule, ${\Theta }_m$ is the surface fraction of group m , Equation 3c, which is further expressed in terms of mole fraction of groups present in the bulk phase, Equation 3d, and lastly, the temperature-dependent parameter as per UNIFAC is ${\Psi }_{mn}$, Equation 3e.[19]

$${\Theta }_m={\displaystyle \frac{Q_m{X}_m}{\sum _m{Q}_m{X}_m}}$$

(3c)

$$X_m={\displaystyle \frac{{\sum }_i^j{\nu }_m^{\left(i\right)}{x}_j}{{\sum }_1^j{\sum }_1^n{\nu }_m^{\left(i\right)}{x}_j}}$$

(3d)

$${\Psi }_{mn}=e^{-{\displaystyle \frac{a_{mn}}{T}}}$$

(3e)

In Equation 3e, $a_{mn}$ is called the group interaction parameter (represents interaction of group m with group n). It is defined as the difference in the energy level between group m and group n or between two of the main group m. In the latter case, the interaction parameter is zero because of the same energy levels of both the groups. These values are obtained from the literature for every group. [19]

In the next subsection, we discuss the development of the optimization problem for solvent selection used in this study along with the parameters important to build the optimization problem.

2.3. Solvent Selection Optimization Problem

This section formulates the solvent-selection optimization problem for extracting a desired compound from a multicomponent mixture. Solvent performance is evaluated using the distribution coefficient ($m$) and solvent selectivity ($\beta$), both estimated using the UNIFAC group contribution method. For a binary mixture containing desired solute B, undesired solute A, and solvent S, these quantities are defined by Equations 4(a) and 4(b), respectively.[5]

$$m_B={\displaystyle \frac{C_{B, S}}{C_{B, A}}}={\displaystyle \frac{{\gamma }_{BA}}{{\gamma }_{BS}}}$$

(4a)

$${\beta }_B={\displaystyle \frac{m_B}{m_A}}={\displaystyle \frac{C_{B,S}}{C_{A,S}}}={\displaystyle \frac{{\gamma }_{SA}}{{\gamma }_{BS}}}$$

(4b)

where C_x,y represents the concentration of x present in y. However, in the present study, the system is a multicomponent. Therefore, a recent study has reported the definitions of $m$ and $\beta$ in a better way to estimate these based on the mole fractions ($x_{A_i}$) of every undesired solute (A_i, i=1,2,.., N) present in the multicomponent system. Equation 5(a) and Equation 5(b) shows the net distribution coefficient ($m_{net}$) and the net solvent selectivity (${\beta }_{net}$), respectively, as presented in the literature.[5]

$$m_{net}={\displaystyle \frac{{\gamma }_{BX}}{{\gamma }_{BS}}}$$

(5a)

$${\beta }_{net}=\sum _{i=1}^N{x}_{A_i}{\beta }_i$$

(5b)

where ${\gamma }_{BX}={\displaystyle \frac{{\sum }_{i=1}^N{x}_{A_i}{\gamma }_{B{A}_i}}{{\sum }_{i=1}^N{x}_{A_i}}}$, and ${\beta }_i$ is the solvent selectivity value obtained between desired solute B and the undesired solute A_i. The system of undesired solutes is represented as X. That is the reason the activity coefficient (${\gamma }_{BX}$) is defined between the desired solute B and every undesired solute A_i by weighing against their corresponding mole fraction.

Next, in this study, the multicomponent system is built based on earlier experimental studies. Sugarcane wax is a complex matrix of compounds. Therefore, we considered a few of the most abundant compounds reported in the literature to be present in this multicomponent mixture, while the desired solute remains to be octacosanol.[13] Figure 2 shows the system of compounds considered in this study along with their corresponding mole fractions as reported in the literature.

Solvent generation requires a defined set of UNIFAC groups. In this study, candidate molecules are constructed by selecting up to 12 groups from the 20 groups listed in Table 1, which were obtained from Fredenslund et al.[19] Each group is indexed from 1 to 20. A molecule may contain multiple instances of the same group, but the final structure must satisfy the octet rule. This feasibility constraint is expressed mathematically in Equation 6.[2]

$${\displaystyle \sum _i^{N_1}}{b}_i=2(N_1-1)$$

(6)

where $b_i$ determines the number of groups that can be attached to the group index i.

Finally, our objective is to maximize the net distribution coefficient of the solvent generated by any of the 20 groups listed in Table 1, so that it can isolate octacosanol with a net solvent selectivity (${\beta }_{net,min}$) greater than or equal to the set value from the ten undesired solutes present in the mixture. Based on these, the optimization problem is defined as shown in Equation 7.

$$\underset{N_1,N_2^{\left(i\right)}}{max} m_{net}$$

(7)

$$\mathrm{s}.\mathrm{t}. {\beta }_{net}\ge {\beta }_{net,min}$$

$$U\left({\gamma }_i\right)=0$$

$${\displaystyle \sum _i^{N_1}}{b}_i=2(N_1-1)$$

$$2\le N_1\le 12$$

$$1\le N_2^{\left(i\right)}\le 20 \text{Ɐ} i \mathsf{ϵ} [1,2,\dots ,N_1]$$

Equation 7 shows that the decision variables, $N_1$ and $N_2^{\left(i\right)}$, are discrete, and the objective function and constraints, which are evaluated using UNIFAC, are nonlinear; therefore, the problem is a mixed-integer nonlinear programming problem. The combinatorial complexity arises because the groups in Table 1 can be connected in many feasible ways. To address this challenge, metaheuristic algorithms such as ACO and SA, and their efficient variants ESA and EACO, are used to identify optimal solvent candidates. The next section presents the four solvers used in this work to solve this optimization problem.

3. Algorithmic Framework

As mentioned, ACO and SA are two metaheuristic approaches for solving combinatorial optimization problems. Firstly, we briefly discussed the two conventional algorithms in this section, followed by a discussion of Hammersley sequence sampling. Then, the implementation of the HSS technique in ACO and SA to develop EACO and ESA, respectively, is discussed. Furthermore, the efficiency of these algorithms compared with their counterparts reported in the literature is discussed. In addition, the frameworks of both EACO and ESA are also presented along with the algorithm parameters used in this study.

3.1. Ant-Colony Optimization

Ant-colony optimization (ACO) is a nature-inspired metaheuristic based on the foraging behaviour of ants, introduced by Dorigo for solving combinatorial optimization problems.[23] In nature, ants deposit pheromones as they return from food sources, allowing other ants to follow promising paths. Because shorter paths are traversed more frequently, pheromone accumulates faster on them, gradually guiding the colony toward the shortest route.[23,24] Similarly, in ACO, artificial ants construct solutions using pheromone information and, when available, heuristic information. At iteration iter, the probability of selecting solution component j from state I, (${\mathbb{P}}_{ij}^{iter}$) is given by Equation 8.[23,24]

$${\mathbb{P}}_{ij}^{iter}={\displaystyle \frac{{\tau }_{ij}^{\alpha }\left(iter\right){\eta }_{ij}^{\beta }}{\sum _{iter}{\tau }_{ij}^{\alpha }\left(iter\right){\eta }_{ij}^{\beta }}}$$

(8)

where and are the pheromone value and edge desirability for edge ij, respectively, at iteration iter, while and determine the influence of pheromone and heuristic information.[24] After all ants construct solutions, pheromone trails are updated according to solution quality, with additional pheromone assigned to edges in better solutions, as shown in Equation 9.[24]

$${\tau }_{ij}\left(iter+1\right)=\rho {\tau }_{ij}\left(iter\right)+∆{\tau }_{ij}\left(iter\right)$$

(9)

where $\rho$ is the pheromone evaporation factor, and $∆{\tau }_{ij}$ is the pheromone addition for edge ij.

The artificial pheromone along the path acts as a communication bridge among artificial ants. In the less favorable solution routes, the artificial pheromone decays gradually, whereas in the better solution routes, the pheromone concentration increases, as observed in Equation 9. Eventually, iterative updating of pheromone values based on information gained from earlier iterations encourages the artificial ants towards the optimal solution in the solution space.

The solution construction in this study is based on incremental manner i.e., an ant, first, probabilistically chooses one of the solutions from the archive based on the Equation 10.[24]

$${\mathbb{P}}_i={\displaystyle \frac{{\omega }_i}{{\displaystyle {\sum }_{j=1}^K}{\omega }_j}}$$

(10)

where ${\omega }_i$ is the weight associated with solution i, i is the rank of the solution in the archive, K is the size of solution archive. The weight is defined as shown in Equation 11.[24]

$${\omega }_i={\displaystyle \frac{\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{{\left(i-1\right)}^2}{2{\left(qK\right)}^2})}{qK\sqrt{2\pi }}}$$

(11)

where q is the algorithmic parameter of ACO, and mean of the Gaussian function is 1. However, in this study, the optimization problem is a mixed integer non-linear programming problem. Therefore, using a set of probability density functions, a multimodal one-dimensional density function is obtained. To estimate the multimodal density function, the following equation, Equation 12, is used.[25]

$$G_j\left(x\right)={\displaystyle \sum _{i=1}^K}{\omega }_i{\displaystyle \frac{\exp \left(-\frac{{\left(x-{\mu }_i\right)}^2}{2{\sigma }_i^2}\right)}{{\sigma }_i\sqrt{2\pi }}}, i,j\in \left\{1,\dots ,K\right\}$$

(12)

where $G_j\left(x\right)$ is the Gaussian kernel, ${\omega }_i$ is the weight associated with individual Gaussian function, ${\mu }_i$ and ${\sigma }_i$ are the vector means of the individual solution components and standard deviations, respectively.

3.2. Simulated Annealing

Simulated annealing (SA) is a probabilistic metaheuristic for finding near-optimal solutions to optimization problems. It is inspired by metallurgical annealing, in which a material is heated and then slowly cooled so that its particles settle into a low-energy state. Similarly, SA gradually reduces a control parameter, analogous to temperature, to guide the search toward lower objective-function values while allowing occasional uphill moves. If cooling is sufficiently slow, the system can approach thermal equilibrium at each temperature level, described by the Boltzmann distribution. At temperature T, the probability of occupying energy level E_i is given by Equation 13.[26]

$${\mathbb{P}}_T={\displaystyle \frac{exp\left(-\frac{E_i}{k_{B}T}\right)}{Z}}$$

(13)

where 1/Z is the normalization factor, and k_B is the Boltzmann’s constant (1.3806 x 10^-23 J/K).

In SA, the objective function represents the system energy. Starting from an initial configuration, the algorithm generates a candidate configuration by perturbing the current one. If the candidate has a lower objective value, the move is accepted. If it has a higher objective value, the move may still be accepted with probability according to the Metropolis criterion, given in Equation 14.[26,27]

$$\mathrm{Accept} \mathrm{the} \mathrm{move} \mathrm{if} {Pr}_{accept}\le A_{ij}=\left\{\begin{array}{c}exp\left(-{\displaystyle \frac{∆E}{T}}\right) if ∆E=E_j-E_i\ge 0\\ 1 otherwise\end{array}\right\} \forall i,j\in \mathrm{move}$$

(14)

where $A_{ij}$ is the probability of accepting a transition from state i to state j. At high temperatures, SA accepts more uphill moves, improving exploration of the solution space. As temperature decreases, the probability of accepting such moves declines, allowing the algorithm to converge while reducing the risk of becoming trapped in local optima.

ACO and SA differ substantially in their search mechanisms. SA follows a single-solution trajectory, exploring the search space through successive perturbations of one candidate solution; therefore, complex problems may require many iterations, and escape from local optima depends on the Metropolis criterion.[26,28] In contrast, ACO uses multiple ants to construct solutions in parallel and share information through pheromone trails, which improves exploration and reduces the likelihood of premature trapping.[23] However, as a population-based heuristic, ACO may still miss some high-quality regions of the search space.

3.3. Hammersley Sequence Sampling

Sampling strongly affects optimization efficiency. Monte Carlo sampling uses independently generated pseudo-random numbers over the sampling domain; these samples are statistically unbiased and useful for estimating probabilities and objective-function values in complex systems.[27] However, pseudo-random samples can cluster in some regions while leaving others unexplored, as shown in Figure 3. This reduces sampling efficiency and often requires many samples to adequately cover the search space.[27,29] Quasi-random methods, such as Hammersley sequence sampling (HSS), reduce this clustering by distributing samples more uniformly. HSS, developed by Kalagnanam and Diwekar based on Hammersley points [30], provides better multidimensional uniformity than Monte Carlo methods and other methods, such as Latin hypercube sampling, as illustrated in Figure 3.[27]

The next section explains how HSS is incorporated into ACO and SA to develop variants that preserve k-dimensional sampling uniformity.

3.4. Efficient Ant-Colony Optimization (EACO)

ACO can suffer from premature convergence when excessive pheromone accumulation reduces population diversity and limits exploration. Several variants have been developed to improve its performance, including the Elitist Ant System, Rank-Based Ant System, Max-Min Ant System, Ant Colony System, hypercube-based ACO, and hybrids with local search methods such as tabu search.[4] Efficient Ant-Colony Optimization (EACO) improves ACO by replacing selected pseudo-random numbers with quasi-random numbers generated using Hammersley sequence sampling (HSS).[24,29] Reported studies show that EACO improves computational efficiency by 3%–71% for selected problems.[24] In this study, HSS is used only during initialization, where k-dimensional uniformity is most important for generating the initial ant population, as shown in Figure 4.

In both ACO and EACO, 30 artificial ants are used, with a pheromone decay factor of 0.75 and a solution archive of size 1500 × 13. The first 12 columns store the group indices used to generate the solvent molecule, and the last column stores the objective-function value. The algorithm terminates when the absolute difference between successive objective-function values falls below 10^-6 or when the maximum limit of 20,000 iterations is reached. Because solvent generation involves multiple constraints, the objective function is formulated as a penalty function. Following the oracle penalty method reported by Kalagnanam and Diwekar,[30] the objective function is transformed into an equality constraint of the form OF−Ω, where Ω is the oracle value and is set to 3. Further details on the oracle penalty method are available in Diwekar and Gebreslassie.[24]

3.5. Efficient Simulated Annealing (ESA)

SA uses two probabilities: the acceptance probability, $A_{ij}$, and the generational probability,$G_{ij}$.[27] The acceptance probability determines whether a candidate configuration is accepted, whereas the generational probability controls how new configurations are sampled near the current solution. SA efficiency depends strongly on the generational probability.[27] In conventional SA, these probabilities are generated using uniformly distributed pseudo-random numbers, which can cluster sampled moves and require longer Markov chains. Several SA variants, including fast SA and hybrid SA, modify the generational probability, cooling schedule, or combine SA with methods such as differential evolution and genetic algorithms.[31,32,33] However, these approaches do not directly address sampling uniformity. Kim and Diwekar improved configurational sampling by incorporating HSS into SA, leading to efficient simulated annealing (ESA).[27] As shown in Figure 5, ESA replaces pseudo-random numbers with HSS to generate configuration moves, thereby improving k-dimensional uniformity and reducing the Markov chain length required at each temperature level. Studies report that ESA is approximately 30%–54% more efficient than conventional SA in terms of total moves.[27]

In this study, the initial SA temperature is 50 and is reduced by a cooling factor of 0.85 until it reaches the freezing temperature of 10^-7. The Markov chain length is 1000 for both SA and ESA. SA uses pseudo-random numbers for all four generational probabilities and the acceptance probability. In ESA, the four generational probabilities are generated using HSS, while the acceptance probability remains pseudo-random. Both algorithms terminate when the temperature falls below the freezing temperature or when the candidate solvent remains unchanged for 10 consecutive moves across different samples.

4. Results and Discussions

This study evaluates the efficiency of the proposed metaheuristic algorithms for solving the CAMD problem of selecting solvents to extract octacosanol from a multicomponent sugarcane wax mixture. In Equation 13, the minimum net solvent selectivity value (β_{net, min}) is set to 1.2, and the top 20 candidate solvents are generated using ACO, SA, EACO, and ESA.

4.1. Candidate Solvents Using ACO and SA

Based on the algorithm parameters mentioned in the methodology section, the candidate solvents generated using ACO and SA are listed in Table 2 and Table 3, respectively.

As shown in Table 2 and Table 3 and Figure 6, the top nine solvents generated by ACO and SA are identical. However, SA identified better candidates among the bottom 11 solvents. Despite this improvement, SA required substantially more objective function evaluations than ACO: 1,271,876 compared with 153,330. ACO takes 12% of the computational cost as compared to SA. This higher computational cost arises because SA produces one solution per run and must be repeated with different initial values to obtain multiple solutions, whereas ACO generates all 20 solutions in a single run.

4.2. Candidate Solvents Using EACO and ESA

The same optimization problem was then solved using EACO and ESA. Table 4 lists the candidate solvents obtained from EACO for the problem defined in Equation 7. Several optimal solvents generated by EACO were also identified by conventional ACO, as shown in Figure 7. The top 15 solvents are identical for EACO and ACO, while EACO produced better candidates among the remaining five. EACO required 131,550 objective function evaluations to obtain these solvents, and reducing the computational requirements by 14.2 %.

Similarly, Table 5 lists the optimal solvents obtained by solving the solvent-selection optimization problem with ESA. Similar to Figure 7, Figure 8 compares the net distribution coefficient values for the solvents generated by SA and ESA, excluding the top five candidates. ESA again outperformed SA by identifying five better solvents among the lower-ranked 20 candidates. ESA required 1,264,327 objective function evaluations, which is 0.59% fewer than SA.

Figure 9 presents solvents obtained by ESA and EACO. This figure shows the same trend as the comparison of SA and ACO in Figure 6. The top 9 solvents found using EACO are same as that of found in ESA case, while the bottom 11 are better solvents found by ESA as compared EACO. However, EACO takes 10.4% of the computational cost as compared to ESA.

5. Conclusion

This study developed and evaluated a CAMD-based framework for selecting solvents to extract octacosanol from multicomponent sugarcane wax. The solvent-selection problem was formulated as a mixed-integer nonlinear programming problem using UNIFAC-based estimates of the net distribution coefficient and net solvent selectivity. Four metaheuristic solvers—ACO, SA, EACO, and ESA—were compared to assess both solvent quality and computational efficiency. Across all methods, the highest-ranked candidates were consistent, indicating that the CAMD formulation reliably identifies strong solvent structures for the target extraction. Chemically, the best candidate solvent, ethane, achieved the highest net distribution coefficient, and other solvents, aldehyde- and carboxylic-acid-containing structures, also satisfied the selectivity constraint and appeared among the top candidates. Incorporating the Hammersley sequence sampling improved search performance by distributing samples more uniformly across the solution space. EACO reduced the number of objective function evaluations compared with ACO while retaining the same top-ranked solvents and improving some lower-ranked candidates. ESA also identified better lower-ranked candidates than SA, although its reduction in computational cost was modest. Overall, ACO- and EACO-based approaches were substantially more computationally efficient than SA-based approaches for generating a list of candidate solvents, whereas SA and ESA provided improved diversity among lower-ranked solutions. These results show that quasi-random sampling can enhance metaheuristic solvent design and that the proposed CAMD framework is robust for complex nutraceutical extraction problems.