Identification of some nonsmooth evolution systems with illustration on adhesive contacts at small strains.
Lukáš Adam, Jiří V. Outrata, Tomáš Roubíček.
Institute of Information Theory and Automation, Academy of Sciences,
Pod vodárenskou věží 4,
CZ–18208 Praha 8, Czech Republic.
Faculty of Mathematics and Physics, Charles University, Sokolovská 83,
CZ–18675 Praha 8, Czech Republic
Institute of Thermomechanics of the ASCR,
Dolejškova 5, CZ–18200 Praha 8, Czech Republic.
Abstract: A class of evolution quasistatic systems which leads, after a suitable time discretization, to recursive nonlinear programs, is considered and optimal control or identification problems governed by such systems are investigated. The resulting problem is an evolutionary Mathematical Programs with Equilibrium Constraints (MPEC). A subgradient information of the (in general nonsmooth) composite objective function is evaluated and the problem is solved by the Implicit programming approach. The abstract theory is illustrated on an identification problem governed by delamination of a unilateral adhesive contact of elastic bodies discretized by finiteelement method in space and a semiimplicit formula in time. Being motivated by practical tasks, an identification problem of the fracture toughness and of the elasticity moduli of the adhesive is computationally implemented and tested numerically on a twodimensional example. Other applications including frictional contacts or bulk damage, plasticity, or phase transformations are outlined.
Keywords: rateindependent systems, optimal control, identification, fractionalstep time discretization, quadratic programming, gradient evaluation, variational analysis, implicit programming approach, limiting subdifferential, coderivative, nonsmooth contact mechanics, delamination.
AMS Classification: 35Q90, 49N10, 65K15, 65N38, 65M32, 74M15, 74P10, 74R20, 90C20.
1 Introduction
Many evolution systems have the structure of the generalized gradient flow
with functionals and . Here is an abstract state of the system and denotes its time derivative. Quite typically, is convex and, making the conjugate of with respect to the “driving force” variable , i.e. , the generalized gradient flow can equivalently be written in the Biotequation form
(1) 
In many cases, the problem is nonsmooth due to a nonsmoothness of or , which is why we wrote inclusion in (1) rather than equality. Ansatz (1) is very general and covers great variety of problems in particular in nonsmooth continuum mechanics. The state variable may involve displacements and various internal parameters (but also various concentrations of some constituents subjected to diffusive processes). In this paper, we focus on a subclass of such problems where the state has the structure
(2) 
for each time instance in a Banach space . In this way, a quasistatic plasticity, or damage, or various phase transformations at small strains can be modelled, and also various problems in contact mechanics like friction or adhesion, together with various combinations of these phenomena.
After a suitable time discretization, (1) gives rise to recursive optimization problems. Often (or, in applications in continuum mechanics we have in mind, rather typically) ranges over an infinitedimensional Banach space and, after a possible “spatial” discretization, these minimization problems have a structure of strictly convex Quadratic Programming problems. It is then relatively easy to use such a discretized evolution problem as a governing system for some optimization problem, e.g., optimal control or identification of parameters. This leads to Mathematical Programs with Evolution Equilibrium Constraints (MPEEC) which have been studied e.g. in [17, 18, 27].
The functionals in (1) depend also on an abstract parameter and have a special form
We then consider an optimalcontrol or an identification problem on a fixed time interval :
(3) 
with some and specified later; here is a closed convex set of a Banach space where lives. In some models, the flow rule for in (3) is purely static, i.e. . In this case, if there is no dissipation in this part, then only but not is decisive when considering . We use the standard notation for Bochner space of Banachspacevalued Bochner measurable functions on .
The main aim of the paper is a deep analysis of a discretized version of MPEEC (3) leading both to sharp necessary optimality conditions as well as to an efficient numerical procedure based on the socalled Implicit Programming approach (ImP), cf. [21, 28]. In particular, on the basis of the subdifferential calculus of B. Mordukhovich [24, 32] we will show that the solution map defined via the discretized equilibrium relations in (3), is singlevalued and locally Lipschitz and satisfies henceforth the basic ImP hypothesis. In the respective proof one has to deal with a difficult multifunction arising in connection with our evolving constraint sets. The application of standard tools of generalized differential calculus provides us in this case only with an upper estimate of the coderivative of the normal cone mapping to the overall constraint set, which could be a substantial drawback both in the optimality conditions as well as in the used numerical approach. To overcome this hurdle, we have employed a new result from [2] which enables us to compute the limiting coderivative of the mentioned normal cone mapping exactly.
The plan of the paper is the following. In Section 2, we briefly introduce a suitable discretization of the identification problem (3) that yields a unique response of the constraint system of (3) for a given and allows for efficient optimization technique. After introducing some notation and notions from variational analysis in Section 3, we formulate in Section 4 firstorder necessary optimality condition for the discrete version of MPEEC (3) and derive a subgradient formula for the composite objective function of the discretized problem. In Section 5, we formulate a specific applicationmotivated identification problem from contact continuum mechanics that fits (and illustrates) the system (3). Eventually, in Section 6, we illustrate the usage of the subgradient evaluation procedure via an adhesive contact problem in a nontrivial twodimensional example involving a spatial discretization by FiniteElementMethod (FEM).
2 Discretization of the identification problem (3)
The natural procedure is to discretize the problem (3) in time. This might be a rather delicate problem especially when the controlled system (3b,c) exhibits responses of different time scales, and in particular with tendencies for jumping, which quite typically happens in rateindependent systems governed by nonconvex potentials .
We consider (for simplicity) an equidistant partition of the time interval with a time step such that and then a fractionalsteptype semiimplicit time discretization of (3b,c). Moreover, if , or is infinitedimensional, the timediscrete problem still remains infinitedimensional, and to implement it on computers, we need to apply also an abstract space discretization controlled by a parameter, let us denote it by . Such an approximation can be considered by replacing , , and in (4) by their finitedimensional subsets , , and . Counting also a possible numerical approximation of , denoted by , altogether, (3) turns into the problem:
(4) 
where is an approximation of the initial condition . Let us note that the controlled system (4b,c) decouples so that, for a given , one is to solve alternating optimization problems
Minimize  (5a)  
and, taking (one of) its solution for , further  
Minimize  (5b) 
which yields as (one of) its solution. Assuming as well as its approximation separately strictly convex (and, of course, coercive with compact level sets) and convex, , both problems in (5) have unique solutions and , respectively, and thus the whole recursive problem in the constraint system (4) has a unique response for a given as well. This allows us to reformulate (4) as a minimization problem for a functional depending on only, cf. (9) below. This will be exactly the situation we will consider in the rest of this article. The fully discretized system (4) can thus be understood as an MPEC for which a developed theory exists.
In what follows, we will confine ourselves to problems with a bit more detailed (but nevertheless still fairly general) structure, namely
(6a)  
(6b)  
(6c) 
where , , and are finite and smooth, , , and are convex closed set, the last one being a cone. We will use as a possible approximation of .
Although, in Section 6, we will illustrate usage of this model on a rather special inverse adhesivecontact problem, most of the considerations can expectedly be applied (after possible modification) to many other problems from continuum mechanics and physics, as (various combination of) damage, phasetransformations, plasticity, etc.
Remark 2.1 (Stability and convergence for and ).
The focus of this article is on the identification of the discrete finitedimensional problem. Nevertheless, the convergence towards the original continuous problem when and is of interest.
Without going into (usually rather technical) details, let us only mention that under certain qualification of , and , a boundedness (= numerical stability) and convergence of a solution to the discrete state problem obtained by interpolation from values towards a weak solution to controlled state system for a fixed can usually be shown at least in terms of subsequences in various situations. A rather simple situations is if , or possibly also , is uniformly convex; this corresponds to some viscosity. In a special fully rateindependent case when and is 1homogeneous and independent of , such convergence was proved in [33]; in this case the weak solutions are called local solutions. The uniqueness is however not guaranteed in general. If is jointly uniformly convex, then this uniqueness and even continuous dependence on holds, cf. [22, 23] for a survey of such situations. This is e.g. the case of linearized rateindependent plasticity with hardening. Sometimes, viscosity can help for this uniqueness. This is the case of frictional normalcompliance contact of viscoelastic bodies which, after a certain algebraic manipulation gets the structure with separately uniformly quadratic with linear constraints in twodimensions, cf. [34], or with conical constraints in threedimensions. The uniqueness of the response of the continuous problems was shown in [13].
As usual, the convergence of solutions to (4) towards solutions to (3) is much more delicate and it is a wellknown fact that it cannot be expected unless the controlled state system in (3) has a unique response or at least any solution to (1) can be attained by the discretized solutions, which is however usually not granted unless the solution to (1) is unique. In any case, one needs to show the continuous convergence of the solution map , i.e. that and and implies . This is usually a relatively simple modification of the convergence for fixed.
3 Notation and selected notions of variational analysis
Having in mind the discrete problem with and , we will use notation
(7a)  
(7b)  
(7c)  
(7d) 
Since problems (5) are convex, necessary optimality conditions are also sufficient and thus, taking into account structure (6), problem (4) can equivalently be written in the form:
(8) 
with . Defining the solution map implicitly via system (8), we may used the socalled implicit programming approach to rewrite problem (8) equivalently into the form
(9) 
In the rest of the paper we will make use of the following standing assumption, which imply in particular the singlevaluedness of the solution map :

and are strictly convex,

and are convex,

and are continuously differentiable mappings, and

and are closed convex sets.
Note that (A1)–(A3) implies that and have a positive definite Jacobian.
In what follows, we will fix time (and, if any, also space) discretization and thus we will omit and in the following sections. The dimension of , , and will be respectively denoted by , , and .
Before devising a (necessarily) quite complicated procedure to evaluate a gradient information for the nonsmooth functional , let us still briefly present basic notions from variational analysis which are essential for this paper. More information can be found in [32] for finitedimensional setting or in [24] and [7] for the general infinitedimensional case.
All objects in this section are finitedimensional. For sequence of sets we define the PainlewéKuratowski upper limit as
Using this construction, we define for any the Bouligand tangent cone, Fréchet normal cone and the limiting normal cone respectively as
where by we understand with . If , then we say that is a regular point of , otherwise it is a nonregular point.
To a function we can define its subdifferential as
If happens to be convex, both normal cones coincide and are equal to the normal cone in the standard sense of convex analysis
and similarly, if is continuously differentiable at , then .
For a setvalued mapping and for any we define the coderivative at this point as
If is singlevalued, we write only instead of . If is singlevalued and smooth, then its coderivative amounts to the adjoint Jacobian.
A setvalued mapping has the Aubin property at if there exist a constant and neighborhoods of and of such that for all the following inclusion holds true
where is the unit ball. If is singlevalued, then the Aubin property reduces exactly to locally Lipschitzian property.
Example 3.1.
For a short illustration of the aforementioned objects, we use a simple example whose extension will be used later in the text. Consider . Since is convex, both normal cones coincide and we have
Fix now and compute
One can see that in this case limiting normal cone is strictly greater than Fréchet one. This means that is a nonregular point. Similarly, one can see that so is also .
4 Evaluation of a subgradient of and firstorder necessary optimality conditions for (8)
To solve problem (8) or equivalently (9) efficiently, we need to compute a subgradient information for the mapping . Unfortunately, we cannot expect that is a differentiable function and thus, we need first to compute some kind of generalized derivative of .
We will work with the generalized differential calculus of Mordukhovich [24, 32] and compute the limiting subdifferential of the objective in (9). To be able to do so, we first have to compute the socalled coderivative , which for continuously differentiable functions amounts to the adjoint Jacobian. First we state a lemma which links these two concepts together.
Lemma 4.1.
Consider the solution mapping being implicitly defined by system (8) and fix some . Assume that is Lipschitz continuous on some neighborhood of and that is continuously differentiable on some neighborhood of . Denoting , we have
To obtain the necessary optimality conditions in the form of original data, we need to compute . This is conducted in the next lemma which will also be the basis for proving the Lipschitzian continuity of later in Corollary 4.4.
Lemma 4.2.
Consider the setting of the solution mapping being implicitly defined by system (8) and fix some . Assuming (A1)–(A4), the upper estimate of is the collection of all quantities
(10) 
such that for the adjoint equations
(11a)  
(11b) 
with the terminal conditions and are fulfilled. For the multipliers we have the relations
(12a)  
(12b) 
where and and where, for and , we have defined
Proof.
Similarly to and , we define
We define the following partially linearized mapping
and show that it is single–valued and locally Lipschitz around . Indeed, the relations defining read for
with and . Since the first inclusion is solved for and the second one for , we obtain that is singlevalued due to (A1)–(A4). By virtue of [10, Corollary 3D.5] we further obtain that is Lipschitz continuous around , so that the system defining is strongly regular (in the sense of Robinson [31]) at .
If is a polyhedral set, then can be computed due to [9, Theorem 2] or [15, Proposition 3.2]. For the computation of , we will consider two cases of , specifically
(13a)  
(13b) 
where in (13b), the inequality is understood componentwise. The former case (13a) corresponds to , while the latter case (13b) corresponds to and . The former case is simple because from (12b) we immediately obtain that and . For the analysis of the more complicated case (13b) we recall the definition of and define its counterpart for a single time instant
This case is more involved than the previous one, because in one has to do with normal cones to moving sets whereby (components of) arise simultaneously both as the arguments as well as parameters specifying the movement of the constraint sets. Such a situation occurs typically in quasivariational inequalities and has been studied, e.g., in [25]. Unfortunately, the results of [25] cannot be directly applied here because the set
does not satisfy even the Mangasarian–Fromovitz constraint qualification.
As we have mentioned in the introduction, it would be possible to use standard calculus rules to obtain a formula for based on multiple computations of . The graph of can easily be visualized in Figure 1 and thus, can be computed by analyzing parts of separately, see definition of below. However, when computing on the basis of and the chain rule from [14, Theorem 4.1] one obtains only an upper estimate and not equality.
That is why we make use of [2] where formulas for Fréchet and limiting normal cones to a finite union of convex polyhedra have been derived and then applied to a special structure arising in time dependent problems. For simplicity, we will show the result only for the simplest case of . However, the generalization to a more–dimensional space is straightforward and can be conducted in componentwise way.
In the following text we will assume that . Define now the following sets
It is not difficult to show that . Moreover, forms the so–called normally admissible partition of as defined in [2]. Now, define the following index sets
where we assume that and all relations are required to hold for all . Further define
As shown in [2], we obtain that and that forms a normally admissible partition of . Now, we may state the result concerning the computation of , which replaces the computation of normal cone to a nonconvex sets by the computation of multiple normal cones to convex sets.
Proposition 4.3.
[2, Section 4] Fix any and denote by the index of the unique component such that . Then
where denotes the common value for any . For any and , this value can be computed as