1 Introduction

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 finite-element method in space and a semi-implicit 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 two-dimensional example. Other applications including frictional contacts or bulk damage, plasticity, or phase transformations are outlined.

Keywords: rate-independent systems, optimal control, identification, fractional-step 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 Biot-equation 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

 q=(u,z) (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 infinite-dimensional 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

 E(t,π,q) =E(t,π,u,z).

We then consider an optimal-control 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 Banach-space-valued 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 so-called 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 single-valued 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 first-order 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 application-motivated 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 two-dimensional example involving a spatial discretization by Finite-Element-Method (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 rate-independent 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 fractional-step-type semi-implicit time discretization of (3b,c). Moreover, if , or is infinite-dimensional, the time-discrete problem still remains infinite-dimensional, 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 finite-dimensional subsets , , and . Counting also a possible numerical approximation of , denoted by , altogether, (3) turns into the problem:

 Minimizeτ\large∑% Kk=1j(ukτh,zkτh)+H(π)subject to∂\large.uR1(π,uk−1τh,zk−1τh;ukτh−uk−1τhτ)+∂uEh(kτ,π,ukτh,zk−1τh)∋0,  u0τh=u0h,∂\large.zR2(π,uk−1τh,zk−1τh;zkτh−zk−1τhτ)+∂zEh(kτ,π,ukτh,zkτh)∋0,    z0τh=z0h,ukτh∈Uh  and  zkτh∈Zh  for k=1,...,K,  π∈Πh,⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (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 τR1(π,uk−1τh,zk−1τh;u−uk−1τhτ)+Eh(kτ,π,u,zk−1τh) subject to u∈Uh (5a) and, taking (one of) its solution for ukτh, further Minimize τR2(π,uk−1τh,zk−1τh;z−zk−1τhτ)+Eh(kτ,π,ukτh,z) subject to z∈Zh (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

 E(t,π,u,z)={E(t,π,u,z)% if u∈Λt0, z∈Kt0,∞otherwise, (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 adhesive-contact 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, phase-transformations, plasticity, etc.

###### Remark 2.1 (Stability and convergence for τ→0 and h→0).

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 rate-independent case when and is 1-homogeneous 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 rate-independent plasticity with hardening. Sometimes, viscosity can help for this uniqueness. This is the case of frictional normal-compliance contact of visco-elastic bodies which, after a certain algebraic manipulation gets the structure with separately uniformly quadratic with linear constraints in two-dimensions, cf. [34], or with conical constraints in three-dimensions. 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 well-known 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

 pkτh(π,\small˜% \normalsizeu,u,% \small˜\normalsizez):=∇uR1(π,\small˜\normalsizeu,\small˜\normalsizez,u−\small˜\normalsizeuτ)+∇uEh(kτ,π,u,\small˜% \normalsizez), (7a) qkτh(π,\small˜% \normalsizeu,u,% \small˜\normalsizez,z):=∇zR2(π,\small˜\normalsizeu,\small˜\normalsizez,z−\small˜\normalsizezτ)+∇zEh(kτ,π,u,z), (7b) \largeKk(~z):=(K1+~z)∩Kkτ0,   and (7c) J(π,^u,^z):=τ\large∑Kk=1j(uk,zk)+H(π)   with  ^u=(u1,...,uK)  and  ^z=(z1,...,zK). (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:

 MinimizeJ(π,uτh,zτh)   with  uτh:=(u1τh,...,uKτh)  and  zτh=(z1τh,...,zKτh)subject to0∈pkτh(π,uk−1τh,ukτh,zk−1τh)+NΛkτ(ukτh),      k=1,…,K,      u0τh=u0h,0∈qkτh(π,uk−1τh,ukτh,zk−1τh,zkτh)+N\normalsizeKk(zk−1τh)(zkτh), k=1,…,K, z0τh=z0h,π∈Πh⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (8)

with . Defining the solution map implicitly via system (8), we may used the so-called implicit programming approach to rewrite problem (8) equivalently into the form

 Minimize  J(π,Sτh(π))  subject to  π∈Πh. (9)

In the rest of the paper we will make use of the following standing assumption, which imply in particular the single-valuedness 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 finite-dimensional setting or in [24] and [7] for the general infinite-dimensional case.

All objects in this section are finite-dimensional. For sequence of sets we define the Painlewé-Kuratowski upper limit as

 Limsupk→∞Ak={x| ∃xk∈Ak, x is an accumulation point of {xk}}.

Using this construction, we define for any the Bouligand tangent cone, Fréchet normal cone and the limiting normal cone respectively as

 TA(¯x) ={v| ∃vk→v, λk↘0, ¯x+λkvk∈A}, ^NA(¯x) =(TA(¯x))∗={x∗| ⟨x∗,v⟩≤0 for all v∈TA(¯x)}, NA(¯x) =Limsupx\lx@stackrelA→¯x^NA(x),

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

 ∂f(¯x)={x∗| (x∗,−1)∈Nepif(¯x,f(¯x))}.

If happens to be convex, both normal cones coincide and are equal to the normal cone in the standard sense of convex analysis

 NA(¯x)={x∗| ⟨x∗,x−¯x⟩≤0 for all x∈A}

and similarly, if is continuously differentiable at , then .

For a set-valued mapping and for any we define the coderivative at this point as

 D∗M(¯x,¯y)(y∗)={x∗| (x∗,−y∗)∈NgphM(¯x,¯y)}

If is single-valued, we write only instead of . If is single-valued and smooth, then its coderivative amounts to the adjoint Jacobian.

A set-valued 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

 M(x)∩V⊂M(x′)+L|x−x′|B(0,1),

where is the unit ball. If is single-valued, then the Aubin property reduces exactly to locally Lipschitzian property.

###### Example 3.1.

For a short illustration of the afore-mentioned 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

 gphNC=gph^NC=({0}×R−)∪([0,1]×{0})∪({1}×R+).

Fix now and compute

 ^NgphNC(¯x,¯y) =R−×R+   and   NgphNC(¯x,¯y)

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 π↦J(π,S(π)) and first-order 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 so-called 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

 ∂~J(¯π)⊂∇πJ(¯π,¯u,¯z)+D∗S(¯π,¯u,¯z)(∇uJ(¯π,¯u,¯z),∇zJ(¯π,¯u,¯z)).
###### Proof.

It follows directly from [26, Theorem 7] and [32, Exercise 8.8]. ∎

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

 −K∑k=1(∇πpk)⊤βk−K∑k=1(∇πqk)⊤δk (10)

such that for the adjoint equations

 −u∗k =αk−(∇upk)⊤βk−(∇uqk)⊤δk−(∇~upk+1)⊤βk+1−(∇~uqk+1)⊤δk+1,  and (11a) −z∗k =γk−(∇zqk)⊤δk−(∇~zpk+1)⊤βk+1−(∇~zqk+1)⊤δk+1 (11b)

with the terminal conditions and are fulfilled. For the multipliers we have the relations

 (αkβk) ∈NgphNΛk(¯uk,−pk(¯π,¯uk−1,¯uk,¯zk−1)),    and (12a) (γδ) ∈NgphQ(¯z,−q(¯π,¯u,¯z)), (12b)

where and and where, for and , we have defined

 q(π,u,z):=⎛⎜⎝q1(π,u0,u1,z0,z1)…qK(π,uK−1,uK,zK−1,zK)⎞⎟⎠:RL×RKN×RKM→RKM,   %andQ(z):=\Large×Kk=1N\normalsizeKk(zk−1)(zk):RKM⇉RKM.
###### Proof.

Similarly to and , we define

 p(π,u,z):=⎛⎜⎝p1(π,u0,u1,z0)…pK(π,uK−1,uK,zK−1)⎞⎟⎠:RL×RKN×RKM→RKN,    andP(u):=\Large×Kk=1NΛk(uk):RKN⇉RKN

We define the following partially linearized mapping

 M(μ,ν):={(π,u,z)∣∣∣μ∈p(¯π,¯u,¯z)+∇up(¯π,¯u,¯z)(u−¯u)+∇zp(¯π,¯u,¯z)(z−¯z)+P(u)ν∈q(¯π,¯u,¯z)+∇uq(¯π,¯u,¯z)(u−¯u)+∇zq(¯π,¯u,¯z)(z−¯z)+Q(z)}

and show that it is single–valued and locally Lipschitz around . Indeed, the relations defining read for

 μk ∈pk(¯π,¯uk−1,¯uk,¯zk−1)+∇upk(¯π,¯uk−1,¯uk,¯zk−1)(uk−¯uk) +∇~upk(¯π,¯uk−1,¯uk,¯zk−1)(uk−1−¯uk−1)+∇~zpk(¯π,¯uk−1,¯uk,¯zk−1)(zk−1−¯zk−1)+NΛk(uk), νk ∈qk(¯π,¯uk−1,¯uk,¯zk−1,¯zk)+∇uqk(¯π,¯uk−1,¯uk,¯zk−1¯zk)(uk−¯uk) +∇~uqk(¯π,¯uk−1,¯uk,¯zk−1,¯zk)(uk−1−¯uk−1)+∇zqk(¯π,¯uk−1,¯uk,¯zk−1)(zk−¯zk) +∇~zqk(¯π,¯uk−1,¯uk,¯zk−1)(zk−1−¯zk−1)+N\normalsizeKk(zk−1)(zk)

with and . Since the first inclusion is solved for and the second one for , we obtain that is single-valued 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 .

This enables us to use [27, Proposition 3.2] and [32, Theorem 6.14] to obtain, with being the identity matrix, that

with some and satisfying

 (αβ)∈NgphP(¯u,−p(¯π,¯u,¯z))    and    (γδ)∈NgphQ(¯z,−q(¯π,¯u,¯z)).

Applying the product rule for normal cones [32, Proposition 6.41] we obtain the statement of the lemma. ∎

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

 \largeKk(zk−1) =RM    or (13a) \largeKk(zk−1) ={z∈RM|0≤z≤zk−1} (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

 Q(z) ={v∈RKM| vk∈N[0,zk−1](zk), k=1,…,K}, ~Q(~z) :={(z,v)∈RM×RM| v∈N[0,~z](z)}.

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

 \largeKk(0)={0}

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

 ~Q1 ={(~z,z,v)∈R3|~z∈(0,∞),z=~z,v∈(0,∞)}, ~Q2 ={(~z,z,v)∈R3|~z∈(0,∞),z=~z,v=0}, ~Q3 ={(~z,z,v)∈R3|~z∈(0,∞),z∈(0,~z),v=0}, ~Q4 =(0,∞)×{0}×{0}, ~Q5 =(0,∞)×{0}×(−∞,0), ~Q6 ={0}×{0}×(−∞,0), ~Q7 ={0}×{0}×{0}, ~Q8 ={0}×{0}×(0,∞).

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

 Θ =⎧⎪ ⎪⎨⎪ ⎪⎩(i1,…,iK)∣∣ ∣ ∣∣ik∈{1,…,8}ik−1∈{1,2,3}⟹ik∈{1,2,3,4,5}ik−1∈{4,5,6,7,8}⟹ik∈{6,7,8}⎫⎪ ⎪⎬⎪ ⎪⎭and I(s) =⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩(i1,…,iK)∈Θ∣∣ ∣ ∣ ∣∣ sk=1⟹ik=1,sk=5⟹ik=5sk=2⟹ik∈{1,2,3},sk=6⟹ik∈{5,6}sk=3⟹ik=3,sk=7⟹ik∈{1,…,8}sk=4⟹ik∈{3,4,5},sk=8⟹ik∈{1,8}⎫⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪⎭,

where we assume that and all relations are required to hold for all . Further define

 Qi:={(z,v)∈R2K∣∣(zk−1,zk,vk)∈~Qik, k=1,…,K}.

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

 NgphQ(¯z,¯v)=⋃s∈I(¯s)⋂i∈I(s)NclQi(Qs),

where denotes the common value for any . For any and , this value can be computed as

 N