On the compatibility of classical multiplier estimates with variable reduction techniques when there are nonlinear inequality constraints

Q¨ UESTII ´ O,vol. 23, 1, p. 61-83, 1999 ON THE COMPATIBILITY OF CLASSICAL MULTIPLIER ESTIMATES WITH VARIABLE REDUCTION TECHNIQUES WHEN THERE ARE NONLINEAR INEQUALITY CONSTRAINTS E. MIJANGOS Zientzi Fakultatea (EHU)∗ N. NABONA Universitat Polit` ecnica de Catalunya∗∗ The minimization of a nonlinear function subject to linear and nonlinear equality constraints and simple bounds can be performed through minimizing a partial augmented Lagrangian function subject only to linear constraints and simple bounds by variable reduction techniques. The first-order procedure for estimating the multiplier of the nonlinear equality constraints through the Kuhn-Tucker conditions is analyzed and compared to that of Hestenes-Powell. There is a method which identifies those major iterations where the procedure based on the Kuhn-Tucker conditions can be safely used and also computes these estimates. This work justifies the extension of the former results to the case of general inequality constraints. To this end two procedures that convert inequalities into equalities are considered. Keywords: Nonlinear programming, general inequality constraints, augmentedlagrangian,variable reduction, Lagrange multiplier estimates. * Department of Applied Mathematic, Statistics and Operations Research, Zientzi Fakultatea (EHU), P.O. Box 644, 48080 Bilbao (Espanya). Corresponding author e-mail: [email protected]. **Department of Statistics and Operations Research, Universitat Polit` ecnica de Catalunya, Pau Gargallo, 5. 08028 Barcelona (Espanya). – Received February 1998. – Accepted October 1998. 61 1. INTRODUCTION AND MOTIVATION Consider the linearly and nonlinearly constrained problem (1) (2) (3) (4) minimize f(x) subject to: Ax =b c(x) = 0 l⩽x⩽u, (EP) where (1) f:RnR.f(x)is nonlinear and twice continuously differentiable on the feasible set defined by constraints (2–4). (2) Ais an m×nmatrix and ban m-vector. (3) c:RnRr,is such that c= [c1,··· , cr]t,ci(x)being linear or nonlinear and twice continuously differentiable on the feasible set defined by constraints (2) and (4) ∀i= 1,··· , r. (4) n≫m+r. To solve this problem one could use, among others, partial augmented Lagrangian techniques [1, 2, 3] as in [9, 11, 6], where only the general constraints (3) are included in the Lagrangian. In the application of these techniques there are two fundamental steps. The first solving (5) (6) (7) minimize xLρ(x, µ) subject to: Ax =b l⩽x⩽u, (ES) where ρ > 0and µare fixed, Lρ(x, µ) = f(x) + µtc(x) + 1 2ρc(x)tc(x). Should the solution exobtained be infeasible with respect to (3), the second step, which is the updating of the estimate µof the Lagrange multipliers of constraints (3), is carried out, also updating, if necessary, the penalty coefficient ρ, then going back to the first step. Should exbe feasible (or the violation of constraints (3) be sufficiently small) 62 the procedure ends. It is of paramount importance that the multipliers estimate µbe as accurate as possible, otherwise the convergence of the algorithm can be severely handicapped, as shown in [1, 2, 3]. In practice there are two first-order procedures to estimate µ. On the one hand the method put forward by Hestenes [8] and Powell [15] eµ=µ+ρc(ex), and on the other hand µLobtained through the classical solution to the system of Kuhn- Tucker necessary conditions (8) ∇f(ex) + ∇c(ex)µ+Atπ+λ= 0 by least squares, as suggested in [7]. However, whereas the first procedure can be always used for any ex, without hampering the convergence [2, 3], this is not the case with the second procedure, as system (8) is only known to be compatible at the optimizer x∗, not being necessarily so at ex, thus possibly giving rise to bad estimates µL, shown up by large residuals for system (8). Section 2 of this work presents a study of the viability of using this multiplier estimation technique within the minimization of a Partial Augmented Lagrangian subject to linear constraints and bounds by the Murtagh and Saunders procedure [13] for problem EP. (An alternative development of the contents of Section 2 can be found in [10].) Sections 3 and 4 consider two ways of extending these results to problem (9) (10) (11) (12) minimize f(x) subject to: Ax =b c≤c(x)≤c l⩽x⩽u, (IP) where cj<cj, for j= 1,...,r. In Section 3, a vector of slacks «y»is used to convert constraints (11) into equalities, and, in Section 4, slacks «y»and artificial variables «w» are used with the same aim. Section 5 contains the conclusions. 2. ANALYSIS OF COMPATIBILITY The compatibility of the multiplier estimate obtained through the classical solution to the system of Kuhn-Tucker necessary conditions with the variable reduction techniques and its relationship with that of Hestenes and Powell are analyzed along this section. 63 Let us consider the first-order conditions associated with a local optimizer x∗of the problem EP, which are ∇f(x∗) + ∇c(x∗)µ∗+Atπ∗+λ∗= 0 Ax∗=b c(x∗) = 0 λ∗ i≤0,if x∗ i=li λ∗ i≥0,if x∗ i=ui λ∗ i= 0,otherwise so that unique vectors µ∗, π∗and λ∗exist, such that the first equation holds. These vectors are denoted Lagrange multipliers; being ∇c(x) = [∇c1(x),· · · ,∇cr(x)]. (The gradient is considered to be a column vector). Throughout this work we assume AS1. x∗is a regular point — i.e., the Jacobian of the active constraints at x∗has full rank. Solving problem EP through a partial augmented Lagrangian techniques consists basically of the following algorithm, where for given ρ > 0and µthe subproblem ES (5–7) is successively solved. Algorithm 2.1. 1. For an initial point x0(not necessarily feasible with respect to c(x) = 0), a given scalar ρ > 0and vector µ, solve subproblem ES and obtain its optimizer ex= x(µ, ρ). 2. Should this exmake c(ex)to be zero or nearly so for a prespecified tolerance, x∗=ex and the problem is solved, otherwise, 3. µis updated by estimating µ∗, for which two possibilities are considered: –UKT . Solving system (13) ∇f(ex) + ∇c(ex)µ+Atπ+λ= 0, (which could be solved through least squares using QR factorization as justified in [7], obtaining eµ, and also eπ). 64 –UHP . Setting (14) eµ=µ+ρc(ex), as established in [2, 3]. 4. Should the solution of ES not reduce kc(x)ksufficiently, ρwould be updated as ρ=νρ, where ν > 1. 5. Make µ=eµand x0=ex, and return to 1. The issue is now the compatibility of system (13) at ex, and in this event the relationship between the procedures UKT and UHP to estimate vector µ∗at that point. Each time subproblem ES is solved exactly by Murtagh and Saunders’s active set method [13], we obtain an optimizer exand an associated partition of matrix A= [BA|SA|NA], and, thus, the variable reduction matrix ZAshown below: (15) ZA=   −B−1 ASA 1l 00   ,which satisfies AZA= 0. Since exis an optimizer of problem ES the necessary first-order optimality conditions must hold; i.e., there exist unique vectors eπand e λsuch that: ∇xLρ(ex, µ) + Ateπ+e λ= 0(16) Aex=b,(17) exi=li, i =t+ 1,··· ,t(18) exi=ui, i =t+ 1,··· , n(19) where tis the number of basic and superbasic variables, and e λi≤0,if exi=li e λi≥0,if exi=ui e λi= 0,otherwise. Expression (16) is equivalent to (20) Ateπ+∇c(ex)[µ+ρc(ex)] + e λ=−∇f(ex), 65 which in matrix form yields    Bt A∇BAc(ex) 00 St A∇SAc(ex) 00 Nt A∇NAc(ex) 1l     eπ µ+ρc(ex) e λ   =−   ∇BAf(ex) ∇SAf(ex) ∇NAf(ex)   , where ∇BAc(ex)stands for the rows of ∇c(ex)associated with the rows of Bt A. Similarly ∇SAc(ex)and ∇NAc(ex), and also the partition of ∇f(ex), are defined. Let us assume that matrix A(ex) = m z}| { s z }| { BASANA}m ∇BAc(ex)t∇SAc(ex)t∇NAc(ex)t}r | {z } n has full row rank. It is now possible to get from A(ex)a full rank basic matrix Bsuch that it contains matrix BA, which was obtained at the end of the first step of Algorithm 2.1. Once the basic matrix Bhas been defined, matrix Scan be established such that [B|S]contains BAand SAas submatrices, thus: (21) B= m+r z}| { BAS′ Ae NA}m ∇BAc(ex)t∇S′ Ac(ex)t∇e NAc(ex)t}r | {z } m| {z } s′|{z } m1 and S= s−s′ z}| { S′′ A ∇S′′ Ac(ex)t | {z } s′′ , with SA= [S′ A|S′′ A]. The rest of columns from A(ex)makes up submatrix N. Let ∇Bv(x)denote the gradient with respect to the variables associated with Bof any differentiable function v(x). Similarly ∇Sv(x)and ∇Nv(x). See [10] for an efficient procedure to build up Bfrom data available at exwhen subproblem ES has been solved. From this partition [B S N]of A(ex)we have the variable reduction matrix Z(ex) =    −B−1S 1l 00   ,which satisfies A(ex)Z(ex) = 0. Premultiplying by Z(ex)both sides of (20) we get the equivalent expression (22) Z(ex)tAteπ+Z(ex)t∇c(ex)[µ+ρc(ex)] + Z(ex)te λ=−Z(ex)t∇f(ex). 66 According to the definition of Z(ex)the following must hold: (23) Z(ex)t∇c(ex) = 0 and Z(ex)tAt= 0. Furthermore, (24) Z(ex)te λ=£−(B−1S)t1l 00¤  e λB e λS e λN  =−(B−1S)te λB+e λS. As we can see, Z(ex)t∇f(ex)will vanish only if the following condition holds (25) e λS= (B−1S)te λB, which in general is not satisfied, as can be easily proved; see [11]. Let σk= [σ1k,··· , σmk]t(with m=m+r) be the column of B−1Sassociated with the superbasic variable xk; then (25) can be recast as e λS k= m X j=1 σjke λB j, k ∈ S, Sbeing the set of indices associated with the columns of S. The basic equivalent path βkof superbasic variable xkis defined as the set of basic variables xlthat have a nonzero entry σlk in the column of B−1Scorresponding to variable xk. Taking into account the expressions (22)-(24) we get (26) Z(ex)t∇f(ex) = α, where αis a vector (whose dimension is the number of columns of S) such that (27) αk=X i∈βk∩e N σike λB i−e λS k, e Nbeing the index set of the columns of matrix (28) e NA ∇e Nc(ex)t selected to make up the basis matrix of A(ex)(see expression (21)) and σik, as previously defined, entry (i, k)of matrix B−1S. Furthermore, since by construction the 67 set of indices associated with the columns of Sis a subset of the set of indices associated with the columns of SA,e λk= 0 holds for all kcorresponding to a column of S, hence (27) becomes (29) αk=X i∈βk∩e N σike λB i. It must be pointed out that vector αturns out to be the nonzero part of the residual vector corresponding to system (13), when, once the partition of matrix A(ex)is fixed, it is solved calculating first (π, µ)through the solution of the compatible system Bt"π µ#=−∇Bf(ex), and then computing: λN=−Nt"π µ#− ∇Nf(ex), as, by definition of Z(ex)and in view of (13) we are led to (30) Z(ex)t∇f(ex) = −(B−1S)t∇Bf(ex) + ∇Sf(ex) =St[−(Bt)−1∇Bf(ex)] + ∇Sf(ex) =St"π µ#+∇Sf(ex). Here a series of propositions are presented to be used later. Let us consider now, for any x, a full-row-rank matrix A(x) = "A ∇c(x)t# partitioned as A(x) = [B S N], where Bis a nonsingular matrix, and Sand N have at least one column. Let (31) b A(x) = "B S N 00 00 1l #, and Z(x) =    −B−1S 1l 00   ,which satisfies b A(x)Z(x) = 0. 68 Let x=xbe a vector such that system (32) Atπ+∇c(x)µ+λ+∇f(x) = 0 is compatible. Suppose that λi= 0 for all iassociated with a column of either Bor S. Then, premultiplying this system by Z(x)twe have Z(x)tAtπ+Z(x)t∇c(x)µ+Z(ex)tλ+Z(x)t∇f(x) = 0, which, since b A(x)Z(x) = 0, implies Z(x)t∇f(x) = 0. Proposition 2.1. System (32) is compatible if and only if (33) Z(x)t∇f(x) = 0 is verified. Proof: The previous result to this proposition proves its first part. To prove the second part we consider the nonsingular (n×n)-matrix W=   1l 00 00 −(B−1S)t1l 00 00 00 1l   , where 00 stands for a zero matrix of suitable dimensions, and such that the unit matrix at the bottom is (n−t)×(n−t),(n−t)being the number of columns of N. Let us consider now that system (32) is not compatible. Hence, the matrix [b A(x)t∇f(x)] of this system has full column rank and we have W[b A(x)t∇f(x)] = W   Bt00 ∇Bf(x) St00 ∇Sf(x) Nt1l ∇Nf(x)   =   Bt00 ∇Bf(x) 00 00 Z(x)t∇f(x) Nt1l ∇Nf(x)   , where ∇Mf(x)is the gradient of f(x)at x=xwith respect to the subset of variables associated with M, for all submatrix Mformed by columns in A(x). Since the product of a nonsingular matrix of order n×nmultiplied by a matrix with full column rank of order n×(m+r+ (n−t) + 1) gives rise to a matrix with these same characteristics, it is shown that the product Z(x)t∇f(x)cannot be the null vector. Therefore, if (33) holds, system (32) is compatible. ¥ 69 technique put forward by Conn et al in [5], which is also employed by Murtagh and Saunders in [14]. Through this section and the following we add to assumption AS1, given in §2, the new one: AS2. µ∗satisfies the strict complementarity condition if cj(x∗) = cj⇒µ∗ j<0, if cj(x∗) = cj⇒µ∗ j>0, otherwise µ∗ j= 0. Solving problem IP is equivalent to solving: (45) (46) (47) (48) (49) minimize f(x) subject to: Ax =b c(x)−y= 0 l≤x≤u c≤y≤c, (EPy) where yrepresents the slacks vector that turns the inequalities (11) into equalities. As regards problem EP (1–4), the difference between this and problem EPy is that in the latter there are constraints (47) and (49) instead of constraints (3). Therefore we must analyze the effect of replacing (3) with (47) and (49) in the results of the former Section to see whether they can still be applied to the current problem. First, it can be observed that the Kuhn-Tucker conditions of an optimal solution (x∗, y∗) to problem EPy are ∇f(x∗) + ∇c(x∗)µ+Atπ+λ= 0(50) −µ+γ= 0(51) Ax∗=b(52) c(x∗)−y∗= 0,(53) such that unique vectors µ∗,π∗,λ∗and γ∗exist that satisfy equations (50) and (51), and with γ∗also satisfying (54) γ∗ i≤0 if y∗ i=ci γ∗ i≥0 if y∗ i=ci γ∗ i= 0 otherwise, 76 and λ∗ (55) λ∗ i≤0 if x∗ i=li λ∗ i≥0 if x∗ i=ui λ∗ i= 0 otherwise. Note that the only consequence on the multipliers of the introduction of the slacks y are expressions (51) and (54), the first being used to determine γ∗once the rest of the variables have been found through (50). Therefore, in order to obtain a first-order estimate of all multipliers at point (x, y)it suffices, as pointed out in [7], to solve (56) ∇f(x) + ∇c(x)µ+Atπ+λ= 0, just as in the case of problem EP, obtaining γafterwards through (51). Furthermore, expression (13) does not change — let us compare it with (56). Here the associated subproblem is defined as (57) (58) (59) (60) minimize x,y Lρ(x, y, µ) subject to: Ax =b l≤x≤u c≤y≤c (ESy) where (61) Lρ(x, y, µ) = f(x) + µt[c(x)−y] + 1 2ρkc(x)−yk2 2 is the augmented Lagrangian function. Furthermore, in Algorithm 2.1, subproblem ES (5–7) is replaced by subproblem ESy and c(x)with c(x)−y, hence the UHP type estimate at the optimizer (ex, ey), obtained from the solution of ESy, can be written as: (62) eµ=µ+ρ[c(ex)−ey], and the UKT type estimate is obtained solving ∇f(ex) + ∇c(ex)µ+Atπ+λ= 0, as in §2 with (13). Thus, as before in §2, we can now define matrices A(ex),Z(ex)and ZA, arriving at the results analogous to those obtained in §2, by propositions 2.1-2.3 and corollary 2.1. 77 Other modifications of the algorithm are associated with the substitution of xwith x, y. From the Kuhn-Tucker conditions for this subproblem we obtain ∇f(ex) + ∇c(ex){µ+ρ[c(ex)−ey]}+Atπ+λ= 0(63) −{µ+ρ[c(ex)−ey]}+γ= 0(64) such that vectors π=eπand λ=e λexist satisfying the first equation and a vector eγcan be obtained through the second one. As in §2 for expresion (20), here (63) is the key expression for studying the compatibility of estimate UKT when using the active set techniques of Murtagh and Saunders [13] to solve ESy, obtaining the results equivalent to those of the proposition 2.4. In conclusion, all the analysis after expression (20) in §2 is also applicable to this case, with the only exception that expression µ+ρc(ex)must be replaced by µ+ρ[c(ex)−ey]. 4. EXTENSION BY USING VECTORS y y yAND w w w Here the results of §2 are extended to the case of problems with general inequality contraints by using vectors yand w. This section describes an alternative way of dealing with problem IP (9–12). Solving problem IP, as put forward in [9] (using slack variables raised to the square, as done by Rockafellar in [16]), is equivalent to solving problem: (65) (66) (67) (68) (69) minimize f(x) subject to: Ax =b (c(x)−c)−y2= 0 l≤x≤u y2+w2=c−c, (EPyw) where y, w ∈Rrare auxiliary vectors of free variables — whithout bounds — that through vectors y2=    y2 1 . . . y2 r    and w2=    w2 1 . . . w2 r     allow the transformation of inequality (11) into an equality. 78 As regards the reference problem EP (1–4), the difference with problem EPyw (65– 69) is that in the latter instead of constraint (3) there are constraints (67) and (69), but contrary to what happens in Section 2, the number of bounds does not increase. Therefore, we must examine the effect of replacing (3) with (67) and (69) in the main calculation steps involved in the results of Section 2 and extend them to the type of problem now considered. It must be noticed first that the Kuhn-Tucker conditions to be satisfied by an optimizer (x∗, y∗, w∗)of problem EPyw are ∇f(x∗) + ∇c(x∗)µ+Atπ+λ= 0(70) (−µ+γ)ty∗= 0(71) γtw∗= 0(72) Ax∗=b(73) (c(x∗)−c)−(y∗)2= 0,(74) (y∗)2+ (w∗)2=c−c,(75) such that unique vectors µ∗,π∗,λand γ∗exist that satisfy equations (70), (71) and (72), and λ∗such that λ∗ i≤0 if x∗ i=li λ∗ i≥0 if x∗ i=ui λ∗ i= 0 otherwise. Note that the only consequence of the introduction of slacks y(apart from the constraints where they appear) are expressions (71) and (72), and γ∗can be obtained from these once the remainig multipliers have been computed from (70), given that in case y∗6= 0,γ∗=µ∗(due to (71)), otherwise, w∗6= 0 yields γ∗= 0 (due to (72)). Therefore, to obtain a first-order estimate of the UKT type at point (x, y, w)of all multipliers it suffices to solve (76) ∇f(x) + ∇c(x)µ+Atπ+λ= 0, as in problem EP (1–4) and then compute γthrough (71) and (72). Furthermore, expression (13), which coincides with (76), does not change. Here the associated subproblem would be minimize x,y Lρ(x, y, µ) subject to: Ax =b l≤x≤u y2+w2=c−c, (ESyw) 79 where (77) Lρ(x, y, µ) = f(x) + µt[(c(x)−c)−y2] + ρ 2k(c(x)−c)−y2k2 2 is the augmented Lagrangian function. Furthermore, see [9], the constraints where variables (y, w)appear can be eliminated by replacing the above augmented Lagrangian by (78) Lρ(x, µ) = f(x) + r X j=1 ½µjϕj[cj(x), µj, ρ] + 1 2ρ|ϕj[cj(x), µ, ρ]|2¾, where ϕj[cj(x), µj, ρ] =      cj(x)−cjif µj+ρ[cj(x)−cj]>0 cj(x)−cjif µj+ρ[cj(x)−cj]<0 −µj/ρ otherwise, the expression in braces that appears in (78) being continuously differentiable with respect to x, for fand cdefined as in §1, and (bearing in mind the strict complementarity assumption AS2 in §2) twice continuously differentiable with respect to xif x∈ X , where X=©x|µj+ρ[cj(x)−cj]6= 0, µj+ρ[cj(x)−cj]6= 0,∀j= 1,··· , rª. Therefore, according to [9], solving problem ESyw is equivalent to solving problem minimize xLρ(x, µ) subject to: Ax =b l≤x≤u. (ESx) Moreover, in Algorithm 2.1, replacing subproblem ES (5–7) with ESx and cj(x)with ϕj[cj(x), µj, ρ], for j= 1,··· , r, we have that the UHP type estimate at the optimizer ex, obtained through the solution of ESx, can be expressed by (see [2, 9]) (79) eµ=µ+ρϕ[c(ex), µ, ρ], such that ϕ≡[ϕ1,··· , ϕr]t, and the UKT type estimate is obtained solving ∇f(ex) + ∇c(ex)µ+Atπ+λ= 0, as in §2 with (13). Thus, as before in §2, we can now define matrices A(ex),Z(ex)and ZA, arriving at results analogous to those obtained in §2, by propositions 2.1-2.3 and corollary 2.1. 80 In addition, from the Kuhn-Tucker conditions associated with subproblem ESx we obtain: (80) ∇f(ex) + ∇c(ex){µ+ρϕ[c(ex), µ, ρ]}+Atπ+λ= 0, where vectors π=eπand λ=e λexist satisfying this equation. Hence, from (80), by an analogous process to that followed from expression (20), we are led to the same results as are obtained in §2 when using the active set techniques of Murtagh and Saunders [13] to solve ESx, obtaining the equivalent results to those of the proposition 2.4. In conclusion, all the analysis after expression (20) in §2 is also applicable to this case, with the only exception that expression µ+ρc(ex)must be replaced by µ+ρϕ[c(ex), µ, ρ]. 5. CONCLUSIONS This work shows when one may compute under certain guarantees the firts-order multiplier estimate based on the Khun-Tucker conditions. In addition, it proves the equivalence between this type of first-order estimate and that obtained through the original multiplier method (Hestenes and Powell’s method) when the exact minimization is used to solve the subproblem and it is not necessary to use columns of the constraint matrix that correspond to strongly active variables (with respect to the subproblem) to obtain a submatrix of the Jacobian (not including the simple bounds) having full row rank. Nevertheless, in practice not even in the latter case both procedures give the same estimate, as usually an inexact minimization of the subproblem is carried out. This work puts forward also a procedure to compute the multiplier estimates by solving a reduced system (41), instead of having to solve the large system (13), if the conditions so permit (see first lines of this paragraph). In the previous two sections specific procedures for transforming problems of type IP (9–12) into problems of type EP (1–4) have been considered. The procedure described in Section 4 as compared with that described in section 3 has the advantage that it does not increase subproblem size with respect to the original problem IP. Results of numerical tests comparing both procedures have not been included yet as, in our view, the construction of appropiate software to exploit the technique put forward by Conn et al in [5] −once fitted to the structure of problem EPy (45-49)−is by no means trivial and its coding is still underway. Using these procedures, inequalities in general constraints are eliminated. The validity for these problems of results and procedures put forward in [10] for type EP problems only has also been proved. 81 6. REFERENCES [1] Bertsekas D.P. (1977). «Approximation procedures based on the method of multipliers».J. Opt. Th. and Appl.,23, 487–510. [2] Bertsekas D.P. (1982). Constrained Optimization and Lagrange Multiplier Methods. Academic Press, New York. [3] Bertsekas D.P. (1995). Nonlinear Programming. Athena Scientific, Belmont, Massachusetts. [4] Bertsekas D.P. (1998). Network Optimization: Continuous and Discrete Models. Athena Scientific, Belmont, Massachusetts. [5] Conn, A.R., Gould, N. and Toint, Ph.L. (1994). «A note on exploiting structure when using slack variables».Mathematical Programming,67, 89–97. [6] Conn, A.R., Gould, N., Saternaer, A. and Toint, Ph.L. (1996). «Convergence properties of an augmented Lagrangian algorithm for optimization with a combination of general equality and linear constraints».SIAM J. on Optimization,6, 3, 674–703. [7] Gill, P.E., Murray, W. and Wright, M.H. (1981). Practical Optimization. Academic Press, New York. [8] Hestenes, M.R. (1969). «Multiplier and gradient methods».J. Opt. Th. and Appl., 4, 303-320. [9] Mijangos, E. and Nabona, N. (1996). «The application of the multipliers method in nonlinear network flows with side constraints».Technical Report 96/10, Dept. of Statistics and Operations Research. Universitat Polit` ecnica de Catalunya, 08028 Barcelona, Spain. [10] Mijangos, E. and Nabona, N. (1997). «On the compatibility of classical first order multiplier estimates with variable reduction techniques.»Technical Report 97/06, Dept. of Statistics and Operations Research. Universitat Polit` ecnica de Catalunya, 08028 Barcelona, Spain. [11] Mijangos, E. (1996). «Optimizaci´ on de flujos no lineales en redes con restricciones laterales mediante t´ ecnicas de multiplicadores».Ph.D. Thesis, Statistics and Operations Research Dept.. Universitat Polit` ecnica de Catalunya, 08028 Barcelona, Spain. [12] Mijangos, E. (1997). «PFNRN03 user’s guide».Technical Report 97/05, Dept. of Statistics and Operations Research. Universitat Polit` ecnica de Catalunya, 08028 Barcelona, Spain. [13] Murtagh, B.A. and Saunders, M.A. (1978). «Large–scale linearly constrained optimization».Mathematical Programming,14, 41–72. 82 [14] Murtagh, B.A. and Saunders, M.A. (1983). MINOS 5.0. User’s guide. Dept. of Operations Research, Standford University, CA 9430, USA. [15] Powell, M.J.D. (1969). «A method for nonlinear constraints in minimization problems».Optimization (R. Fletcher, ed.), 283–298. Academic Press New, York. [16] Rockafellar, R.T. (1971). «New applications of duality in convex programming». Proc. Confer. Probab., 4th, Brasov, Romania, 73–81. 83