Two-Phase Generalized Reduced Gradient Method for Constrained Global Optimization

Abstract

The random perturbation of generalized reduced gradient method for optimization under nonlinear differentiable constraints is proposed. Generally speaking, a particular iteration of this method proceeds in two phases. In the Restoration Phase, feasibility is restored by means of the resolution of an auxiliary nonlinear problem, a generally nonlinear system of equations. In the Optimization Phase, optimality is improved by means of the consideration of the objective function, on the tangent subspace to the constraints. In this paper, optimal assumptions are stated on the Restoration Phase and the Optimization Phase that establish the global convergence of the algorithm. Some numerical examples are also given by mixture problem and octagon problem.

1. Introduction

We consider the problem (𝐗)(𝐗)Maximize∶𝐟subjectto∶𝐡=0,𝐗∈Ω, (1.1) where 𝐟∶ℝ^𝑛→ℝ and 𝐡∶ℝ^𝑛→ℝ^𝑚 are continuously differentiable, and Ω⊂ℝ^𝑛 is closed and convex set (e.g., Ω=∏^𝑛_𝑖=1[𝐚_𝑖,𝐛_𝑖]).

Notice that any nonlinear programming can be put into standard form (1.1), by introduction of nonnegative slack variables, if there are inequalities (other than bounds on the variables) among the constraints, and by allowing some of the bounds to be +∞ or −∞ if necessary. The standard form is adopted here for ease in notations and discussion.

Feasible methods for optimization with constraints like generalized reduced gradient play an important role in practicing and are still widely used in technological applications. The main techniques proposed for solving constrained optimization problems in this work are generalized reduced gradient method via random perturbation. We are mainly interested in the situation where, on one hand, 𝐟 is not concave and, on the other hand, the constraints are in general not linear. It is worth noting that some variant of the generalized gradient method reduces, in the case where all the constraints are linear, to the reduced gradient method [1], and some other variant, in the case of linear programming, to the Dantzig simplex method.

The problem (1.1) can be numerically approached by using sequential quadratic programming in [2], other methods for nonlinear programming in [3, 4], and generalized reduced gradient method [5], which generates a sequence {𝐗^𝑘}_𝑘≥0, where 𝐗⁰ is an initial feasible point and, for each 𝑘>0, a new feasible point 𝐗^𝑘+1 is generated from 𝐗^𝑘 by using an operator 𝐐_𝑘 (see Section 3). Thus the iterations are given by𝑘≥0∶𝐗^𝑘+1=𝐐_𝑘𝐗^𝑘. (1.2)

A fundamental difficulty arises due to the lack of concavity: the convergence of the sequence {𝐗^𝑘}_𝑘≥0 to a global maximum point is not ensured for the general situation considered. In order to prevent from converging to local maximum, various modifications of these basic methods have been introduced in the literature. For instance, we can find in the literature modifications of the basic descent methods [6–9], stochastic methods combined to penalty functions [10], evolutionary methods [11], and simulated annealing [12]. We introduce in this paper a different approach, inspired from the method of random perturbations introduced in [13] for unconstrained minimization of continuously differentiable functions and adapted to linearly constrained problems in [14].

In such a method, the sequence {𝐗^𝑘}_𝑘≥0 is replaced by a random vectors sequence {𝐗^𝑘}_𝑘≥0 and the iterations are modified as follows:𝑘≥0∶𝐗^𝑘+1=𝐐_𝑘𝐗^𝑘+𝒫_𝑘, (1.3) where 𝒫_𝑘 is a suitable random variable, called the stochastic perturbation. The sequence {𝒫_𝑘}_𝑘≥0 goes to zero slowly enough in order to prevent convergence to a local maximum (see Section 4), the generalized reduced gradient method is recalled in Section 3, the notations are introduced in Section 2, and the results of some numerical experiments are given in Section 5.

2. Notations and Assumptions

We have the following. (i)𝐄=ℝ^𝑛 is the 𝑛-dimensional real Euclidean space.(ii)𝐗 stands for (𝐱₁,…,𝐱_𝑛)^𝑇∈𝐄.(iii)𝐡(𝐗) is the column vector whose components are 𝐡₁(𝐗),…,𝐡_𝑚(𝐗).(iv)‖𝐗‖=(𝐱²₁+⋯+𝐱²_𝑛)^1/2 is the Euclidean norm of 𝐗: ‖𝐀‖_∞=sup{‖𝐀𝐗‖∶‖𝐗‖=1}, (2.1)(v)𝐀^𝑇 is the transpose matrix associated to 𝐀.

Definition 2.1. A point 𝐗 is said to be feasible if |𝐡_𝑖(𝐗)|≤𝜀₀, for 𝑖=1,…,𝑚, where 𝜀₀ is same preassigned small positive constant.

Assume that we know some feasible point 𝐗⁰. We make the following nonedegeneracy assumption. The vector 𝐗 can be split into two components: an 𝑚-dimensional component 𝐲 (the basic part), and 𝐱, a component of dimension 𝑛−𝑚 (the nonbasic part) such that the following two properties hold: (H1)𝐲⁰ is strictly between bounds; (H2) the square 𝑚×𝑚 matrix 𝜕𝐡/𝜕𝐲, computed at 𝐗⁰, is nonsingular.

If the property (H2) does not hold or 𝑛<𝑚, we add artificial variables to the constraints, so the property (H2) holds and 𝑛≥𝑚 in associate problem (1.1).

Let 𝐂={𝐗∈𝐄∣𝐡(𝐗)=0,𝐗∈Ω}. (2.2) The objective function is 𝐟∶𝐄→ℝ, and its upper bound on 𝐂 is denoted by 𝐥^∗: 𝐥^∗=max_𝐂𝐟. (2.3) Let us introduce 𝐂_𝛼=𝐒_𝛼∩𝐂;where𝐒_𝛼={𝐗∈𝐄∣𝐟(𝐗)≥𝛼}. (2.4) We assume that𝐟istwicecontinuouslydifferentiableon𝐄,∀𝛼<𝐥^∗∶𝐂_𝛼isnotempty,closedandbounded,∀𝛼<𝐥^⋆∶meas𝐂_𝛼>0, (2.5)(2.6)(2.7) where meas(𝐂_𝛼) is the measure of 𝐂_𝛼.

Since 𝐄 is a finite dimensional space, assumption (2.6) is verified when 𝐂 is bounded or 𝐟 is coercive, that is, lim_{‖𝐗‖→+∞}𝐟(𝐗)=−∞. Assumption (2.6) is verified when 𝐂 contains a sequence of neighborhoods of a point of optimum 𝐗^⋆ having strictly positive measure, that is, when 𝐗^⋆ can be approximated by a sequence of points of the interior of 𝐂. We observe that the assumptions (2.5)-(2.6) yield that 𝐂=_𝛼<𝐥⋆𝐂_𝛼,thatis,forall𝐗∈𝐂∶∃𝛼<𝐥^∗suchthat𝐗∈𝐂_𝛼. (2.8) From (2.5)-(2.6),𝛾₁=sup‖∇𝐟(𝐗)‖∶𝐗∈𝐂_𝛼<+∞. (2.9) Then 𝛾₂=sup‖𝐝‖∶𝐗∈𝐂_𝛼<+∞, (2.10) where 𝐝 is the direction determined on the basis of the gradient ∇𝐟(𝐗).

Thus,𝛽(𝛼,𝜀)=sup𝐘−𝐗+𝜂𝐝∶(𝐗,𝐘)∈𝐂_𝛼×𝐂_𝛼,0≤𝜂≤𝜀<+∞, (2.11) where 𝜀,𝜂 are positive real numbers.

3. Generalized Reduced Gradient Method

By the implicit function theorem, there exists, in some neighborhood 𝐕 of 𝐱⁰, a unique continuous function (mapping), say 𝐲(𝐱), such that 𝐟(𝐱,𝐲(𝐱) is identically zero in 𝐕. In addition, 𝐲(𝐱) has a continuous derivative 𝐝𝐲/𝐝𝐱, which can be computed by the chain rule: 𝜕𝐡+𝜕𝐱𝜕𝐡𝜕𝐲𝐝𝐲𝐝𝐱=0, (3.1) or, more conveniently,𝐝𝐲𝐝𝐱=−𝜕𝐡𝜕𝐲⁻¹𝜕𝐡.𝜕𝐱 (3.2) In what follows, we call 𝐀 the Jacobian of 𝐡(𝐗) computed at 𝐗⁰. Similarly, we set 𝐁=𝜕𝐡𝜕𝐲,𝐍=𝜕𝐡𝜕𝐱computedat𝐗⁰. (3.3) Substituting 𝐲(𝐱) into the objective function 𝐟(𝐱,𝐲), we obtain the reduced function:𝜙(𝐱)=𝐟(𝐱,𝐲(𝐱)), (3.4) the gradient of which at 𝐱⁰ is, by the chain rule again,𝐠=𝜕𝐟+𝜕𝐱𝜕𝐟𝜕𝐲𝐝𝐲𝐝𝐱 (3.5) (all derivatives computed at 𝐗⁰). Setting 𝐯=𝜕𝐟𝜕𝐱,𝐰=𝜕𝐟𝜕𝐲computedat𝐗⁰, (3.6) we have the following formula for the reduced gradient 𝑔: 𝐠=𝐯−𝐰𝐁⁻¹𝐍. (3.7) The generalized reduced gradient method tries to extend the methods of linear optimization to the nonlinear case. These methods are close to, or equivalent to, projected gradient methods [15]; just the presentation of the methods is frequently quite different.

Let us define the projected reduced gradient 𝐩 (in the space of 𝐠) by its components: (i)𝐩_𝑗=0 if 𝐱⁰_𝑗=𝐚_𝑗 and 𝐠_𝑗<0; (ii)𝐩_𝑗=0 if 𝐱⁰_𝑗=𝐛_𝑗 and 𝐠_𝑗>0; (iii)𝐩_𝑗=𝐠_𝑗, otherwise.

It is convenient to set 𝐮=𝐰𝐁⁻¹, (3.8) and hence 𝐰−𝐮𝐁=0,𝐯−𝐮𝐍=𝐠. (3.9) It is a simple matter to verify that Kuhn-Tucker condition for the problem in standard form reduces to 𝐩=0, and that 𝐮 is the row-vector of multipliers corresponding to the equations 𝐡(𝐗)=0. We assume from now on that 𝐩≠0.

The following relations 𝐱_𝑗=𝐚_𝑗 if 𝐱⁰_𝑗=𝐚_𝑗 and 𝐠_𝑗<0, 𝐱_𝑗=𝐛_𝑗 if 𝐱⁰_𝑗=𝐛_𝑗 and 𝐠_𝑗>0

define what we call the face (at 𝐗⁰), denoted by 𝐅. The row-vector 𝐩 (the projected reduced gradient) is also named the projection of 𝐠 onto 𝐅.

Let 𝐝 be any nonzero vector column-vector in 𝐅 such that 𝐩𝐝>0; the vector 𝐝 is an ascent direction for the reduced function 𝜙(𝐱).

There is a striking analogy with what is usually done in linear programming, where 𝐲(𝐱) can be computed in close form. This is generally not the case if the constraints 𝐡(𝐗) are nonlinear. Even if close form is available, actual substitution may very well be undesirable.

The generalized reduced gradient algorithm consists of the following steps.

Step 1. Assume that some that feasible 𝐗⁰ is known. Set 𝑘=0 and go to the next step.

Step 2. Step 2 is conveniently divided into substeps. (1.1) Compute the Jacobian 𝐀 and the gradient of the objective function.(1.2) Determine a splitting of 𝐗 into (𝐱,𝐲) and corresponding of 𝐀 into (𝐍,𝐁), such that 𝐲^𝑘 is strictly between bound and 𝐁 is nonsingular; invert 𝐁.(1.3) Compute the Lagrange multipliers 𝐮 and the reduced gradient 𝐠. (1.4) Determine the face 𝐅 and the projection 𝐩 of 𝐠 onto 𝐅 [16]. (1.5) If 𝐩 is zero (or almost zero in some sense), then terminate. 𝐗^𝑘 is a KKT point [17, 18]. Otherwise, go to the next step.

Step 3. Choose ascent direction 𝐝=𝐩^𝑡 [18], that is, 𝑑_𝑗=0if𝑔_𝑗>0,𝑢_𝑗=𝑎_𝑗,0if𝑔_𝑗<0,𝑢_𝑗=𝑏_𝑗,𝑔_𝑗,otherwise. (3.10)

Step 4. Choose a first stepsize 𝜂₁.

Step 5. Maximize, with respect to 𝜂, the function 𝜓𝜂=𝐟𝐱^𝑘+𝜂𝐝,𝐲𝐱^𝑘+𝜂𝐝=𝜙𝐱^𝑘+𝜂𝐝 (3.11) with more or less accuracy (the linear search). For each value of 𝜂 under consideration, this step requires solving the following system of 𝑚 equations: 𝐡𝐱^𝑘+𝜂𝐝,𝐲=0, (3.12) where 𝐲 is the 𝑚-dimensional vector of unknowns.

Step 6. Assuming that Step 5 succeeds, an improved feasible point is obtained, which replaced 𝐗^𝑘. Replace 𝑘 by 𝑘+1 and go back to Step 2.

Let us recall briefly the essential points of the generalized reduced gradient method: an initial feasible guess 𝐗⁰∈𝐂 is given and a sequence {𝐗^𝑘}_𝑘≥0⊂𝐂 is generated by using iterations of the general form:∀𝑘≥0𝐗^𝑘+1=𝐐_𝑘𝐗^𝑘=𝐗^𝑘+𝜂_𝑘𝐝^𝑘. (3.13)

Remark 3.1. The theoretical convergence of generalized reduced gradient method has been proved by [17, 19, 20].

4. Two-Phase Generalized Reduced Gradient (TPGRG) Method

The main difficulty remains the lack of concaveity: if 𝐟 is not concave, the Kuhn-Tucker points may not correspond to a global maximum (see e.g., [7, 8]). In the next, this point is improved by using an appropriate random perturbation.

The sequence of real numbers {𝐗^𝑘}_𝑘≥0 is replaced by a sequence of random variables {𝐗^𝑘}_𝑘≥0 involving a random perturbation 𝒫_𝑘 of the deterministic iteration (3.13); then we have 𝐗⁰=𝐗⁰;∀𝑘≥0𝐗^𝑘+1=Q_𝑘𝐗^𝑘+𝒫_𝑘=𝐗^𝑘+𝜂_𝑘𝐝^𝑘+𝒫_𝑘=𝐗^𝑘+𝜂_𝑘𝐝^𝑘+𝒫_𝑘𝜂_𝑘, (4.1) where 𝜂_𝑘≠0 satisfied the Step 5 in GRG algorithm, and∀𝑘≥1∶𝒫_𝑘isindependentfrom𝐗^𝑘−1,…,𝐗⁰, (4.2) see [9], and 𝐗∈𝐂⟹𝐐_𝑘𝐗+𝒫_𝑘∈𝐂. (4.3) Equation (4.1) can be viewed as perturbation of the ascent direction 𝐝^𝑘, which is replaced by a new direction 𝐃_𝑘=𝐝^𝑘+𝒫_𝑘/𝜂_𝑘, and the iterations (4.1) become 𝐗^𝑘+1=𝐗^𝑘+𝜂_𝑘𝐃_𝑘. (4.4) General properties defining convenient sequences of perturbation {𝒫_𝑘}_𝑘≥𝑂 can be found in the literature [13, 14]: usually, sequence of Gaussian laws may be used in order to produce elements satisfying these properties.

We introduce a random vector 𝒵_𝑘, and we denote by Φ_𝑘 and 𝜙_𝑘 the cumulative distribution function and the probability density of 𝒵_𝑘, respectively.

We denote by 𝐅_𝑘+1(𝐘∣𝐗^𝑘=𝐗) the conditional cumulative distribution function:𝐅_𝑘+1𝐘∣𝐗^𝑘=𝐗=𝑃𝐗^𝑘+1<𝐘∣𝐗^𝑘=𝐗, (4.5) and the condition probability density of 𝐗^𝑘+1 is denoted by 𝐟_𝑘+1.

Let us introduce a sequence of 𝑛-dimensional random vectors {𝐙_𝑘}_𝑘≥0∈𝐂. We consider also {𝜉_𝑘}_𝑘≥0, a suitable decreasing sequence of strictly positive real numbers converging to 0 and such that 𝜉₀≤1.

The optimal choice for 𝜂_𝑘 is determined by Step 5. Let 𝒫_𝑘=𝜉_𝑘𝐙_𝑘:𝐅_𝑘+1𝐘∣𝐗^𝑘=𝐗=𝑃𝐗^𝑘+1<𝐘∣𝐗^𝑘=𝐗. (4.6) It follows that 𝐅_𝑘+1𝐘∣𝐗^𝑘=𝐗=𝑃𝐙_𝑘<𝐘−𝐐_𝑘(𝐗)𝜉_𝑘=Φ_𝑘𝐘−𝐐_𝑘(𝐗)𝜉_𝑘. (4.7) So, we have𝐟_𝑘+1𝐘∣𝐗^𝑘=1=𝐗𝜉^𝑛_𝑘𝜙_𝑘𝐘−𝐐_𝑘(𝐗)𝜉_𝑘𝐘∈𝐂. (4.8) The relation (2.11) shows that 𝐘−𝐐_𝑘(𝐗)≤𝛽(𝛼,𝜀)for(𝐗,𝐘)∈𝐂_𝛼×𝐂_𝛼. (4.9) We assume that there exists a decreasing function 𝑡↦𝐠_𝑘(𝑡), 𝐠_𝑘(𝑡)>0 on ℝ⁺ such that𝐘∈𝐂_𝛼⟹𝜙_𝑘𝐘−𝐐_𝑘(𝐗)𝜉_𝑘≥𝐠_𝑘𝛽(𝛼,𝜀)𝜉_𝑘. (4.10) For simplicity, let𝐙_𝑘=1_𝐂𝒵_𝑘𝒵_𝑘, (4.11) where 𝒵 is a random variable; for simplicity let 𝒵∼𝐍(𝟎,1).

The procedure generates a sequence 𝐔_𝑘=𝑓(𝐗^𝑘). By construction this sequence is increasing and upper bounded by 𝐥^∗: ∀𝑘≥0∶𝐥^∗≥𝐔_𝑘+1≥𝐔_𝑘. (4.12) Thus, there exists 𝐔≤𝐥^∗ such that 𝐔_𝑘⟶𝐔for𝑘⟶+∞. (4.13)

Lemma 4.1. Let 𝒫_𝑘=𝜉_𝑘𝐙_𝑘 and 𝛾=𝐟(𝐗⁰) if 𝐙_𝑘 is given by (4.11); then there exists 𝜈>0 such that 𝑃𝐔_𝑘+1>𝜃∣𝐔_𝑘≥𝐂≤𝜃meas_𝛾−𝐂_𝜃𝜉^𝑛_𝑘𝐠_𝑘𝛽𝛾,𝜀𝜉_𝑘>0∀𝜃∈𝐥^∗−𝜈,𝐥^∗, (4.14) where 𝐧=dim(𝐄).

Proof. Let 𝐂_𝜽={𝐗∈𝐂∣𝐟(𝐗)>𝜃}, for 𝜃∈[𝐥^∗−𝜈,𝐥^∗).

Since 𝐂_𝛼⊂𝐂_𝜽,𝐥^∗>𝛼>𝜃, it follows from (2.7) that 𝐂_𝜃 is not empty and has a strictly positive measure.

If meas(𝐂−𝐂_𝜃)=0 for any 𝜃∈[𝐥^∗−𝜈,𝐥^∗), the result is immediate, since we have 𝑓(𝐗)=𝐥^∗ on 𝐂.

Let us assume that there exists 𝜀>0 such that meas(𝐂−𝐂_𝜀)>0. For 𝜃∈[𝐥^∗−𝜀,𝐥^∗), we have 𝐂_𝜽⊂𝐂_𝜀 and meas(𝐂−𝐂_𝜽)>0.

𝑃(𝐗^𝑘∉𝐂_𝜽)=𝑃(𝐗^𝑘∈𝐂−𝐂_𝜽)=∫_𝐂−𝐂_𝜽𝑃(𝐗^𝑘∈𝐝𝐗)>0 for any 𝜃∈[𝐥^∗−𝜀,𝐥^∗), since the sequence {𝐔_𝑖}_𝑖≥0 is increasing, and we have also𝐗^𝑖_𝑖≥0⊂𝐂_𝛾. (4.15) Thus 𝑃𝐗^𝑘∉𝐂_𝜽=𝑃𝐗^𝑘∈𝐂_𝛾−𝐂_𝜽=_𝐂𝜸−𝐂_𝜽𝑃𝐗^𝑘∈𝐝𝐗>0forany𝜃∈𝐥^∗−𝜀,𝐥^∗. (4.16) Let 𝜃∈[𝐥^∗−𝜀,𝐥^∗); we have from (4.12) 𝑃𝐔_𝑘+1>𝜃∣𝐔_𝑘≤𝜃=𝑃𝐗^𝑘+1∈𝐂_𝜽∣𝐗^𝑖∉𝐂_𝜽,𝑖=0,…,𝑘. (4.17) But Markov chain yields that 𝑃𝐗^𝑘+1∈𝐂_𝜽∣𝐗^𝑖∉𝐂_𝜽,𝑖=0,…,𝑘=𝑃𝐗^𝑘+1∈𝐂_𝜽∣𝐗^𝑘∉𝐂_𝜽. (4.18) By the conditional probability rule 𝑃𝐗^𝑘+1∈𝐂_𝜽∣𝐗^𝑘∉𝐂_𝜽=𝑃𝐗^𝑘+1∈𝐂_𝜽,𝐗^𝑘∉𝐂_𝜽𝑃𝐗^𝑘∉𝐂_𝜽. (4.19) Moreover 𝑃𝐗^𝑘+1∈𝐂_𝜽,𝐗^𝑘∉𝐂_𝜽=_𝐂−𝐂_𝜽𝑃𝐗^𝑘𝐂∈𝐝𝐗_𝜽𝐟_𝑘+1𝐘∣𝐗^𝑘=𝐗𝐝𝐘. (4.20) From (4.15),we have 𝑃𝐗^𝑘+1∈𝐂_𝜽,𝐗^𝑘∉𝐂_𝜽=_𝐂𝛾−𝐂_𝜽𝑃𝐗^𝑘𝐂∈𝐝𝐗_𝜽𝐟_𝑘+1𝐘∣𝐗^𝑘𝑃𝐗=𝐗𝐝𝐘,_𝑘+1∈𝐂_𝜽,𝐗^𝑘∉𝐂_𝜽≥inf_{𝐗∈𝐂𝛾−}𝐂_𝜽𝐂_𝜽𝐟_𝑘+1𝐘∣𝐗^𝑘=𝐗𝐝𝐘_𝐂𝛾−𝐂_𝜽𝑃𝐗^𝑘∈𝐝𝐗. (4.21)(4.22) Thus 𝑃𝐗^𝑘+1∈𝐂_𝜽∣𝐗^𝑘∉𝐂_𝜽≥inf_{𝐗∈𝐂𝛾−}𝐂_𝜽𝐂_𝜽𝐟_𝑘+1𝐘∣𝐗^𝑘=𝐗𝐝𝐘. (4.23) Taking (4.8) into account, we have 𝑃𝐗^𝑘+1∈𝐂_𝜽∣𝐗^𝑘∉𝐂_𝜽≥1𝜉^𝑛_𝑘inf_{𝐗∈𝐂𝛾−}𝐂_𝜽𝐂_𝜽𝜙_𝑘𝐘−𝐐_𝑘(𝐗)𝜉_𝑘𝐝𝐘. (4.24) The relation (2.11) shows that 𝐘−𝐐_𝑘(𝐗)≤𝛽𝛾,𝜀, (4.25) and (4.10) yields that 𝜙_𝑘𝐘−𝐐_𝑘(𝐗)𝜉_𝑘≥𝐠_𝑘𝛽𝛾,𝜀𝜉_𝑘. (4.26) Hence 𝑃𝐗^𝑘+1∈𝐂_𝜽∣𝐗^𝑘∉𝐂_𝜽≥1𝜉^𝑛_𝑘inf_{𝐗∈𝐂𝛾−}𝐂_𝜽𝐂_𝜽𝐠_𝑘𝛽𝛾,𝜀𝜉_𝑘𝑃𝐗𝐝𝐘,_𝑘+1∈𝐂_𝜽∣𝐗^𝑘∉𝐂_𝜽≥meas𝐂_𝛾−𝐂_𝜃𝜉^𝑛_𝑘𝐠_𝑘𝛽𝛾,𝜀𝜉_𝑘. (4.27)

4.1. Global Convergence

The global convergence is a consequence of the following result, which yields from the Borel-Catelli's lemma (e.g., see [13]):

Lemma 4.2. Let {𝐔_𝑘}_𝑘≥0 be an increasing sequence, upper bounded by 𝐥^∗. Then, there exists 𝐔 such that 𝐔_𝑘→𝐔 for 𝑘→+∞. Assume that there exists 𝜈>0 such that for any 𝜃∈[𝐥^∗−𝜈,𝐥^∗), there is a sequence of strictly positive real numbers {𝐜_𝑘(𝜃)}_𝑘≥0 such that ∀𝑘≥0∶𝑃𝐔_𝑘+1>𝜃∣𝐔_𝑘≤𝜃≥𝐜_𝑘(𝜃)>0,^+∞_𝑘=0𝐜_𝑘(𝜃)=+∞. (4.28) Then 𝐔=𝐥^∗ almost surely.

Proof. For instance, see [13, 21].

Theorem 4.3. Let 𝛾=𝐟(𝐗⁰); assume that 𝐗⁰∈𝐂, the sequence 𝜉_𝑘 is nonincreasing, and ^+∞_𝑘=0𝐠_𝑘𝛽𝛾,𝜀𝜉_𝑘=+∞. (4.29) Then 𝐔=𝐥^∗ almost surely.

Proof. Let 𝐜_𝑘(𝜃)=meas𝐂_𝛾−𝐂_𝜃𝜉^𝑛_𝑘𝐠_𝑘𝛽𝛾,𝜀𝜉_𝑘>0. (4.30) Since the sequence {𝜉_𝑘}_𝑘≥0 is nonincreasing, 𝐜_𝑘(𝜃)≥meas𝐂_𝛾−𝐂_𝜃𝜉^𝑛₀𝐠_𝑘𝛽𝛾,𝜀𝜉_𝑘>0. (4.31) Thus, (4.29) shows that ^+∞_𝑘=0𝐜_𝑘(𝜃)≥meas𝐂_𝛾−𝐂_𝜃𝜉^𝑛₀^+∞_𝑘=0𝐠_𝑘𝛽𝛾,𝜀𝜉_𝑘=+∞. (4.32) Using Lemmas 4.1 and 4.2 we have 𝐔=𝐥^∗ almost surely.

Theorem 4.4. Let 𝐙_𝑘 be defined by (4.11), and let 𝜉_𝑘=𝐚𝐝,log𝑘+ (4.33) where 𝐚>0, 𝐝>0, and 𝑘 is the iteration number. If 𝐗⁰∈𝐶, then, for 𝐚 large enough, 𝐔=𝑙^∗ almost surely.

Proof. We have 𝜙_𝑘(𝐙)=1√2𝜋^𝑛−1exp2‖𝐙‖²=𝐠_𝑘(‖𝐙‖)>0. (4.34) So, 𝐠_𝑘𝛽𝛾,𝜀𝜉_𝑘=1√2𝜋^𝑛𝑘+𝐝^{𝛽(𝛾,𝜀)2/(2}𝐚₎. (4.35) For 𝐚 such that 𝛽0<𝛾,𝜀²2𝐚<1, (4.36) we have ^∞_𝑘=0𝐠_𝑘𝛽𝛾,𝜀𝜉_𝑘=+∞, (4.37) and, from Theorem 4.4, we have 𝐔=𝑙^∗ almost surely.

4.2. Practical Implementation

The above results suggest the following numerical algorithm.(1)An initial guess 𝐗⁰∈𝐶 is given. (2)At the iteration number 𝑘≥0, 𝐗^𝑘 is known and 𝐗^𝑘+1 is determined by performing the following three substeps.(2.1) Unperturbed ascent: we determine the descent direction 𝐝^𝑘 and the step 𝜂^𝑘 using ascent method (3.13). This generates the first trial point: 𝐓⁰_𝑘+1=𝐐_𝑘𝐗^𝑘. (4.38)(2.2) Perturbation: we determine a sample (𝒫¹_𝑘,…,𝒫^𝑘sto_𝑘) of 𝑘_sto new trial points: 𝐓^𝑖_𝑘+1=𝐓⁰_𝑘+1+𝒫^𝑖_𝑘,𝑖=1,…,𝑘_sto. (4.39)(2.3) Dynamics: we determine 𝐗^𝑘+1 by selecting it from the set of available points: 𝒜_𝑘=𝐗^𝑘,𝐓⁰_𝑘+1,…,𝐓^𝑘sto_𝑘+1. (4.40)

The computation of 𝐗^𝑘+1 is performed in two phases. First phase: we determine a trial point (unperturbed ascent step and perturbation step). Second phase: we determine 𝐗^𝑘+1 by selection it from 𝒜_𝑘 (dynamics step).

As was shown in Theorem 4.4, Substep (2.2) may use 𝒫^𝑖_𝑘=𝜉_𝑘𝐙_𝑘, where 𝒵=(𝐙₁,…,𝐙_𝑘+1) is a sample of 𝑁(0,1) and 𝜉_𝑘 is given by (4.33).

For instance, we can consider elitistic dynamics:𝐗^𝑘+1=argmin_{𝐗∈𝒜𝑘}𝐟(𝐗). (4.41)

5. Numerical Results

In order to apply the method presented in (4.41), we start at the initial value 𝐗⁰=𝐗⁰∈𝐶. At step 𝑘≥0, 𝐗^𝑘 is known and 𝐗^𝑘+1 is determined.

We generate 𝑘_sto the number of perturbation; the iterations are stopped when 𝐗^𝑘 is a Kuhn-Tucker point. We denote by 𝑘_end the value of 𝑘 when the iterations are stopped (it corresponds to the number of evaluations of the gradient of 𝐟). The optimal value and optimal point are 𝐟_opt and 𝐗^opt, respectively. The perturbation is normally distributed and samples are generated by using the log-trigonometric generator and the standard random number generator of the FORTRAN library. We use 𝜉_𝑘=𝐚/log(𝑘+2), where 𝐚>0.

Concern experiments performed on a workstation HP Intel(R) Celeron(R) M processor 1.30 GHz, 224 Mo RAM. The row cpu gives the mean CPU time in seconds for one run.

5.1. Octagon Problem

Consider polygon in the plane with 5 sides (5-gons for short) and unit diameter. Which of them have maximum area ? (see e.g., [22]).

Graham's conjecture states that the optimal octagon can be illustrated as in Figure 1, in which a solid line between two vertices indicates that the distance between these points is one.

Figure 1 

Definition of variables for the configuration conjectured by Graham to have maximum area.

This question can be formulated to the quadratically constrained quadratic optimization problem defining this configuration that appears as follows:1Maximize2𝐱₂+𝐱₃−4𝐱₁𝐲₁+3𝐱₁−2𝐱₃+𝐱₅𝐲₂+3𝐱₁−2𝐱₂+𝐱₄𝐲₃+𝐱₃−2𝐱₁𝐲₄+𝐱₂−2𝐱₁𝐲₅−𝐱₁𝐀Subjectto₀−𝐀₁≤1∶𝐱₁−𝐱₂²+𝐲₁−𝐲₂²+𝐳₁𝐀=1,₀−𝐀₂≤1∶−𝐱₁+𝐱₂−𝐱₅²+𝐲₁−𝐲₃+𝐲₅²+𝐳₂𝐀=1,₀−𝐀₆≤1∶𝐱₁−𝐱₂+𝐱₄²+𝐲₁−𝐲₂+𝐲₄²+𝐳₃𝐀=1,₀−𝐀₇≤1∶−𝐱₁+𝐱₃²+𝐲₁−𝐲₃²+𝐳₄𝐀=1,₁−𝐀₂≤1∶2𝐱₁−𝐱₂−𝐱₃+𝐱₅²+𝐲₂−𝐲₃+𝐲₅²+𝐳₅𝐀=1,₁−𝐀₃≤1∶2𝐱₁+𝐱₂²+𝐲²₂+𝐳₆𝐀=1,₁−𝐀₄≤1∶𝐱₁−𝐱₂²+𝐲₁−𝐲₂−1²+𝐳₇𝐀=1,₁−𝐀₇≤1∶2𝐱₁−𝐱₂−𝐱₃²+−𝐲₂+𝐲₃²+𝐳₈𝐀=1,₂−𝐀₃≤1∶𝐱₃−𝐱₅²+−𝐲₃+𝐲₅²+𝐳₉𝐀=1,₂−𝐀₄≤1∶−𝐱₁+𝐱₃−𝐱₅²+𝐲₁−𝐲₃+𝐲₅−1²+𝐳₁₀𝐀=1,₂−𝐀₅≤1∶2𝐱₁+𝐱₃−𝐱₅²+−𝐲₃+𝐲₅²+𝐳₁₁𝐀=1,₂−𝐀₆=1∶2𝐱₁−𝐱₂−𝐱₃+𝐱₄+𝐱₅²+−𝐲₂+𝐲₃+𝐲₄−𝐲₅²𝐀=1,₃−𝐀₆≤1∶−2𝐱₁+𝐱₂−𝐱₄²+𝐲₂−𝐲₄²+𝐳₁₂𝐀=1,₄−𝐀₆≤1∶𝐱₁−𝐱₂+𝐱₄²+𝐲₁−𝐲₂+𝐲₄−1²+𝐳₁₃𝐀=1,₄−𝐀₇≤1∶𝐱₁−𝐱₃²+1−𝐲₁+𝐲₃²+𝐳₁₄𝐀=1,₅−𝐀₆≤1∶𝐱₂−𝐱₄²+𝐲₂−𝐲₄²+𝐳₁₅𝐀=1,₅−𝐀₇≤1∶2𝐱₁−𝐱₃²+𝐲²₃+𝐳₁₆𝐀=1,₆−𝐀₁≤1∶2𝐱₁−𝐱₂−𝐱₃+𝐱₄²+−𝐲₂+𝐲₃+𝐲₄²+𝐳₁₇𝐱=1,₂−𝐱₃−𝐳₁₈𝐱=0,_2𝑖+𝐲²_𝑖𝐱=1,𝑖=1,2,3,4,5,₁+𝐳₁₉𝐱=0.5,_𝑖+𝐳_18+𝑖=1,𝑖=2,3,4,5,0≤𝐱_𝑖,𝐲_𝑖𝑖=1,2,…,5,0≤𝐳_𝑖𝑖=1,2,…,23. (5.1)

There are 33 variables and 34 constraints.

We use 𝐚=10𝑑+6 and 𝑘_sto=15000 and used regular octagon as a starting point (area of regular octagon ≈0.7071 see Figure 2), and if we let the error less than 10⁻⁵.

Figure 2 

Regular octagon.

The Fortran code of TPGRG furnishs the following optimal solutions, 𝑓_opt≈0.72687:𝐗^opt=(),𝐲0.25949,0.66664,0.66664,0.90703,0.90730_opt𝑖=1−𝐱^opt_𝑖,𝑘_end=626,CPUtime=4684s(1.5hours). (5.2)

Branch and cut method [23] was solving the octagon problem with using regular octagon as a starting point, the maximum area found is ≈0.72687, but the cpu time is more than 30 hours.

5.2. Mixture Problem

In this example of petrochemical mixture, we have four reservoirs: 𝐑₁, 𝐑₂, 𝐑₃, and 𝐑₄. The two first receive three distinct source products, and then their content is combined in the two other in order to create the wanted miscellanies. The question is to determine the quantity of every product to buy in order to maximize the profits (see e.g., [24, 25]).

The 𝐑₁ reservoir receives two quality products 3 and 1 in quantities 𝐪₀ and 𝐪₁, respectively. The 𝐑₂ reservoir contains a product of the third source, of quality 2 in quantity 𝐪₂. One wants to get in the reservoirs 𝐑₃ and 𝐑₄ of capacities 10 and 20, respectively.

Figure 3, where the variable 𝐱₁ to 𝐱₂ represents some quantities, illustrates this situation. The unit prices of the products bought are, respectively, 60, 160, and 100; those of the products finished are 90 and 150. Therefore the difference between the costs of purchase and sale is 60𝐪₀+160𝐪₁+100𝐪₂−90𝐱₁+𝐱₃−150𝐱₂−𝐱₄. (5.3) The quality of the mixture done contained in 𝐑₁ is𝐱₅=3𝐪₀+𝐪₁𝐪₀+𝐪₁∈[1,3]. (5.4) The qualities of the final miscellanies are𝐱₁𝐱₅+2𝐱₃𝐱₁+𝐱₃𝐱≤2.5,²𝐱₅+2𝐱₄𝐱₂+𝐱₄≤1.5. (5.5) The addition of the constraints of volume conservation 𝐪₀+𝐪₁=𝐱₁+𝐱₂ and 𝐪₂=𝐱₃+𝐱₄ does permit the elimination of the variables 𝐪₀,𝐪₁, and 𝐪₂: 𝐪₀=𝐱₅−1𝐱₁+𝐱₂2,𝐪₁=3−𝐱₅𝐱₁+𝐱₂2,𝐪₂=𝐱₃+𝐱₄. (5.6)

Figure 3 

A Problem of mixture.

Remark 5.1. One has 1≤𝐱₅≤1.5. (5.7) By the relation (5.5), 𝐱₂𝐱₅≤1.5𝐱₂. (5.8) Suppose that 𝐱₅>1.5, then 𝐱₂𝐱₅>1.5𝐱₂, and so 𝐱₂𝐱₅>1.5𝐱₂−0.5𝐱₄ that contradict with inequality (5.8).

The mathematical model transformed is the following: Maximize−120𝐱₁−60𝐱₂−10𝐱₃+50𝐱₄+50𝐱₁𝐱₅+50𝐱₂𝐱₅Subjectto2.5𝐱₁+0.5𝐱₃−𝐱₁𝐱₅−𝐳₁=0,1.5𝐱₂−0.5𝐱₄−𝐱₂𝐱₅−𝐳₂𝐱=0,₁+𝐱₃+𝐳₃𝐱=10,₂+𝐱₄+𝐳₄𝐱=20,₅+𝐳₅𝐱=1.5,₅−𝐳₆=1,0≤𝐱_𝑖𝑖=1,…,5;0≤𝐳_𝑖𝑖=1,…,6. (5.9) There are 11 variables and 10 constraints.

We use 𝐚=10𝑑10, 𝑘_sto=20000, and 𝐗⁰=(2,9,0,8,1),𝐟𝐗⁰=170. (5.10) The Fortran code of TPGRG furnishes the following optimal solutions: cpu = 100 second, and 𝑘_end=1390,𝐗_opt=(0,10,0,10,1),𝐟_opt=400, (5.11) since the mixture problem is known by this global solution.

6. Concluding Remarks

A two-phase generalized reduced gradient method is presented for nonlinear constraints, involving the adjunction of a stochastic perturbation. This approach leads to a stochastic ascent method where the deterministic sequence generated by the two-phase generalized gradient method is replaced by a sequence of random variables.

TPGRG method converges to global maximum for all differential objective function, but GRG method converges to local maximum.

The numerical experiments show that the method is effective to calculate for global optimization problems. Here yet, we observe that the adjunction of the stochastic perturbation improves the result, with a larger number of evaluations of the objective function. The main difficulty in the practical use of the stochastic perturbation is connected to the tuning of the parameters 𝐚 and 𝑘_sto.