An Improved Interval Global Optimization Method and Its Application to Price Management Problem Bartlomiej Jacek Kubica1 and Ewa Niewiadomska-Szynkiewicz2 1 Institute of Control and Computation Engineering, Warsaw University of Technology, Nowowiejska 15/19, PL-00-665 Warsaw, Poland
[email protected] http://www.ia.pw.edu.pl 2 Research and Academic Computer Network (NASK), Wawozowa 18, PL-02-796 Warsaw, Poland
[email protected] http://www.nask.pl
Abstract. We present an interval global optimization algorithm using a modified monotonicity test. The improvement applies to constrained problems and can result in significant speedup, when constraints are sparse, i.e. they “bind” a few of the variables, not all of them. A theorem that ensures the correctness of the new tool, is given and proved. The improved method is applied to an economic problem of setting optimal prices on a couple of products.
1
Introduction
Many issues related to the practical problems require the solution of the following very general optimization problem: min f (x) ,
x∈X
(1)
where the function f (x) can be nonconvex and X is defined by the set of (possibly nonconvex) constraints: gj (x) ≤ 0, j = 1, . . . , m and bounds on each component of x. The problem of designing algorithms to compute global solutions is very difficult. In general there are no local criteria in deciding whether the local solution is the global one. Many of methods apply heuristics able to find an approximate solution only, see [2], [8]. Interval methods aim to find the global solution(s). Unfortunately, they are usually slow and memory demanding so, the acceleration is worthwhile. The general framework of the optimization algorithm is the branch–and–bound schema. Many different tools are used to reject or at least narrow subboxes that do not contain the global optimum, e.g. [3], [4], [9], [2]. The objective of this paper is to identify some of the advantages and disadvantages of existing algorithms and provide suitable modifications to increase their efficiency. We propose a modified monotonicity test for constrained optimization problems. It allows to remove from further considerations far more boxes than its traditional relative. The presented algorithm was applied to a practical marketing problem. B. Kragstr¨ om et al. (Eds.): PARA 2006, LNCS 4699, pp. 1055–1064, 2007. c Springer-Verlag Berlin Heidelberg 2007
1056
2
B.J. Kubica and E. Niewiadomska-Szynkiewicz
Interval Branch–and–Bound Method
The general scheme of interval global optimization methods is branch–and– bound (see e.g. [1], [3], [4]). It can be expressed by the following pseudocode: IBB (x(0) ; f, ∇f, ∇2 f; g1 , ∇g1 , ∇2 g1 , . . . , gm , ∇gm , ∇2 gm , . . .) // x(0) is the initial box // f(·) is the interval extension of the objective function f (·) // ∇f(·) and ∇2 f(·) are interval extensions of gradient and Hessian of f (·) // gi (·) are interval extensions of the constraints, etc. // Lsol is the list of solutions [y (0) , y (0) ] = f(x(0) ) ; compute fmin = the upper bound on the global minimum (e.g. objective value in a feasible point) L = {(x(0) , y (0) )} ; // the list of boxes Lsol = ∅ ; while (L = ∅) do x = the element of L with the lowest function value underestimation ; compute the values of interval extensions of the constraint functions ; if (x is infeasible) then discard x ; update fmin if possible ; perform other rejection/reduction tests on x ; if (x is verified to contain a unique critical point or x is small and not infeasible) then add x to Lsol ; else bisect x to subboxes x(1) and x(2) ; compute lower bounds y (1) and y (2) on the function value in the obtained boxes ; delete x from L ; for i = 1, 2 do put (x(i) , y (i) ) on the list L preserving the increasing order of the lower bounds ; end for delete from L boxes with y(i) > fmin ; end if end while delete from Lsol the boxes with y(i) > fmin ; return Lsol ; end IBB
The above pseudocode mentions some “rejection/reduction tests” that may be used in conjunction with the IBB algorithm. There are several such tests. Most important of them are: several kinds of Newton operators, constraint propagation steps and, probably oldest of them all, monotonicity tests. We do not describe them all, as they are widely available in literature, e.g. [1], [4], etc. Most of the current research concentrate on several types of interval Newton operators and tools related to them e.g. [6]. It is reasonable – interval Newton methods are one of the most powerful interval techniques, indeed. They can not
An Improved Interval Global Optimization Method and Its Application
1057
only reject or narrow subboxes during the branch–and–bound process, but they can also verify that a subbox contains a root (or a feasible point) certainly; sometimes, they can even verify the uniqueness of this root. Nevertheless, in this paper, we shall concentrate on monotonicity tests. This tool is less investigated recently, but – as we shall see – it can also be improved, resulting in a significant speedup of the overall algorithm.
3
The Classical Monotonicity Test
Monotonicity test ([2]) allows to discover boxes where the objective function is strictly monotone – such boxes cannot contain a global minimum. In the case of a general constrained optimization problem, monotonicity test can be applied only to boxes where all constraints are satisfied, otherwise it does not guarantee the non-existence of an optimum. By x ⊆ X we denote the box on X and the inclusion functions of f , ∇f and g, respectively – by f(x), ∇f(x) and g(x). The classical monotonicity test (CMT) can be formulated by the following pseudocode: for i = 1, . . . , n do compute [y, y] = ∇fi (x) ; if (y < 0) then // the function is decreasing compute [gj , g j ] = gj (x) for all j = 1, . . . , m ; if (g j < 0 ∀j = 1, . . . , m) then discard x ; else if (g j ≤ 0 ∀j = 1, . . . , m) then update xi to xi = xi ; end if end if if (y > 0) then // the function is increasing ... end if end for The simplest way to remove infeasible solutions is the flag system, proposed by Ratschek and Rokne [2]. The idea is as follows. We associate with each box x the binary vector r = (r1 , . . . , rm ), where m is the number of constraints. So, rj = 1 indicates that the j-th constraint is satisfied in the box x, otherwise rj = 0. It is obvious that if a considered constraint is satisfied in a box x, it would also be satisfied in all subboxes, resulting from its bisections. The algorithm, determining the flags for a new box xk , obtained from the bisection of x, may be expressed as follows: for j = 1, . . . , m do if (rj == 1) then rjk = 1 ; else compute [g, g] = gj (xk ) ;
1058
B.J. Kubica and E. Niewiadomska-Szynkiewicz
if (g ≤ 0) then rjk = 1 ; else if (g > 0) then discard xk ; else rjk = 0 ; end if end for In the case of problems in which the constraints bind only a group of variables the efficiency of the monotonicity test can be increased. Next section proposes a suitable modification.
4
Main Results
The proposed modified version of the monotonicity test allows to discard a box x in the case when the constraint g(x) ≤ 0, x ∈ x is not satisfied but does not “bind” the variable, with respect to which the objective is monotonous. We say that gj (x1 , . . . , xn ) does not “bind” the variable xi when the variable xi is not present in gj (x) formulation. The modified algorithm (MMT) is as follows: for i = 1, . . . , n do compute [y, y] = ∇fi (x) ; if (y < 0) then // the function is decreasing ∂gj (x) = [0, 0] and j = 1, . . . , m ; compute [g j , g j ] = gj (x) for all j such that ∂x i if (all computed g j ’s satisfy g j < 0) then discard x ; else if (all computed g j ’s satisfy g j ≤ 0) then update xi to xi = xi ; end if end if if (y > 0) then // the function is increasing ... end if end for
It can be proved that the modified test preserves all the solutions (see Theorem 1, below). A minor improvement, associated with the modified monotonicity test is a new variant of the flag vector. We can distinguish four possibilities: a) the j-th constraint is satisfied and inactive in the box x (g j < 0), b) the constraint is satisfied, but may be active (g j ≤ 0), c) the constraint may be violated for some points of the box (0 ∈ [gj , g j ]), d) the constraint is violated in the whole box (g j > 0) and the box should be discarded from further processing. Well then, to be consistent with the monotonicity test it is better to use “flags” of three possible values, instead of binary ones: for j = 1, . . . , m do if (rj == 2) then rjk = 2 ; else compute [g, g] = gj (xk ) ;
An Improved Interval Global Optimization Method and Its Application
1059
if (g < 0) then rjk = 2 ; else if (g ≤ 0) then rjk = 1 ; else if (g > 0) then discard xk ; else rjk = 0 ; end if end for And now let us present the main theorem: Theorem 1. Consider the optimization problem (1). Consider a box x = (x1 , . . . , xn )T , contained in the interior of the feasible set X. Assume for some i ∈ {1, . . . , n} the following two conditions are both fulfilled: – ∀x ∈ x either
∂f (x) ∂xi
< 0 or
∂f (x) ∂xi
> 0,
– ∀j = 1, . . . , m either we have ∀x ∈ x
gj (x) < 0 or
∂gj (x) ∂xi
≡ 0.
Then there is no optimum in the box x. Proof Assume a point x0 ∈ x is the (at least local) optimum of the problem (1). The < 0 or ∂f∂x(x) > 0, i.e. suppositions of Theorem 1 say that ∀x ∈ x either ∂f∂x(x) i i the objective f (·) is monotone on the box x. If f (·) is monotone at x0 that the only case when x0 can be the optimum of the problem (1) is when some constraints are active at this point and all improvement directions are infeasible. Nevertheless, we can simply show that at least one improvement direction is feasible at x0 . Let us consider the direction pointed by the vector d = (d1 , . . . , dn )T , where: 0, for k = i dk = if f (·) is decreasing, or 1, for k = i 0, for k = i if f (·) is increasing. dk = −1, for k = i From now on let us assume that f (·) is decreasing, i.e. ∂f∂x(x) < 0. The proof for i the second case would be analogous. Obviously, d is an improvement direction, because: ∂f (x) ∇f (x0 )T · d = 0
∀α
0 < α < σ → (x0 + αd) is feasible.
(3)
Indeed, consider a constraint gj (x) < 0, j ∈ {1, . . . , m}. From the Taylor expansion we get: gj (x0 + αd) = gj (x0 ) + ∇gj (ξ) · (αd) = gj (x0 ) = α
∂gj (ξ) , ∂xi
(4)
1060
B.J. Kubica and E. Niewiadomska-Szynkiewicz
where ξ ∈ (x0 , x0 + α · d) ⊆ x. But this implies that to: gj (x0 + α · d) = gj (x0 ) .
∂gj (ξ) ∂xi
= 0 and (4) reduces (5)
But x0 is feasible by supposition, so gj (x0 + α · d) = gj (x0 ) < 0. Consequently, the direction d is feasible. Thus, there is a feasible improvement direction in any x0 ∈ x, so no optimum can be in x. Interpretation. The essence of the assumptions in Theorem 1 is quite simple, actually. The first supposition says that the objective is monotone with respect to the variable xi . The second one says that all of the constraints are either inactive or do not contain the variable xi in their formulae. This simply means that all of the active constraints are invariant with respect to xi and the feasible set has the form of some kind of a pipe along xi . An example of such a set is shown on Figure 1. It should be obvious that in such case the constraints cannot block the direction along xi , which is an improvement direction.
Fig. 1. When can the modified monotonicity test be applied ?
5
How to Treat Bound Constraints ? Peeling Process
To investigate the usefulness of the developed rejection/reduction test it is reasonable to check how it cooperates with other tools, like Newton methods, etc. We cannot present all of the results here. Investigated variants of the IBB algorithm use the interval Newton step (precisely: interval Gauss-Seidel step). Whereas we examine how does the modified monotonicity test cooperate with the “peeling” process (see e.g. [3], [4]), so we apply four variants of the algorithm: with classical or modified monotonicity test and with or without peeling. What is the peeling process ? It is one of the ways to deal with bound constraints. Actually there are two main approaches to solve bound–constrained problems:
An Improved Interval Global Optimization Method and Its Application
1061
– treat bound constraints as any other inequality constraints (or, roughly speaking, similarly to them), – consider the “boundary boxes” and “interior boxes” separately. In this second variant, we use reduced gradients and Hessians in boundary boxes (because they have a reduced dimension) and in interior boxes bounds constraints are ignored (in monotonicity tests or Newton steps; only other constraints are used then). The boundary boxes are created at the beginning of the IBB algorithm, by the peeling process. Let us formulate the peeling procedure in a pseudocode. We present it in a slightly simpler form than the one in [3] or [4]: peeling (x, i) if (i == n) then let y = x ; let yi = [xi , xi ] ; enqueue (y) ; let yi = [xi , xi ] ; enqueue (y) ; let yi = xi ; enqueue (y) ; else let y = x ; let yi = [xi , xi ] ; peeling (y, i + 1) ; let yi = [xi , xi ] ; peeling (y, i + 1) ; let yi = xi ; peeling (y, i + 1) ; end if end peeling According to [3], complexity of the peeling process is of the degree 3k , where k is the number of bound constraints. This leads to inefficiency of this process if the number of bound constraints is large. All the same, not all of the bounds may actually be “bound constraints” – some of them may be known to be redundant or non–significant. Unfortunately, that is not the case for the problem we are trying to solve.
6
Numerical Experiments
Numerical experiments were performed for both testing problems and real-life applications. Results obtained for the selected test problem defined as: min x
n
(−1)i · x2i ,
(6)
i=1
s.t. ni=1 xi ≤ 1, xi−1 + xi ≤ 0.5, −1 ≤ xi ≤ 2, i = 1, . . . , n. are given in Fig. 2. One of the considered case study was concerned with computing of the optimal prices for products that are sold in a shop. The goal was to maximize the total n xi profit defined as: P R = i=1 1+vi − di · Si , where n denotes the number of
1062
B.J. Kubica and E. Niewiadomska-Szynkiewicz 30000000
function evaluations
25000000 20000000 CMT
15000000
MMT
10000000 5000000 0 2
3
4
5
6
7
8
problem dimension
Fig. 2. Number of function evaluations; CMT – classical and MMT – modified monotonicity tests
products, vi and di given constants, corresponding to the market entities of VAT and cost per product i, Si are expected sales of products within the considered period and xi prices, we are trying to set. Two models describing the market response Si on the price were considered: Cobb-Douglas model described in [10]. This model implements the cross-effects ij n with other substitute or complementary own products: Si = αi · j=1 xβj , where β ij is the elasticity of sales of the i–th product with respect to the price of the j–th product, β ii is referred to as the direct elasticity and β ij , where i = j – the cross elasticity. This function is widely used, but it does not capture some important effects, such as different market sensitivities to small and large price changes. These features are expressed by the so-called s-shape models. Hybrid model formulated in [7] that exhibits an s-shape and includes cross ij n co effects: Si = ai + αi · j=1 xβj + c1i · sinh c2i · (xi − xco i ) , where xi is the average competitive price and a, c1 and c2 are model parameters. The model combines descriptions proposed by Cobb-Douglas and Gutenberg (see [10]). The following constraints for price, sale and cash of each product were considered: bounds on prices and bounds on the cash flows; sales, n constraints for total sale n and cash: ST min ≤ i=1 Si ≤ ST max, CT min ≤ i=1 xi · Si ≤ CT max , linear constraints for price differences of substitute or complementary products dmin ij ≤ xi − xj ≤ dmax (i, j) ∈ D, where D is the set of pairs of correlated products. ij , The overall optimization problem was as follows: n xi − di · Si (x) , max P R = x 1 + vi i=1
(7)
max ≤ xi ≤ xmax , Simin ≤ Si (x) ≤ Simax , Cimin ≤ x , for i = s.t. xmin i · Si (x) ≤ Ci i i n n min max min 1, . . . , n and ST ≤ i=1 Si (x) ≤ ST , CT ≤ i=1 xi · Si (x) ≤ CT max , max dmin (i, j) ∈ D. ij ≤ xi − xj ≤ dij ,
An Improved Interval Global Optimization Method and Its Application
1063
Tables 1 and 2 present the results of application of the flag-based IBB algorithm with classical (CMT) and modified (MMT) monotonicity tests to price management problem formulated for 8 products. The first one applies to the case of Cobb-Douglas model and the second one – to Hybrid model. The accuracy for all eight cases was set to: ε = 0.05 and there were two pairs of correlated products in the set D. Meaning of the columns is as follows: the first one describes type of the monotonicity test, the other ones – the number of bisections, objective evaluations, objective’s gradient and Hessian evaluations, constraints evaluations, respectively, and the last one – the number of boxes, enclosing the global optimum, resulting form the IBB algorithm. Table 1. Price management problem – results for Cobb-Douglas model IBB without peeling mon. test exec.time bisec. f.evals. grad.evals. Hess.evals. constr.evals. res.boxes CMT 25917.8 s 262430 776208 516851 514803 777351 252373 MMT 28.92 s 940 2285 1848 1377 2435 437 IBB with peeling mon. test exec.time bisec. f.evals. grad.evals. Hess.evals. constr.evals. res.boxes CMT 41561.7 s 317514 856786 574011 564052 1684852 223292 MMT 273.48 s 2603 9587 8356 3217 271166 263
Table 2. Price management problem – results for Hybrid model IBB without peeling mon. test exec.time bisec. f.evals. grad.evals. Hess.evals. constr.evals. res.boxes CMT 1772.78 s 33810 100203 66848 66578 100554 32768 MMT 11.75 s 435 850 836 447 1040 12 IBB with peeling mon. test exec.time bisec. f.evals. grad.evals. Hess.evals. constr.evals. res.boxes CMT 8289.94 s 109380 249401 184489 175911 699398 30205 MMT 395.64 s 4023 11890 8879 4757 286847 7
Results. It can be observed that the modified monotonicity test results – at least for investigated problems – in dramatical speedup. The improvement is significant in both cases – when we use peeling process and when we do not. Actually, for the price management problem, the variant of the IBB algorithm using peeling was much slower. This was caused by the fact that we had significant bound constraints on all of the eight decision variables. Peeling is usually inefficient in such cases. Nevertheless, the monotonicity test gave a large speedup also for this variant of the algorithm. Obviously, the improvement was that large, because the examined problem had several constraints biding only few variables. Otherwise the speedup would not be achieved.
1064
7
B.J. Kubica and E. Niewiadomska-Szynkiewicz
Conclusions
As a final conclusion we can say that the proposed modifications to the flagbased IBB algorithm and monotonicity test increase the speed of convergence with respect to classical ones. It seems to be a powerful accelerating device, which gives approximately exponential improvement for constraints binding few variables.
References 1. Hansen, E.: Global Optimization Using Interval Analysis. Marcel Dekker, New York (1992) 2. Horst, R., Pardalos, P.M. (eds.): Handbook of Global Optimization. Kluwer, Dordrecht (1995) 3. Kearfott, R.B.: A Review of Techniques in the Verified Solution of Constrained Global Optimization Problems. In: Kearfott, R.B., Kreinovich, V. (eds.) Applications of Interval Computations, Kluwer, Dordrecht (1996) 4. Kearfott, R.B.: Rigorous Global Search: Continuous Problems. Kluwer, Dordrecht (1996) 5. Kearfott, R.B., Nakao, M.T., Neumaier, A., Rump, S.M., Shary, S.P., van Hentenryck, P.: Standardized notation in interval analysis, available at http://www.mat.univie.ac.at/∼ neum/software/int/notation.ps.gz 6. Kubica, B.J., Malinowski, K.: An Interval Global Optimization Algorithm Combining Symbolic Rewriting and Componentwise Newton Method Applied to Control a Class of Queueing Systems. Reliable Computing 11(5), 393–411 (2005) 7. Malinowski, K.: PriceStrat 4.0 Initial Research Paper. KSS Internal Document, Manchester (2000) 8. Michalewicz, Z., Fogel, D.B.: How to Solve It: Modern Heuristics. Springer, Heidelberg (2004) 9. Neumaier, A.: Complete Search in Continuous Global Optimization and Constraint Satisfaction. In: Acta Numerica, pp. 271–369. Cambridge University Press, Cambridge (2004) 10. Simon, H.: Price Management. North-Holland (1989)