Operations Research Letters 15 (1994) 73-79
ELSEVIER
An algorithm for solving general D.C. programming problems Thai Quynh
Phong
LM1 - INSA Rouen. B.P. 08, 76131, Mt St Aignan Cedex. France
(Received 20 January 1993; revised 31 March 1993)
Abstract We consider the problem of globally minimizing a d.c. (difference of two convex functions) function under d.c. constraints. Under a certain condition the problem can be reduced to an unconstrained d. c. minimization problem. A simple outer approximation algorithm is then proposed. Key words." Nonlinear; d.c. programming; Global optimization; Outer approximation
1. Problem statement In this paper we consider the following global optimization problem (DCP)
min
f(x)
s.t.
.~(x) 4 0 ,
i = 1. . . . . m,
x~, where f~ is a closed convex set in •" a n d f , f l . . . . . fm are d.c. functions. Recall that a function f : ~" ---, ~ is called a d.c. function if it can be written as the difference of two convex functions. P r o b l e m (DCP) is of interest from a practical as well as from a theoretical viewpoint. A detailed discussion on d.c. functions, d.c. p r o g r a m m i n g and its applications can be found in Horst and T u y [2]. Specially, if fl = Pl - ql a n d f 2 = P2 - q2 are d.c. functions then max{f1 ,f2} = m a x { p l + q2, ql + P2} -- (ql + q2),
(1)
m i n { f l , J 2 } = (Pl + P2) -- m a x { p , + q2, qx + P2},
(2)
i.e. m a x { f l , f 2 } and m i n { f l , f 2 } are also d.c. functions. It follows, in particular, that any system of d.c. constraints f/(x)~< 0, i = 1. . . . . m, (fi(x) is a d.c. function) is equivalent to a single d.c. constraint max f . ( x ) <~ O. Thus, problem (DCP) can be rewritten as (Q)
min s.t.
h(x) - g(x), p(x) -- q(x) ~ O,
where h, g, P, q are finite convex functions on R". 0167-6377/94/$07.00 © 1994 Elsevier Science B.V. All rights reserved SSDI 0 1 6 7 - 6 3 7 7 ( 9 3 ) E 0 0 8 1 - 2
T.Q. Phong / Operations Research Letters 15 (1994) 73-79
74
The aim of this paper is to propose an algorithm for solving general d.c. programming problems of the form (Q). Our method is an analogue of the Relief indicator one by Thach and Tuy for globally minimizing a continuous function f over a compact set S in R". Their approach is centered on the concept of relief indicator O(~, x) such that a feasible point x* is a global optimal solution if and only if 0 = min{¢P(~*, x): x • •"} where ~* =f(x*). Thus the considered problem is reduced to an unconstrained optimization problem. An outer approximation algorithm (see [7]) and a branch and bound algorithm (see [8]) were then proposed for solving the later problem. In general, a relief indicator is defined by mean of the so-called separator for the functionfon the set S. In [7], it has been taken as a d.c. function of the form O(e, x) = h(~, x) - ]l x Ir2
(3)
where h(ct, x) is a certain, implicitly defined, convex function. It should be noted that by definition the function ~(~, x) is not uniquely defined. Moreover, the decomposition of a d.c. into the difference of two convex functions is not unique. This gives a motivation to us to construct a specific relief indicator for a special class of optimization problems. More precisely, we will show in our case that i f f a n d S have the d.c. structure as in (Q) then such a function ~(~, x) can be directly and explicitly derived from the functions h, g, P, q. In fact, the obtained expression is more natural than (3) and well suited to the problem under consideration. This, of course, simplifies the implementation of the algorithm. In Section 2 we will show how problem (Q) can be reduced to an unconstrained d.c. minimization problem under an assumption of robustness. In Section 3 we develop an outer approximation procedure for solving the specific unconstrained d.c. minimization problem. The convergence of the method is established in Section 4. Finally, an illustrative example will be given in the Section 5.
2. Reduction to unconstrained d.c. minimization problem
We assume that the set S = {x e R": p(x) - q(x) <<.0}
is bounded so S is a compact subset of R". If we denote S, = { x • S: h(x) - g(x) < ~},
(4)
S, = {x • S: h(x) - g(x) ~ ~},
(5)
then, it is evident that ~* = min(Q) ¢~ if,. 4:0 and S,. = 0, where min(Q) denotes the optimal value of problem (Q). Now we define the function F(a, x) = max {h(x) - g(x) - ~, p(x) - q(x)}. It is easy to see that F(a, x) = l(e, x) - k ( x ) where l(a, x) = max{g(x) + p(x), h(x) + q(x) - ~}, k(x) = q(x) + g(x),
(6)
T.Q. Phong / Operations Research Letters 15 (1994) 73 79
75
are convex functions and :t > ~' ~ F(~, x) ~< F(c(, x),
(7)
Throughout the paper we assume the robustness of the feasible region, i.e. c l { x e R " : p ( x ) - q(x) < 0} = S,
(9)
where cl A denotes the closure of A. It is not difficult to verify in this case that
S~ :~ 0 ¢¢- Int S, ¢ 0.
(10)
It can be verified that F(~, x) is indeed a relief indicator in the sense of Definition 4.1 ([7]). The following fundamental theorem is similar to that of Tuy [8] (see also [7]). For the sake of convenience we give below a direct proof of the theorem. Theorem 1. 1. ct* --- min(Q)c~min F(c(*, x) = O for all x e ~ " , and any point x* satisfyino 0 = F(c(*, x*) = min{F(ct*, x): x E ~"} is a 9lobal optimal solution of(Q). 2. l f i n f F ( ~ , x) < 0 then any point x such that F(~, x) < 0 satisfies x E S , . 3. l f i n f F ( ~ , x) >~ 0 then h(x) - g(x) > orfor all x ~ S (so ~ = + ~
implies S = 0).
Proof. Let ~* = min(Q), then ;~,, # 0, or equivalently, 3x*: F(~*, x*) ~< 0. Suppose the contrary, i.e., there exists a point x' E ~" such that F(c(*, x') < 0. Then h(x') - 9(x') - ~* < O,
p(x') - q(x') < O,
hence S,. 4: 0, that contradicts (6). Therefore min F(~*, x) = F(a*, x*) = 0. Conversely, if min F(:(*, x) = 0, then there exists a solution x*: F(~*, x*) = 0, which implies x* ~ S,., i.e., g~, :~ O.
(11)
Moreover, Vx E [~": F(~*, x) >~ 0 so Int S~, must be empty. By virtue of (I0), one has S~, = 0.
(12)
This, together with (11), (6) implies min(Q) = ~*. It is easily seen for a point x* satisfying 0 = F(a*, x*) = min{F(~*, x): x e ~"}, one has p(x*) - q(x*) <<,0 and h(x*) - g(x*) = c~*. This means x* is a global optimal solution of (Q). The assertions 2 and 3 are immediate from the definition of function F. The importance of this theorem is that it reduces problem (Q) to the unconstrained minimization problem: min {F(~*, x) = l(~*, x) - k(x): x ~ ~"} where the objective function is a d.c. function. This problem is equivalent to the following concave minimization problem with convex constraint: min {t - k(x): l(~*, x) <~ t, x ~ ~"}, which can be solved by several available methods (cf. Horst and Tuy [2]). Of course, the main difficulty is that e* is not known in advance. In the next section we shall overcome this difficulty by using outer approximation technique.
T.Q. Phong / Operations Research Letters 15 (1994) 73 79
76
3. Description of the method In this section we shall develop an algorithm for solving the following problem: min{t - k(x): 1(5", x) <~ t, x e R"},
(13)
where 5" is the optimal value of problem (Q). From Theorem 1, if 5" < + oo and (x*, t*) is the solution of (13) then x* is a solution of (Q). Problem (13) is a concave minimization problem, but not in conventional form, since the convex constraint /) = {(x, t)~l~ "+1" l(5*,x) ~< t}
(14)
is not known explicitly. In fact, 5" is unknown as long as the problem has not been solved yet. We shall proceed similarly as in Tuy [8]. Namely, at iteration k, the best value 5k SO far obtained is used as an upper estimation of 5", so that the unknown function l(5", x) is approximated by its underestimator
l(Sk, X) = max{h(x) + p(x) - 5k, q(x) + g(x)}.
(15)
Furthermore, if a couple (x k, t k) ~ ~n × ~ is generated, which does not belong to D--~= {(x, t)~ R "+~" l(sk, x)~< t}
(16)
(so does not belong t o / ) ) , one can construct a hyperplane
n k = {(x, t): ( qk, X -- X k) + l(Sk, X k) -- t = 0},
(17)
with qk s M(Sk, Xk), which separates (x k, t k) and Dk. The algorithm can be summarized as follows.
Algorithm • Initialization: Select an n-simplex Mo in R" that is known to contain a global optimal solution. Take a point x ° s Mo and set ~o = x o, 5o = h(x °) - g(x °) if xCe S. Otherwise set 5o = + oo. Compute a subgradient qo of function 1(5o, x) at point x ° and let PI = {(x, t)E ~,+1: x ~ M o , Set k = l. • Iteration k = 1, 2 .... k.1. Solve subproblem (Qk) min{t -- k(x): (x,
(qo, x -
x ° ) + l(5o, x °) <<. t}.
t)ePk}.
Let (x k, t k) be a global optimal solution of (Qk). k.2. If t k >1 k(xk), stop: If 5k- ~ < + ~ then £k- ~ solves (Q), if 5k- 1 = + ~ then (Q) is infeasible. k.3. If t k < k(x k) set 'Sk_
1
5k = ~ m i n ( s k - l, h(x k) - g(xk)}
if p(x*) -- q(x*) > O, otherwise,
(18)
and define ~k = arg (Sk). k.4. Compute qk+l 66~l(5k, xR) - Let
Pk+l = Pk C~ ((X, t)elt~"+l: (qk+l, X -- X k) + 1(5k, X k) <<.t). k.5. S e t k = k +
1 and go tok.1.
(19)
T.Q. Phong / Operations Research Letters 15 (1994) 73 79
77
Remark 1. (i) In step k.3, it follows from the construction o f t x k by (18) that x k ¢ S,~, hence Xk ~ lnt S,~ (see (6)). Therefore F(Ctk, x k) >1 0 and, since t k < k(xk), we must have t k < l(Ctk, xk), i.e. (t k, x k) ¢ Ok. A cutting plane is then added to PR t o obtain Pk+l (see (19)). (ii) If t k >1 k(x k) then min F(~ k_ 1, x)
=
min{/(ak_ 1, x) -- k(x): x E R"}
XER n
= min{t - k(x): l ( e k - l , x ) < ~ t} >/min{t -- k(x): X e P k } = t k -- k ( x k) >10.
If~k- 1 = h( yck - 1) _ 9( yk - 1) < + Oe then, by the construction, p( yck - 1) _ q( yk - ~) <<.0 SO F (~k - 1, ~k - 1) <~ O. Hence min F ( ~ k - 1, x) = F ( ~ k - 1, ~k- 1) = 0 and, by Theorem 1, yk- 1 solves (Q). If~k- 1 = + ~ then, as seen above, rain F(otk_ 1, x) >~ O. Suppose that 3~: min F ( ~ k - 1, x) = F ( ~ k - 1,2) = 0 then, by the definition of F ( ~ k - l , x ) (taking into account the fact ~k-~ = + oo), one must have p(~) - q(2) = 0 and p(x) - q(x) >~ O, Vx E R", that means the robustness assumption is not verified. Hence min F(~k- 1, x) > 0 and, by Theorem 1, S = 0. (iii) The main source of computational effort of the algorithm is the work on subproblems (Qk). By virtue of the concavity of function t - k(x) the problem is done if the vertex set of Pk is available. By the construction, Pk + 1 is derived from Pk by adding a cutting plan. The vertex set of Pk + ~ then can be found via knowledge of vertex set of Pk by, e.g., method of Horst-Thoai-Vries [1] or by the recent modification proposed by the author in [6] which exploits the fact that Pk is unbounded polyhedral set with only one infinite direction along axis t.
4. Convergence of the algorithm The following properties follow from the construction of cutting planes. Lemma 1. The hyperplane constructed by (19) separates (x k, t k) and Dk. Lemma 2. I f (x k, tk) ~ (Y, t-) and ~k ~ ~t then t = l(~, if). We are in position to establish the convergence theorem.
Theorem 2. I f S ~ 0 and the algorithm does not terminate after a finite number of iterations then ~k ~ ~* = min(Q) and any accumulation point o f the sequence {2k} solves (Q). Proof. It is easily seen, from the construction, that {(xk, tk)} is bounded. Therefore, one may select a subsequence, if necessary, {(xkin, tk")} that converges to (if, r). Let ~ = limk~Ctk. By Lemma 2 we have { = l(~, Y). Since t k < k(Xk) then ?~< k(2) so we have F(& £) = l(8, £) -- k(£) = ~ - k(£) ~< 0.
(20)
T.Q. Phong / Operations Research Letters 15 (1994) 73-79
78
As it was shown in Remark (i) F(~k, x k) >>-0 and, since {~k} is nonincreasing, we have F(o~g, x k) >~ 0 for all j/> k. By passing to the limit, one has F(~, x k) >~ O, for k = 1, 2,..., that follows F(0~,£) >~ 0. This, together with (20) implies F(0~,£) = 0. Indeed, min F(0~, x) = 0, since for all k = 0, 1,..., and for all x: F(~k- 1, X) >>.t k -- k(x k) so F(~, x) >t F(~k- 1, x) >~ t k -- k(x k)
and again, by passing to the limit. F(cT, x) ~> [ - k(~) = l(~, ~) - k(2) = F(~, ~) = 0. Using the same argument as in Remark (2), we see that 0~< + oo, and, by Theorem 1, or* = ~ = min(Q). Since ~k = h(x k) - g(£k) __.}~, any cluster point of the sequence {£k) is a global optimal solution of (Q).
It can be seen from the proof that any accumulation point of the sequence {x k} (may be infeasible) is also a global optimal solution for (Q). The following also follows from Theorem 2. Corollary 1. If S # 0 then the algorithm generates a feasible point after finitely many steps. Remark 2. Our method may be applied to the general class of indefinite quadratic programming problems, i.e. those of minimizing a quadratic function under quadratic constraints. Indeed, as is well known, any indefinite quadratic function is a d.c. function. However, the expression of a d.c. function as the difference of two convex functions is not unique in general. Practically, the choice of such an expression is very important. In [61, a specific decomposition of an indefinite quadratic function into the difference of two convex functions was proposed and successfully used for dealing with a feasible region given by several indefinite quadratic constraints. The reader should be warned that such a feasible region needs not necessarily satisfy the assumption of robustness defined in Section 2. Remark 3. It should be noticed that our most important contribution is dealing with the d.c. structure of a feasible region which is still little investigated in the literature. Indefinite quadratic programming problems considered in [5] (see also [4]) have as feasible region a polytope given by a system of linear inequations. Therefore, the methods developed in [5] should be more efficient than the method herein which is designed for a more general framework.
5. Example We illustrate the Algorithm by the following example: min
( - x 2 + x 2 + x 2)
s.t.
xl + x2 + x3 <%2, xl-4x
2-x
2<~ - 1,
x l , X2, x3 >10.
The objective function can be rewritten as h - g where h ( x , , x2, x3) = x 2 + x32,g(xl, x2, x3)
=
X2
T.Q. Phong / Operations Research Letters 15 (1994) 73-79
79
a n d the c o n s t r a i n t can be expressed, using Eq. (1), as p(x) - q(x) ~ 0 where p(xl,xz,x3)
= m a x { m a x { - x l , - x2, - x 3 , x ,
+ x2 + x 3 -
2} + 4x z + x~,x~ + 1},
q(x1, x2, x3) = 4x~ + x3z. Initialization:
Mo = {(0, 0, 0), (2, 0, 0), (0, 2, 0), (0, 0, 2)}, ~Zo = (0, 0, 2), Vo = {(0, 0, 0, - 5), ( 2 , 0 , 0 , - 5), (0, 2, 0, - 5), (0,0, 2,4)}. Iteration 1: T h e cutting h y p e r p l a n e
0.1x I + 16.5x2 -4- 0 . 0 X
3 --
t --
16.999 = 0
a n d the vertex set of p o l y h e d r a l set P~ V1 = {(0, 0, 0, - 5), (2,0,0, - 5), (0, 2, 0, 16), (0, 0, 2, 4), (0, 0.727, 0, - 5), (1.273, 0.727, 0, - 5),
(0, 1, 1, - 0.5)} were determined. T h e a l g o r i t h m s t o p p e d after 22 iterations a n d p r o v i d e d the exact global o p t i m a l solution = (1.250, 0.750, 0.000) with the o p t i m a l value ~ = - 1.
Acknowledgement
T h e a u t h o r w o u l d like to t h a n k the referee for his helpful c o m m e n t s .
References [1] R. Horst, N.V. Thoai and J. de Vries 0988), "On finding new vertices and redundant constraints in cutting plane algorithms for global optimization", Oper. Res. Lett. 7, 85-90. [2] R. Horst and H. Tuy (1990), Global Optimization: Deterministic Approaches, Springer, Berlin. [3] R. Horst, T.Q. Phong, N.V. Thoai and J. de Vries (1991), "On solving a d.c. programming problem by a sequence of linear programs", J. Global Optim. 1, 183-203. [4] P.M. Pardalos (ed.) (1993), Complexity in Numerical Optimization, World Scientific. [5] P.M. Pardalos and J.B. Rosen (1987), Constrained Global Optimization: Algorithm and Applications, Springer Verlag, Berlin, Lecture Notes in Computer Science 268. [6] T.Q. Phong, D.T. Pham and T.H.A. Le, "A New Method for Solving d.c. Programming Problem, Application to Fuel Mixture Nonconvex Optimization Problem", J. Global Optim., (submitted). [7] P.T. Thach and H. Tuy (1990), "'The relief indicator method for constrained global optimization", Naval Res. Logist. 37, 473-497. [8] H. Tuy (1990), "The relief indicator method as a new approach to constrained global optimization", Prep. Inst. of Math. Hanoi.