minimize j(z)
,
s u b j e c t z:
5
zi5 z y
fori=
1,...,
n ,can be solved in NOA2 by a modification of the method described in the preceding section (Kiwiel, 1985c,1986c,1987).
The presence of finite upper and lower bounds ensures the existence of a solution and prevents divergence of the algorithm, which must occur when there is no solution (then (z k
I
tends, in theory, t o infinity; in practice-
until an arithmetic overflow terminates the calculation). It is always advisable t o place bounds of the form - 1 0 0 0 ~ ~ , ~ 1 0 0 0 , which should not be active when the solution lies inside the box.L U The objective j and its subgradient g j will be evaluated only inside the box [z , Z
I.
This ma be used t o eliminate regions where j is undefined. For example, if j(z)=z:Y2
+
e z p ( z 2 ) , it is essential t o place bounds of the form z l > ~ ~ 5 2 0 .If the user specifies an infeasible initial point zl, i t is projected on the box (by
L U k
replacing
21
withmax{z,
, m i n ( z i , z ,))
) . Successive z remain in the box.At the k-th iteration, an approximate feasible descent direction d is found t o k -k k
minimize j ( z +d)
+
u(dI2/2, ( 1 9 4L U
s u b j e c t z,
5
zf+d,5
z,,
for i = l ,...,
n.This subproblem is a natural extension of (16). Consequently, the preceding remarks on the choice of parameters remain in force.
We may add t h a t t h e introduction of box constraints only slightly increases the work a t the search direction finding.
For the problem with general linear constraints
minimize j(z), s u b j e c t A z
5
b, (20)K . Kiwiel, A . Stachurski
-
52 - Theoretical guide to N O A 2 the search direction finding subproblem becomes^k k
minimize
f
( z + d )+
uldI2/2,s u b j e c t A ( x k + d )
5
b. (21b)Due t o rounding errors, the calculated direction dk need not be "strictly" feasible. T o measure the infeasibility of a direction d we use the constraint violation function
v 2 ( d ) =max{ h ( x k + d ) ,0) zk
defined in terms of
h ( x ) = m a x
{
A,x-b': i = l ,...,
m ~ ) , (22) where A, denotes the i-th row of A . Subproblem (21) is equivalent t o t h e unconstrained problem^k k
minimize f ( x +d)
+
uldI2/2+
c v 2 ( d ) -all d (23) when the penalty parameter c is sufficiently large. Hence we may test increasing values of c until the solution of (23) is feasible, and hence solves (21). Starting from c = p , where p>O may be provided by the user, each successive c is multiplied by 10 until the solutiondk of (23) passes the feasibility test
h ( x k + d k )
<
e p , (24)where c~ is a positive absolute feasibility tolerance. If this test is failed by even "very large" c , the calculation terminates. This occurs if c > l / c M , where E M is the relative machine accuracy (the smallest positive c for which l + c > l in the computer's arithmetic).
No computational difficulties should arise if the linear constraints are well-scaled and the feasibility tolerance cF is large enough. In particular, it may be necessary t o ensure t h a t the coefficients of A are of order 1 and eF>cM2. For instance, if the coefficients of A result from measurements corrupted by errors of magnitude 10- 6
,
one should set c ~ = I o - ~ .If the initial point specified by the user is not feasible to within the tolerance EF, t h e algorithm tries t o project it onto the feasible set (by using a version of (23)). If the projec- tion is successful, each successive z k satisfies the linear constraints t o within c ~ . More- over, f(y) and g,( y) are calculated only a t eF - feasible points with
h(
y) <eF.A combination of the preceding techniques is used for the problem minimize
f
( x ) -all zsatisfying A x
5
b, x L5
x<
x u .In this case, all trial points satisfy the simple bounds exactly, and the general linear con- straints t o within e ~ .
K . Kiwiel, A . Stachurski - 53 - Theoretical guide to N O A 2 2.3. Exact penalty methods for convex constrained problems
The convex minimization problem
minimize I ( z ) -all z ( 2 5 4
satis/ying F , ( z )
<_
0 for j = l ,...,my, F , ( z ) = 0 for j = m I + l ,...,
m I + m E ,where the functions
I
andF ,
j = l ,...,
ml, are convex and the functionsF,,
j = m I + l ,...,
m I + m E , are affine (linear), may be solved in NOA2 by the uncon- I strained minimization of the ezact penalty functione ( z ; p ) = I ( z )
+
P F + ( z ) , (26)where p>O is a fized penalty coefficient, and the constraint violation is measured by F + ( z ) = m a x { F ( z ) , O ) ,
F ( z ) = m a x { F , ( z ) : j= l ,
...,
my, IF,(z)I: j = m I + l,...,
m I + m E } . Each solution z p t o the problemminimize e ( z ; p ) -all z in
R n
(27)solves (25) if it is feasible ( F ( z p ) < O ) . This holds if p is sufficiently large, (25) has a solu- tion and its constraints satisfy the generalized Slater constraint qualification, i.e. for some
F , ( z S ) < 0 , j= 1
,...,
my, Fj(zs)=O, j = m I + l,...,
m I + m E .The methods with a fixed penalty coefficient require the user t o specify a sufficiently large p . For well-scaled problems one may usually choose p in the interval I10,100]. If p is too small, (27) need not be equivalent t o (25)) and the algorithm may diverge when the penalty function has no finite minimum. On the other hand, too large p hinders the minimization of the penalty function, which becomes ill-conditioned. (If p is large, the algorithm must hug the boundary of the feasible set.)
The first method in NOA2 solves (26) by one of the algorithms for unconstrained minimization. A t the k-th iteration, a polyhedral approximation ~ ~ ( e t o ; ~e ( . ; p ) ) is con- structed from the past linearizations of e ( . ; ~ ) (see (5) and (6)). (These linearizations are calculated as in (5) from subgradients of the functions of (25)) which are evaluated by the user's subroutine.) The k-th search direction dk is chosen t o
-k k
minimize e ( z + d ; p )
+
u ldI2/2 -all d (28)(see (16)). Termination occurs if IukI
5 ss
( l + l e ( z k ; p ) l ) andF ( z k )
5
( F ,where E S and E F are positive final accuracy and feasibility tolerances, provided by the
k k k k
user, whereas uk is a dual estimate of the predicted descent e ( z + d ; p ) - e ( z ; p ) , which
K . Kiwiel, A . Stachurski - 54 - Theoretical guide t o N O A 2 satisfies the optimality estimate
where X* is a solution t o (25). This method does not exploit the specific structure of e ( . ; p ) .
T h e second method exploits the additive structure of e ( . ; p ) by constructing separate
^ k
polyhedral approximations f k and F t o the objective
f
and constraint function F . Thus the method may use a more accurate polyhedral approximation t o e ( . ; p )t k ( z ; p ) = j"(z)
+
p r n a x { F k ( z ) , o ) (31) in the search direction finding subproblem (28), which usually enhances faster conver- gence.Both methods may be allowed to choose the penalty coefficient automatically during the calculations (Kiwiel, 1985d). Then a t the k-th iteration we set p = p k in (28) and (31).
The initial p 1 may be specified by the user. The penalty coefficient is increased only if z k
!
k kis an a proximate solution t o (27) (i.e. z minimizes e ( . ; p
)
to within some positive toler- ance6
), but it is significantly infeasible (i.e. F ( z k ) is "large"). T h e specific rule for updating p k is- p and bk+'=bk;
if
-uk2 $
or F ( z k ) j - u k set pi+'- ( 3 2 4otherwise set pk+'=pk and bk+'= c,@, (32b)
where c p > l and c v € ( O , l ) are parameters t h a t increase the penalty and decrease the accu- racy tolerance of unconstrained minimization 6 k ; 6'=(u1(. Usually one may use pl=lO, c =2 or c = l o , and c,=0.1. Larger values of c and c, enable a faster
P P 1
growth of the penalty coefficient a t earlier iterations, if the initial p was too small. On the other hand, very large values of penalty coefficients slow down convergence.
When employing the exact penalty methods, the user should place sensible upper and lower bounds on all variables. If the box defined by such bounds is not too large, the penalty coefficient will quickly reach a suitable value and then will stay constant. More- over, box constraints ensure the existence of a solution and prevent the algorithm from diverging.
We may add t h a t the automatic choice of the penalty coefficient may produce a very large value of p k . The methods terminate a t the k-th iteration if p k + ' > l / f M , where E M is the relative machine precision. Such abnormal termination may indicate t h a t the con- straints are not regular (e.g. are inconsistent), or that they are ill-scaled.
In the current version of NOA2 additional general linear constraints A z j b can be handled only by the first method t h a t does not exploit the structure of the penalty func- tion.
2.4. The constraint linearization method