an algorithm for linearly constrained programs with a ... - Science Direct

0 downloads 0 Views 782KB Size Report
We present an algorithm for solving the following optimization problem: Find a (local) minimum of the objective ..... P.E. Gill, W. Murray and M.M. Wright, Practical.
0895-7177192 $5.00 + 0.00 Copyright@ 1992 Pergamon Press plc

Mathl. bornput. Modelling Vol. 16, No. 4, pp. 41-50, 1992 Printed in Great Britain. All rights reserved

AN ALGORITHM FOR LINEARLY CONSTRAINED PROGRAMS WITH A PARTLY LINEAR OBJECTIVE FUNCTION H. GREINER Philips Forschungslabor GmbH Weisshausstra&, 51 Aachen, Germany (Received

January

1991)

Abstract-A Newton-type algorithm for finding a minimum of a function, subject to linear constraints, which is linear in some of its arguments is presented. It employs an active set strategy especially adapted to this hind of problem. Compared with other methods, it is particularly efficient for problems in which the number of nonlinear variables greatly exceeds the number of linear constraints, and for which the Hessian matrix of the nonlinear part can be easily inverted. This is exemplified for the calculation of complex chemical equilibria involving an ideal gas phase with many species and pure condensed phases.

1. INTRODUCTION We present an algorithm for solving the following optimization problem: Find a (local) minimum of the objective function G(x, y) = F(x) + L(y) = MIN, (la) where x = (~i,...,z,) E 54” are “nonlinear” variables and y = (~1, . . . , yr) E Rr are “linear” variables. F is assumed to be a twice continuously differentiable function defined on I?‘, which grows like IX]” for x --+ co with Q > 1, and to have a positive definite Hessian. The linear part of the objective function is defined by

L(Y) =

2

Pi

(lb)

Yi,

i=l

with real coefficients p = (~1, . . . , ,up) E W. The problem variables (x, y) are subject to the following constraints: The linear variables y are required to be non-negative j= l,...,p, Yj L 09 (Ic) and the variables (x, y) are furthermore restricted by the linear equalities Ax+By=f,

(14

where A=

(al,...

B=

(bl,...,bp)

,a,) is a real m x n matrix with column vectors aj E R”, is a real m x p matrix with column vectors bj E R”,

and

f is a column vector from R”. Linearly constrained problems having the form of Problem (1) arise in many applications [l]. They can be computed by Newton-type methods as described in [l-3]. Depending on how the approximating quadratic programs are solved, these algorithms are qualified as “null-space” or “range-space” methods [2,3]. If one considers large-scale problems, null-space methods are Typeset by AA&TEX 41

42

H. GREINER

particularly suited for problems where the number of nonlinear variables n is about equal to the number of linear constraints m. In the case that n is much greater than m, the application of a range-space method is recommended. This is explained in more detail in Section 4. Computational problems may arise from the indefiniteness of the Hessian with respect to the linear variables y (in [4, Chapter 9.21, this is discussed for the chemical equilibrium problem). Also, for problems with a large number of nonlinear variables and few linear constraints, nullspace methods are impractical as discussed in [2, Chapter 5.6.21. Our algorithm uses an active set strategy [2,3] for the constraints yj 1 0, which allows one to eliminate the free linear variables from the subproblems defined by the current active set. In this way, problems with an indefinite Hessian can be avoided. Furthermore, in problems where the application of a range-space method is indicated, the solution of the subproblems can be considerably simplified. We also present a version of the algorithm for the case when the objective function is convex. Its underlying idea is due to Zangwill [5]. Finally, we discuss the application of the algorithm to the calculation of complex chemical equilibria. 2. PRELIMINARIES In this section, we give some definitions and facts from optimization theory concerning the solution of Problem (1). A necessary condition for a feasible point (x, y) to be an optimal point for Problem (1) is the existence of a row vector of multipliers

(v,w) satisfying the Kuhn-Tucker tion, we recommend [S]):

E

IFP x RP,

conditions (as a general reference on linear and non-linear optimiza-

pi + vbj - wj = 0 YjWj =

3

0,

Wj

10,

forj=

l,...,p,

(2b)

forj=

l,...,p.

(2c)

= (WI,... , wp) E OWP are the multipliers corresponding to the inequality constraints (lc), v = (Wl,...,Vm) E OWm are the multipliers corresponding to the linear equality constraints (Id). To solve Problem (l), one has to identify the active and inactive inequality constraints at the solution. An inequality constraint yj > 0 is called active, or binding, at (x, y) if the equality Otherwise, it is called inactive. An active set strategy is a method for % = 0 is satisfied. finding the correct active set [2]. Starting with an assumed set of active inequality constraints, a “subproblem” of Problem (1) ( in which the active constraints are treated as equality constraints) is solved. If in the course of the optimization, one of the inactive constraints becomes binding, it is added to the active set and the minimization is restarted. If an optimal solution of a subproblem has been identified, it is tested to see whether the deletion of any constraints from the active set leads to a lower minimum, and the algorithm is continued with an accordingly modified active set. It terminates when a point satisfying the Kuhn-Tucker conditions has been found. A thorough discussion and comparison of various active set strategies is given in [7]. A fast and reliable algorithm for the solution of Problem (1) consequently consists of w

(1) an efficient method for solving the subproblems defined by the various active sets; (2) an efficient active set strategy, which defines the new active set after the solution of each subproblem and guarantees that the optimal solution will be found after solving a finite number of subproblems, i.e., that the same active set cannot recur. Clearly, at an optimal solution (x0, y”) of Problem (l), the vectors bj, with yj > 0, are linearly independent (consider the linear program L(y) = MIN, with constraints B y = f - A x0

43

Linearly constrained programs

and y 2 0). Therefore, the subproblems are defined in our case by subsets D of (1,. . . ,p}, such that the vectors D = {bj 1j E D} are linearly independent, and are formulated as follows. SUBPROBLEM

(D).

(x, y) E R” x R,

G(x, y) = F(x) + L(y) = MIN, subject to Ax+By=f,

for j E (1,. . . ,p} - D.

Yj = 0,

Provided that Subproblem (D) is feasible, it can be solved by the method described in Section 4. The solution (x, y, v) has to fulfill the Kuhn-Tucker conditions dF(x) a21,+vak pj +vbj Imposing the additional constraints will refer in the sequel: SUBPROBLEM

fork=l,...,n,

(3a)

= 0,

for j E D.

(3b)

yj > 0 for j E D defines a related subproblem to which we

(D+).

G(x, subject

= 0,

y) = F(x) + L(y) = MIN,

(x, y) E Rfl x HP,

to Ax+By=f,

The Kuhn-Tucker

Yj=O,

fOrjE{l,...,p}-D,

yj>O,

fOrjED.

conditions for this subproblem are: Wx) z+vak pj+vbj-wj Yjwj -- 0,

Wj

=O,

for k=l,...,n,

(4a)

~0,

for j E D,

(4b)

for j E D.

(4c)

>

0,

The conditions imposed on F ensures that the Subproblems (D) and (D+) are bounded below, and possess a solution. To establish the convergence of the algorithms we are going to formulate in the next section, we need the following regularity condition on Problem (1). PROPERTY. If (x,y,v) fulfills the optimality the multipliers (u, v) are uniquely determined.

UNIQUENESS

conditions (3) for Subproblem

(D),

This condition is not very restrictive and is clearly fulfilled if the linear constraint matrix (A, B) has maximal rank. 3. THE

ALGORITHM

In this section, we formulate our algorithm for Problem (1). We present two versions: a general one and one which applies only if the function F in (1) is convex. But first we introduce some notation. Notation For y E BP, define P(Y) = {j I Yj > 01, and the corresponding set of vectors P(Y) = {bj

Ij

E P(Y)}.

IDI designates the number of elements in a finite set D.

H. GREINER

44 Statement

of the Algorithm

(General

Case)

Initialization

Choose a feasible point (x0, y”) for Problem (l), such that D1 = P(y”) is linearly independent, as a starting point for the algorithm. Set k = 1 and go to (I). (I) Solve the Subproblem

(D$)

Apply a Newton type algorithm for the solution of Subproblem (Dt), starting point (see Section 4). Consider the following possibilities:

with (x’-l,

y”-I)

as a

(Ia) For all iterates, we have yj > 0 and the algorithm terminates with a solution (Xk,jrk,v”) with yk > 0 for Subproblem (Dk) . Define Gk = G(TZk,yk) and proceed to (II). (Ib) One of the constraints yj 2 0, j E D k, becomes binding during the iteration. Define (x’, yk) as that iterate, G” as the corresponding function value and Dk+’ = D” - {j}. Set k = k + 1 and return to (I). (II) Test Optimality If the optimality conditions Wj

=

pj

+

v”bj

2 0,

j=

l,...,p,

are fulfilled, (rrl”,yk) is an optimal solution for Problem (1) and the algorithm stops. Otherwise, proceed to (III). (III) Modify Active Set If the optimality conditions are not fulfilled, choose a j(k) E (1,. . . ,p} satisfying pj(k)

+

Vkbj(k) = min{&

+ v’bj

1j = 1,. . . ,p} < 0.

Construct a new set Dk+l according to the following two possibilities: (IIIa) If bj(k) is linearly dependent on the set P(yk) write bj(k) =

C

cj bj.

jWP) Now, let t,,

since $ > 0 for j E P(p). stops. Otherwise, define

= max{t

If t,,,

Y:

I$

- t cj 2 0,

> 0,

is unbounded, Problem (1) is unbounded and the algorithm

=$'ttmmCj,

~7;~~) = t,,

for all j E P@‘)}

and yj” = 0,

forj E P(P), for all other j.

Define xk = T”, Dk+l = P(yk), set k = k + 1 and proceed to (I). (IIIb) If bj(k)k not linearly dependent on the set P(yk), define (xk,yk) P(y”) Uj(k),k = k + 1, and go to (I).

PROPOSITION. The algorithm converges in a finite number holds for Problem (1).

= (Tk,yk),

D”+l

=

of steps, if the uniqueness condition

PROOF. At each iteration, the strict inequality G”+l < Gk holds (for Case (Ia), this follows from Lemma Al in the Appendix, and for Case (Ib) from the fact that G decreases with each iterate for Subproblem (D:)). Let k(i), i = 1,2,. . . , denote the indices for which Case (Ia) occurs. If

45

Linearly constrained programs k(i) < II < 12 < lc(i + l), D’z must be strictly contained in D’l. occur, at most, m times in succession. From

Gk(“+l)

= min

~k(“+l)

< G’(j)

Consequently,

Case (Ib) can

= min D’E(‘),

we have that the linearly independent sets D”ci), i = 1,2,. . . , are distinct. terminate after a finite number of steps.

The iteration must I

If the function F is convex, the following variant of the algorithm can be established (see The only difference with the formulation of the general also [5,6], for a related formulation). algorithm arises in Step (I). Instead of the Subproblems (D$), the Subproblems (D”) are solved without regard to the constraints yj, j E D k. If the solution should turn out to be infeasible, feasibility is restored by taking the intersection of the line segment defined by the solution, and the previous iterate with the boundary of the feasible set. Algorithm

for Convex F

(I) Solve Subproblem

(D”)

Obtain a solution (Kk, yk, v”) for the Subproblem (D”). Distinguish the following cases: (Ia) If yk >_ 0, set Gk = G(zk,yk) (Ib) Otherwise, define

and go to (II).

s=max{tEIO,l]~(l-t)$l+t$~O,

j=l,...,

p},

(Xk,yk)=(l-S)(Xk-l,yk-l)+S(Xk,~k), Gk =

G(x",y”),

Dk+’ = P(yk),

set k = k + 1 and go to (I). The initialization

and Steps (II) and (III) are as in the general case.

The algorithm for convex F converges in a finite number of steps, if Problem PROPOSITION. fulfills the uniqueness condition.

(1)

The proof is similar to the one given in [5,6]. Depending on whether Case (Ia) or (Ib) occurs at the kth iteration, one obtains the strict inequality Gk+’ < Gk, or the inequality G”+l < G” (Lemma A2 and A3 in the Appendix). The rest of the argument follows the proof for thegeneral case given above. I PROOF.

Although, in general, one would apply the general algorithm even to problems with convex F (as the optimization procedure stays in the feasible set), there are special problems for which the algorithm for convex F may present an advantage. This may be the case if the non-linear equations (3), characterizing an optimal point for Subproblems (D”) can be solved directly in an efficient way (see [8, Chapter 5.41). A n example is the method described by Schnedler [9] for the calculation of complex chemical equilibria we will discuss in Section 6. We finally remark that, at the expense of a little more notational clutter, a version of the algorithm with prescribed lower and upper bounds on the nonlinear and linear variables can be formulated. 4. NUMERICAL

IMPLEMENTATION

We describe numerical procedures for the implementation

of the algorithm.

Initialization Starting with an arbitrary feasible point (x, y), a feasible point (x0, y’) such that P(y’) linearly independent can be obtained as the solution of the linear program L(y) = MIN,

subjecttoBy=f-Ax,

Y 20.

is

46

H. GREINER of the Subproblems

(D) and (D+)

To solve the Subproblems from the linear constraints

(Dk) or (D:),

Solution

we first eliminate

Ax+C

the free linear

variables

yj, j E D

bjgj =f jE:O

(the columns bi corresponding to these variables problems are transformed to the form

E(x)

are linearly

= F(X) + e

independent!).

As a result,

ei Zi = MIN,

the

(5a)

j=l

subject

to the constraints

Ax=&

(5b)

Axxo.

(5c)

Here, the second term of the objective function results from the elimination of the free variables. d is an (m - ]D],n)-matrix, J% a (jD],n)-matrix and $ is an (m - ID])-vector. The transformed Problems (5) are linearly constrained optimization problems of the form H(x) with a positive-definite Hessian An iteration of a Newton-type (i) For a given feasible

= MIN, matrix. method

subject

to Rx

= b,

for their solution

point x k, find a descent

direction

Cx > 0,

consists

steps

dk by solving the quadratic

Q(d) = g”d + f dTGkd = MIN, subject

of the following

[2,3,8]: program

(f-w

to Rd=O.

(6b)

Here, the notation gk = VH(xk),

Gk = ‘0%(x”)

is used. (ii) Determine cYkax = max{a

I C(x” + adk) 10).

(iii)

Find an ok, 0 < ok 5 crmax, such that xk+r = xk + okdk and H(xk+‘) < H(x’) by an appropriate line search algorithm. (iv) If the current iterate is optimal, or if an inequality constraint becomes binding (for Subproblem (D+)), stop. Otherwise, set L = k + 1 and resume at (i).

The necessary

optimality

conditions

for Problem

(5) are

gk + Gkd = RTw, Rd=O.

and

(7a) (7b)

Here w E fpm designates the Lagrangian multipliers of the problem. Obviously, the main computational effort of our algorithm is spent on the quadratic subproblems (6). Their solution is determined from Equations (7) according to one of the following methods [2,3].

Linearly constrained programs

Null-Space

47

Methods

Find an (n, n - m)-matrix Z, whose columns span the null-space of R. As d belongs to the null-space of R, it has the form Z v. Introducing d = Zv into Equation (7a), and multiplying with ZT from the left gives a linear system of equations with (n - m) unknowns: (ZTG@)Z)

v = -ZTgk.

(8)

(N.B. ZTWT = (RZ)T = O!) From the solution v, one obtains dk = Zv. Range-Space

Methods

Multiplying Equations

(7a) with Gck)-l

from the left, we obtain

dk = dk-(RTw Multiplying

(9)

(9) with R from the left leads to a linear system of equations with m unknowns

(RG(k)-l~T) Substituting

-g”).

w

R

=

($k)-lgk.

the solution w into Equations (9), then defines dk. 5. COMPUTATIONAL

EFFICIENCY

Obviously, the main computational effort for Problem (1) is invested into the solution of the quadratic subproblems (6). For n large compared to m, null-space methods become impractical as they involve the solution of a dense system of linear equations with (n - m) unknowns. In this case, range-space methods are to be preferred (provided the G” have a sparse structure, such that the linear equations with coefficient matrix G” can be easily solved). They require only the solution of a system of linear equations with m unknowns. It is a particular feature of our algorithm that it reduces the dimension of this sytem of linear equations even further by the number of the free linear variables. In this way, the computational effort for solving the Subproblems (Dk) and (D$) can be diminished considerably. 6. APPLICATION: THE CALCULATION COMPLEX CHEMICAL EQUILIBRIA

OF

In the following, we want to illustrate the application of the algorithm to the so-called “gas solids” problem. We consider a closed chemical reaction system, consisting of an ideal gas phase and pure condensed species at fixed temperature T and volume V (the case of constant temperature and pressure is a little bit more complicated, but can be treated along similar lines). These systems are important in many applications. The system’s equilibrium is described by the minimum of its total free energy, subject to mass balance constraints. We assume that the volume contains k= l,... , m different elements with mole numbers f = (fi , . . . , fm). The elements form various gaseous and pure condensed species. The mole numbers of the gaseous species will be denotedbyzi, i= l,..., n, and the mole numbers of the pure condensed species by yj, j = 1,. . . ,p. For a gaseous species i, the kth component of the vector A’ gives the number of moles of element k contained in one mole of that species. Similar notation applies to the pure condensed species with vectors Bj. The chemical potentials of the species at the given temperature are denoted by ,$ and ~5. The free energy function of the gas phase is given by (energy is expressed in units of RT): Fgas = g

xi (d

+ ln(zi) + In V - 1).

(19)

i=l Chemical equilibrium is then defined by the solution of the following minimization

G(x, Y) = Fgas(X) +

2 P: j=l

MCM 16:4-D

Yj = MIN,

problem

(114

H. GREINER

48

subject

to

2 A”Zi+f: Bjyj

i=l

WI

xLOandy20.

=f,

j=l

In dealing with complex systems with many pure condensed phases, one of the difficulties is to determine the condensed phases actually present at equilibrium, out of a multitude of possible phase combinations, in an effective way. Here we show that our algorithm can be directly applied to this problem (for a different approach to the phase finding problem, see [lo]). It is not difficult to verify that the function Fgas defined on the open set R; is differentiable and convex, and that it can be extended continuously to @+. As the feasible set for Problem (11) is bounded (the linear constraint matrix (A, B) has only non-negative entries!), the corresponding Subproblems (D) and (D+) possess solutions with unique nonvanishing mole numbers zi > 0, i= 1 ,..., n, for the gaseous species. (This is due to the logarithmic terms in Fgas.) Therefore, Problem (11) can be handled by our algorithm for the convex case. The Subproblems (D) become, after the elimination of the free linear variables (see Equations (5)) = Fgas(x) + g

Fsas(x)

ei xi = MIN,

(12a)

i=l

subject

to Ax==.

(I2b)

In the sequel, we shall write f(z) for the vector with components f(~i) where f is a real valued function, and z = (tr , . . . , zk) a number tuple. D(z) will designate the diagonal matrix with entries zi . Finally, c stands for the vector with components

With

these designations,

ci = pi(T)

+ ei + In(V),

the necessary

optimality

In(x)

(i = 1,. . . ,n). conditions

= dTw

for Problem

(12) can be written

- c, and

as

(I3a)

Ax=I.

(13b)

To solve Equations (13), Schnedler [9] introduces the logarithmic variables p = In(x) and applies Newton’s method to the transformed equations. By changing to logarithmic variables, the numerical problems associated with exceedingly small mole numbers are avoided. Linearizing about a given point (pa, w”) leads to the linear equations (we write x0 = exp(p’)) : Ap=~T~o+~TA~-po-c,

and

(14a)

iiD(~~)Ap=~-ii.~~. Equations (14) are solved by substituting m - IDI linear equations (A D(x’) AT) Aw = i - ix0

(14b) (14a) into (14b) and solving

the resulting

system

of

- k D(x’) iTwo + i D(x’) (p” + G),

for Aw. Inserting Aw into (14a) then gives Ap. In the literature on the computation of chemical equilibria algorithms [4,11,12], they are classified as stoichiometric or nonstoichiometric, depending on how the mass-balance constraints are handled. Examples for either class are discussed in [4,11]. In our terminology, stoichiometric algorithms correspond to null-space and nonstoichiometric ones to range-space methods. According to the literature, stoichiometric algorithms are more efficient than the nonstoichiometric ones for the gas-solids problem [11,13]. Nonstoichiometric algorithms encounter problems of convergence (see [4, Chapter 9.21 and [ll]) w h’ic h are related to singular coefficient matrices (A, B) (V’G)-I@,

B)T,

49

Linearly constrained programs

caused by an indefinite Hessian and the attendant difficulty to find the correct pure condensed phases at equilibrium. Our range-space algorithm avoids these problems by eliminating the pure condensed species from the calculation and by using a consistent strategy for finding the correct pure condensed species. Numerical singularities of the matrix ii

(V2FgJ1iiT

can be escaped by choosing the most abundant species as the components (basis species) of the system (for more details on this point, see [9]). Thus, we have demonstrated that an efficient and reliable nonstoichiometric algorithm for the gas-solids problem can be formulated (cf. also [14]). 7. CONCLUSIONS We have presented an algorithm for solving linearly constrained programs whose objective function is partly linear. Compared to other established methods, the algorithm should be particularly efficient for problems in which the number of nonlinear variables is much larger than the number of linear constraints (provided the Hessian matrix of the nonlinear part can be easily inverted). This is confirmed by practical experience with an implementation of the algorithm for the calculation of complex chemical equilibria. REFERENCES 1. B.A. Murtagh and M.A. Saunders, Largescale linearly constrained optimization, Mathematical Programming 14, 41-72 (1978). 2. P.E. Gill, W. Murray and M.M. Wright, Practical Optimization, Academic Press, San Diego, CA, (1981). 3. P.E. Gill and W. Murray, Editors, Numerical Methods for Constnzined Optimization, Academic Press, San Diego, CA, (1974). Equilibrium Analysis: Theory and Applications, John 4. W.R. Smith and R.W. Missen, Chemical Reaction Wiley, New York, NY, (1982). 5. W.I. Zangwill, A decomposable nonlinear programming approach, Operations Research XV (6), 1068-1087 (1967). 6. W.I. Zangwill, Nonlinear Progmmming, A Unified Approach, Prentice Hall, Englewood Cliffs, NJ, (1969). 7. M.L. Lenard, A computational study of active set strategies in nonlinear programming with linear constraints, Mathematical Programming 16, 81-97 (1979). 8. E. Minoux, Mathematical Programming: Theory and Algokthms (Chapter 5.4), John Wiley, New York, NY, (1986). 9. E. Schnedler, The calculation of complex chemical equilibria, CALPHAD 8 (3), 265-279 (1984). 10. H. Greiner, Computing complex chemical equilibria by generalized linear programming, Mathematical and Computer Modelling 10 (7), 529-550 (1988). 11. W.R. Smith, The computation of chemical equilibria in complex systems, Ind. Eng. Chem. Fundam. 19, l-10 (1980). 12. W.R. Smith, Computational aspects of chemical equilibrium in complex systems, In Theoretical Chemistry: Advances and Perspectives, Vol. 5 (Edited by H. Eyring and D. Henderson), Academic Press, San Diego, CA, (1980). 13. W.R. Smith and R.W. Missen, Calculating complex chemical equilibria by an improved reaction-adjustment method, Canadian Journal of Chemical Engineering 46, 269-272 (1968). 14. H. Greiner, An efficient implementation of Newton’s method for complex non-ideal chemical equilibria, Computers Chem. Engrg. 15 (2), 115-123 (1991).

APPENDIX We prove the lemmas needed to establish the convergence of the algorithms described in Section 3. LEMMA Al.

Consider

the

algorithm

the general case. If Case

for

(III)

occurs

at the kth iteration,

the strict

inequality Gk+l

< Gk

holds.

PROOF. To prove the assertion, we treat the Cases (IIIa) and (IIIb) separately. Case (I&): Assume that Dk+’ has been constructed according to the prescription given in (I&).

by Yj(t)=$-ttj yj(k)(t)

=

t ad

for ?/j(t)

=

0,

j E p(Vk),

for all other j.

Define y(t)

E W’

H. GREINER

50 Clearly, yk = y(tmaw). Using the identity

one obtains L(Y(t)) Hence, we have L(yk) that the point Case (IIIb): was

= L(Y(0))

= L(y (tmax )) < L(jrk).

+

t(Vkbj(k)

+

(x”, yk) is feasible for Subproblem (D$+l).

We prove that (xk, yk) camrot be optimal for Subproblem (Dttl)

optimal for Subproblem (D:+‘),

the optima&y

one could llnd multipliers (ck,wk),

conditions (4) for Subproblem (D:+‘).

But, asyJ > Oforj E 03, (xk,yk,ck) by the uniqueness condition. Hence,

the algorithm for the convex

Gktl

(xk,

yk)

would fulflh

would meet the optimality conditions (3) for Subproblem (Dk) and ckbj(kl

Consider

by redzlctio ab absurdurn: if

such that (xk,yk,jk,Gk)

In particular, one would have

+

pJ(k)


0 must hold, or equivalently, jr;&; > 0. But -‘+* yJck) _< 0 would contradict the optimality of (xk,yk) for the problem G(x, Y) = Mm,

subject to the constraints (Id) and $!j(k) < 0.

(It is easily checked that (x”, yk,vk) fulftllsthe Kuhn-Tuclcer conditions for this problem. This and the convexity of G imply that (xk,yk) must be optimal for this problem. Because of the inequality G(Zktl ,yktl) < G(xk,Yk) I proved above, (Zktl, yk*+l ) cannot be optimal for the problem.)