D.C. programming approach for solving an applied ore-processing problem

This paper was motivated by a practical optimization problem formulated at the Erdenet Mining Corporation (Mongolia). By solving an identification problem for a chosen design of experiment we developed a quadratic model that quite adequately represents the experimental data. The problem obtained turned out to be the indefinite quadratic program, which we solved by applying the global search theory for a d.c. programming developed by A.S. Strekalovsky [ 13 ]-[ 15 ]. According to this d.c. optimization theory, we performed a local search that takes into account the structure of the problem in question, and constructed procedures of escaping critical points provided by the local search. The algorithms proposed for d.c. programming were verified using a set of test problems as well as a copper content maximization problem arising at the mining factory.

1. Introduction. We consider optimization problems that arise at the ore mining corporation of Mongolia and can be formulated as a quadratic programming with d.c. reduction form. The criteria or objective functions in these problems are the copper content and the copper recovery in a rougher concentrate.
Quadratic programming frequently arises in optimization of technological processes and design of experiments. It is assumed that the experimenter is concerned with a technological process involving some response f which depends on the input variables x 1 , x 2 , . . . , x n from a given experimental region. The standard assumptions on f are that f is a twice differentiable function on the experimental region and the independent variables x 1 , x 2 , . . . , x n are controlled in the experimental process and measured with a negligible error. The feasible region of variables can be a nonconvex set but in most cases, for simplicity, the experimenter usually restricts himself to the box type of regions. As a rule, the researcher has the following second-order regression or a quadratic model expressed by a quadratic function that adequately represents the experimental data: where coefficients a ij , b j , c, i, j = 1, . . . , n, are assumed to be found by solving an identification problem for a chosen design of the experiment, for example, the orthogonal central composite design [8]. It is required to find the global extremum of the function f (·) over an experimental region or to equivalently reduce the problem to the indefinite quadratic programming over a box constraint [1,2].
As well known [4,7,14], any quadratic program can be formulated as a problem of the d.c. maximization [12]- [15]: where g(·), h(·) are convex functions, and the set Π ∈ IR n is also convex. From now on, let us suppose that the following assumption on the function f (·) over IR n is satisfied: It is known [7,12,13,14] that the A.D. Alexandov's functions (or the (d.c.) functions that can be represented in a form of the difference of two convex functions), form a linear space, which is dense in the space of continuous functions (in the topology of homogeneous convergence on the compacts). Thus, problems of the d.c. programming represent a rather large and, besides, very attractive class of optimization problems, for which A.S. Strekalovsky developed the theory of the global search [12]- [15]. According to the theory based on the global optimality conditions, the process of finding a global solution in nonconvex optimization problems (see [3,14,16,17]) consists of the two principal stages: (i) a special local search, which takes into account the structure of the problem under scrutiny, and (ii) the procedures of escaping from critical points (provided by the local search) based on the global optimality conditions [13,14]. There are many works [5,6,7,10,11,18,19,20] devoted to theory and algorithms for d.c. optimization.
The paper is organized as follows. In section 2, we formulate the identification problem, construct the function f (·) and examine its properties. In section 3, we recall the theoretical basis for solving the d.c. maximization problems and the global search strategy. In the final section, we address applied optimization problems that arise in ore-processing and demonstrate numerical results obtained by the algorithms proposed.
2. Identification problem and model formulations. In general case, the coefficients of the quadratic model can be estimated by the least-squares method based on the following experimental observation data: Then, according to the least-squares method, in order to find the coefficients A = {a ij }, b j , i, j = 1, . . . n, and c, we need to solve the following unconstrained minimization problem: Further, for simplicity, we omit the index l and write down the function F l (A, b, c) as follows Proof.
Then, it is clear that Now let us show that the functions φ t , t = 1, . . . , m, are convex. Further, we omit the index t and denote by y and u the following vectors: y = (a 11 , . . . , a 1n , a 21 , . . . , a 2n , . . . , a n1 , . . . , a nn ) T ∈ IR n×n , u = (x 1 x 1 , x 1 x 2 , . . . , x n x 1 , . . . , x n x n ) T ∈ IR n×n , Compute the Hessian of the function φ(·): Therefore, the Hessian has the form which is obviously a positive semidefinite matrix due to (2). It means that φ t is convex for all t = 1, . . . , m, and, consequently, the function F (A, b, c) is convex. The lemma is proved. Therefore, problem (1) can be solved by the suitable classical convex optimization methods [9] and respective software packages (FICO Xpress, IBM CPLEX, etc.).
In our case as industrial application, we consider the process of collective flotation of copper and molybdenium minerals at the Erdenet Mining Corporation (Mongolia). That is why it is important to maximize the copper content and the copper recovery in the rougher concentrate subject to technological constraints.
To find the parameters a ij , b i , i, j = 1, . . . , n, and c, we processed m = 5000 of the real industrial data, solved two (for l = 1 and l = 2) identification problems (1) with n = 7 employing the software package IBM CPLEX, and constructed two functions The function f 1 (·) approximates the copper content in the rougher concentrate (measured in % of mass); meanwhile, the function f 2 (·) represents the copper recovery in the rougher concentrate (in % of mass). These functions depend on the following variables: x 1 -consumption of collector agent AeroMix, in grams per ton; x 2 -consumption of collector agent VK-901, in grams per ton; x 3 -consumption of foaming agent MIBK, in grams per ton; x 4 -content of -74 micrometer grain class in the hydrocyclone overflow, in % of mass; x 5 -total content of copper in the ore prior to treatment, in % of mass; x 6 -total content of primary copper in the ore prior to treatment, in % of mass; x 7 -total content of oxidized copper in the ore prior to treatment, in % of mass.
Since the matrices A 1 and A 2 have both positive and negative eigenvalues, the corresponding problems (l = 1, 2) turn out to be the nonconvex optimization problems, where The box constraint (4) is the normalized form of the technological requirements for the variables (3). Below we will show how Problems (P 1 ) and (P 2 ) can be represented as the d.c. programs and solved by an algorithm based on the global search theory [12]- [15]. 3. D.C. representation. In order to solve Problems (P l ), we need an explicit d.c. representation of the nonconvex functions f l (·) = h l (x) − g l (x), l = 1, 2. For the sake of simplicity, we omit the index l and describe a simple method [14] of constructing functions h(x) and g(x).
It is well known that any quadratic matrix may be transformed into a symmetrical matrix Q which, in turn, can be represented as the difference of two symmetric positive definite matrices Q = Q 1 −Q 2 . This allows us to get the d.c. representation of the quadratic function f (·) where h(·) and g(·) are strongly convex functions (Q 1 , Q 2 > 0 are positive definite). For example, it can be done in the following way [14]. First, we represent the matrix Q via the difference of two matrices with nonnegative components: Second, we construct the matrices ij is the sum of nondiagonal elements of the row i in the matrix D 1 , and the number ε > 0. Thus, Γ 1 is a positive definite matrix.
Similarly, we obtain Q 1 = Γ 1 +Λ 2 , Q 2 = Γ 2 +Λ 2 , where Λ 2 is a diagonal matrix: ij , the sum of nondiagonal elements of the row i in the matrix Γ 2 . Hence, the matrix Q is represented as the difference Q = Q 1 − Q 2 of matrices Q 1 and Q 2 with non-negative components and dominant diagonals, and we received the d.c. representation (5).
In what follows, according to the global search theory we should focus on finding a local maximizer to the Problem (P ).

Local search.
In order to find a local solution to the Problem (P ), we apply the well-known DC-Algorithm [12,13,14,18]. As known, it consists of linearizing, at a current point, the function h(·) which defines the basic non-convexity of Problem (P ), and minimizing the convex approximation of the goal function f (·) obtained by replacing the non-convex part with its linearization. It is easy to see that the algorithm constructed in this way provides critical points by employing only tools of the convex analysis. Therefore, we start with an initial point x 0 ∈ IR n . Suppose a point x s ∈ Π is provided. Then, we find x s+1 ∈ Π as an approximate solution to the linearized problem It means that the next iteration x s+1 satisfies the following inequality where the sequence {δ s } fulfils the following conditions δ s ≥ 0, s = 0, 1, 2, . . . ; ∞ s=0 δ s < ∞.
Note that the linearized Problem (P L s ) turned out to be convex, meanwhile Problem (P ) was a nonconvex one.
As it was suggested in [12,13,14], one of the following inequalities can be employed as a stopping criterion for the local search method: where τ is a given accuracy. If one of the inequalities (7) is fulfilled, it can be easily shown that the point x s turns out to be a critical point to Problem (P ) with the accuracy τ and under the condition δ s ≤ τ 2 . Indeed, (7) together with the inequality (6) imply that Therefore, if δ s ≤ τ 2 , the point x s is a τ -solution to Problem (P L s ).
In the next section we show how to escape from critical points provided by the local search method. 5. Optimality conditions and the global search strategy. Let us recall the fundamental result of the Global Search Theory.
Theorem 5.1. [13,14,15] Suppose that ∃v ∈ Π : Then, a point z ∈ Π is a global solution to Problem (P ) if and only if As we can see, the verifying condition (E) for a given y requires solving the convex program (P L(y)): depending on 'perturbation' parameters (y, β) satisfying h(y) = β + ζ. According to Theorem 5.1, in order to conclude whether a given point z ∈ Π is a global solution to Problem (P ) or not, we need to solve a family of linearized problems (8) by one of the well known convex optimization methods [9].
On the other hand, we can see that if the condition (E) is violated at a given tuple (ỹ,β, u), u ∈ Π g(u) −β < ∇h(ỹ), u −ỹ , due to convexity of h(·), then we get g(u) <β + h(u) − h(ỹ) and conclude that z ∈ Π is not optimal. Moreover, on each level ζ k = f (z k ), it is not necessary to investigate all pairs of (y, β) satisfying (E), ζ k = h(y) − β, but it is sufficient to discover the violation of the variational inequality (E) only for one pair (ỹ,β) and u ∈ Π.
The properties of the Optimality Conditions (E) allow one to construct an algorithm for solving the d.c. maximization problems. The algorithm comprises two principal stages: a) the local search, which provides for an approximately critical point z k with the value corresponding to the goal function ζ k = f (z k ); b) procedures of escaping from critical points, which are based on the Optimality Conditions (E).
Global Search Scheme.
Step 1. By using the local search method find a critical point z k in Problem (P ).
Step 2. Choose a number β : inf(g, Π) ≤ β ≤ sup(g, Π). Choose an initial Step 3. Construct a finite approximation Step 4. Find a δ k -solutionū i of the following Linearized Problem: Step 5. Starting from the pointū i , find a local maximizer u i by the local search method.
Step 6. Find a δ k -solution w i (h(w i ) = β − ζ k ) of the level problem, i.e.
Step 7. Compute the value η k (β) : Step 8. If η k (β) > 0 go to Step 1 with u j , a new starting point for the local search.
Step 9. Otherwise, choose a new value of β and go to Step 3.

6.
Approximation of a level surface. One of the principal features of the Global Search Algorithm is an approximation of the level surface of the convex function h(·) which generates the basic nonconvexity in Problem (P ) (Step 3). In particular, an approximation R k (β) of the level surface h(x) = β + ζ for each pair (β, ζ k ), ζ k = f (z k ) can be constructed by the following rule where e i is the unit vector from the Euclidean basis of IR n . The search of µ i turns out to be rather simple and, moreover, analytical (i.e. it reduces itself to the solution of a quadratic equation of one variable) for the quadratic function. When h(x) = Q 1 x, x the number µ i for each i = 1, . . . , n, is computed by the following formula The set (approximation) (9) has proven to be rather competitive [14,16,17] during the computational simulations.   ). 0.408, 1.000, 0.572, 1.000, 0.628, 1.000, 0.167 The experimental results showed that beginning with the various starting points, the LSM delivered quite a varying set of local solutions to Problems (P 1 ) and (P 2 ). At the same time, there exist the values f 1 (z 5 ) = 1.365 and f 2 (z 1 ) = f 2 (z 3 ) = 1.1013 (bold font in the tables) that might be the best known solution to the Problems (P 1 ) and (P 2 ), respectively.

7.2.
Testing the global search method. On the basis of the Global Search Scheme from Subsection 5, we developed a global search algorithm (GSA) for searching a global solution to the Problems (P 1 ) and (P 2 ).
Further, starting from the same points (see Tables 1 and 2), the GSA has found the best known solution to Problems (P 1 ) and (P 2 ). The results of computational simulations are presented in Tables 3 and 4. In addition to the denotations employed in Tables 1 and 2, we denoted by f * l ( l = 1, 2) the value of the function at the point provided by the GSA; loc denotes a number of the LSA runs during the GSA; it stands for a number of the GSA iterations. Thus, the global search algorithm provided us with the global (best-known) solutions to Problems (P l ), l = 1, 2, in a normalized form: 8. Conclusions. We considered the real life problem of maximizing the copper content in the concentrate which arises at the Erdenet Mining Corporation (Mongolia). This problem has been originally formulated as an indefinite quadratic programming over a box constraint [1,2] of technological variables. In general, it is known that such problems are NP-hard. Further, we reduced this problem to a d.c. programming problem, so that we could apply the global optimality conditions developed by A.S. Strekalovsky [12]- [15] and a corresponding algorithm proposed in [13]- [17]. The global (best-known) solution provided by the algorithm meets the technological requirements given by the Erdenet Mining Corporation.