# Two-Phase Generalized Reduced Gradient Method for Constrained Global Optimization

Authors: Abdelkrim El Mouatasim
Journal: Journal of Applied Mathematics
 (2010)

---

## 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.


---

## Body

## 1. Introduction

We consider the problem(1.1)Maximize:f(X)subjectto:h(X)=0,X∈Ω,
where f:ℝn→ℝ and h:ℝn→ℝm are continuously differentiable, and Ω⊂ℝn is closed and convex set (e.g., Ω=∏i=1n[ai,bi]).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,f 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 {Xk}k≥0, where X0 is an initial feasible point and, for each k>0, a new feasible point Xk+1 is generated from Xk by using an operator Qk (see Section 3). Thus the iterations are given by(1.2)k≥0:Xk+1=Qk(Xk).A fundamental difficulty arises due to the lack of concavity: the convergence of the sequence{Xk}k≥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{Xk}k≥0 is replaced by a random vectors sequence {X̂k}k≥0 and the iterations are modified as follows:(1.3)k≥0:X̂k+1=Qk(X̂k)+Pk,
where 𝒫k is a suitable random variable, called the stochastic perturbation. The sequence {𝒫k}k≥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)
E=ℝn is the n-dimensional real Euclidean space.(ii)
X stands for (x1,…,xn)T∈E.(iii)
h(X) is the column vector whose components are h1(X),…,hm(X).(iv)
∥X∥=(x12+⋯+xn2)1/2 is the Euclidean norm of X:
(2.1)‖A‖∞=sup⁡⁡{‖AX‖:‖X‖=1},(v)
AT is the transpose matrix associated to A.Definition 2.1.
A pointXis said to be feasible if|hi(X)|≤ɛ0, fori=1,…,m, whereɛ0is same preassigned small positive constant.Assume that we know some feasible pointX0. We make the following nonedegeneracy assumption. The vector X can be split into two components: an m-dimensional component y (the basic part), and x, a component of dimension n-m (the nonbasic part) such that the following two properties hold:  (H1)
y0 is strictly between bounds;(H2)
the squarem×m matrix ∂h/∂y, computed at X0, is nonsingular.If the property (H2) does not hold orn<m, we add artificial variables to the constraints, so the property (H2) holds and n≥m in associate problem (1.1).Let(2.2)C={X∈E∣h(X)=0,X∈Ω}.
The objective function is f:E→ℝ, and its upper bound on C is denoted by l*: (2.3)l*=max⁡⁡Cf.
Let us introduce (2.4)Cα=Sα∩C;whereSα={X∈E∣f(X)≥α}.
We assume that(2.5)fistwicecontinuouslydifferentiableonE,(2.6)∀α<l*:Cαisnotempty,closedandbounded,(2.7)∀α<l⋆:meas⁡(Cα)>0,
where meas(Cα) is the measure of Cα.SinceE is a finite dimensional space, assumption (2.6) is verified when C is bounded or f is coercive, that is, lim∥X∥→+∞f(X)=-∞. Assumption (2.6) is verified when C contains a sequence of neighborhoods of a point of optimum X⋆ having strictly positive measure, that is, when X⋆ can be approximated by a sequence of points of the interior of C. We observe that the assumptions (2.5)-(2.6) yield that (2.8)C=⋃α<l⋆Cα,thatis,forallX∈C:∃α<l*suchthatX∈Cα.
From (2.5)-(2.6),(2.9)γ1=sup⁡⁡{‖∇f(X)‖:X∈Cα}<+∞.
Then (2.10)γ2=sup⁡⁡{‖d‖:X∈Cα}<+∞,
where d is the direction determined on the basis of the gradient ∇f(X).Thus,(2.11)β(α,ɛ)=sup⁡⁡{‖Y-(X+ηd)‖:(X,Y)∈Cα×Cα,0≤η≤ɛ}<+∞,
where ɛ,η are positive real numbers.


## 3. Generalized Reduced Gradient Method

By the implicit function theorem, there exists, in some neighborhoodV of x0, a unique continuous function (mapping), say y(x), such that f(x,y(x) is identically zero in V. In addition, y(x) has a continuous derivative dy/dx, which can be computed by the chain rule: (3.1)∂h∂x+(∂h∂y)(dydx)=0,
or, more conveniently,(3.2)dydx=-(∂h∂y)-1(∂h∂x).
In what follows, we call A the Jacobian of h(X) computed at X0. Similarly, we set (3.3)B=∂h∂y,N=∂h∂x(computedatX0).
Substituting y(x) into the objective function f(x,y), we obtain the reduced function:(3.4)ϕ(x)=f(x,y(x)),
the gradient of which at x0 is, by the chain rule again,(3.5)g=∂f∂x+(∂f∂y)(dydx)
(all derivatives computed at X0). Setting (3.6)v=∂f∂x,w=∂f∂y(computedatX0),
we have the following formula for the reduced gradientg: (3.7)g=v-wB-1N.
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 gradientp (in the space of g) by its components:  (i)
pj=0 if xj0=aj and gj<0;(ii)
pj=0 if xj0=bj and gj>0;(iii)
pj=gj, otherwise.It is convenient to set(3.8)u=wB-1,
and hence (3.9)w-uB=0,v-uN=g.
It is a simple matter to verify that Kuhn-Tucker condition for the problem in standard form reduces to p=0, and that u is the row-vector of multipliers corresponding to the equations h(X)=0. We assume from now on that p≠0.The following relationsxj=aj if xj0=aj and gj<0,xj=bj if xj0=bj and gj>0define what we call theface (at X0), denoted by F. The row-vector p (the projected reduced gradient) is also named the projection of g onto F.Letd be any nonzero vector column-vector in F such that pd>0; the vector d is an ascent direction for the reduced function ϕ(x).There is a striking analogy with what is usually done in linear programming, wherey(x) can be computed in close form. This is generally not the case if the constraints h(X) 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 feasibleX0 is known. Set k=0 and go to the next step.Step 2.
Step2 is conveniently divided into substeps.  (1.1)
Compute the JacobianA and the gradient of the objective function.(1.2)
Determine a splitting ofX into (x,y) and corresponding of A into (N,B), such that yk is strictly between bound and B is nonsingular;invertB.(1.3)
Compute the Lagrange multipliersu and the reduced gradient g.(1.4)
Determine the faceF and the projection p of g onto F [16].(1.5)
Ifp is zero (or almost zero in some sense), then terminate. Xk is a KKT point [17, 18]. Otherwise, go to the next step.Step 3.
Choose ascent directiond=pt [18], that is,
(3.10)dj={0ifgj>0,uj=aj,0ifgj<0,uj=bj,gj,otherwise.Step 4.
Choose a first stepsizeη1.Step 5.
Maximize, with respect toη, the function
(3.11)ψ(η)=f(xk+ηd,y(xk+ηd))=ϕ(xk+ηd)
with more or less accuracy (the linear search). For each value of η under consideration, this step requires solving the following system of m equations:
(3.12)h(xk+ηd,y)=0,
where y is the m-dimensional vector of unknowns.Step 6.
Assuming that Step5 succeeds, an improved feasible point is obtained, which replaced Xk. Replace k by k+1 and go back to Step 2.Let us recall briefly the essential points of the generalized reduced gradient method: an initial feasible guessX0∈C is given and a sequence {Xk}k≥0⊂C is generated by using iterations of the general form:(3.13)∀k≥0Xk+1=Qk(Xk)=Xk+ηkdk.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: iff 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{Xk}k≥0 is replaced by a sequence of random variables {X̂k}k≥0 involving a random perturbation 𝒫k of the deterministic iteration (3.13); then we have X̂0=X0;(4.1)∀k≥0X̂k+1=Qk(X̂k)+Pk=X̂k+ηkdk+Pk=X̂k+ηk(dk+Pkηk),
where ηk≠0 satisfied the Step 5 in GRG algorithm, and(4.2)∀k≥1:Pkisindependentfrom(X̂k-1,…,X̂0),
see [9], and (4.3)X̂∈C⟹Qk(X̂)+Pk∈C.
Equation (4.1) can be viewed as perturbation of the ascent direction dk, which is replaced by a new direction Dk=dk+𝒫k/ηk, and the iterations (4.1) become (4.4)X̂k+1=X̂k+ηkDk.
General properties defining convenient sequences of perturbation {𝒫k}k≥O 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𝒵k, and we denote by Φk and ϕk the cumulative distribution function and the probability density of 𝒵k, respectively.We denote byFk+1(Y∣X̂k=X) the conditional cumulative distribution function:(4.5)Fk+1(Y∣X̂k=X)=P(X̂k+1<Y∣X̂k=X),
and the condition probability density of X̂k+1 is denoted by fk+1.Let us introduce a sequence ofn-dimensional random vectors {Zk}k≥0∈C. We consider also {ξk}k≥0, a suitable decreasing sequence of strictly positive real numbers converging to 0 and such that ξ0≤1.The optimal choice forηk is determined by Step 5. Let 𝒫k=ξkZk:(4.6)Fk+1(Y∣X̂k=X)=P(X̂k+1<Y∣X̂k=X).
It follows that (4.7)Fk+1(Y∣X̂k=X)=P(Zk<Y-Qk(X)ξk)=Φk(Y-Qk(X)ξk).
So, we have(4.8)fk+1(Y∣X̂k=X)=1ξknϕk(Y-Qk(X)ξk)Y∈C.
The relation (2.11) shows that (4.9)‖Y-Qk(X)‖≤β(α,ɛ)for(X,Y)∈Cα×Cα.
We assume that there exists a decreasing  functiont↦gk(t), gk(t)>0 on ℝ+ such that(4.10)Y∈Cα⟹ϕk(Y-Qk(X)ξk)≥gk(β(α,ɛ)ξk).
For simplicity, let(4.11)Zk=1C(Zk)Zk,
where 𝒵 is a random variable; for simplicity let 𝒵~N(0,1).The procedure generates a sequenceUk=f(X̂k). By construction this sequence is increasing  and upper bounded by l*: (4.12)∀k≥0:l*≥Uk+1≥Uk.
Thus, there exists U≤l* such that (4.13)Uk⟶Ufork⟶+∞.Lemma 4.1.
Let𝒫k=ξkZk and γ=f(X0) if Zk is given by (4.11); then there exists ν>0 such that
(4.14)P(Uk+1>θ∣Uk≤θ)≥meas⁡(Cγ-Cθ)ξkngk(β(γ,ɛ)ξk)>0∀θ∈[l*-ν,l*),
where n=dim(E).Proof.
LetĈθ={X∈C∣f(X)>θ}, for θ∈[l*-ν,l*).SinceCα⊂Ĉθ,l*>α>θ, it follows from (2.7) that Ĉθ is not empty and has a strictly positive measure.Ifmeas(C-Ĉθ)=0 for any θ∈[l*-ν,l*), the result is immediate, since we have f(X)=l* on C.Let us assume that there existsɛ>0 such that meas(C-Ĉɛ)>0. For θ∈[l*-ɛ,l*), we have Ĉθ⊂Ĉɛ and meas(C-Ĉθ)>0.P(X̂k∉Ĉθ)=P(X̂k∈C-Ĉθ)=∫C-ĈθP(X̂k∈dX)>0 for any θ∈[l*-ɛ,l*), since the sequence {Ui}i≥0 is increasing, and we have also(4.15){X̂i}i≥0⊂Cγ.
Thus (4.16)P(X̂k∉Ĉθ)=P(X̂k∈Cγ-Ĉθ)=∫Cγ-ĈθP(X̂k∈dX)>0foranyθ∈[l*-ɛ,l*).
Let θ∈[l*-ɛ,l*); we have from (4.12) (4.17)P(Uk+1>θ∣Uk≤θ)=P(X̂k+1∈Ĉθ∣X̂i∉Ĉθ,i=0,…,k).
But Markov chain yields that (4.18)P(X̂k+1∈Ĉθ∣X̂i∉Ĉθ,i=0,…,k)=P(X̂k+1∈Ĉθ∣X̂k∉Ĉθ).
By the conditional probability rule (4.19)P(X̂k+1∈Ĉθ∣X̂k∉Ĉθ)=P(X̂k+1∈Ĉθ,X̂k∉Ĉθ)P(X̂k∉Ĉθ).
Moreover (4.20)P(X̂k+1∈Ĉθ,X̂k∉Ĉθ)=∫C-ĈθP(X̂k∈dX)∫Ĉθfk+1(Y∣X̂k=X)dY.
From (4.15),we have (4.21)P(X̂k+1∈Ĉθ,X̂k∉Ĉθ)=∫Cγ-ĈθP(X̂k∈dX)∫Ĉθfk+1(Y∣X̂k=X)dY,(4.22)P(X̂k+1∈Ĉθ,X̂k∉Ĉθ)≥inf⁡X∈Cγ-Ĉθ{∫Ĉθfk+1(Y∣X̂k=X)dY}∫Cγ-ĈθP(X̂k∈dX).
Thus (4.23)P(X̂k+1∈Ĉθ∣X̂k∉Ĉθ)≥inf⁡X∈Cγ-Ĉθ{∫Ĉθfk+1(Y∣X̂k=X)dY}.
Taking (4.8) into account, we have (4.24)P(X̂k+1∈Ĉθ∣X̂k∉Ĉθ)≥1ξkninf⁡X∈Cγ-Ĉθ{∫Ĉθϕk(Y-Qk(X)ξk)dY}.
The relation (2.11) shows that (4.25)‖Y-Qk(X)‖≤β(γ,ɛ),
and (4.10) yields that (4.26)ϕk(Y-Qk(X)ξk)≥gk(β(γ,ɛ)ξk).
Hence (4.27)P(X̂k+1∈Ĉθ∣X̂k∉Ĉθ)≥1ξkninf⁡X∈Cγ-Ĉθ∫Ĉθgk((β(γ,ɛ)ξk))dY,P(X̂k+1∈Ĉθ∣X̂k∉Ĉθ)≥meas⁡(Cγ-Cθ)ξkngk(β(γ,ɛ)ξk).


### 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{Uk}k≥0 be an increasing sequence, upper bounded by l*. Then, there exists U such that Uk→U for k→+∞. Assume that there exists ν>0 such that for any θ∈[l*-ν,l*), there is a sequence of strictly positive real numbers {ck(θ)}k≥0 such that
(4.28)∀k≥0:P(Uk+1>θ∣Uk≤θ)≥ck(θ)>0,∑k=0+∞ck(θ)=+∞.
Then U=l* almost surely.Proof.
For instance, see [13, 21].Theorem 4.3.
Letγ=f(X0); assume that X0∈C, the sequence ξk is nonincreasing, and
(4.29)∑k=0+∞gk(β(γ,ɛ)ξk)=+∞.
Then U=l* almost surely.Proof.
Let(4.30)ck(θ)=meas⁡(Cγ-Cθ)ξkngk(β(γ,ɛ)ξk)>0.
Since the sequence {ξk}k≥0 is nonincreasing,
(4.31)ck(θ)≥meas⁡(Cγ-Cθ)ξ0ngk(β(γ,ɛ)ξk)>0.
Thus, (4.29) shows that
(4.32)∑k=0+∞ck(θ)≥meas⁡(Cγ-Cθ)ξ0n∑k=0+∞gk(β(γ,ɛ)ξk)=+∞.
Using Lemmas 4.1 and 4.2 we have U=l* almost surely.Theorem 4.4.
LetZk be defined by (4.11), and let
(4.33)ξk=âlog⁡⁡(k+d̂),
where â>0, d̂>0, and k is the iteration number. If X0∈C, then, for â  large enough, U=l* almost surely.Proof.
We have(4.34)ϕk(Z)=1(2π)nexp⁡⁡(-12‖Z‖2)=gk(‖Z‖)>0.
So,
(4.35)gk(β(γ,ɛ)ξk)=1(2π)n(k+d̂)β(γ,ɛ)2/(2â).
For â such that
(4.36)0<β(γ,ɛ)22â<1,
we have
(4.37)∑k=0∞gk(β(γ,ɛ)ξk)=+∞,
and, from Theorem 4.4, we have U=l* almost surely.


### 4.2. Practical Implementation

The above results suggest the following numerical algorithm.(1)
An initial guessX0∈C is given.(2)
At the iteration numberk≥0, Xk is known and Xk+1 is determined by performing the following three substeps.
(2.1)
Unperturbed ascent: we determine the descent directiondk and the step ηk using ascent method (3.13). This generates the first trial point:
(4.38)Tk+10=Qk(Xk).(2.2)
Perturbation: we determine a sample(𝒫k1,…,𝒫kksto) of ksto new trial points:
(4.39)Tk+1i=Tk+10+Pki,i=1,…,ksto.(2.3)
Dynamics: we determineXk+1 by selecting it from the set of available points:
(4.40)Ak={Xk,Tk+10,…,Tk+1ksto}.The computation ofXk+1 is performed in two phases.First phase: we determine a trial point (unperturbed ascent step and perturbation step).Second phase: we determineXk+1 by selection it from 𝒜k (dynamics step).As was shown in Theorem4.4, Substep (2.2) may use 𝒫ki=ξkZk, where 𝒵=(Z1,…,Zk+1) is a sample of N(0,1) and ξk is given by (4.33).For instance, we can consider elitistic dynamics:(4.41)Xk+1=arg⁡min⁡⁡X∈Akf(X).


## 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{Uk}k≥0 be an increasing sequence, upper bounded by l*. Then, there exists U such that Uk→U for k→+∞. Assume that there exists ν>0 such that for any θ∈[l*-ν,l*), there is a sequence of strictly positive real numbers {ck(θ)}k≥0 such that
(4.28)∀k≥0:P(Uk+1>θ∣Uk≤θ)≥ck(θ)>0,∑k=0+∞ck(θ)=+∞.
Then U=l* almost surely.Proof.
For instance, see [13, 21].Theorem 4.3.
Letγ=f(X0); assume that X0∈C, the sequence ξk is nonincreasing, and
(4.29)∑k=0+∞gk(β(γ,ɛ)ξk)=+∞.
Then U=l* almost surely.Proof.
Let(4.30)ck(θ)=meas⁡(Cγ-Cθ)ξkngk(β(γ,ɛ)ξk)>0.
Since the sequence {ξk}k≥0 is nonincreasing,
(4.31)ck(θ)≥meas⁡(Cγ-Cθ)ξ0ngk(β(γ,ɛ)ξk)>0.
Thus, (4.29) shows that
(4.32)∑k=0+∞ck(θ)≥meas⁡(Cγ-Cθ)ξ0n∑k=0+∞gk(β(γ,ɛ)ξk)=+∞.
Using Lemmas 4.1 and 4.2 we have U=l* almost surely.Theorem 4.4.
LetZk be defined by (4.11), and let
(4.33)ξk=âlog⁡⁡(k+d̂),
where â>0, d̂>0, and k is the iteration number. If X0∈C, then, for â  large enough, U=l* almost surely.Proof.
We have(4.34)ϕk(Z)=1(2π)nexp⁡⁡(-12‖Z‖2)=gk(‖Z‖)>0.
So,
(4.35)gk(β(γ,ɛ)ξk)=1(2π)n(k+d̂)β(γ,ɛ)2/(2â).
For â such that
(4.36)0<β(γ,ɛ)22â<1,
we have
(4.37)∑k=0∞gk(β(γ,ɛ)ξk)=+∞,
and, from Theorem 4.4, we have U=l* almost surely.


## 4.2. Practical Implementation

The above results suggest the following numerical algorithm.(1)
An initial guessX0∈C is given.(2)
At the iteration numberk≥0, Xk is known and Xk+1 is determined by performing the following three substeps.
(2.1)
Unperturbed ascent: we determine the descent directiondk and the step ηk using ascent method (3.13). This generates the first trial point:
(4.38)Tk+10=Qk(Xk).(2.2)
Perturbation: we determine a sample(𝒫k1,…,𝒫kksto) of ksto new trial points:
(4.39)Tk+1i=Tk+10+Pki,i=1,…,ksto.(2.3)
Dynamics: we determineXk+1 by selecting it from the set of available points:
(4.40)Ak={Xk,Tk+10,…,Tk+1ksto}.The computation ofXk+1 is performed in two phases.First phase: we determine a trial point (unperturbed ascent step and perturbation step).Second phase: we determineXk+1 by selection it from 𝒜k (dynamics step).As was shown in Theorem4.4, Substep (2.2) may use 𝒫ki=ξkZk, where 𝒵=(Z1,…,Zk+1) is a sample of N(0,1) and ξk is given by (4.33).For instance, we can consider elitistic dynamics:(4.41)Xk+1=arg⁡min⁡⁡X∈Akf(X).


## 5. Numerical Results

In order to apply the method presented in (4.41), we start at the initial value X̂0=X0∈C. At step k≥0, X̂k is known and X̂k+1 is determined.We generateksto the number of perturbation; the iterations are stopped when Xk is a Kuhn-Tucker point. We denote by kend the value of k when the iterations are stopped (it corresponds to the number of evaluations of the gradient of f). The optimal value and optimal point are fopt and Xopt, 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 ξk=â/log(k+2), where â>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 Figure1, 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:(5.1)Maximize12{(x2+x3-4x1)y1+(3x1-2x3+x5)y2+(3x1-2x2+x4)y3+(x3-2x1)y4+(x2-2x1)y5}-x1Subjectto‖A0-A1‖≤1:(x1-x2)2+(y1-y2)2+z1=1,‖A0-A2‖≤1:(-x1+x2-x5)2+(y1-y3+y5)2+z2=1,‖A0-A6‖≤1:(x1-x2+x4)2+(y1-y2+y4)2+z3=1,‖A0-A7‖≤1:(-x1+x3)2+(y1-y3)2+z4=1,‖A1-A2‖≤1:(2x1-x2-x3+x5)2+(y2-y3+y5)2+z5=1,‖A1-A3‖≤1:(2x1+x2)2+y22+z6=1,‖A1-A4‖≤1:(x1-x2)2+(y1-y2-1)2+z7=1,‖A1-A7‖≤1:(2x1-x2-x3)2+(-y2+y3)2+z8=1,‖A2-A3‖≤1:(x3-x5)2+(-y3+y5)2+z9=1,‖A2-A4‖≤1:(-x1+x3-x5)2+(y1-y3+y5-1)2+z10=1,‖A2-A5‖≤1:(2x1+x3-x5)2+(-y3+y5)2+z11=1,‖A2-A6‖=1:(2x1-x2-x3+x4+x5)2+(-y2+y3+y4-y5)2=1,‖A3-A6‖≤1:(-2x1+x2-x4)2+(y2-y4)2+z12=1,‖A4-A6‖≤1:(x1-x2+x4)2+(y1-y2+y4-1)2+z13=1,‖A4-A7‖≤1:(x1-x3)2+(1-y1+y3)2+z14=1,‖A5-A6‖≤1:(x2-x4)2+(y2-y4)2+z15=1,‖A5-A7‖≤1:(2x1-x3)2+y32+z16=1,‖A6-A1‖≤1:(2x1-x2-x3+x4)2+(-y2+y3+y4)2+z17=1,x2-x3-z18=0,xi2+yi2=1,i=1,2,3,4,5,x1+z19=0.5,xi+z18+i=1,i=2,3,4,5,0≤xi,yii=1,2,…,5,0≤zii=1,2,…,23.There are 33 variables and 34 constraints.We useâ=10d+6 and ksto=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-5.Figure 2
Regular octagon.The Fortran code of TPGRG furnishs the following optimal solutions,fopt≈0.72687:(5.2)Xopt=(0.25949,0.66664,0.66664,0.90703,0.90730),yiopt=1-xiopt,kend=626,CPUtime=4684s(1.5hours).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:R1, R2, R3, and R4. 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]).TheR1 reservoir receives two quality products 3 and 1 in quantities q0 and q1, respectively. The R2 reservoir contains a product of the third source, of quality 2 in quantity q2. One wants to get in the reservoirs R3 and R4 of capacities 10 and 20, respectively.Figure3, where the variable x1 to x2 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 (5.3)60q0+160q1+100q2-90(x1+x3)-150(x2-x4).
The quality of the mixture done contained in R1 is(5.4)x5=3q0+q1q0+q1∈[1,3].
The qualities of the final miscellanies are(5.5)x1x5+2x3x1+x3≤2.5,x2x5+2x4x2+x4≤1.5.
The addition of the constraints of volume conservation q0+q1=x1+x2 and q2=x3+x4 does permit the elimination of the variables q0,q1, and q2: (5.6)q0=(x5-1)(x1+x2)2,q1=(3-x5)(x1+x2)2,q2=x3+x4.Figure 3
A Problem of mixture.Remark 5.1.
One has(5.7)1≤x5≤1.5.
By the relation (5.5),
(5.8)x2x5≤1.5x2.
Suppose that x5>1.5, then x2x5>1.5x2, and so x2x5>1.5x2-0.5x4 that contradict with inequality (5.8).The mathematical model transformed is the following:(5.9)Maximize-120x1-60x2-10x3+50x4+50x1x5+50x2x5Subjectto2.5x1+0.5x3-x1x5-z1=0,1.5x2-0.5x4-x2x5-z2=0,x1+x3+z3=10,x2+x4+z4=20,x5+z5=1.5,x5-z6=1,0≤xii=1,…,5;0≤zii=1,…,6.
There are 11 variables and 10 constraints.We useâ=10d10, ksto=20000, and (5.10)X0=(2,9,0,8,1),f(X0)=170.
The Fortran code of TPGRG furnishes the following optimal solutions: cpu = 100 second, and kend=1390,(5.11)Xopt=(0,10,0,10,1),fopt=400,
since the mixture problem is known by this global solution.


## 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 Figure1, 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:(5.1)Maximize12{(x2+x3-4x1)y1+(3x1-2x3+x5)y2+(3x1-2x2+x4)y3+(x3-2x1)y4+(x2-2x1)y5}-x1Subjectto‖A0-A1‖≤1:(x1-x2)2+(y1-y2)2+z1=1,‖A0-A2‖≤1:(-x1+x2-x5)2+(y1-y3+y5)2+z2=1,‖A0-A6‖≤1:(x1-x2+x4)2+(y1-y2+y4)2+z3=1,‖A0-A7‖≤1:(-x1+x3)2+(y1-y3)2+z4=1,‖A1-A2‖≤1:(2x1-x2-x3+x5)2+(y2-y3+y5)2+z5=1,‖A1-A3‖≤1:(2x1+x2)2+y22+z6=1,‖A1-A4‖≤1:(x1-x2)2+(y1-y2-1)2+z7=1,‖A1-A7‖≤1:(2x1-x2-x3)2+(-y2+y3)2+z8=1,‖A2-A3‖≤1:(x3-x5)2+(-y3+y5)2+z9=1,‖A2-A4‖≤1:(-x1+x3-x5)2+(y1-y3+y5-1)2+z10=1,‖A2-A5‖≤1:(2x1+x3-x5)2+(-y3+y5)2+z11=1,‖A2-A6‖=1:(2x1-x2-x3+x4+x5)2+(-y2+y3+y4-y5)2=1,‖A3-A6‖≤1:(-2x1+x2-x4)2+(y2-y4)2+z12=1,‖A4-A6‖≤1:(x1-x2+x4)2+(y1-y2+y4-1)2+z13=1,‖A4-A7‖≤1:(x1-x3)2+(1-y1+y3)2+z14=1,‖A5-A6‖≤1:(x2-x4)2+(y2-y4)2+z15=1,‖A5-A7‖≤1:(2x1-x3)2+y32+z16=1,‖A6-A1‖≤1:(2x1-x2-x3+x4)2+(-y2+y3+y4)2+z17=1,x2-x3-z18=0,xi2+yi2=1,i=1,2,3,4,5,x1+z19=0.5,xi+z18+i=1,i=2,3,4,5,0≤xi,yii=1,2,…,5,0≤zii=1,2,…,23.There are 33 variables and 34 constraints.We useâ=10d+6 and ksto=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-5.Figure 2
Regular octagon.The Fortran code of TPGRG furnishs the following optimal solutions,fopt≈0.72687:(5.2)Xopt=(0.25949,0.66664,0.66664,0.90703,0.90730),yiopt=1-xiopt,kend=626,CPUtime=4684s(1.5hours).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:R1, R2, R3, and R4. 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]).TheR1 reservoir receives two quality products 3 and 1 in quantities q0 and q1, respectively. The R2 reservoir contains a product of the third source, of quality 2 in quantity q2. One wants to get in the reservoirs R3 and R4 of capacities 10 and 20, respectively.Figure3, where the variable x1 to x2 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 (5.3)60q0+160q1+100q2-90(x1+x3)-150(x2-x4).
The quality of the mixture done contained in R1 is(5.4)x5=3q0+q1q0+q1∈[1,3].
The qualities of the final miscellanies are(5.5)x1x5+2x3x1+x3≤2.5,x2x5+2x4x2+x4≤1.5.
The addition of the constraints of volume conservation q0+q1=x1+x2 and q2=x3+x4 does permit the elimination of the variables q0,q1, and q2: (5.6)q0=(x5-1)(x1+x2)2,q1=(3-x5)(x1+x2)2,q2=x3+x4.Figure 3
A Problem of mixture.Remark 5.1.
One has(5.7)1≤x5≤1.5.
By the relation (5.5),
(5.8)x2x5≤1.5x2.
Suppose that x5>1.5, then x2x5>1.5x2, and so x2x5>1.5x2-0.5x4 that contradict with inequality (5.8).The mathematical model transformed is the following:(5.9)Maximize-120x1-60x2-10x3+50x4+50x1x5+50x2x5Subjectto2.5x1+0.5x3-x1x5-z1=0,1.5x2-0.5x4-x2x5-z2=0,x1+x3+z3=10,x2+x4+z4=20,x5+z5=1.5,x5-z6=1,0≤xii=1,…,5;0≤zii=1,…,6.
There are 11 variables and 10 constraints.We useâ=10d10, ksto=20000, and (5.10)X0=(2,9,0,8,1),f(X0)=170.
The Fortran code of TPGRG furnishes the following optimal solutions: cpu = 100 second, and kend=1390,(5.11)Xopt=(0,10,0,10,1),fopt=400,
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 parametersâ and ksto.



---

*Source: 976529-2010-12-08.xml*