Comparing Nonsmooth Nonconvex Bundle Methods in Solving Hemivariational Inequalities.
Hemivariational inequalities can be considered as a generalization of variational inequalities. Their origin is in nonsmooth mechanics of solid, especially in nonmonotone contact problems. The solution of a hemivariational inequality proves to be a substationary point of some functional, and thus can be found by the nonsmooth and nonconvex optimization methods. We consider two type of bundle methods in order to solve hemivariational inequalities numerically: proximal bundle and bundle-Newton methods. Proximal bundle method is based on first order polyhedral approximation of the locally Lipschitz continuous objective function. To obtain better convergence rate bundle-Newton method contains also some second order information of the objective function in the form of approximate Hessian. Since the optimization problem arising in the hemivariational inequalities has a dominated quadratic part the second order method should be a good choice. The main question in the functioning of the methods is how remarkable is the advantage of the possible better convergence rate of bundle-Newton method when compared to the increased calculation demand.
Introduction
Hemivariational inequalities introduced by Panagiotopoulos are generalizations
of variational inequalities. By means of them problems involving
nonmonotone and multivalued constitutive laws and boundary conditions
can be defined mathematically. In many cases hemivariational
inequalities can be reformulated as substationary point problems of the
corresponding nonsmooth nonconvex energy functionals. For mathematical
theory and the applications of hemivariational inequalities we
refer to [23],[22].
The aim of this paper is to apply nonsmooth nonconvex optimization
methods for the numerical solution of hemivariational inequalities
and analyse their efficiency. As a typical example of hemivariational
inequalities we consider a laminated composite structure under loading
when the binding material between laminae obeys a nonmonotone multivalued
law. This kind of a mechanical problem has been investigated
numerically in [21],[26]. There it has been replaced by a sequence of convex
subproblems which have been solved by using convex minimization
methods.
The discretization of the considered problem is realized by the finite
element method scheme for nonmonotone multivalued differential inclusions
presented in [17],[18],[19]. This scheme has been proved to be
mathematically well-posed: stable and convergent.
Bundle methods are at the moment the most promising methods
for nonsmooth optimization. Their origin is the classical cutting plane
method of [4] and [7] and they base on the piecewise linear approximation
of the objective function. Due to the numerical experiments
the proximal bundle methods (see [10], [25] and [16]) seem to work in
the most efficient and reliable way. They can be called also to diagonal
variable metric methods, since a stabilizing quadratic term in form of
diagonal matrix was added to the polyhedral approximation in order
to accumulate some second order information about the curvature of
the objective function.
However, the development of "real" second order method has been
fascinating the researchers of nonsmooth optimization during its whole
history. Several attempts has been done in order to exploit the second
order subderivative information. Already in his pioneering work [11]
Lemar'echal derived a version of variable metric bundle method utilizing
the classical BFGS secant updating formula from smooth optimization.
Due to the disappointing numerical results in [12] this idea was buried
nearly for two decades. Several modifications of the variable metric concepts
has been proposed for example in [6], [13], [3] and [20]. According
to very limited numerical experiments (see for example [6]) it seems
that the variable metric bundle methods works fairly well. However,
when proportion the results to the extra computational efforts needed
with the full matrix algebra they do not offer substantial advancement
in numerical solution process.
More recently a new second order approach has been proposed in
[15], where the bundle-Newton method was introduced. The main difference
compared to the earlier methods was include second order information
directly to the model, in other words the piecewise linear model
was replaced by piecewise quadratic model. In numerical tests of [15]
it turns out to be very effective for quadratic type of problems.
The aim of this paper is to compare the proximal bundle method
and the bundle-Newton method in solving hemivariational inequalities.
Since the optimization problem arising in the solution of hemivaria-
tional inequalities has a dominated quadratic part the second order
method should be a good choice. The main question in the functioning
of the method is how remarkable is the advantage of the possible better
convergence rate of the bundle-Newton method when compared to
the increased calculation demand. The numerical results indicate the
superiority of the bundle-Newton method for the problems having a
quadratic nature.
The paper is organized as follows. In chapter 2 we introduce the
considered optimization problem and give some definitions needed in
continuation. The chapters 3 and 4 are devoted to the proximal bundle
method and bundle-Newton method, respectively. In chapter 5 the test
hemivariational inequality problem is introduced and finally in chapter
6 we give some numerical results.
2. Nonsmooth and Nonconvex Optimization Problem
We consider the following nonsmooth and nonconvex optimization problem
ae minimize f(x)
subject to x 2 K;
where the objective function f from R n to R is supposed to be locally
Lipschitz continuous function. The feasible set K has a more specific
structure, i.e.
where x l and x u are the lower and upper bounds for variables, respec-
tively. We suppose that at each x 2 K we can evaluate the function
value f(x) and an arbitrary subgradient g(x) from the subdifferential
of Clarke (see [5])
The normal cone of K at x 2 K can be defined (see [16]) by
cl
ae [
oe
R is the distance function of K, i.e.
ck
Definition 1. Suppose that f : R n ! R is locally Lipschitz continu-
ous. Then x   is called a substationary point of the problem (P) if
Now we can formulate the following necessary optimality condition.
Theorem 1. Every local minimizer of the problem (P) is substation-
ary.
For the proof we refer to [16].
In the following chapters we describe two methods: proximal bundle
method and bundle-Newton method for finding local minimizers (and
thus substationary) points of the problem (P). For the convergence of
the methods we need the following semismoothness assumption due to
[2].
Definition 2. The function f : R n ! R is said to be upper semis-
mooth, if for any x 2 R n , d 2 R n and sequences g i ae R n and t i ae (0; 1)
satisfying d) and t i # 0, one has
lim sup
3. Proximal Bundle Method
In this chapter we shortly describe the ideas of the proximal bundle
method for nonsmooth and nonconvex minimization. For more details
we refer to [10], [25] and [16].
3.1. Direction finding
Our aim is to produce a sequence fx k g 1
converging to some
local minimum of the problem (P). Suppose that the starting point x 1
is feasible and at the k-th iteration of the algorithm we have the current
iteration point x k and some trial points y j 2 R n (from past iterations)
and subgradients where the index set J k is a
nonempty subset of
The idea behind the proximal bundle method is to approximate the
objective function below by a piecewise linear function, in other words,
we replace f by so called cutting-plane model
which equivalently can be written in the form
with the linearization error
Note, that in convex case the subdifferential can be rewritten as (see
Then it is easy to proof that
In other words, if f is convex, then the cutting-plane model -
f k is an
under estimate for f and the nonnegative linearization error ff k
how good an approximation the model is to the original problem.
In nonconvex case these facts are not valid anymore: ff k
may have a tiny
(or even negative) value, although the trial point y j lies far away from
the current iteration point x k and thus the corresponding subgradient
j is useless. For these reasons the linearization error ff k
j is replaced by
so called subgradient locality measure (cf. [8])
where is the distance measure parameter
and
is the distance measure estimating
without the need to store the trial points y j . Then obviously fi k
for all j 2 J k and min x2K
min
In order to calculate the search direction d k 2 R n we replace the
original problem (P) by the cutting plane model
subject to x k
(CP)
where the regularizing quadratic penalty term 1u k d T d is added to guarantee
the existence of the solution d k and keep the approximation local
enough. The weighting parameter added to improve the
convergence rate and to accumulate some second order information
about the curvature of f around x k . It was adapted from the proximal
point algorithm by [24] and [1] and was first time used in [10] and [25].
Notice, that the problem (CP) still is a nonsmooth optimization
problem. However, due to piecewise linear nature it can be rewritten
as a (smooth) quadratic programming subproblem finding the solution
subject to \Gamma fi k
and
3.2. Line search
In previous section we calculated the search direction d k . Next we
consider the problem of determining the step size into that direction.
We assume that mL 2 (0; 1), mR 2 (mL ; 1) and - t 2 (0; 1] are fixed
line search parameters. First we shall search for the largest number
is the predicted amount of descent and it holds
due to (11). If such a parameter exists we take a long serious step
Otherwise, if (11) holds but
then a short serious step
R d k
is taken and if t k
we take a null step
L is such that
In long serious step there occurs a significant decrease in the value
of the objective function. Thus there is no need for detecting discontinuities
in the gradient of f , and so we set g k+1 2 @f(x k+1 ). In short
serious steps and null steps there exists discontinuity in the gradient of
f . Then the requirement (12) ensures that x k and y k+1 lie on the opposite
sides of this discontinuity and the new subgradient g k+1 2 @f(y k+1 )
will force a remarkable modification of the next search direction finding
problem. In what follows we are using the line search algorithm presented
in [16]. The convergence proof of the algorithm is assumed f to
be upper semismooth (see (2)).
The iteration is terminated if
final accuracy tolerance supplied by the user.
3.3. Weight updating
One of the most important questions concerning proximal bundle method
is the choice of the weight u k . The simplest strategy might be to keep
it constant u k j u fix . This, however, leads to several difficulties. Due
to Theorem 4 we observe the following:
If u fix is very large, we shall have small jv k j and kd k k, almost all steps
are serious and we have slow descent.
If u fix is very small, we shall have large jv k j and kd k k, and each serious
step will be followed by many null steps.
Therefore, we keep the weight as a variable and update it when neces-
sary. For updating u k we use the safeguarded quadratic interpolation
algorithm due to [10].
4. Bundle-Newton Method
Next we describe the main ideas of the second order bundle-Newton
method. For more details we refer to [15].
4.1. Direction finding
We suppose that at each x 2 K we can evaluate, in addition to the
function value f(x) and an arbitrary subgradient g(x) 2 @f(x), also
an n \Theta n symmetric matrix G(x) approximating the Hessian matrix
f(x). For example, at the kink point y of piecewise twice differentiable
function we can take "infinitely close"
to y.
Instead of piecewise linear cutting-pane model (3) we introduce a
piecewise quadratic model of the form
~
damping parameter. The
model (14) can again equivalently be written as
~
and for all j 2 J k the linearization error takes now the form
Note that now even in the convex case ff k
j might be negative. Therefore
we replace the linearization error (16) again by the subgradient locality
measure (7) and we remain the property (see [15])
min
~
The search direction d k 2 R n is now calculated as the solution of
minimize ~
d)
subject to x k
(CN)
Note, that since the model already has second order information no
regularizing quadratic terms are needed like in (CP). The problem (CN)
is transformed to a nonlinear programming problem, which is then
solved by a recursive quadratic programming method (see [15]). If we
denote
this procedure leads to a quadratic programming subproblem finding
the solution (d
subject to \Gamma fi k
and
where
and -
are the Lagrange multipliers of (QN) from the
previous iteration k \Gamma 1. In calculations W k is replaced by its positive
definite modification, if necessary.
4.2. Line search
The line search operation of the bundle-Newton method is following
the same principles than in section 3.2 for proximal bundle method.
The only remarkable difference occurs in the termination condition for
short and null steps, in other words (12) is replaced by two conditions
and
where supplied by the user.
The bundle-Newton method is using the line search algorithm presented
in [15]. The convergence proof of the algorithm is assumed f to
be upper semismooth (see (2)).
5. Formulation of the problem
5.1. Continuous problem
We consider a two-dimensional laminated composite structure consisting
of two elastic laminae and a interlayer binding material under loading
(see Fig.1). The mechanical behaviour of the binding material, the
interlaminar bonding forces versus the corresponding relative displacements
of the laminae, is depicted by Fig.2. This relation is typically
nonmonotone and multivalued. Therefore, this kind of mechanical
problems lead to the hemivariational inequalities or in the energy formulation
to the substationary point problems of nonconvex, nonsmooth
functionals. The similar mechanical structures have been investigated
in [21],[26].
Due to the symmetry of the structure and by assuming that the
forces applied to the upper and lower part of the structure are the
same it is enough to study the upper lamina. Next, we formulate
mathematically the problem. We denote
by\Omega ae R 2 the upper lamina
in its undeformed state. The Lipschitz boundary \Gamma
of\Omega consists
of four nonoverlapping open subsets We denote
by
i;j=1 the stress tensor,
i;j=1 the strain tensor,
F
F
G
G
G
G24

Figure

1. Laminated composite structure under loading
Relative normal displacement of laminae in mm
Interlaminar
binding
force
in
kN/m
F
G2000040000
A
Figure

2. Nonmonotone adhesive force between laminae
the displacement,
the outward unit normal
vector to \Gamma and
the boundary force.
Assuming that the deformations are small the lamina obeys the
Hooke's law of form
@x l
and
i;j;k;l=1 is the elasticity tensor satisfying the usual symmetry
and ellipticity conditions. The equation of the equilibrium state
of\Omega is as follows:
since there are no volume forces.
Next, we define the boundary conditions. On \Gamma 1 we have given displacement

and on both \Gamma 2 and \Gamma 3 given boundary forces
where F is a constant force.
Furthermore, on \Gamma 4 we have in the tangential direction a given
boundary force
On the other hand, due to the binding interlayer material it holds a
nonmonotone, multivalued boundary condition (see Fig.2) expressed
by means of the subdifferential @j of a locally Lipschitz function j
and because of the nonpenetration of the laminae a unilateral boundary
condition
in the normal direction. Note that we have to scale by 1the x 1 -axis of
the graph in Fig.2 in order to get the nonmonotone law of (26), because
in Fig.2 we have relative displacements of the laminae.
By employing Green-Gauss theorem and taking into account both
the equilibrium equation (21) and the boundary conditions (22)-(27)
we obtain the following hemivariational inequality: Find u 2 K such
that
Z
Z
where K is a set of kinematical admissible displacements defined by
G
G
F
GG

Figure

3. 32\Theta4 triangulation Th
of\Omega
(\Omega\Gamma denotes the Sobolev space), j ffi is the generalized directional
derivative of j and a the bilinear form of the linear elasticity defined
by
Z
Due to the symmetry of the bilinear form a the potential energy of
the considered mechanical system has the form:
Z
Z
Then, we can formulate the following problem: Find u 2 K such
that u is a substationary point (cf. (1)) of f on K, i.e.
where NK (u) is the normal cone to the nonempty, closed and convex
set K at u. It can be shown that the problem (28) has at least one
solution (see [22]), and that the solutions of the inclusion (33) are also
the solutions of the inequality (28) (see [19]).
5.2. Discrete problem
We apply the finite element approximation scheme developed in [17],[18].
Let h be a discretization parameter related to the mesh size of the triangulation
of\Omega (see Fig.3). The set K is approximated by the set
of piecewise linear functions over the triangulation T h defined by
For the approximation of the bilinear form a and the boundary inte-
gral
R
F (v 2 \Gammau 2 )d\Gamma we use appropriate numerical integration formulae
(this is standard in the finite element method). The nonmonotone term
in (28) is approximated by the following numerical integration formula
Z
i2I
where fx i g i2I is the set of the nodal points of the triangulation T h on
are the coefficients of the integration formula.
Let us make the identification X h with the corresponding subset of
R 2n , i.e. the displacement function v is identified with the displacement
vector of the nodal points (n is the number of the nodal points). Then
the discrete problem in the matrix form is as follows: Find u 2 K h such
that
i2I
where A is the stiffness matrix of the structure, F the (discrete) load
vector and I the set of the indices corresponding to the x 2 -displacements
of the nodal points on \Gamma 4 .
Then we can formulate the discrete substationary point problem on
where
i2I
Because of Corallary 1 of Proposition 2.3.3 in [5] and the terms j(u i ),
are independent of each others we have that
i2I
Now it holds that every solution of (38) is a solution of the problem (37)
(see [19]). Also due to the results in [18],[19] we know that the discrete
problem (38) is solvable and its solutions converge in subsequences to
the solutions of the continuous problem (33).
Remark 1. The functional f h is upper semismooth. Indeed, let u; d 2
R 2n be given. ?From Definition 2 of the upper semismoothness we see
that it is enough to study the restriction of f h onto the line
Rg, on which f h is smooth except finite number of points (cf.
and Fig.2). If f h j L is smooth at u, it is smooth also on some small
(one-dimensional) neighborhood of u and the condition
lim sup
is trivially satisfied. On the other hand, if f h j L is nonsmooth in u, there
exists small neighborhood of u in which u is the only nonsmooth point.
Therefore, f h j L is continuously differentiable in this set except at the
point u and the classical one-sided derivatives exist at u implying (41).
Since the nonlinear behaviour and the constraints of the problem
have an effect only on the x 2 -displacements on the nodes on \Gamma 4 , we can
make our computation much more effective by applying the method of
condensation of unknowns. Assuming that m components corresponding
to the index set I are listed first, we have the following decomposition
of the matrix A and the vectors F and u
A 21 A 22
where A 11 is an m \Theta m matrix and F 1 ; . Then the elimination
of u 2 (2 R 2n\Gammam ) from (37) leads to the following discrete potential
functional
~
where ~
called the Schur complement and
~
. Hence, the eliminated substationary point
problem is formulated as follows: Find u 1 2 K h such that
where K h is interpreted as a set of R m .
6. Numerical Results
The optimization algorithms has been implemented in Fortran 77 and
the test runs have been performed on an HP9000/J280 (180 MHz)
computer. The tested optimization codes are presented in Table I.
All the codes are utilizing the subgradient aggregation strategy of
[8] to keep the storage requirements bounded. The code PB is employing
the quadratic solver QPDF4, which is based on the dual active

Table

I. Tested codes and software.
Code Software Package Author(s) Method
PB NSOLIB M-akel-a Proximal bundle
PBL UFO Luk-san &Vl-cek Proximal bundle
BNL UFO Luk-san &Vl-cek Bundle-Newton

Table

II. Eliminated with discretization 16 \Theta 2.
Code Ni Nf CPU f Ni Nf CPU f
set method derived in [9], while the codes PBL and BNL are using
the solver ULQDF1 implementing the dual projected gradient method
proposed in [14].
The above optimization codes have been applied to the mechanical
structure of Fig.1 under eleven constant loadings (F=20000, 21000,
In calculation we have used the incremental procedure: the loading of
the structure is increased uniformly and as an initial guess for the
solution of the next load it is used the solution of the previous load.
The bigger increment from 25000 kN/m 2 to 26200 kN/m 2 is due to the
fact that the structure is very sensitive to the increase of the loading
between 25000-27000 kN/m 2 where the partial (branches B-G in Fig.2)
and the complete delamination (branch G-H in Fig.2) take place. With
the load 26200 kN/m 2 we can illustrate the case in which only the
partial delamination occurs. We have applied the plane stress model

Table

III. Eliminated with discretization 4.
Code Ni Nf CPU f Ni Nf CPU f
PBL 426 434 3.15 -0.871201 3007 3063 14.16 -15.21929

Table

IV. Eliminated with discretization 64 \Theta 8.
Code Ni Nf CPU f Ni Nf CPU f

Table

V. Eliminated with discretization 128 \Theta 16.
Code Ni Nf CPU f Ni Nf CPU f
with the elasticity modulus the Poissons ratio
0:3 and the thickness
The calculated results with different discretizations are presented in

Tables

II-VI, in which Ni denotes the number of iterations, Nf denotes
the number of objective function (and also subgradient) evaluations,
CPU the used total computer time in seconds and f is the objective
function value at termination. Note that the elimination, the computation
of the Schur complement ~
A and the new load vector ~
F , for the
discretization 4, 64 \Theta 8 and 128 \Theta 16 is taking about 0.67,
1.35, 23.1 and 545.0 seconds, respectively. The stopping criterions are
chosen to be comparable, in other words the relative accuracy in the
optimal objective function value is about 6 digits.

Table

VI. Not eliminated with discretization 4.
Code Ni Nf CPU f Ni Nf CPU f
PB 4458 4474 156.97 -0.880013 40112 40128 1458.51 -15.22886

Table

VII. Average error of the normal displacements on \Gamma 4 by
BNL.
Discretization Error in mm Error in mm
128 \Theta 16 *
From the numerical results in Tables II-VI we can conclude the
superiority of the bundle-Newton code BNL. In all cases it used less
computing resources and found the local minimum in the most reliable
way. Due to the matrix operations the individual iteration is more
costly than in proximal bundle methods. However, the total amount of
iterations stays very low, and thus the used CPU time of BNL is always
less than with the other codes.
Note, that the iteration number and the function evaluations of BNL
do not depend on the dimension of the problem. On the other hand,
when the size of problem is doubled, BPL needed about two times and
PB about 1.6 times more iterations and function evaluations.
When comparing the CPU times versus discretization the behaviour
is nonlinear and the multiplier is growing exponentially. We have to
remember, that especially in the large problems the elimination is taking
most of the time and the differences do not seem to be so remark-
able. For example for the discretization 128 \Theta 16 the optimization time
of BNL is only about 2 % of the total time.
By comparing the Tables III and VI we can see the influence of the
elimination; the solution of the problem without elimination is taking
nearly ten times more CPU time than the eliminated one. Note that
especially for the proximal bundle methods the elimination is essential.
Although the codes PB and PBL are realizations of the same method,
the difference in their functioning is remarkable. PBL is clearly more
effective when comparing the function evaluations and CPU times.
However, it has slight difficulties to reach the desired accuracy in larger
problems (see the function values denoted by    in Tables V-VI). Notice
also, that PB converges to the different local minimum with the load
(see the function values denoted by   in Tables III-V).
When we changed the starting point to be zero in all the components it
found the same optimum as the other codes. This behaviour is due the
fact that the delamination takes place between the loads 25000-27000
Normal
displacements
in
mm
Interface in mm
load 23000 kN/m^2
load 26200 kN/m^2
load 27000 kN/m^2
load 30000 kN/m^2

Figure

4. Normal displacements of the interface of the upper lamina
making the structure very unstable and sensitive to the increase
of the load.
In

Table

VII we have listed the average errors of the normal displacements
with the loads 26200 kN/m 2 and 27000 kN/m 2 obtained
by BNL. As an exact solution it is used the solution obtained by the
discretization 128 \Theta 16. Table VII indicates that the applied approximation
scheme has a very good convergence rate and the results calculated
by BNL are reliable.
In Fig.4 and Fig.5 it is presented the normal displacements of the
upper lamina (the relative displacements of laminae are twice larger)
and the binding forces between laminae. Further, Fig.6 illustrates the
progress of delamination: With the load 23000 kN/m 2 there is no damage
of the binding material, all the nodes are on the branch A-B. The
partial delamination occurs when the load is increased to 26200 kN/m 2 ,
some of the nodes are on the branch C-D, and the complete delamination
when the load is greater than 27000 kN/m 2 , most of the nodes are
on the branch G-H.
7. Conclusions
The solution of a hemivariational inequality can be found as a sub-
stationary point of some nonconvex functional being composed of a
Interlaminar
binding
force
in
Interface in mm
load 23000 kN/m^2
load 26200 kN/m^2
load 27000 kN/m^2
load 30000 kN/m^2

Figure

5. Nonmonotone adhesive force between laminae
branch G-H
F=23000 kN/m
F=27000 kN/m
F=30000 kN/m22
branch A-B branch C-D branch E-F

Figure

6. Progress of delamination
dominating quadratic part and a nonsmooth piecewise quadratic part.
We have tested the functioning of different nonconvex and nonsmooth
optimization methods in the solution of a laminated composite structure
under loading when the binding material between laminae obeys
a nonmonotone multivalued law.
Due to the strong quadratic nature of this kind of problem, the
bundle-Newton method based on the second order piecewise quadratic
model proves to be superior when compared to proximal bundle
method based on the first order polyhedral approximation. Bundle-
Newton method was clearly faster and more reliable, and what is best,
the iteration number and the function evaluations of bundle-Newton
method do not depend on the dimension of the problem. The same
trend can be seen also in the academic tests of [15].

Acknowledgements

This work was supported by the grants SA-34063 and SA-8583 given by
the Academy of Finland and the grant No. 201/96/0918 given by the
Czech Republic Grant Agency. The research was partially done during
the third author's visit at the Laboratory of Scientific Computing



--R

Numerical Methods for Nondifferentiable Convex Opti- mization
Optimization of Upper Semidifferentiable Functions
Gilbert J.
Golstein A.
Optimization and Nonsmooth Analysis
Variants to the Cutting Plane Approach for Convex Nondifferentiable Optimization
The Cutting Plane Method for Solving Convex Programs
Methods of Descent for Nondifferentiable Optimization
A Method for Solving Certain Quadratic Programming Problems Arising in Nonsmooth Optimization
Proximity Control in Bundle Methods for Convex Non-differentiable Optimization
Nonsmooth Optimization and Descent Methods
Numerical Experiments in Nonsmooth Optimization
An Approach to Variable Metric Bundle Methods
Dual Method for Solving a Special Problem of Quadratic Programming as a Subproblem at Linearly Constrained Nonlinear Minmax Approximation
A Bundle-Newton Method for Nonsmooth Unconstrained Minimization

Approximation of Non-Monotone Multivalued Differential Inclusions
Finite Element Approximation of Vector-Valued Hemivariational Inequalities

A Quasi-second-order Proximal Bundle Algorithm
Baniotopoulos C.
Mathematical Theory of Hemi- variational Inequalities and Applications
Hemivariational Inequalities
Monotone Operators and the Proximal Point Algo- rithm
A Version of the Bundle Idea for Minimizing a Nonsmooth Functions: Conceptual Idea

--TR
bundle methods;hemivariational inequalities;nondifferentiable programming;nonmonotone contact problems;substationary points