Design of Flotation Circuits Using Tabu-Search Algorithms: Multispecies, Equipment Design, and Profitability Parameters

The design of a flotation circuit based on optimization techniques requires a superstructure for representing a set of alternatives, a mathematical model for modeling the alternatives, and an optimization technique for solving the problem. The optimization techniques are classified into exact and approximate methods. The first has been widely used. However, the probability of finding an optimal solution decreases when the problem size increases. Genetic algorithms have been the approximate method used for designing flotation circuits when the studied problems were small. The Tabu-search algorithm (TSA) is an approximate method used for solving combinatorial optimization problems. This algorithm is an adaptive procedure that has the ability to employ many other methods. The TSA uses short-term memory to prevent the algorithm from being trapped in cycles. The TSA has many practical advantages but has not been used for designing flotation circuits. We propose using the TSA for solving the flotation circuit design problem. The TSA implemented in this work applies diversification and intensification strategies: diversification is used for exploring new regions, and intensification for exploring regions close to a good solution. Four cases were analyzed to demonstrate the applicability of the algorithm: different objective function, different mathematical models, and a benchmarking between TSA and Baron solver. The results indicate that the developed algorithm presents the ability to converge to a solution optimal or near optimal for a complex combination of requirements and constraints, whereas other methods do not. TSA and the Baron solver provide similar designs, but TSA is faster. We conclude that the developed TSA could be useful in the design of full-scale concentration circuits.


Introduction
Froth flotation is a process used in mining, based on the different surface properties of ore, for separating the valuable mineral from gangue. This process is performed in an aerated tank where the ore is mixed with water and reagents to render the valuable mineral hydrophobic [1]. Due to the complexity of the process in practice, multiple interconnected stages are used to form a flotation circuit.
The design of flotation circuits using optimization techniques has been widely studied in the literature [2,3]. The design requires three elements: (1) defining superstructures for representing a set of alternatives for design; (2) a mathematical model for modeling the different alternatives of the design, here considered goals and constraints, and determining the objectives to be optimized; and (3) an optimization technique for solving the problem. size increases [5]. In this context, according to Acosta-Flores et al. [34] and Cisternas et al. [2], the genetic algorithms are the only metaheuristic algorithms that have been used for designing flotation circuits [1,[35][36][37][38]. However, GA has been used only for small problems because its convergence is slow.
The methodologies proposed in the literature considered different assumptions for simplifying the design problem. For example, many of these methodologies considered a maximum of six flotation banks and a maximum of eight cells; however, these studies usually considered only two species in the fed ore to the flotation circuit. The exceptions are the works of Calisaya et al. [33] and Acosta-Flores et al. [34]. These authors examined several species in the fed ore, but their works were based on the fact that there are few structures that are optimal for a given problem [39]. These few structures were identified before applying the optimization search, which reduces the size of the optimization problem. Therefore, the studied problem is complex and difficult to solve when obtaining a global solution.
In this study, we propose using the TSA for designing flotation circuits. The developed algorithm incorporates diversification and intensification, the first of which is used for exploring new regions, and the second, for exploring regions close to a good solution. The algorithm was implemented for designing circuits that process several species. The superstructure implemented involves five flotation banks and each bank could use 3-15 cells, and the objective functions are economical. Four case studies were developed for illustrating the applicability of the algorithm.

Superstructure
The proposed superstructures in the literature usually correspond to an equipment superstructure, which allows the streams of concentrate and tail of one equipment to be sent to any other equipment. The superstructures are important because they define the alternatives and the size of the problem [2]. According to Sepúlveda et al. [40], the circuits generally use between three and five flotation stages. Then, the equipment superstructure ( Figure 1) used in this work considers the following stages: rougher stage (R), cleaner stage (C1), re-cleaner stage (C2), scavenger stage (S1), and re-scavenger stage (S2).

Mathematical Model
The model includes mass balances in flotation stages, splitters, and mixers, and goals, constraints, and objective function are considered. The mass balance in flotation stages is determined with: Many of the proposed superstructures in the literature considered nonsensical alternatives and/or presented degeneracy, i.e., alternatives that are equivalent [2]. All the design alternatives shown in Figure 1 have sense and do not present degeneracy. In Figure 1, the triangles with label M represent mixers of the feed that arrive at stage i, where i ∈ {1, 2, 3, 4, 5} = I. Note that the numbers 1, 2, 3, 4, and 5 are related to R, C 1 , C 2 , S 1 , and S 2 , respectively. The triangles with label D represent splitters that allow sending the concentrate and tail streams from stage i to the other stages of the superstructure.

Mathematical Model
The model includes mass balances in flotation stages, splitters, and mixers, and goals, constraints, and objective function are considered. The mass balance in flotation stages is determined with: where F ik is the mass flow of the species k fed to the flotation stage i, C ik is the mass flow of the species k in the concentrate stream C i of the flotation stage i, T ik is the mass flow of the species k in the tail stream T i of the flotation stage i, and R ik is the recovery of the species k in the flotation stage i, where k = 1, 2, . . . , m, i ∈ I, and m is the number of species. The flotation model used for representing the recovery in the flotation stages was proposed by Yianatos and Henríquez [41]: where R ik is the recovery of the species k in the flotation stage i, k max,i,k is the maximum rate constant of the species k in the flotation stage i, τ i is the cell residence time in the flotation stage i, N i is the number of flotation cells used in the flotation stage i, and R max,i,k is the maximum recovery at the infinite time of the species k in the flotation stage i. In practice, the flotation circuits do not use stream branching because a large number of pumps and junction boxes are necessary for carrying out the concentration process [37]. Then, the mass balances in the splitters of the concentrate streams are expressed as: where α ij ∈ {0, 1} are decision variables indicating the destination of the concentrate stream ( Figure 1), C ijk is the mass flow of species k in the concentrate stream from stage i to stage j, and C ik is the mass flow of species k in the concentrate stream of stage i. Similarly, for the tail streams: where β ij ∈ {0, 1} are decision variables indicating the destination of the tail stream (Figure 1), T ijk is the mass flow of species k in the tail stream from stage i to stage j, and T ik is the mass flow of species k in the tail stream of stage i. The mass balances in mixers are expressed as: where M 1k is the mass flow of the species k fed to the flotation circuit. Note that mass balance for each species k processed in the circuit can be rewritten as a matrix system. These systems are solved numerically for C ik and T ik using a linear programming algorithm. The TSA has the ability to make use of this type of algorithm. The market for copper concentrate establishes a minimum grade (g LO ); then, the following equation is included in the mathematical model: where g k,cu is the copper grade of the species k, and gradeCu is the copper grade in the final concentrate of flotation circuit. In this work, the value of g LO is 0.25. Several objective functions have been used in the literature. Mehrotra and Kapur [25], Green [42], and Pirouzan et al. [38] implemented technical expressions for designing flotation circuits. The technical functions are difficult to define. For example, if the recovery is maximized, it is necessary to restrict the concentration grade to a value that is not known a priori. This difficulty was overcome by some authors using a multi-objective function; however, it is difficult to define the relative weight of each objective. Schena et al. [29] and Cisternas et al. [32,39,43] implemented economic expressions for designing circuits. These expressions highlight the maximization of revenues, maximization of profits, and the maximization of the net present worth. These last two expressions require estimating both the equipment costs and the operating costs. Cisternas et al. [43] found that the objective function has a significant effect on both the solution and circuit structure obtained. In this work, we used the maximization of revenues and maximization of the net present worth as the objective functions.
The revenue can be calculated using different models depending on the type of product and its market. In the case of copper concentrate, the net-smelter-return formula can be used [32,43]: where p is the fraction of metal paid, µ is the grade deduction, Trc is the treatment charge, q is the metal price, R f c is the refinery charge, CF k is the mass flow of the species k in the final concentrate, and g k,cu is the copper grade of the species k. For determining net present worth, the capital costs and total costs of the process must first be estimated. The capital cost considers the fixed capital and working capital. The first is estimated with the following equation: where [44]: with where V i is the cell volume (m 3 ) of flotation stage i, F i is the feed stream to stage i, E g is the gas factor, ρ p is the pulp density, I F,i is the fixed capital cost to stage i, I F is the fixed capital cost of circuit, and FL is the Lang factor [45]. Equation (10) is valid for volume between 5 m 3 and 200 m 3 . The working capital costs are estimated with the following equation: where FL w is the Lang factor for working capital and is assumed to be 0.9. The total costs of the process are estimated with the following equation: where C op,i is the operating cost of flotation stage i, P k is the kilowatt-hours cost, E g is the gas factor, and MCM are the costs of mine-crushing-grinding per ton of fed ore to the flotation plant. The profits generated by the project are estimated with: The annual cash flows are estimated using the following expression: with where P B are the profits before taxes, r t is the tax rate, D is the annual depreciation, and n is the life time of the project (the value of n used here is 15). The net present worth is calculated using: where W NP is the net present worth, I cap = I F + I w , and r d is the discount rate. It is assumed that cash flows are equal in all years of the project, so Equation (19) is simplified: with where γ denotes the penalty parameter [46]. For each violation of the constraints of the mathematical model, a value v i > 0 is considered, defined by the user. This penalty parameter must also be considered when the objective function is the maximization of revenues.

Optimization Technique: Tabu-Search Algorithm
TSA is a local search methodology that was proposed by Glover and Laguna in 1998 [15]. TSA is a strategy for solving combinatorial optimization problems ranging from graph theory to mixed integer programming problems. It is an adaptive procedure with the ability to making use of other methods, such as linear programming algorithms, which help to overcome the limitations of local optimality [16]. Gogna and Tayal [47] stated that the TSA uses the information gathered during the iterations to produce a more efficient search process. Here, the search space is simply the space of all possible solutions that can be considered during the search. For example, in vehicle routing problems, the search space considers both binary and continuous variables [48]. The TSA accepts non-improving solutions to the global solution to move out of local optima. The distinguishing feature of the TSA is the use of memory structures.
The main memory structure used by the TSA is the Tabu list (TL) short-term memory, which has a record of previously visited solutions. This key idea can be linked to artificial intelligence concepts [48]. The TL should be carefully formulated for an effective search while minimizing the computation time and the memory requirements. Medium-and long-term memories can be used for improving the intensification and diversification of the TSA [4].
Usually, the TSA starts with an initial solution, randomly selected, which is entered to TL. Then the algorithm uses a local search procedure or neighborhoods to move iteratively from a potential solution x to an improved solution x , also called the best neighbor, in the neighborhood of x (N(x)). Local search procedures often become stuck in local optima. In order to avoid these pitfalls and to explore regions of search space that would be left unexplored by other local search procedures, TSA carefully explores the neighborhood of each solution as the search progresses [47]. The solutions admitted to the new neighborhood are determined through the use of memory structures. The best neighbor (x best ) of N(x) is accepted as a global solution if x best / ∈ TL and if it maximizes the objective function. If the best neighbor x best ∈ TL, then the next best neighbor of N(x) is the new postulant to enter into TL. This procedure is repeated until x best is entered into T. Next, the new neighborhood of the best neighbor is generated, and the procedure described earlier is repeated until some stopping criterion is satisfied [49].
The search space of the design problem studied in this work considers binary, discrete, and continuous variables, which are related to circuit structure, the number of cells in flotation stages, and the operating conditions in flotation stages, respectively. The developed TSA implements short-term and long-term memories: the TL and the frequency matrix (FM), respectively. The latter was proposed by De los Cobos [50], which allows the exploration of new regions of the search space.
In our case, the algorithm starts with an initial solution ( Here, the variables α ij , β ij i,j are associated with the circuit structure, and (N t , τ t ) t are associated with number of equipment and operating conditions of circuit. Then, the neighborhood N(x) of the initial solution x is generated to create natural permutations of the structural variables, i.e., a set of structures is generated. Subsequently, the operating conditions and equipment number are assigned to each structure through a uniform distribution function. The uniform distribution function is defined using the operating conditions and the equipment number of the initial solution. This method of generating neighborhoods is based on the structure being more influential on the objective function than the operating conditions and equipment number [43]. Subsequently, the equipment size, copper grade of the concentrate, and profitability parameters of each flotation circuit, neighbor of N(x), are determined via mass balances. If the constraints of the mathematical model are violated for some neighbor of N(x), then its objective function ( f ) is penalized (γ). Then, the best neighbor (x best ) of N(x) is determined, i.e., the neighbor maximizing the objective function. The best solution (x best ) is accepted as a global solution if x best / ∈ TL and f x global < f (x best ). If x best ∈ TL, the next best neighbor of N(x) is taken as x best . This step is repeated until x best is entered to TL. The structural variables of x best are entered into FM. Subsequently, the neighborhood of x best is generated and the search procedure described earlier is repeated Iteration max times. Notably, TL has a determined number of rows (n r ), i.e., once TL is full, the update of its information is carried out at each iteration of the algorithm. Next, we explain how diversification and intensification strategies are included in the search procedure. Each time that x global does not improve before D max iterations of algorithm, then the diversification in the search procedure is incorporated. The diversification allows the generation of a new best neighbor, whose structural variables are obtained from the gathered information in FM, and operating conditions and equipment number are assigned through a uniform distribution function. Note that FM records all the structural variables of x best from the first iteration of algorithm. This information allows us to determine the structures more often visited (high frequency in FM) by the algorithm. A structure not visited, or rarely visited (low frequency in FM), is assigned to the new best neighbor. The intensification is incorporated in the search procedure each time passing I max iterations of the algorithm. The intensification allows the exploration of regions close to a good solution. In this work, this is carried out using the gathered information in TL. The new best neighbor is aleatorily selected among the better neighbors recorded in TL. The full procedure is shown in Figure 2.

Applications
Four case studies were developed for illustrating the proposed algorithm. The first and second cases involve designing a copper ore concentrator plant considering the maximization of revenues and maximization of the net present worth as the objective function, respectively. The third case analyses a benchmarking between the Tabu-search algorithm and the Baron solver. Finally, the fourth case provides the comparison between our approach and the methodology proposed by Acosta-Flores et al. [34].

Maximization of Revenues
The equipment superstructure, the mathematical model, included maximization of revenues as the objective function, and the TSA were used for solving the design problem. The developed algorithm was capable of solving the design problem despite the number of species processed, and the requirements and constraints established in the design procedure. The obtained structure is shown in Figure 3, and the revenue, the net present worth, profit before taxes, total capital investment, and total costs of the circuit were USD $130,960,079/year, USD $595,475,455, USD $91,716,855/year, USD $53,265,087, and USD $36,358,032/year, respectively. The The number of cells and residence time in rougher, scavenger, and re-scavenger are the maxima available. These results are logical because the metallurgical aim of these stages is increasing the recovery of the valuable ore. Also, the objective function does not consider either the equipment costs or the operating costs. The number of cells and residence time in the cleaner and recleaner stages are the minima available. Again, these results are logical because the aim of these stages is increasing the copper grade in the concentrate. These results were obtained by previously evaluating different combinations of the algorithm parameters. In this case, the TSA used 170 neighbors in each neighborhood, 2000 iterations, diversification was applied each time the global solution did not improve after 30 iterations, the intensification was applied each time passing 50 iterations of the algorithm, and the number of rows of the TL was 50. The runtime of the TSA was 601.63 s. The number of cells and residence time in rougher, scavenger, and re-scavenger are the maxima available. These results are logical because the metallurgical aim of these stages is increasing the recovery of the valuable ore. Also, the objective function does not consider either the equipment costs or the operating costs. The number of cells and residence time in the cleaner and recleaner stages are the minima available. Again, these results are logical because the aim of these stages is increasing the copper grade in the concentrate. These results were obtained by previously evaluating different combinations of the algorithm parameters. In this case, the TSA used 170 neighbors in each neighborhood, 2000 iterations, diversification was applied each time the global solution did not improve after 30 iterations, the intensification was applied each time passing 50 iterations of the algorithm, and the number of rows of the TL was 50. The runtime of the TSA was 601.63 s.

Maximization of the Net Present Worth
The equipment superstructure and mathematical model included maximization of the net present worth as the objective function, and the TSA was used for solving the design problem. Again, the developed algorithm was capable of solving the design problem despite the number of species processed, and the requirements and constraints established in the design procedure. The obtained structure is shown in Figure 4, and the net present worth, revenue, profit before taxes, total capital investment, and total costs of the circuit were USD $724,828,320

Maximization of the Net Present Worth
The equipment superstructure and mathematical model included maximization of the net present worth as the objective function, and the TSA was used for solving the design problem. Again, the developed algorithm was capable of solving the design problem despite the number of species processed, and the requirements and constraints established in the design procedure. The obtained structure is shown in Figure 4, and the net present worth, revenue, profit before taxes, total capital investment, and total costs of the circuit were USD $724,828,320, USD $129,830,340/year, USD $107,977,400/year, USD $ 8,386,189, and USD $21,398,690/year, respectively. The values of τ R , τ C1 , τ C2 , τ S1 , τ S2 , N R , N C1 , N C2 , N S1 , N S2 , V R , V C1 , V C2 , V S1 , and V S2 were 3 min, 3

Benchmarking between the Tabu-Search Algorithm and the Baron Solver
Many authors have used the Baron solver to solve design problems. As such, we completed benchmarking between the TSA and the Baron solver. The mathematical model for the Baron solver appears in Appendix A. Note that the Baron solver is based on the branch and cut algorithm, i.e., belongs to the family of exact methods.
Initially, we completed benchmarking using the maximization of the revenues as the objective function; the results are shown in Table 3 and Figures 3-5. Then, we performed benchmarking using the maximization of the net present worth as the objective function, and the results are shown in Table 4 and Figures 3-5. Many authors only used the revenues for designing the circuit, i.e., they included neither equipment design nor operational costs in the mathematical model [1,31,32,34,51]. Thus, we performed benchmarking using a simplified mathematical model, and the results are shown in Table 5 and Figures 3-5. Note that the superstructure of Figure 1 represents a total of 144 flotation circuit configurations and that only 11 configurations were selected in all the analyzed cases. These structures correspond to 4 for the case maximization of the revenues as the objective function, 4 for the case the maximization of the net present worth as the objective, and 3 when a simplified mathematical model is used. The structures of Figures 3 and 4 were the most frequently selected circuits. Table 3 shows the benchmarking between the TSA and the Baron solver when the maximization of revenues was the objective function. This table depicts five cases. The first considered the species Cpf, Cps, S, and G; the second considered the species Cpf, Cps, P, S, and G; the third considered the species Cpf, Cps, Cf, P, S, and G; the fourth considered the species Cpf, Cps, Cf, Cs, P, S, and G; and the fifth considered the species Cpf, Cps, Cf, Cs, Bof (bornite fast), P, S, and G. In each case, the following output variables were considered: the revenues, net present worth, profit before taxes, total capital investment, total annual cost, equipment size, operating conditions, number of equipment, copper grade in concentrate, circuit structure, runtime of algorithm, number of iterations of TSA, neighborhood size, number of iterations of each diversification, number of iterations of each intensification, and number of rows of TL. Table 3 shows that TSA is capable to converge independently of the number of species processed in the flotation circuit. The results provided by the TSA and the Baron solver were similar in all analysed cases, except in the fifth case. In this last case, both optimization techniques provided similar revenue; however, the other output variables were different. This result could be related to the simultaneous imposition of all goals and constraints of the mathematical model. Perhaps a relaxation of constraints could help with the exploration of regions, offering a better quality solution. The runtime of TLA in the first, second, third, fourth, and fifth cases was 21.3, 183.4, 287.7, 344.7, and 386.2 s, respectively. The runtime for the Baron solver in the first, second, third, fourth, and fifth cases was 233.0, 423.0, 9026.2, 14,640.0, and 10,257.0 s, respectively. The TSA provided results faster than the Baron solver, but only provides a good quality solution.                        Table 3 shows the benchmarking between the TSA and the Baron solver when the maximization of revenues was the objective function. This table depicts five cases. The first considered the species Cpf, Cps, S, and G; the second considered the species Cpf, Cps, P, S, and G; the third considered the species Cpf, Cps, Cf, P, S, and G; the fourth considered the species Cpf, Cps, Cf, Cs, P, S, and G; and C2 S2 S1  Table 4 shows the benchmarking between TSA and the Baron solver when the maximization of the net present worth was the objective function. This table demonstrates six cases. The first considered the species Cpf, Cps, S, and G; the second considered the species Cpf, Cps, P, S, and G; the third considered the species Cpf, Cps, Cf, P, S, and G; the fourth considered the species Cpf, Cps, Cf, Cs, P, S, and G; the fifth considered the species Cpf, Cps, Cf, Cs, Bof, P, S, and G; and the sixth considered the species Cpf, Cps, Cf, Cs, Bof, Bos (bornite slow), P, S, and G. In each case, the following output variables were considered: net present worth, revenues, profit before taxes, total capital investment, total annual cost, equipment size, operating conditions, number of equipment, copper grade in concentrate, circuit structure, runtime of algorithm, number of iterations of TSA, neighborhood size, number of iterations of each diversification, number of iterations of each intensification, and number of rows of TL. Table 4 shows that TSA is capable to converge independently of the number of species processed in the flotation circuit and objective function used. The results provided by both optimization techniques were similar in all analysed cases. The runtime of TLA in the first, second, third, fourth, fifth, and sixth cases was 106.260 secs, 424.72, 464.99, 527.61, 710.04, and 875.98 s, respectively. The runtime for the Baron solver in the first, second, third, fourth, fifth, and sixth cases was 94.2, 294.92, 533.40, 355.20, 1082.64, and 4299.86 s, respectively. In the first, second, and fourth cases, the Baron solver provided results before the TSA. In the third, fifth, and sixth cases, the TSA provided results before the Baron solver. Table 5 shows the benchmarking between the TSA and the Baron solver when the mathematical model is simplified. This table demonstrates six cases. The first considered the species Cpf, S and G; the second considered the species Cpf, Cps, S, and G; the third considered the species Cpf, Cps, Pf (pyrite fast), S, and G; the fourth considered the species Cpf, Cps, Pf, Ps (pyrite slow), S, and G; the fifth considered the species Cpf, Cps, Cf, Pf, Ps, S, and G; and the sixth considered the species Cpf, Cps, Cf, Cs, Pf, Ps, S, and G. In each case, are the following output variables were considered: revenue, operating conditions, number of equipment, copper grade in concentrate, circuit structure, runtime of algorithm, number of iterations of TSA, neighborhood size, number of iterations of each diversification, number of iterations of each intensification, and number of rows of TL. The runtime of TLA in the first, second, third, fourth, and fifth cases was 162.36, 344.06, 368.83, 480.36, and 1098.57 s, respectively. The runtime for the Baron solver in the first case was 1961.8 s, and in the second to fifth cases, the Baron solver did not converge after five days. Both optimization techniques do not provide a solution for the sixth case because the minimal copper grade constraint in the final concentrate cannot be satisfied. Note, the results provided by the TSA could be used for reducing the runtime of the Baron solver. When we delimited the search space of the Baron solver in the fifth case in Table 5, based on the TSA results, runtime of the Baron solver was 101.22 s, and the output variables of mathematical model were similar to those provided by the TSA.

Comparison with Another Approach
Section 3.3 showed that the Baron solver could converge after a long time period depending on the mathematical model and number of species used in the design procedure. Maybe, due to these results, Acosta-Flores et al. [34] proposed a design methodology based on two phases. In the first phase, they identified a set of optimal structures using discrete values for the flotation stages, then, in the second phase, they determined the optimal design of each structure obtained in the previous phase. Note that they used a superstructure of six flotation stages. However, the optimal structures only used five stages (R, C1, C2, S1, and S2). These authors modeled the flotation stages using the expressions proposed by Yianatos et al. [41]. When all parameters used by Acosta-Flores et al. [34] were included in the proposed methodology, the design shown in Table 6 was obtained. Considering that versions of the Baron solver in the General Algebraic Modeling System (GAMS) environment improve over time, we determined the final design of the structures proposed by Acosta-Flores et al. [34] using our version of GAMS. The results are shown in Table 6.   Table 6 shows that TSA is capable of converging despite changing the parameters in the methodology proposed. The best design provided by both approaches was similar. However, our method required fewer computational resources. The TSA used 320 neighbors at each neighborhood, 1500 iterations, the diversification was applied each time that the global solution was not improved after 20 iterations; the intensification was applied each time passing 40 iterations of the algorithm, and the number of rows of the TL was 70. Note that the TSA provided secondary designs (Table 7). However, neither approach guarantees finding a global solution.   In general, Tables 3-6 show that the TSA converges faster than the Baron solver, and, the TSA always provided a solution in a reasonable amount of time when the mathematical model was well-defined. Notably, the TSA parameters must be adjusted. Evidently, this procedure took time. For example, Lin et al. [52] proposed determining the neighborhood size using the number of free variables of a mathematical model. A similar strategy was applied in this work; however, in general, the TSA did not provide good results. This behavior could be related to the considered multispecies in the design procedure. Then, the neighborhood size was tuned by trial and error.
Cisternas et al. [2] and Lin et al. [5] mentioned that both exact and approximate methods have a low probability of finding the global solution in a reasonable amount of time for large. However, our results contradict these observations. Tables 3 and 4 show that both optimization techniques converge and provided similar designs despite the design procedure considering several species, operating costs, capital costs, size of equipment, number of equipment, and metallurgical parameters, among other variables. In the case of the TSA, these results could be explained by a minimally dynamic and noncomplex superstructure, which did reduce the number of alternatives to 144. In the case of the Baron solver, these results could be explained by the delimitation of the search space due to large number of equations used by the mathematical model, which helped the solver to converge in a reasonable amount of time. This finding is corroborated by the results shown in Table 5. In this case, the mathematical model included neither equipment design nor operational cost, and generally, the solver did not converge.
The results indicate that the diversification and neighborhood have a strong influence on the results of TSA. When the diversification was not incorporated into the search procedure, the results of the TSA were distinct compared to the results provided by the Baron solver. Although a large neighborhood did not help with obtaining good solutions, better solutions were obtained using small neighborhoods and many iterations of the TSA, i.e., achieving a gradual improvement in the best neighbor.
The TSA developed in this work is flexible, i.e., would allow including into the design procedure grinding models, cell models, or the geological uncertainty (non-linear algorithm). The latter is related to the variation of copper grade that occurs during real operation of the circuit. Initially, the consideration of geological uncertainty or degree of liberation of valuable minerals (grinding) drives optimization under uncertainty. There are several approaches for addressing this type of problem such as stochastic optimization [53]. A possible strategy for solving the stochastic optimization problem is scenario trees based on stochastic parameters of the model. The solution of each scenario would be determined via TSA. A summary of the advantages and disadvantages of the optimization techniques used in this work is shown in Table 8. Table 8. Comparison between the Tabu-search algorithm and the Baron solver.

Algorithm
Tabu Search Baron Solver

Advantages
The convergence is fast.
The algorithm always provides a solution when the mathematical model is well defined. The algorithm is flexible, i.e., it allows the use of other methods, such as linear programming algorithms.
The algorithm provides good quality solutions.
The algorithm provides secondary designs.
The solver provides a global optimal design when converged. The solver does not need to adjust parameters for providing a solution.
The obtained designs do not change if the solver is run again.

Disadvantages
Some algorithm parameters, such as neighborhood size and number of iterations of the algorithm, among others, must be adjusted for finding a good quality solution. Penalty parameters must be used for satisfying the constraints of the mathematical model. The obtained designs could change if the algorithm is run again. The algorithm must incorporate diversification and intensification for finding a good quality solution.
Depending on the mathematical model and the number of species, the convergence is slow or the algorithm does not converge.
The variables of the model must be bounded for guaranteeing the finding of global optimal design, i.e., experience in circuit design is required.
The solver provides a single design. The obtained designs depend on the version of the solver.

Discussion and Conclusions
An algorithm based on Tabu-search was developed and used for designing flotation circuits for several species. The algorithm incorporates diversification for exploring new regions and intensification for exploring regions close to a good solution. The circuits designed using the Tabu-search algorithm were logical and allowed incorporating the influence of the objective function into the design of concentration plants. A comparison between the Tabu-search algorithm and the Baron solver in the GAMS environment was performed. In general, the solutions provided by both optimization techniques were similar. The Tabu-search algorithm always quickly provided a solution when the mathematical model was well-defined, whereas the Baron solver, in some cases, converged slowly or did not converge at all. Finally, we compared our approach and the methodology proposed by Acosta-Flores et al. [34]. The results indicated the best design provided by both proposals was similar. Both approaches provided secondary designs, but our proposal required fewer computational resources. Thus, the developed algorithm can converge to an optimal solution or near optimal for a complex combination of requirements and constraints in a design problem, whereas other methods do not. The developed algorithm represents progress in the design of flotation circuits, which could be useful in the design of full-scale mineral concentration circuits.
In future work, we will seek to include geological uncertainty, operating uncertainty, and grinding in the design procedure. Furthermore, we will analyze the hybridization between TSA and genetic algorithms.

Conflicts of Interest:
The authors declare no conflict of interest.

C1
Cleaner stage C2 Re-cleaner stage C i Concentrate stream of the flotation stage i C ik Mass flow of species k in concentrate C i C ijk Mass flow of species k in the concentrate stream from stage j to stage i C op,i Operating cost of flotation stage i D Annual  For the mass balance in splitters, the following equations are considered: CC(r, k1) = ∑ s F c (r, s, k 1 ), (r, s) ∈ LC, k 1 ∈ K 1 , r ∈ M2 CT(r, k1) = ∑ s F w (r, s, k 1 ), (r, s) ∈ LT, k 1 ∈ K 1 , r ∈ M2 where F c (r, s, k 1 ) is the concentrate flows of species k 1 in the stream from r to s, F w (r, s, k 1 ) is tail flow of species k 1 from r to s, CC(r, k 1 ) and CT(r, k 1 ) are the mass flow of species k 1 in the concentrate and tail streams of stage r, respectively. Stream branching is not allowed, so the following equations are included: where FS(s, k 1 ) is the mass flow of species k 1 in feed streams to stage s. The final concentrate of the flotation circuits is calculated using: CF(k 1 ) = ∑ r F c (r, P, k 1 ), (r, P) ∈ LC, k 1 ∈ K 1 where CF(k 1 ) is the mass flow of species k 1 in final concentrate. The mass balance in the flotation stages is determined using: CC(r, k 1 ) = T(r, k 1 )·FS(s, k 1 ), k 1 ∈ K 1 , r ∈ M2 CT(r, k 1 ) = (1 − T(r, k 1 ))·FS(s, k 1 ), k 1 ∈ K 1 , r ∈ M2 where T(r, k 1 ) is the flotation model proposed by Yianatos and Henríquez [41]. The objective function is defined using the Equations (8)- (20), without considering the penalty parameter.