René Henrion

A critical note on empirical (sample average, Monte Carlo) approximation of solutions to chance constrained programs

Keywords: chance constrained programming, probabilistic constraints, empirical approximation, sample average approximation, convergence of solution sets, sample size

The solution of chance constrained optimization problems by means of empirical approximation of the underlying multivariate distribution has recently become a popular alternative to conventional methods due to the efficient application of appropriate mixed integer programming techniques. As the complexity of required computations depends on the sample size used for approximation, exponential estimates for the precision of optimal solutions or optimal values have become a key argument for controlling the sample size. However, these exponential estimates may involve unknown constants such that the required sample size to approximate the solution of a problem may become arbitrarily large. We will illustrate this effect for Gaussian distributions.

Introduction

A chance constrained optimization problem has the general form min{g(x)|P(h(x, ξ) ≥ 0) ≥ p, x ∈ C},

where x ∈ R n is a decision vector, g : R n → R is an objective function, ξ is an s-dimensional random vector defined on some probability space (Ω, A, P), h : R n × R s → R m is a Borel measurable with respect to the second argument mapping, C ⊆ R n represents some abstract deterministic constraint and p ∈ [0, 1] is a fixed probability level. The random inequality system h(x, ξ) ≥ 0 may reflect some technological constraints in engineering problems which are affected by uncertainty. Since usually a decision x has to be taken before the uncertain parameter ξ is observed, it has become a standard approach of robust modeling to define x as feasible, whenever the probability of satisfying the random inequality is at least p. This is expressed in the so-called chance constraint P(h(x, ξ) ≥ 0) ≥ p. For a standard introduction to chance constrained programming we refer to the classical monograph [START_REF] Prékopa | Stochastic Programming[END_REF] and to the more recent treatise in [START_REF] Shapiro | Lectures on Stochastic Programming[END_REF].

Recent progress in mixed integer programming techniques tailored to chance constraints has led to the idea of solving (1) by empirical approximation (or sample average approximation) of the original random vector ξ (e.g., [START_REF] Luedtke | A sample approximation approach for optimization with probabilistic constraints[END_REF][START_REF] Pagnoncelli | Sample Average Approximation Method for Chance Constrained Programming: Theory and Applications[END_REF]). This means that an i.i.d. sample ξ 1 , . . . , ξ N of size N is drawn from the distribution of ξ and that the law P•ξ -1 of ξ in (1) is replaced by the empirical measure

N -1 N i=1 δ(ξ i ),
where δ(z) is the Dirac measure centered on z. Doing so, the original chance constraint P(h(x, ξ) ≥ 0) ≥ p turns into its empirical counterpart # i|h(x, ξ i ) ≥ 0 ≥ pN.

Of course this change of constraint leads to another optimization problem whose solution deviates from the solution of ( 1) and one has to answer two questions: do the approximating solutions converge with N → ∞ and if so, how large should N be chosen in order to guarantee a given precision of the solution obtained? The first question is answered by the stability theory of chance constrained programming mainly developed in [START_REF] Henrion | Metric regularity and quantitative stability in stochastic programs with probabilistic constraints[END_REF][START_REF] Henrion | Hölder and Lipschitz Stability of Solution Sets in Programs with Probabilistic Constraints[END_REF] which applies to arbitrary approximations of the original distribution and, in particular, to empirical ones. Under some explicitly verifiable conditions, not only qualitative convergence of approximating solutions can be guaranteed but also rates for this convergence can be derived. The latter allow us to obtain exponential bounds (in terms of sample size) for the precision of solutions in case of empirical approximations. However, one has to take into account that the exponential term involves apart from the sample size N also some other constants which may depend on the conditioning of the problem and may be hard to estimate. Thus, exponential estimates of solutions do not exclude the need for a large sample size even in small dimension in order to arrive at a reasonable precision of the solution. This situation occurs in particuar if the law of the original random vector ξ has unbounded support. We will illustrate and explain this effect for a multivariate Gaussian distribution (but similar observations could be made for other classes of multivariate distributions such as log-normal or t-). In order to keep the presentation as simple as possible we restrict ourselves to the simplest yet meaningful instance of problem (1):

min{c T x | P(ξ ≤ x) ≥ p}.
This means that we consider just linear objective functions, we forget about additional abstract deterministic constraints and we assume the chance constraint being in elementary separated form. Recalling the definition of the distribution function F ξ (x) := P(ξ ≤ x) of a random vector ξ, we may rewrite this problem as min{c

T x | F ξ (x) ≥ p} (P c,ξ,p ). (2)
In order to emphasize the dependence on the problem data c, ξ and p, we label problem (2) as (P c,ξ,p ). Before coming back to the issue of empirical distributions discussed above, we derive in the next section our main result demonstrating the difficulty of approximating a chance constrained program with unbounded support of the underlying distribution by means of distributions with bounded support. Note that, in particular, empirical measures have bounded support.

Main Result

We start by recalling the following well-known representation for partial derivatives of Gaussian distribution functions. We make use of the familiar notation ξ ∼ N (µ, Σ) to designate a Gaussian random vector with expectation µ and covariance matrix Σ.

Theorem 1 ([6], p. 204). Let ξ ∼ N (µ, Σ) with some positive definite covariance matrix Σ = (σ ij ) of order (s, s). Then, the distribution function F ξ is continuously differentiable at any z ∈ R s and

∂F ξ ∂z j (z) = f ξj (z j ) • F ξ(zj ) (z 1 , . . . , z j-1 , z j+1 . . . , z s ) (j = 1, . . . , s) .
Here, f ξj denotes the one-dimensional Gaussian density of the component ξ j , ξ(z j ) is an (s-1)-dimensional Gaussian random vector distributed according to ξ(z j ) ∼ N (μ, Σ), μ results from the vector µ + σ -1 jj (z jµ j ) σ j by deleting component j and Σ results from the matrix Σσ -1 jj σ j σ T j by deleting row j and column j, where σ j refers to column j of Σ.

Corollary 1. In the context of the previous Theorem, one has that ∂F ξ ∂z j (z) > 0 ∀z ∈ R s , ∀j ∈ {1, . . . , s} .

Proof. This follows immediately from the formula in Theorem 1 and the fact that both the density and the distribution function of a regular Gaussian distribution are strictly positive.

With each problem (P c,ξ,p ) in (2) we associate its (possibly empty) solution set

Ψ c,ξ,p := arg min{c T x | P(ξ ≤ x) ≥ p}. Lemma 1. For problem (P c,ξ,p ) in (2) assume that c i > 0 for i = 1, . . . , s. Then, Ψ c,ξ,p ⊆ [a, b] := {x ∈ R s | a i ≤ x i ≤ b i (i = 1, . . . , s)},
where, with 'supp' denoting the support of a random vector,

a i := inf{z i |z ∈ supp ξ} b i := sup{z i |z ∈ supp ξ} (i = 1, . . . , s).
Proof. Assume that there is some x * ∈ Ψ c,ξ,p and some i such that x * i > b i . Then, b i < ∞ and we may define x by

xi := (x * i + b i ) /2, xj := x * j (j = i) .
From x * being feasible for problem (P c,ξ,p ), we conclude that

p ≤ P(ξ ≤ x * ) = P(ξ ≤ x) + P(x i ≤ ξ i ≤ x * i , ξ j ≤ xj (j = i)). Now, since xi > b i , it follows that {x ∈ R n | xi ≤ x i ≤ x * i } ∩ supp ξ = ∅ , whence P(x i ≤ ξ i ≤ x * i , ξ j ≤ xj (j = i)) = 0
and P(ξ ≤ x) ≥ p. Therefore, x too is feasible for problem (P c,ξ,p ). On the other hand, c T x < c T x * due to c i > 0, xi < x * i and xj = x * j for j = i. This contradicts the assumption x * ∈ Ψ c,ξ,p . Consequently, x * ≤ b for any x * ∈ Ψ c,ξ,p . Similarly one shows that x * ≥ a for any

x * ∈ Ψ c,ξ,p . It follows that Ψ c,ξ,p ⊆ [a, b].
Now, we are in a position to state our main result:

Theorem 2. Let s > 1.
Assume that ξ has a regular normal distribution according to ξ ∼ N (µ, Σ) and that η is a random vector with compact support. Then, for any p ∈ (0, 1) there exists a sequence c

(n) ∈ R n with c (n) i > 0 for i = 1, . . . , s such that Ψ c (n) ,ξ,p = ∅ and inf{ x -y | x ∈ Ψ c (n) ,ξ,p , y ∈ Ψ c (n) ,η,p } > n ∀n ∈ N. (3)
Proof. Fix an arbitrary p ∈ (0, 1) and an arbitrary n ∈ N. Since supp η is compact, we may apply Lemma 1 to η in order to derive the existence of some compact(!) rectangle [a, b] such that

Ψ c,η,p ⊆ [a, b] ∀c ∈ R s : c i > 0 (i = 1, . . . , s). (4)
With [a, b] being compact, we may choose L n > 0 such that

y -z ≥ n ∀y ∈ [a, b], ∀z : z ≥ L n . (5)
Since s > 1 by assumption, we may define the ratio

κ(z) := ∂F ξ ∂z 1 (z) ∂F ξ ∂z 2 (z) (z ∈ R s ) . (6)
Note that the partial derivatives of F ξ are continuous (see Theorem 1) and strictly positive (see Corollary 1), hence κ is correctly defined and continuous. Consequently the quantity

κ := sup {κ(z)|z ∈ B (0, L n )} (7)
is finite. Next, let q p be the p-quantile of the first marginal distribution of ξ, i.e., of the distribution of the first component ξ 1 . Since ξ 1 has a one-dimensional normal distribution and p ∈ (0, 1), q p ∈ R is uniquely defined by P (ξ 1 ≤ q p ) = p. We claim that for each k ∈ N there exists some t k ∈ R such that

F ξ q p + k -1 , t k , . . . , t k = p. (8)
Indeed, for arbitrarily fixed k one has that

lim τ →-∞ F ξ q p + k -1 , τ, . . . , τ = 0
as a general property of distribution functions and that

p = P (ξ 1 ≤ q p ) < P ξ 1 ≤ q p + k -1 = lim τ →∞ F ξ q p + k -1 , τ, . . . , τ .
Now, the existence of t k with the desired property (8) follows from continuity of F ξ and from p > 0. Next we claim that t k → k ∞. If there existed some subsequence t k l and some r ∈ R such that t k l ≤ r for all l ∈ N, then

F ξ q p + k -1 l , r, . . . , r ≥ F ξ q p + k -1 l , t k l , . . . , t k l = p ∀l ∈ N
which again by continuity of F ξ as a function of each of its components yields the contradiction

p ≤ F ξ (q p , r, . . . , r) = P (ξ 1 ≤ q p , ξ 2 ≤ r, . . . , ξ s ≤ r) < P (ξ 1 ≤ q p ) = p.
Here, the strict inequality relies on the fact that a regular Gaussian distribution has a density which is strictly positive everywhere. Hence, we have shown that

t k → k ∞. Consider the sequence z (k) := q p + k -1 , t k , . . . , t k (k ∈ N) .
Then, by [START_REF] Talagrand | Sharper bounds for Gaussian and empirical processes[END_REF],

F ξ z (k) = p (k ∈ N) . (9)
Moreover, Theorem 1 yields that

∂F ξ ∂z 1 (z (k) ) = f ξ1 (q p + k -1 ) • F ξ z (k) 1 (t k , . . . , t k ) , (10)
where f ξ1 denotes the one-dimensional Gaussian density of the component ξ 1 and ξ z

(k) 1
is an (s-1)-dimensional Gaussian random vector distributed according to ξ z

(k) 1
∼ N (μ (k) , Σ) where μ(k) and Σ result from the original parameters µ and Σ, respectively, of ξ as detailed in Theorem 1. In particular,

μ(k) = (µ 2 , . . . , µ s ) + σ -1 11 q p + k -1 -µ 1 (σ 21 , . . . , σ s1 )
and we observe that μ(k) → k μ := (µ 2 , . . . , µ s ) + σ -1 11 (q pµ 1 ) (σ 21 , . . . , σ s1 ) .

Note also that in contrast to μ(k) , the covariance matrix Σ does not depend on the index k. Now we define the centered random vector

ξ := ξ z (k) 1 -μ(k) ∼ N (0, Σ),
whose distribution does no longer depend on the index k. Exploiting the relation

F ξ z (k) 1 (t k , . . . , t k ) = F ξ t k - μ(k) 1 , . . . , t k - μ(k) s-1
and noting that all components of the argument

t k - μ(k) 1 , . . . , t k - μ(k)
s-1 tend to infinity due to (11) and t k → k ∞, we conclude that

F ξ t k - μ(k) 1 , . . . , t k - μ(k) s-1 → k 1
because the values of the (fixed) distribution function F ξ tend to one if all its components tend to infinity. This implies

F ξ z (k) 1 (t k , . . . , t k ) → k 1, whence (10) leads to ∂F ξ ∂z 1 (z (k) ) → k f ξ1 (q p ) > 0 (12)
by continuity and positivity of the density f ξ1 . Similarly, the second partial derivative of F ξ calculates from Theorem 1 as

∂F ξ ∂z 2 (z (k) ) = f ξ2 (t k ) • F ξ z (k) 2 q p + k -1 , t k , . . . , t k , (13)
where f ξ2 denotes the one-dimensional Gaussian density of the component ξ 2 and ξ z

(k) 2
is a certain (s-1)-dimensional Gaussian random vector. From f ξ2 (t k ) → k 0 (due to t k → k ∞) and from the fact that distribution functions are bounded between zero and one, we infer that

∂F ξ ∂z 2 (z (k) ) → k 0,
which along with (12) and ( 6) provides that κ(z (k) ) → k ∞. Therefore, with our arbitrarily fixed number n ∈ N we may associate an index k n ∈ N such that κ(z (kn) ) > κ where κ is defined in [START_REF] Shapiro | Lectures on Stochastic Programming[END_REF]. Now, we assign to n the cost vector c (n) := ∇F ξ z (kn) for the linear objective function in problem (P c (n) ,ξ,p ) in [START_REF] Henrion | Hölder and Lipschitz Stability of Solution Sets in Programs with Probabilistic Constraints[END_REF].

Then, by Corollary 1, we have that c

(n) i > 0 for i = 1, . . . , s as required in the statement of our theorem. Knowing that log F ξ is a concave function (see [START_REF] Prékopa | Stochastic Programming[END_REF]), the problem (P c (n) ,ξ,p ) in ( 2) may be written equivalently as a convex optimization problem min{c (n)T x |log F ξ (x) ≤log p}.

(P c (n) ,ξ,p ) (14)
With c (n) = 0, a solution x * to this problem is equivalently characterized by the conditions log F ξ (x * ) =log p and c (n) + λ∇ (-log F ξ ) (x * ) = 0 for some λ > 0.

Simplifying these leads to the equivalent conditions

F ξ (x * ) = p and c (n) = λ∇ (F ξ ) (x * ) for some λ > 0. ( 15
)
Now, since c (n) = ∇F ξ z (kn) and F ξ (z (kn) ) = p by (9), we conclude that z (kn) is a solution to (P c (n) ,ξ,p ). This shows that Ψ c (n) ,ξ,p = ∅ as asserted in our theorem. Finally, we show that Ψ c (n) ,ξ,p ∩ B (0, L n ) = ∅ with L n defined in [START_REF] Pagnoncelli | Sample Average Approximation Method for Chance Constrained Programming: Theory and Applications[END_REF]. Assume the contrary and choose some x * ∈ Ψ c (n) ,ξ,p with x * ≤ L n . From ( 6) and (15) we derive

κ (x * ) = ∂F ξ ∂z 1 (x * ) ∂F ξ ∂z 2 (x * ) = c (n) 1 c (n) 2 = ∂F ξ ∂z 1 (z (kn) ) ∂F ξ ∂z 2 (z (kn) ) = κ(z (kn) ) > κ.
which is a contradiction with [START_REF] Shapiro | Lectures on Stochastic Programming[END_REF]. Consequently, 4), it follows from (5) that xy ≥ n. Since x and y were arbitrarily chosen, we end up at the final assertion (3) of our theorem.

Ψ c (n) ,ξ,p ∩ B (0, L n ) = ∅. Now, select arbitrary x ∈ Ψ c (n) ,ξ,p and y ∈ Ψ c (n) ,η,p . Then, x > L n . Since also Ψ c (n) ,η,p ⊆ [a, b] by (
Theorem 2 can be interpreted as follows in the context of empirical approximation upon observing that the support of empirical measures is finite, hence compact: no matter how large the sample size N for the empirical approximation of the original random vector ξ is chosen, there is always an instance of problem (2) (by choosing an appropriate cost vector c) such that the solutions between the original problem and its empirical approximation are arbitrarily far from each other. Note that relation (3) implies (and actually is much stronger than) the Hausdorff distance between both solution sets being larger than any prescribed n. Moreover, this effect of ill-conditioning is not caused by letting the probability level tend to one, because the result of the theorem holds true for any fixed p. In the following section we look at the same phenomenon from a slightly different viewpoint.

Exponential estimates with ill-conditioned constants

The recent literature on empirical or sample average approximation on chance constraints [START_REF] Luedtke | A sample approximation approach for optimization with probabilistic constraints[END_REF][START_REF] Pagnoncelli | Sample Average Approximation Method for Chance Constrained Programming: Theory and Applications[END_REF] compiles several convergence results for feasible sets, optimal values and solutions, most of them of qualitative nature (continuity, upper semicontinuity), some of them providing exponential estimates (convergence of feasible sets, lower bounds for optimal values). For the general stability theory of chance constrained programming with arbitrary approximations (not just empirical ones) and quantitative convergence results even for solution sets, we refer to [START_REF] Henrion | Metric regularity and quantitative stability in stochastic programs with probabilistic constraints[END_REF][START_REF] Henrion | Hölder and Lipschitz Stability of Solution Sets in Programs with Probabilistic Constraints[END_REF]. In this section we will appy these results to the special case of empirical approximations in order to obtain an exponential bound for the convergence (in the sense of Hausdorff distance!) of solution sets. Despite this positive result we will show then, that the existence of exponential estimates does not exclude the need for a possibly excessive sample size in the empirical approximation. We start by citing the following stability result for chance constraints whch we present here in a simplified form sufficient for our purposes Theorem 3 ([2], Corollary 3). In problem (P c,ξ,p ) in (2) let p ∈ (0, 1) and the following assumptions be satisfied:

1. log F ξ is a strongly concave function.

2. Ψ c,ξ,p is nonempty and compact.

Then, there exist L, δ > 0 such that

d H (Ψ c,ξ,p , Ψ c,η,p ) ≤ L sup z∈R s |F ξ (z) -F η (z)| ∀η : sup z∈R s |F ξ (z) -F η (z)| < δ.
(16) Here, d H refers to the Hausdorff distance.

A prototype example for a problem (2) which automatically satisfies all assumptions of Theorem 3 is given by a random vector ξ having a standard Gaussian distribution. As a preparation we show the following property which is of independent interest: Proposition 1. In problem (P c,ξ,p ) in (2) let p ∈ (0, 1), c i > 0 for i = 1, . . . , s and ξ ∼ N (0, I s ). Then, the problem has a solution.

Proof. Referring back to the proof of Theorem 2, a solution of problem (P c,ξ,p ) is equivalently characterized by the conditions (15) applied to c rather than c (n) . The distribution assumption implies that for all z ∈ R s and all i = 1, . . . , s,

F ξ (z) = Φ (z 1 ) • • • Φ (z s ) ; ∂F ξ ∂z i (z) = f (z i ) Φ (z 1 ) • • • Φ (z i-1 ) Φ (z i+1 ) • • • Φ (z s ) ,
where f and Φ denote the one-dimensional standard normal density and distribution function, respectively. From (15) we derive that x * is a solution to (P c,ξ,p ) if there exists some λ > 0 such that

Φ(x * 1 ) • • • Φ(x * s ) = p; λ f (x * i ) Φ(x * i ) = c i (i = 1, . . . , s) . (17)
(recall that Φ is strictly positive). Defining α := f /Φ, we observe that due to c i > 0 the second relation of (17) amounts to the fact that α(x * i )/c i is constant for all i. We may assume, without loss of generality, that c 1 is the largest of the coefficients c i . Then, x * is a solution to (P c,ξ,p ) if there exist coefficients ρ 1 , . . . , ρ s-1 ≥ 1 such that

Φ(x * 1 ) • • • Φ(x * s ) = p; α(x * 1 ) = ρ 1 α(x * 2 ) = • • • = ρ s-1 α(x * s ). ( 18
)
We recall from the properties of the one-dimensional standard normal density and distribution function that α(t) → 0 for t → ∞. Consequently, given any t ∈ R, there exist values β i (t) for i = 1, . . . , s -1 such that β i (t) ≥ t and

α(t) = ρ 1 α(β 1 (t)) = • • • = ρ s-1 α(β s-1 (t)). ( 19
)
Taking into account that lim t→-∞ Φ (t) = 0 and Φ (β i (t)) ≤ 1 on the one hand and lim t→∞

Φ (t) = lim t→∞ Φ (β i (t)) = 1
due to β i (t) ≥ t, on the other hand, we conclude that lim t→-∞(+∞)

Φ (t) Φ (β 1 (t)) • • • Φ (β s-1 (t)) = 0(1).
For continuity reasons, there exists some t * such that

Φ (t * ) Φ (β 1 (t * )) • • • Φ (β s-1 (t * )) = p.
Setting x * 1 := t * and x * i := β i-1 (t * ) for i = 2, . . . , s, one verifies via (19) that ( 18) is satisfied and, hence, x * is a solution to (P c,ξ,p ).

Corollary 2. Under the assumptions of Proposition 1 the estimate (16) holds true.

Proof. We have to check that the assumptions of Theorem 3 are satisfied. The strong concavity of the log of Gaussian distribution functions with independent components is easy to verify (see [START_REF] Henrion | Hölder and Lipschitz Stability of Solution Sets in Programs with Probabilistic Constraints[END_REF]Prop. 14]). As shown in Proposition 1, the solution set Ψ c,ξ,p is nonempty. On the other hand, there may not exist more than one solution to problem (P c,ξ,p ) because in its equivalent description min{c T x |log F ξ (x) ≤log p} (P c,ξ,p ) the inequality constraint is strongly convex according to what we have mentioned in the begining of this proof.

We emphasize that in the result of Theorem 3 the approximating random vector η can be arbitrary. In the special case that η is an empirical approximation, one may exploit exponential bounds from empirical process theory (e.g., [START_REF] Talagrand | Sharper bounds for Gaussian and empirical processes[END_REF]) to further interpret the obtained stability result. In order to keep the presentation simple, we refer here to a classical inequality by Kiefer [START_REF] Kiefer | On large deviations of the empiric d.f. of vector chance variables and a law of the iterated logarithm[END_REF] stating in any dimension s the existence of constants k 1 and k 2 < 2 (where k 2 may be chosen arbitrarily close to 2) such that for all ε > 0 and all η N having an empirical distribution of an i.i.d. sample of ξ with size N the following estimate applies:

P sup z∈R s |F ξ (z) -F η N (z)| ≥ ε ≤ k 1 exp -k 2 ε2 N . (20)
Since by ( 16) one has for all ε > 0 the implication sup

z∈R s |F ξ (z) -F η (z)| < min{δ, (ε/L) 2 } =⇒ d H (Ψ c,ξ,p , Ψ c,η,p ) < ε, it follows from (20) with ε := min{δ, (ε/L) 2 } that P d H Ψ c,ξ,p , Ψ c,η N ,p ≥ ε ≤ k 1 exp -k 2 min{δ, (ε/L) 2 } 2 N . (21)
This last relation establishes an exponential bound for the convergence of Hausdorff distance between the solution sets of the original problem and the problem approximated by a sample of size N . Obviously, the quantity

k 2 min{δ, (ε/L) 2 } 2
determines the exponential decay of the required sample size. However, for practical use, one would have to know the values or at least estimates for δ and L which is difficult or impossible in general. Then, the availability of an exponential convergence result does not exclude excessive sample sizes even in order to give a sense to (21), i.e., to ensure that the right-hand side is smaller than one as an upper probability estimate. Revisiting Theorem 3, one observes that the couple (δ, L) in ( 16) is not uniquely determined. Therefore, let us define the best possible coefficient of exponential decay by ϑ (c, ξ, p, ε) := sup k 2 min{δ, (ε/L) 2 } 2 (δ, L) satisfy ( 16) for (P c,ξ,p ) .

Then, (21) can be formally improved to

P d H Ψ c,ξ,p , Ψ c,η N ,p ≥ ε ≤ k 1 exp (-ϑ (c, ξ, p, ε) • N ) .
The following result demonstrates that for a given significant problem class, this coefficient of exponential decay may be arbitrarily close to zero, thus driving the required sample size to infinity. From (20) we infer that, for any N ,

P sup z∈R s |F ξ (z) -F η N (z)| < δ/2 ≥ 1 -k 1 exp -k 2 δ2 N/4 ,
where η N has an empirical distribution of an i.i.d. sample of ξ with size N . For N → ∞, the right-hand side tends to one such that the probability on the lefthand side is at least strictly positive. As a consequence, for some fixed N large enough, there exists a discrete random vector η N with N atoms such that sup

Theorem 4 . 2 ≥ τ / 2 .

 422 In[START_REF] Henrion | Hölder and Lipschitz Stability of Solution Sets in Programs with Probabilistic Constraints[END_REF], let s > 1, ξ ∼ N (0, I s ) and p ∈ (0, 1) be arbitrarily given. Then, for any ε > 0 one has that inf {ϑ (c, ξ, p, ε) |c ≻ 0} = 0, where 'c ≻ 0' means c i > 0 for i = 1, . . . , s.Proof. Denote by τ the infimum above and assume that τ > 0. Then, ϑ (c, ξ, p, ε) ≥ τ for all c ≻ 0. By definition of ϑ (c, ξ, p, ε), we infer that ∀c ≻ 0 ∃ (δ, L) : (δ, L) satisfy (16) for (P c,ξ,p ) and k 2 min{δ, (ε/L) 2 } and L do not depend on c. Consequently, we have shown that there exist δ > 0 and L such that ∀c ≻ 0 ∃δ ≥ δ, L ≤ L : (16) holds true for (P c,ξ,p ). This statement can be evidently reduced to:∀c ≻ 0 : d H (Ψ c,ξ,p , Ψ c,η,p ) ≤ L sup z∈R s |F ξ (z) -F η (z)| (22)∀η : supz∈R s |F ξ (z) -F η (z)| < δ.

  z∈R s |F ξ (z) -F η N (z)| < δ/2. Then, (22) implies that ∀c ≻ 0 : d H Ψ c,ξ,p , Ψ c,η N ,p ≤ L δ/2. (23) Now, since supp η N is compact, Theorem 2 yields the existence of some c ≻ 0 such that inf{ xy | x ∈ Ψ c,ξ,p , y ∈ Ψ c,η N ,p } > L δ/2.This, however, is a contradiction with (23) becaused H (A, B) ≥ inf{ xy | x ∈ A, y ∈ B}for any closed sets A, B. Hence, τ = 0, as was to be shown.

⋆ This work was supported by the DFG Research Center Matheon "Mathematics for key technologies" in Berlin

The effect of the previous theorem can already be illustrated in two dimensions: Example 1. Consider the following 2-dimensional problem:

By independence of components, we may rewrite the chance constraint as

where Φ denotes the one-dimensional standard normal distribution function.

Referring to the optimality conditions as in (15), the solution of this problem is equivalently characterized by the following three nonlinear equations in the three variables x 1 , x 2 , λ:

Here, ϕ is the density of the one-dimensional standard normal distribution. This system is easily solved numerically, providing a unique optimal solution x * 1 = 2.33, x * 2 = 4.88. Now, suppose that we want to approximate this solution empirically up to a precision of ε = 0.1. Then, in particular, the second component of the solution to the problem with empirical approximation has to exceed the value 4.78. Now, the second components ξ

For instance, for N = 10 6 and t = 4.78, one has Φ N (t) ≈ 0.42. This means that the probability of obtaining a one-digit precise solution by empirical approximation with a huge sample size like one million is less than 1 -0.42 = 0.58.

Certainly, the effect of the example relies on the highly unbalanced cost vector c = 1, 10 -4 . Making it more reasonably balanced like c = (1, 0.1), one would still need a sample size of N ≈ 6.300 for estimating the solution of the problem with a precision of 0.1 at a reasonably high probability of 0.99. Taking into account that N corresponds to the number of binary variables required in the discrete optimization problem, this is already a considerable quantity given the trivial dimension s = 2 of the problem. Of course, things may be expected to become much worse in larger dimension.


Barbara Kaltenbacher
email: barbara.kaltenbacher@aau.at

Convergence rates for the iteratively regularized Landweber iteration in Banach space

Keywords: regularization, nonlinear inverse problems, Banach space, Landweber iteration

In this paper we provide a convergence rates result for a modified version of Landweber iteration with a priori regularization parameter choice in a Banach space setting.

An increasing number of inverse problems is nowadays posed in a Banach space rather than a Hilbert space setting, cf., e.g., [START_REF] Burger | Convergence rates of convex variational regularization[END_REF][START_REF] Hofmann | A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators[END_REF][START_REF] Schuster | Regularization Methods in Banach Spaces[END_REF] and the references therein.

An Example of a model problem, where the use of non-Hilbert Banach spaces is useful, is the identification of the space-dependent coefficient function c in the elliptic boundary value problem

-∆u + cu = f in Ω (1)
u = 0 on ∂Ω (2)
from measurements of u in Ω ⊆ R d , d ∈ {1, 2, 3}, where f is assumed to be known. Here e.g., the choices p = 1 for recovering sparse solutions, q = ∞ for modelling uniformly bounded noise, or q = 1 for dealing with impulsive noise are particulary promising, see, e.g., [START_REF] Clason | A semi-smooth Newton method for nonlinear parameter identification problems with impulsive noise[END_REF] and the numerical experiments in Section 7.3.3 of [START_REF] Schuster | Regularization Methods in Banach Spaces[END_REF].

Motivated by this fact we consider nonlinear ill-posed operator equations

F (x) = y (3)
where F maps between Banach spaces X and Y . In the example above, the forward operator F maps the coefficient function c to the solution of the boundary value problem (1), [START_REF] Burger | Convergence rates of convex variational regularization[END_REF], and is well-defined as an operator F : D(F ) ⊆ L p (Ω) → L q (Ω) ,

⋆ Support by the German Science Foundation DFG under grant KA 1778/5-1 and within the Cluster of Excellence in Simulation Technology (EXC 310/1) at the University of Stuttgart is gratefully acknowledged

where

D(F ) = {c ∈ X | ∃ ĉ ∈ L ∞ (Ω)
, ĉ ≥ 0 a.e. : c -ĉ X ≤ r}, r sufficiently small, for any

p, q ∈ [1, ∞], f ∈ L 1 (Ω) if d ∈ {1, 2} p ∈ [1, ∞], q ∈ ( d 2 , ∞], f ∈ L s (Ω), s > d 2 if d ≥ 3 ,
see Section 1.3 in [START_REF] Schuster | Regularization Methods in Banach Spaces[END_REF].

Since the given data y δ are typically contaminated by noise, regularization has to be applied. We are going to assume that the noise level δ in

y -y δ ≤ δ (4)
is known and provide convergence results in the sense of regularization methods, i.e., as δ tends to zero. In the following, x 0 is some initial guess and we will assume that a solution x † to (3) exists. Variational methods in Banach space have been extensively studied in the literature, see, e.g., [START_REF] Burger | Convergence rates of convex variational regularization[END_REF][START_REF] Neubauer | Improved and extended results for enhanced convergence rates of Tikhonov regularization in Banach spaces[END_REF][START_REF] Hofmann | A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators[END_REF] and the references therein.

Since these generalizations of Tikhonov regularization require computation of a global minimizer, iterative methods are an attractive alternative especially for large scale problems. After convergence results on iterative methods for nonlinear ill-posed operator equations in Banach spaces had already been obtained in the 1990's (cf. the references in [START_REF] Bakushinsky | Iterative methods for approximate solution of inverse problems[END_REF]) in the special case X = Y , the general case X = Y has only been treaten quite recently, see e.g. [START_REF] Hein | Accelerated Landweber iteration in Banach spaces[END_REF], [START_REF] Kaltenbacher | Convergence rates for the iteratively regularized Gauss-Newton method in Banach spaces[END_REF], and [START_REF] Kaltenbacher | Convergence of some iterative methods for the regularization of nonlinear ill-posed problems in Banach spaces[END_REF] for an analysis of gradient and Newton type iterations. While convergence rates have already been established for the iteratively regularized Gauss-Newton iteration in [START_REF] Kaltenbacher | Convergence rates for the iteratively regularized Gauss-Newton method in Banach spaces[END_REF], the question of convergence rates is still open for gradient type, i.e. Landweber methods. It is the aim of this paper to provide such a result.

In order to formulate and later on analyze the method, we have to introduce some basic notations and concepts.

Consider, for some q ∈ (1, ∞), the duality mapping J X q (x) := ∂ 1 q x q , which maps from X to its dual X * . To analyze convergence rates we employ the Bregman distance

D jq (x, x) = 1 q x q - 1 q x q -j X q (x), x -x X * .X
(where j X q (x) denotes a single valued selection of J X q (x)) or its shifted version

D x0 q (x, x) := D jq (x -x 0 , x -x 0 ) .
Throughout this paper we will assume that X is smooth, which means that the duality mapping is single-valued, and moreover, that X is q-convex, i.e., D jq (x, y) ≥ c q xy q (5)

for some constant c q > 0. As a consequence, X is reflexive and we also have

D j q * (x * , y * ) ≤ C q * x * -y * q * , (6)
for some C q * . Here q * denotes the dual index q * = q q-1 . Moreover, the duality mapping is bijective and J -1 q = J X * q * , the latter denoting the (by q-convexity also single-valued) duality mapping on X * . We will also make use of the identities

D jq (x, y) = D jq (x, z) + D jq (z, y) + J X q (z) -J X q (y), x -z X * ,X (7)
and

D jq (y, x) = D j q * (J X q (x), J X q (y)) . (8)
For more details on the geometry of Banach spaces we refer, e.g., to [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF] and the references therein.

We here consider the iteratively regularized Landweber iteration

J X q (x δ n+1 -x 0 ) = (1 -α n )J X q (x δ n -x 0 ) -µ n A * n j Y p (F (x δ n ) -y δ ) , (9)
x δ n+1 = x 0 + J X * q * (J X q (x δ n+1 -x 0 )) , n = 0, 1, . . .
where we abbreviate

A n = F ′ (x δ n )
, which, for an appropriate choice of the sequence {α n } n∈IN ∈ (0, 1], has been shown to be convergent with rates under a source condition

x † -x 0 ∈ R(F ′ (x † ) * F ′ (x † )) ν/2 ), (10)
with ν = 1 in a Hilbert space setting in [START_REF] Scherzer | A modified Landweber iteration for solving parameter estimation problems[END_REF]. Since the linearized forward operator F ′ (x) typically has some smoothing property (reflecting the ill-posedness of the inverse problems) condition [START_REF] Neubauer | Improved and extended results for enhanced convergence rates of Tikhonov regularization in Banach spaces[END_REF] can often be interpreted as a regularity assumption on the initial error x †x 0 , which is stronger for larger ν.

In the Hilbert space case the proof of convergence rates for the plain Landweber iteration (i.e., [START_REF] Kaltenbacher | Convergence of some iterative methods for the regularization of nonlinear ill-posed problems in Banach spaces[END_REF] with α n = 0) under source conditions [START_REF] Neubauer | Improved and extended results for enhanced convergence rates of Tikhonov regularization in Banach spaces[END_REF] relies on the fact that the iteration errors x δ nx † remain in the range of (F ′ (x † ) * F ′ (x † )) ν/2 and their preimages under (F ′ (x † ) * F ′ (x † )) ν/2 form a bounded sequence (cf., Proposition 2.11 in [START_REF] Kaltenbacher | Iterative Regularization Methods for Nonlinear Ill-posed Problems[END_REF]). Since carrying over this approach to the Banach space setting would require more restrictive assumptions on the structure of the spaces even in the special case ν = 1, we here consider the modified version with an appropriate choice of {α n } n∈IN ∈ (0, 1].

In place of the Hilbert space source condition [START_REF] Neubauer | Improved and extended results for enhanced convergence rates of Tikhonov regularization in Banach spaces[END_REF], we consider variational inequalities

∃ β > 0 ∀x ∈ B D ρ (x † ) : | J X q (x † -x 0 ), x -x † X * ×X | ≤ βD x0 q (x † , x) 1-ν 2 F ′ (x † )(x -x † ) ν , (11)
cf., e.g., [START_REF] Hein | Approximate source conditions for nonlinear ill-posed problems -chances and limitations[END_REF], where

B D ρ (x † ) = {x ∈ X | D x0 q (x † , x) ≤ ρ 2 } with ρ > 0 such that x 0 ∈ B D ρ (x † ).
Using the interpolation and the Cauchy-Schwarz inequality, it is readily checked that in the Hilbert space case [START_REF] Neubauer | Improved and extended results for enhanced convergence rates of Tikhonov regularization in Banach spaces[END_REF] implies [START_REF] Scherzer | A modified Landweber iteration for solving parameter estimation problems[END_REF]. For more details on such variational inequalities we refer to Section 3.2.3 in [START_REF] Schuster | Regularization Methods in Banach Spaces[END_REF] and the references therein.

The assumptions on the forward operator besides a condition on the domain

B D ρ (x † ) ⊆ D(F ) (12)
include a structural condition on its degree of nonlinearity (cf. [START_REF] Hein | Approximate source conditions for nonlinear ill-posed problems -chances and limitations[END_REF] (

F ′ (x † + v) -F ′ (x † ))v ≤ K F ′ (x † )v c1 D x0 q (x † , v + x † ) c2 , v ∈ X, x † + v ∈ B D ρ (x † ) , (13)
whose strength depends on the smoothness index in [START_REF] Scherzer | A modified Landweber iteration for solving parameter estimation problems[END_REF]. Namely, we assume that

c 1 = 1 or c 1 + c 2 p > 1 or (c 1 + c 2 p ≥ 1 and K is sufficiently small) ( 14
)
c 1 + c 2 2ν ν+1 ≥ 1 , (15)
so that in case ν = 1, a Lipschitz condition on F ′ , corresponding to (c 1 , c 2 ) = (0, 1) is sufficient.

Here F ′ denotes the Gateaux derivative of F , hence a Taylor remainder estimate

F (x δ n ) -F (x † ) -F ′ (x † )(x δ n -x † ) (16) = g(1) -g(0) -F ′ (x † )(x δ n -x † ) = 1 0 g ′ (t) dt -F ′ (x † )(x δ n -x † ) = 1 0 F ′ (x † + t(x δ n -x † ))(x δ n -x † ) dt -F ′ (x † )(x δ n -x † ) ≤ K F ′ (x † )(x δ n -x † ) c1 D x0 q (x † , x δ n ) c2 (17)
where g : t → F (x † + t(x δ nx † )), follows from (13). We will assume that in each step the step size µ n > 0 in ( 9) is chosen such that

µ n 1 -3C(c 1 )K 3(1 -C(c 1 )K) F (x δ n ) -y δ p -2 q * +q-2 C q * µ q * n A * n j Y p (F (x δ n ) -y δ ) q * ≥ 0 (18) where C(c 1 ) = c c1 1 (1 -c 1 ) 1-c1
, and c 1 , K are as in [START_REF] Schuster | Regularization Methods in Banach Spaces[END_REF], which is possible, e.g., by a choice 0

< µ n ≤ C µ F (x δ n )-y δ q-p q-1 An q * =: µ n with C µ := 2 2-q * -q 3 1-3C(c1)K (1-C(c1)K)C q * If p ≥ q (19)
and F , F ′ are bounded on B D ρ (x † ), it is possible to bound µ n away from zero

µ n ≥ C µ sup x∈B D ρ (x † ) ( F (x) -y + δ) p-q F ′ (x) q -1/(q-1) =: µ (20)
for δ ∈ [0, δ], provided the iterates remain in B D ρ (x † ) (which we will show by induction in the proof of Theorem 1). Hence, there exist µ, µ > 0 independent of n and δ such that we can choose

0 < µ ≤ µ n ≤ µ , (21)
(e.g., by simply setting µ n ≡ µ). Moreover, we will use an a priori choice of the stopping index n * according to

n * (δ) = min{n ∈ IN : α ν+1 p(ν+1)-2ν n ≤ τ δ} , (22)
and of {α n } n∈IN such that

α n+1 α n 2ν p(ν+1)-2ν + 1 3 α n -1 ≥ cα n (23)
for some c ∈ (0, 1 3 ) independent of n, where ν ∈ [0, 1] is the exponent in the variational inequality [START_REF] Scherzer | A modified Landweber iteration for solving parameter estimation problems[END_REF].

Remark 1. A possible choice of {α n } n∈IN satisfying (23) and smallness of α max is given by α n = α 0 (n + 1) x with x ∈ (0, 1] such that 3xθ < α 0 sufficiently small, since then with c :=

α n+1 α n θ + ( 1 3 -c)α n -1 = α n α 0 α 0 ( 1 3 -c) - (n + 2) xθ -(n + 1) xθ (n + 2) xθ (n + 1) x = α n α 0 α 0 ( 1 3 -c) - xθ(n + 1 + t) xθ-1 (n + 2) xθ (n + 1) x ≥ α n α 0 α 0 ( 1 3 -c) -xθ (n + 1) x n + 1 + t ≥ 0 , for some t ∈ [0, 1].
Theorem 1. Assume that X is smooth and q-convex, that x 0 is sufficiently close to x † , i.e., x 0 ∈ B D ρ (x † ), (which by [START_REF] Hein | Accelerated Landweber iteration in Banach spaces[END_REF] implies that x †x 0 is also small), that a variational inequality [START_REF] Scherzer | A modified Landweber iteration for solving parameter estimation problems[END_REF] with ν ∈ (0, 1] and β sufficiently small is satisfied, that F satisfies ( 13) with ( 14), (15), that F and F ′ are continuous and uniformly bounded in B D ρ (x † ), that (12) holds and that

q * ≥ 2ν p(ν + 1) -2ν + 1 . ( 24
)
Let n * (δ) be chosen according to (22) with τ sufficiently large. Moreover assume that (19) holds and the sequence {µ n } n∈IN is chosen such that (21) holds for 0 < µ < µ according to (20), and let the sequence {α n } n∈IN ⊆ [0, 1] be chosen such that (23) holds, and α max = max n∈IN α n is sufficiently small. Then, the iterates x δ n+1 remain in B D ρ (x † ) for all n ≤ n * (δ) -1 with n * according to (22). Moreover, we obtain optimal convergence rates

D x0 q (x † , x n * ) = O(δ 2ν ν+1 ) , as δ → 0 ( 25
)
as well as in the noise free case δ = 0

D x0 q (x † , x n ) = O α 2ν p(ν+1)-2ν n (26) for all n IN.
Remark 2. Note that the rate exponent in (26

) 2ν p(ν+1)-2ν = 2ν ν+1 (p -2ν ν+1 ) -1 , always lies in the interval [0, 1 p-1 ], since 2ν ν+1 ∈ [0, 1]
. Moreover, note that Theorem 1 provides a results on rates only, but no convergence result without variational inequality. This corresponds to the situation from [START_REF] Scherzer | A modified Landweber iteration for solving parameter estimation problems[END_REF] in a Hilbert space setting.

Proof. First of all, for x δ n ∈ B D ρ (x † ), ( 13) allows us to estimate as follows (see also ( 16)) in case c 1 ∈ [0, 1):

F (x δ n ) -F (x † ) -A(x δ n -x † ) ≤ K A(x δ n -x † ) c1 D x0 q (x † , x δ n ) c2 ≤ C(c 1 )K A(x δ n -x † ) + D x0 q (x † , x δ n ) c 2 1-c 1 , (27)
where we have used the abbreviation A = F ′ (x † ) and the elementary estimate

a 1-λ b λ ≤ C(λ)(a + b) with C(λ) = λ λ (1 -λ) 1-λ for a, b ≥ 0 , λ ∈ (0, 1) , (28)
and therewith, by the second triangle inequality,

A(x δ n -x † ) ≤ 1 1 -C(c 1 )K F (x δ n ) -F (x † ) + C(c 1 )KD x0 q (x † , x δ n ) c 2 1-c 1
(29) as well as analogously

F (x δ n ) -F (x † ) -A n (x δ n -x † ) ≤ 2C(c 1 )K A(x δ n -x † ) + D x0 q (x † , x δ n ) c 2 1-c 1 ≤ 2C(c 1 )K 1 -C(c 1 )K F (x δ n ) -F (x † ) + D x0 q (x † , x δ n ) c 2 1-c 1 . ( 30
) ≥ 1 -3C(c 1 )K 1 -C(c 1 )K F (x δ n ) -y δ p -C( p-1 p ) 1 -3C(c 1 )K 3C( p-1 p )(1 -C(c 1 )K) F (x δ n ) -y δ p + (3C( p-1 p ) p-1 (1 -C(c 1 )K) 2 p-1 (2C(c 1 )K) p D x0 q (x † , x δ n ) c 2 p 1-c 1 + (1 + C(c 1 )K) p δ p . ( 35
)
Using the variational inequality ( 11), (29), and

(a + b) λ ≤ (a λ + b λ ) for a, b ≥ 0 , λ ∈ [0, 1] , (36)
we get

|α n J X q (x † -x 0 ), x † -x δ n X * ×X | ≤ βα n D x0 q (x † , x δ n ) 1-ν 2 F ′ (x † )(x δ n -x † ) ν ≤ βα n D x0 q (x † , x δ n ) 1-ν 2 1 (1 -C(c 1 )K) ν F (x δ n ) -y δ + δ + C(c 1 )KD x0 q (x † , x δ n ) c 2 1-c 1 ν ≤ βα n D x0 q (x † , x δ n ) 1-ν 2 ǫ -ν ǫ 1 (1 -C(c 1 )K) ν ( F (x δ n ) -y δ + δ) ν +βα n C(c 1 )K (1 -C(c 1 )K) ν D x0 q (x † , x δ n ) 1-ν 2 + νc 2 1-c 1 ≤ C( ν p ) βα n D x0 q (x † , x δ n ) 1-ν 2 ǫ -ν p p-ν + ǫ 1 (1 -C(c 1 )K) ν ( F (x δ n ) -y δ + δ) p +βα n C(c 1 )K (1 -C(c 1 )K) ν D x0 q (x † , x δ n ) 1-ν 2 + νc 2 1-c 1 = C( ν p )      (βǫ -ν ) p p-ν (3C( ν p )C( p(1-ν) 2(p-ν) )) p(1-ν) 2(p-ν) α p(1+ν) 2(p-ν) n   α n D x0 q (x † , x δ n ) 3C( ν p )C( p(1-ν) 2(p-ν) )   p(1-ν) 2(p-ν)      + ǫ 1 (1 -C(c 1 )K) ν ( F (x δ n ) -y δ + δ) p +βα n C(c 1 )K (1 -C(c 1 )K) ν D x0 q (x † , x δ n ) 1-ν-c 1 +νc 1 +2νc 2 2(1-c 1 ) ≤ C( ν p ) C( p(1-ν) 2(p-ν) ) βǫ -ν (3C( ν p )C( p(1-ν) 2(p-ν) )) 1-ν 2 2p p(ν+1)-2ν α p(1+ν) p(ν+1)-2ν n +   α n D x0 q (x † , x δ n ) 3C( ν p )C( p(1-ν) 2(p-ν) )    
Here the last term is nonpositive due to the choice (18) of µ n , so that we arrive at

d 2 n+1 ≤ C 0 d 2c 2 p 1-c 1 n + (1 -1 3 α n + C 1 α q * n )d 2 n + (C 2 + C 3 + C 4 τ -p
C 0 α θω-1 n γ ω - α n+1 α n θ + 1 3 α n -1 -C 1 α q * n α -1 n γ + C 5 ≤ 0 , (42)
where by q * > 1, (15) the factors C 0 α θω-1 n , C 1 α q * -1 n and C 5 can be made small for small α max , β, x †x 0 and large τ . We use this fact to achieve

C 0 α θω-1 n ρ ω-1 + C 1 α q * -1 n ≤ c < c
with c independent of n, which together with (23) yields sufficiency of

C 5 c - c ≤ γ ≤ ρ 2
for (41), (42), which for any (even small) prescribed ρ is indeed enabled by possibly decreasing β, x †x 0 , τ -1 , and therewith C 5 .

1 3 -

 3 xθ α0 > 0, using the abbreviation θ = 2ν p(ν+1)-2ν ∈ [0, 1 p-1 ] we get by the Mean Value Theorem

)

  n + C 1 α q * n )γ n + C 5 (15). Hence as sufficient conditions for uniform boundedness of {γ n } n≤n * by γ and for x δ n+1 ∈ B D ρ (x † ) we get γ ≤ ρ 2 (41)

For any n ≤ n * according to (22), by [START_REF] Kaltenbacher | Convergence rates for the iteratively regularized Gauss-Newton method in Banach spaces[END_REF] we have

where the terms on the right hand side can be estimated as follows. By ( 6) and ( 8) we have

where we have used the triangle inequality in X * and X, the inequality

and [START_REF] Hein | Accelerated Landweber iteration in Banach spaces[END_REF].

For the second term on the right hand side of (31) we get, using (30), ( 28), (34),

where we have used (28) two times and ǫ > 0 will be chosen as a sufficiently small number below. Moreover, by (15), the exponent

) is larger or equal to one and β is sufficiently small so that β

-1 < 1 3 . Finally, we have that

Inserting estimates (32)-( 38) with ǫ = 2 p-1 µ

into (31) and using boundedness away from zero of µ n and the abbreviations

we obtain

In case c 1 = 1, estimates (29), (30) simplify to

and

Therewith, the terms containing

1-c 1 are removed and C(c 1 ) is replaced by ρ 2c2 in (32)-( 38), so that we end up with a recursion of the form (40) (with C 0 replace by zero) as before. Hence the remainder of the proof of uniform boundedness of γ n can be done in the same way in case c 1 < 1.

In case δ = 0, i.e., n * = ∞, uniform boundedness of {γ n } n∈IN implies (26). For δ > 0 we get (25) by using (22) in

Remark 3. In view of estimate (39), an optimal choice of α n would be one that minimizes the right hand side. At least in the special case that the same power of α n appears in the last two terms, i.e., p(1+ν) p(ν+1)-2ν = q * , elementary calculus yields

, which shows that the obtained relation D x0 q (x † , x δ n ) ∼ α 2ν p(ν+1)-2ν n is indeed reasonable and probably even optimal.


Saheed Akindeinde

Daniel Wachsmuth

Adaptive methods for control problems with finite-dimensional control space ⋆

Keywords: optimal control, numerical approximation, a-posteriori error estimates, adaptive refinement

We investigate adaptive methods for optimal control problems with finitely many control parameters. We analyze a-posteriori error estimates based on verification of second-order sufficient optimality conditions. Reliability and efficiency of the error estimator is shown. The estimator is used in numerical tests to guide adaptive mesh refinement.

Introduction

We study optimal control problems of the following type: Minimize the functional J given by J(y, u) = g(y) + j(u) (P) over all (y, u) ∈ Y × U that satisfy the non-linear elliptic partial differential equation E(y, u) = 0 and the control constraints u ∈ U ad .

Here, Y is a real Banach space, U = R n . The set U ad ⊂ U is a non-empty, convex and closed set given by U ad = {u ∈ U : u a ≤ u ≤ u b }, where the inequalities are to be understood component-wise. Here, the cases u a = -∞ and u b = +∞ are allowed, such that problems with one-sided constraints or without control constraints are included in the analysis as well. Examples that are covered by this framework include parameter identification and optimization problems with finitely many parameters, see for instance our previous work [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF].

Adaptive mesh refinement remains a valuable tool in scientific computation. The main objective of an adaptive procedure is to find a discrete solution to a problem while maintaining as few as possible numbers of unknowns with respect to a desired error estimate. As the solution and hence the error distributions on the mesh are unknown a-priori, one has to rely on a-posteriori error estimates.

A-posteriori error estimates for non-linear control and identification problems can be found for instance in [START_REF] Becker | Estimating the control error in discretized pde-constrained optimization[END_REF][START_REF] Griesbaum | Efficient computation of the Tikhonov regularization parameter by goal-oriented adaptive discretization[END_REF][START_REF] Kunisch | Adaptive finite element approximation for a class of parameter estimation problems[END_REF][START_REF] Vexler | Adaptive finite elements for elliptic optimization problems with control constraints[END_REF]. However, they depend on two crucial apriori assumptions: the first is that a second-order sufficient condition (SSC) has to hold at the solution of the continuous problem. With this assumption, error estimates of the type ūu h U ≤ c η + R can be derived, where η is a computable error indicator and R is a second-order remainder term. Here, the second a-priori assumption comes into play: one has to assume that R is small enough, in order to guarantee that mesh refinement solely based on η is meaningful. A different approach with respect to mesh refinement was followed in [START_REF] Ziems | Adaptive multilevel inexact SQP methods for PDEconstrained optimization[END_REF]. There the residuals in the first-order necessary optimality condition were used to derive an adaptive procedure. However, smallness of residuals does not imply smallness of errors without any further assumption. Here again, SSC as well as smallness of remainder terms is essential to draw this conclusion.

In our previous work [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF], we applied a different strategy: There the sufficient optimality condition as well as smallness of remainders is checked a-posteriori. If both conditions are fulfilled, an error-estimator of the form

u -u h U ≤ 2 α (ω y r y + ω p r p )
is available, see [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF]Thm 3.22]. This error estimator is localizable if r y and r p are localizable error estimates for the norm of the residual in the state and adjoint equations, respectively. For earlier and related work on a special problem calls with infinite-dimensional control space, we refer to [START_REF] Rösch | Numerical verification of optimality conditions[END_REF][START_REF] Rösch | How to check numerically the sufficient optimality conditions for infinite-dimensional optimization problems[END_REF].

In this article, we will prove a lower bound of the error estimator. For the setting Y = H 1 0 (Ω), we obtain

r y + r p ≤ c( u -u h U + y -y h Y + ∇y -σ h L 2 (Ω) + p -p h Y + ∇p -τ h L 2 (Ω) + õ),
where y, p and y h , p h are solutions of continuous and discrete state and adjoint equations, respectively, and σ h and τ h are approximations of ∇y h and ∇p h in H(div). The term õ is a higher-order oscillation term. In addition, we have localized lower bounds for the residuals in the state and adjoint equations, respectively. This justifies the use of the a-posteriori estimator above in an adaptive mesh-refinement procedure.

The abstract framework

Let Ω be a polygonal domain in R m , m = 2, 3. The function space for the states of the optimal control problem is chosen as Y := H 1 0 (Ω). Let us now specify assumptions on the abstract problem (P).

Assumption 1

The mapping E : Y × U → Y * , g : Y → R, and j : U → R are twice continuously Fréchet-differentiable with locally Lipschitz continuous second-derivatives. Furthermore, we assume that the mapping E is strongly monotone with respect to the first variable at all points u ∈ U ad .

The assumptions on E are met for instance for semilinear elliptic equations with monotone nonlinearities. Under Assumption (1), the state equation E(y, u) = 0 is uniquely solvable for each admissible control u ∈ U ad , [12, Theorem 26.A, p. 557]. We remark that the differentiability assumption on E can be relaxed to accomodate a more general class of problems, see [

1, Remark 1.2], e.g. differentiability of E from (Y ∩ L ∞ (Ω)) × U to Y * is sufficient.
Let us define the Lagrange functional for the abstract problem:

L(u, y, p) := g(y) + j(u) -E(y, u), p Y * ,Y .
Let (ȳ, ū) be locally optimal for (P). Then the first-order necessary optimality conditions can be expressed as L ′ y (ȳ, ū, p) = 0 and L u (ȳ, ū, p)(u -ū) ≥ 0 for all u ∈ U ad , which is equivalent to

E y (ū, ȳ) * p = g ′ (ȳ) j ′ (ū) -E u (ū, ȳ) * p, u -ū U * ,U ≥ 0 ∀u ∈ U ad .
Since the problem (P) is in general non-convex, the fulfillment of these necessary conditions does not imply optimality. In order to guarantee this, one needs additional sufficient optimality conditions of the type: There exists α > 0 such that

L ′′ (ū, ȳ, p)[(z, v) 2 ] ≥ α v 2 U (1)
holds for all v = u -ū, u ∈ U ad , and z solves the linearized equation E y (ū, ȳ)z + E u (ū, ȳ)v = 0. This condition can be weakened taking strongly active inequality constraints into account, see e.g. [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF][START_REF] Casas | Second-order sufficient optimality conditions for a nonlinear elliptic control problem[END_REF]. For simplicity, we chose to work with this stronger condition. The results of this article hold also under the weakened sufficient condition. Although, the sufficient condition is of high interest, it is difficult to check numerically even when (ū, ȳ, p) are given, see e.g. [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF][START_REF] Rösch | Numerical verification of optimality conditions[END_REF][START_REF] Rösch | How to check numerically the sufficient optimality conditions for infinite-dimensional optimization problems[END_REF]. The main difficulty here is that the function z appearing in (1) is given as solution of a partial differential equation, which cannot be solved explicitly. Any discretization of this equation introduces another error that has to be analyzed.

Discretization

In order to solve (P) numerically, we discretize the problem. Let Y h be a finitedimensional subspace of Y . Here and in the following, the index h denotes a discrete quantity. Then a discretization of the state equation can be obtained in the following way: A function y h ∈ Y h is a solution of the discretized equation for given u ∈ U ad if and only if

E(y h , u), φ h Y * ,Y = 0 ∀φ h ∈ Y h . (2)
The discrete optimization problem is then given by: Minimize the functional J(y h , u h ) over all (y h , u h ) ∈ Y h × U ad , where y h solves the discrete equation.

Let (ȳ h , ūh ) be a local solution of the discrete problem. Then it fulfills the discrete first-order necessary optimality condition, which is given as: there exists a uniquely determined discrete adjoint state ph ∈ Y h such that it holds

E y (ȳ h , ūh ) * ph , φ h Y * ,Y = g ′ (ȳ h ), φ h Y * ,Y ∀φ h ∈ Y h j ′ (ū h ) -E u (ȳ h , ūh ) * ph , u -ūh U * ,U ≥ 0 ∀u ∈ U ad . (3)
Throughout this work, we will assume that errors in discretizing the optimality system are controllable in the following sense. We will not make any further assumptions on the discretization, in particular, we do not assume a sufficient fine discretization.

Assumption 2 For a fixed finite-dimensional subspace Y h , let (u h , y h , p h ) be approximations of the discrete optimal control and the corresponding state and adjoint. There are positive constants r y , r p such that the following holds

E(y h , u h ) Y * ≤ r y , (4)
g ′ (y h ) -E y (y h , u h ) * p h Y * ≤ r p , (5)
j ′ (u h ) -E u (y h , u h ) * p h , u -u h U * ,U ≥ 0 ∀u ∈ U ad . (6)
Here, r y and r p are dual norms of residuals in the state and adjoint equation, respectively, and hence reflect the discretization error. We report on the computation of these residuals in Section 2.

As already mentioned, without any further assumption, smallness of the residuals in ( 4)-( 6) does not imply smallness of the error uu h U in the control. In order to establish such a bound, it is essential to check that a secondorder sufficient optimality condition is satisfied.

Here it is important to recognize that sufficient optimality conditions for the discrete problem alone are still not enough. The sufficient optimality condition for the discrete problem is given by: There exists α h > 0 such that

L ′′ (ū h , ȳh , ph )[(z h , v) 2 ] ≥ α h v 2 U (7)
holds for all v = u -ū, u ∈ U ad , and z h solves the linearized discrete equation

E y (ū h , ȳh )z h + E u (ū h , ȳh )v, φ h Y * ,Y = 0 ∀φ h ∈ Y h . ( 8
)
This condition is equivalent to positive definiteness of a certain computable matrix, see [1, Section 3.5]. Moreover we have the following estimate relating the coercivity constants α and α h appearing in ( 1) and ( 7):

α ≥ α h -E 2 ,
where E 2 is the norm of an error matrix taking the discretization error in the linearized equation E y (ū, ȳ)z + E u (ū, ȳ)v = 0 into account, see [1, Section 3.5] for the details. If the computable lower bound α h -E 2 of α is positive, then it follows that (1) is satisfied. Moreover, we have the following result:

Theorem 1 (Upper bound of the error). Let Assumptions 1 and 2 be satisfied. Let (y h , u h , p h ) be a solution of the discrete optimal control problem. If α h -E 2 > 0 holds and the residuals r y and r p are small enough, then there exists a local solution ū of (P) that satisfies the error bound

ū -u h U ≤ 2 α h -E 2 (ω y r y + ω p r p ) , (9)
where the weights ω y , ω p depend on the discrete solution (y h , u h , p h ). If for different discretizations the discrete solutions {(y h , u h , p h )} h>0 are uniformly bounded in Y × U × Y then the weights ω y , ω p are bounded as well.

For the proof, we refer to [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF]Thm. 3.22]. There, precise estimates of the weights ω y , ω p are given. Moreover, a quantification of the smallness assumption on r y and r p is given, which makes this assumption verifiable a-posteriori.

Corollary 1. Let the assumptions of Theorem 1 be satisfied. Let ȳ, p denote the solutions of the state and adjoint equations to ū, respectively. Then it holds

ȳ -y h Y ≤ υ yu ū -u h U + δ -1 r y , p -p h Y ≤ υ pu ū -u h U + δ -1 r p + υ py r y , with δ -1 being the global bound of E -1 y (y, u) L(Y * ,Y )
, and weights υ yu , υ pu , and υ py depending on (y h , u h , p h ) in the same way as the weights ω y , ω p in Theorem 1.

Proof. The result is a consequence of Theorem 1 and [1, Lemma 3.1, 3.3].

Lower error bounds

In this section, we assume that the general operator E can be written as E(y, u) = -∆y + d(y, u) with d being a superposition operator induced by a smooth function d : R 2 → R. We remark that the subsequent analysis can be easily extended to operators in divergence form with bounded coefficients possibly depending on u [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF]. We will work with a classical finite-element discretization: The discrete space Y h is the classical space of piecewise quadratic and continuous elements (P2) on a given conforming triangulation T h of Ω. The diameter of an element T ∈ T h is denoted by h T . We denote by Σ h ⊂ H(div), a conforming Raviart-Thomas (RT 1 ) discretization of the space H(div).

Let us endow Y = H 1 0 (Ω) with the norm y 2 Y := ∇y 2 L 2 (Ω) + y 2 L 2 (Ω) . In the sequel, let us denote the norm of the embedding H 1 0 (Ω) ֒→ L 2 (Ω) by I 2 . Now let us report on the computation of the residual r y in the state equation. As required by Assumption 2, we are interested in constant-free error estimates, i.e. all constants appearing in the a-posteriori error estimate must be computable. Here, we apply the results of Vohralík [START_REF] Vohralík | Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods[END_REF].

Theorem 2. Let y h ∈ Y h ⊂ H 1 0 (Ω), u h ∈ U ad satisfy the discrete equation (2). Let σ h ∈ Σ h ⊂ H(div) be given such that (div σ h , 1) L 2 (T ) = (d(y h , u), 1) L 2 (T ) for all cells T ∈ T h . ( 10
)
Let us define the cell-wise indicator η y,T , T ∈ T h ,

η y,T := 2 ∇y h -σ h L 2 (T ) + π -1 h T d(y h , u h ) -div σ h L 2 (T ) (11)
Then it holds

-∆y h + d(y h , u h ) 2 H -1 (Ω) ≤ (1 -I 2 2 ) -1 K∈T h η 2 y,T =: r 2 y . ( 12
)
If moreover, T h is shape-regular, then it holds

η y,T ≤ C ∇(ỹ-y h ) L 2 (T ) +c ∇ỹ-σ h L 2 (T ) +ch T d(y h , u h )-Π h d(y h , u h ) L 2 (T ) , (13
) where ỹ := ∆ -1 d(y h , u h ) and Π h denotes the orthogonal L 2 -projection onto Y h . The constants C, c depend only on the spatial dimension m and the shape regularity of the triangulation.

Proof. The upper bound ( 12) is a consequence of [START_REF] Vohralík | Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods[END_REF]Thm. 6.8,6.12] taking [START_REF] Vohralík | Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods[END_REF]Remark 6.3] into account for σ h satisfying (div σ h , 1) L 2 (T ) = (d(y h , u), 1) L 2 (T ) , T ∈ T h . The lower bound [START_REF] Ziems | Adaptive multilevel inexact SQP methods for PDEconstrained optimization[END_REF] follows from [START_REF] Vohralík | Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods[END_REF]Thm. 6.16], see also [START_REF] Vohralík | A posteriori error estimates for lowest-order mixed finite element discretizations of convection-diffusion-reaction equations[END_REF]Lemma 7.6]. Since d(y h , u) is in general not in the discrete space Y h , we obtain the additional oscillation term h T d(y h , u)π h d(y h , u) L 2 (T ) by a standard argument.

Estimates of the residual in the adjoint equation can be obtained after obvious modifications: for τ h ∈ Σ h satisfying (div τ h , 1) L 2 (T ) = (d ′ (y h , u h )p hg ′ (y h ), 1) L 2 (T ) for all cells T ∈ T h (14) and with the local error indicators defined by

η p,T := 2 ∇p h -τ h L 2 (T ) + π -1 h T d ′ (y h , u h )p h -g ′ (y h ) -div τ h L 2 (T ) (15)
we obtain the upper bound

-∆p h + d ′ (y h , u h )p h -g ′ (y h ) 2 H -1 (Ω) ≤ (1 -I 2 2 ) -1 K∈T h η 2 p,T =: r 2 p . ( 16
)
as well as the lower bound

η p,T ≤ C ∇(p -p h ) L 2 (T ) + c ∇p -τ h L 2 (T ) + ch T (I -Π h )(d ′ (y h , u h )p h -g ′ (y h )) L 2 (T ) , (17)
where p :

= ∆ -1 (d ′ (y h , u h )p h -g ′ (y h )).
We remark that the upper bounds ( 12) and ( 16) are constant-free, making them explicitly computable. In our computations, we computed the functions σ h and τ h as a minimizer of the right-hand side in ( 12) and ( 16), respectively, using Raviart-Thomas elements for discretization of H(div). This shows that the requirements of Assumption 2 on the computability of upper bounds on the residuals can be fulfilled. Now let us argue that under the assumptions of Theorem 1 we also obtain lower bounds for the error, which proves efficiency of the error bound. Theorem 3. Let the assumptions of Theorem 1 be fulfilled. Let r y and r p be computed according to [START_REF] Zeidler | Nonlinear functional analysis and its applications[END_REF] and (16). Let (ȳ, ū, p) be the local solution of (P) provided by Theorem 1. Then it holds

T ∈T h r y,T ≤ C ū -u h U + ȳ -y h Y + ∇ȳ -σ h L 2 (Ω) + h T (I -Π h )d(y h , u h ) L 2 (Ω) , (18) T ∈T h r p,T ≤ C ū -u h U + ȳ -y h Y + p -p h Y + ∇p -τ h L 2 (Ω) + h T (I -Π h )(d ′ (y h , u h )p h -g ′ (y h )) L 2 (Ω) , (19)
where C > 0 depends only on the spatial dimension m, the shape regularity of the triangulation, and global bounds of derivatives d y , d u , d yy , and

d yu of d : Y × U → Y * near (y h , u h ).
Proof. Let ỹ be given by ỹ := ∆ -1 d(y h , u h ). Then we can estimate

∇(ỹ -y h ) L 2 (Ω) + ∇ỹ -σ h L 2 (Ω) ≤ 2 ∇(ỹ -ȳ) L 2 (Ω) + ∇(ȳ -y h ) L 2 (Ω) + ∇ȳ -σ h L 2 (Ω) .
Using Lipschitz continuity of d, we find

∇(ỹ -ȳ) L 2 (Ω) ≤ C( ū -u h + ȳ -y h Y ),
with C depending on bounds of d ′ L(Y ×U,Y * ) near (y h , u h ). This together with (13) proves (18). The estimate (19) can be obtained analogously.

These lower bounds together with [START_REF] Vexler | Adaptive finite elements for elliptic optimization problems with control constraints[END_REF] and the local lower bounds ( 13) and (17) justify the use of the error indicators in an adaptive mesh-refinement procedure.

Remark 1. Another possibility of constant-free a-posteriori error estimators based on H(div)-functions is described in [START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF]. There, fluxes across edges in a dual mesh are prescribed instead of the integrals on elements as in [START_REF] Vohralík | A posteriori error estimates for lowest-order mixed finite element discretizations of convection-diffusion-reaction equations[END_REF] and (14). In [START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF] it is proven that the resulting error estimate is reliable and efficient. Moreover, the terms ∇y hσ h L 2 (Ω) and ∇p hτ h L 2 (Ω) do not appear in the lower error bound when compared to (18) and (19), respectively.

Adaptivity

We will compare the performance of adaptive mesh refinement using different strategies to mark elements for refinement. The first one, referred to as 'verified adaptive', is implemented as follows: in each step the verification procedure of [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF] is carried out. If it confirms that the assumptions of Theorem 1 are satisfied, then the error indicator ω y r y + ω p r p given by ( 9) is used to guide the meshrefinement. If the requirements of Theorem 1 cannot be verified, then a uniform refinement step is carried out. Here, we expect that after a small number of uniform refinement steps the requirements of Theorem 1 are confirmed a-posteriori, which coincides with the numerical experiments done in earlier work [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF]. After these initial uniform refinements steps, we expect that the method proceeds with adaptive steps.

A second strategy, called 'fully adaptive', omits the verification step, and simply uses ω y r y + ω p r p from (9) without checking the validity of this bound.

Numerical results

Let us report about the outcome of the above described adaptive methods for a selected example, taken from [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF]. The functional J was chosen as

J(y, u) := 1 2 y -y d 2 L 2 (Ω) + κ 2 u 2 R n .
The nonlinear mapping E represents a semi-linear elliptic equation given by

E(y, u) := -∆y + n k=1 u k d k (y) -g, (20)
where the functions d k are chosen as d 1 (y) = 1, d j (y) = y|y| j-2 for j = 2 . . . n. This example is motivated by parameter identification: given a state y d and source term g, find the set of coefficients u such that the resulting solution y of E(y, u) = 0 is as close as possible to y d . In order to make the operator E strongly monotone, we require positivity of the coefficients u k , i.e. we set U ad = {u ∈ R n : u k ≥ 0 ∀k = 1 . . . n}. For the computations we used the following data: the source term g = 10.0001 and

Ω = (0, 1) 2 , u a = 0, u b = 0.5, κ = 10 -2 , y d (x 1 , x 2 ) = 0.5 sin(2πx 1 x 2 ), n = 4.
Let us remark, that the function d 3 is not of class C 2 globally. Since g is nonnegative, every solution y of (20) to u ∈ U ad will be non-negative. For nonnegative functions y it holds d 3 (y) = y 2 , which is C 2 , so the assumptions on E are satisfied. See also the discussion in [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF]Section 4.3].

We employed a discretization scheme as described in Section 2. After the resulting non-linear optimization problem is solved, the error indicators according to the chosen strategy are computed. For an adaptive refinement, a subset T ⊂ T of elements T with large local error contributions η T were selected for refinement that satisfies T ∈ T η 2 T ≥ θ 2 T ∈T η 2 T with θ = 0.8. Let us report on the outcome of the different adaptive strategies as described in Section 3. For all the methods, we compare the residual norms as given in [START_REF] Vohralík | Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods[END_REF], i.e. with the notation of that section ǫ residual := ω y r y + ω p r p .

Moreover, we employed the verification procedure of Theorem 1, see e.g. [START_REF] Akindeinde | A-posteriori verification of optimality conditions for control problems with finite-dimensional control space[END_REF], and report about the upper error bound

ǫ bound := 2 α h -E 2
(ω y r y + ω p r p ) . As can be expected, the assumptions of Theorem 1 are only fulfilled on a sufficiently fine discretization. This is reflected by our numerical results.

Plots of ǫ residual and ǫ bound versus the number of degrees of freedom can be seen in Figure 4. For reference, we provided the numerical values in Tables 1 and2. In the tables, L refers to refinement level, where L = 0 is the initial mesh, which is the same for all the different adaptive methods. Moreover, dof denotes the number of degrees of freedom.

Let us comment on the observed behavior of the verified adaptive methods. The conditions of Theorem 1 are fulfilled for the first time after three uniform refinement steps. The fourth and all further refinement levels were reached by using adaptive refinement according to the error indicator based on [START_REF] Vohralík | Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods[END_REF]. The fully adaptive scheme, which refines according to the residuals in the optimality system, obtains verified error bounds as of level 7. After the verified adaptive methods actually start adaptive refinement, they quickly reach the same ratio of error bound versus number of degrees of freedom as the full adaptive method. That

Fig. 1 .

 1 Fig. 1. (a) upper bound of residuals versus number of unknowns, (b) verified error bound versus number of unknowns

⋆ This work was funded by Austrian Science Fund (FWF) grant P21564-N18.

means, the early (unverified) adaptive refinements of the full adaptive methods does not seem to give this method an advantage over the verified method. The same observation also applies to the residual error versus number of degrees of freedom ratio, as can be seen in Figure 4.


Chantal Landry
email: chantal.landry@wias-berlin.de

Matthias Gerdts
email: matthias.gerdts@unibw.de

René Henrion
email: rene.henrion@wias-berlin.de

Dietmar Hömberg
email: dietmar.hoemberg@wias-berlin.de

Path-Planning with Collision Avoidance in Automotive Industry

Keywords: Optimal control, collision avoidance, backface culling, active set strategy

An optimal control problem to find the fastest collision-free trajectory of a robot is presented. The dynamics of the robot is governed by ordinary differential equations. The collision avoidance criterion is a consequence of Farkas's lemma and is included in the model as state constraints. Finally an active set strategy based on backface culling is added to the sequential quadratic programming which solves the optimal control problem.

Collision Avoidance

In automotive industry robots have to work simultaneously on the same workpiece. The challenge is to find for each robot the fastest trajectory that avoids any collision with the surrounding obstacles and the other robots. We start with the establishment of the collision avoidance criterion.

For simplicity, we suppose that only one obstacle is present in the workspace. As in [START_REF] Gilbert | A New Algorithm for Detecting the Collision of Moving Objects[END_REF][START_REF] Gilbert | Distance Functions and their Application to Robot Path Planning in the Presence of Obstacles[END_REF] a collision detection can be obtained when the robot is approximated by a union of convex polyhedra. This union is called P and it given by P = nP i=1 P (i) , with

P (i) = {y ∈ R 3 |A (i) y ≤ b (i) }
where n P is the number of polyhedra in P . If p i denotes the number of faces in P (i) , then A (i) ∈ R pi×3 and b (i) ∈ R pi for i = 1, . . . , n P .

Similarly, the obstacle is approximated by the following union of convex polyhedra, called

Q Q = nQ j=1 Q (j) , with Q (j) = {y ∈ R 3 |C (j) y ≤ d (j) }
where n Q is the number of polyhedra in Q. If q j is the number of faces in Q (j) , then C (j) ∈ R qj ×3 and d (j) ∈ R qj for j = 1, . . . , n Q . In the following, n P , A, b and i are always associated with the robot, and n Q , C, d, and j are related to the obstacle. Furthermore, the robot will be identified with its approximation P and the obstacle with Q.

A first characterization of the collision-freeness between P and Q is given by P (i) ∩ Q (j) = ∅, ∀i = 1, . . . , n P and ∀j = 1, . . . , n Q .

The definition of the polyhedra P (i) and Q (j) implies that P (i) does not collide with Q (j) if and only if there does not exist any point y (i,j) ∈ R 3 satisfying

A (i) C (j) y (i,j) ≤ b (i) d (j) .
According to Farkas's lemma [START_REF] Berkovitz | Convexity and Optimization in R n[END_REF], this linear system does not have any solution if and only if there exists a vector w (i,j) ∈ R pi+qj such that

w (i,j) ≥ 0, A (i) C (j) ⊤ w (i,j) = 0 and b (i) d (j) ⊤ w (i,j) < 0.
In conclusion, the pair of polyhedra (P (i) , Q (j) ) is collision-free if and only if such a vector w (i,j) exists. This forms the collision avoidance characterization between a pair of polyhedra. Between the robot and the obstacle, the characterization reads as follows:

Proposition 1. Two unions of convex polyhedra P = nP i=1 P (i) and Q = nQ j=1 Q (j)
do not collide if and only if for each pair of polyhedra (P (i) , Q (j) ), i = 1, . . . , n P , j = 1, . . . , n Q , there exists a vector w (i,j) ∈ R pi+qj such that w (i,j) ≥ 0,

A (i) C (j) ⊤ w (i,j) = 0 and b (i) d (j) ⊤ w (i,j) < 0.

Optimal Control Problem

To express the dynamics of the robot, we need to describe P differently from the previous section. As an industrial robot, P is composed by m links and is asked to move from its current position to a desired point [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF]. Let q = (q 1 , . . . , q m ) denote the vector of joint angles at the joints of the robot. The vector v = (v 1 , . . . , v m ) contains the joint angle velocities and u = (u 1 , . . . , u m ) describes the torques applied at the center of gravity of each link. The Lagrangian form of the dynamics of the robot depends on these three vectors as follows

q ′ (t) = v(t) and M(q(t)) v ′ (t) = G(q(t), v(t)) + F (q(t), u(t)), (1)
where M(q) is the symmetric and positive definite mass matrix, G(q, v) contains the generalized Coriolis forces and F (q, u) is the vector of applied joint torques and gravity forces [START_REF] Lavalle | Planning Algorithms[END_REF][START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF].

For the remainder of the paper, let us define the vector x = (q, v) ∈ R nx with n x := 2m. With the definition of x and the non-singularity of the matrix M, we can define the function f : R nx × R m → R nx as follows

f (x, u) = v M -1 (q) (G(q, v) + F (q, u))
.

The fastest trajectory of a robot is the solution of an optimal control problem where the system of ordinary differential equations (ODE) are given by [START_REF] Berkovitz | Convexity and Optimization in R n[END_REF], see [START_REF] Diehl | Fast Direct Multiple Shooting Algorithms for Optimal Robot Control[END_REF][START_REF] Gerdts | Path Planning and Collision Avoidance for Robots[END_REF]. If an obstacle is present in the workspace, the collision-freeness is assured as soon as the vector w (i,j) of Proposition 1 is found at each time t and for all pairs of polyhedra. However, to be written as state constraints, the strict inequality in Proposition 1 has to be relaxed. Furthermore, since the robot moves, the matrices A (i) and the vectors b (i) evolves in time. Their evolution depends explicitly on x(t). A complete formulation of A (i) (x(t)) and b (i) (x(t)) is given in [START_REF] Gerdts | Path Planning and Collision Avoidance for Robots[END_REF].

Before writing down the model, let us define the index transformation I = (i -1)n Q + j. Hence, to each pair (i, j) ∈ {1, . . . , n P } × {1, . . . , n Q } there corresponds an index I in {1, . . . , n P n Q }, and reciprocally. In the sequel, the index I is used instead of the pair (i, j). The variable w (i,j) is then numbered as w I . Let us also define the functions G I : R nx → R (pi+qj )×3 and g I : R nx → R pi+qj for I = 1, . . . , n P n Q as follows

G I (x) = A (i) (x) C (j) and g I (x) = b (i) (x) d (j) .
G I and g I allow us to write Proposition 1 as a function of time. Finally let set M := n P n Q the number of indices I and n I := p i + q j the size of w I for I = 1, . . . , M , and let t f denote the travel time. Then, after transformation onto the fixed time interval T := [0, 1] the optimal control problem reads as follows:

(OCP): minimize ϕ(x(0), x(1), t f )

with respect to the state variable x ∈ W nx 1,∞ (T ), the control variables u ∈ L m ∞ (T ) and w I ∈ L nI ∞ (T ), I = 1, . . . , M , and t f ≥ 0, subject to:

-ODE: x ′ (t) -t f f (x(t), u(t)) = 0, a.e. in T, -state constraints: G I (x(t)) ⊤ w I (t) = 0, I = 1, . . . , M, a.e. in T, (2) g
I (x(t)) ⊤ w I (t) ≤ -ε, I = 1, . . . , M, a.e. in T, (3)
-boundary conditions: ψ(x(0), x(1)) = 0, -box constraints:

w I (t) ≥ 0, I = 1, . . . , M, a.e. in T, u(t) ∈ U := {u ∈ R m | u min ≤ u ≤ u max }.
where the function ψ :

R nx × R nx → R 2nx is given as follows ψ(x(t 0 ), x(t f )) = (R(q(t 0 ))-R 0 , 0, R(q(t f ))-R f , 0)
where R(q) denotes the position of the barycenter of the last link of the robot and R 0 , R f ∈ R m are given. The vectors u min and u max are also given. The relaxation parameter ε is positive and small. As usual L m ∞ (T ) denotes the Banach space of essentially bounded functions mapping from T into R m and W nx 1,∞ (T ) denotes the Banach space of absolutely continuous functions with essentially bounded derivative that map from T into R nx .

If multiple obstacles are present in the workspace, the anti-collision constraints (2)-( 3) and the associated control variables w I are defined for each obstacle.

Depending on the number M of anti-collision constraints, the problem is inherently sparse since the artificial control variables w I , I = 1, . . . , M , do not enter the dynamics, the boundary conditions, and the objective function of the problem, but only appear linearly in the anti-collision constraints with one-sided coupling through the state.

We attempt to solve the problem (OCP) numerically with a reduced discretization approach [START_REF] Gerdts | Direct Shooting Method for the Numerical Solution of Higher-Index DAE Optimal Control Problems[END_REF]. Let us consider the control grid G N := {t k = kh | k = 0, 1, . . . , N }, which, for simplicity, is chosen equidistantly with the fixed step-size h = 1/N . We use B-spline of second order to approximate the control variables:

u h (t; u 0 , . . . , u N ) := N i=0 u i B i (t), w I,h (t; w I,0 , . . . , w I,N ) := N i=0 w I,i B i (t), I = 1, . . . , M
where (u 0 , . . . , u N ) ⊤ ∈ R m(N +1) and (w I,0 , . . . , w I,N ) ⊤ ∈ R nI (N +1) are the vector of de Boor points, and B i , i = 0, . . . , N , denote elementary B-splines. As the elementary B-splines sum up to one for all t ∈ T , the box constraints u h (t) ∈ U are satisfied, if u i ∈ U, i = 0, . . . , N. The choice of B-splines is convenient as it is easy to create approximations with prescribed smoothness properties and, even more important, the elementary B-splines B i have a local support only.

We solve the differential equations for the initial value x 0 and a given t f by the classical explicit Runge-Kutta method of order 4. The state approximations at the grid points t k , k = 0, . . . , N depend on the vector z := (x 0 , u 0 , . . . , u N , t f ) ⊤ ∈ R nz with n z = n x + (N + 1)m + 1.

Let us define J(z) := ϕ(x 0 , x N (z), t f ), h(z) := ψ(x 0 , x N (z)), as well as ḠI,k (z) := G I (x k (z)) and ḡI,k (z) := g I (x k (z)) for I = 1, . . . , M, k = 0, . . . , N . With these new notations the discretized form of (OCP) can be formulated as follows (DOCP):

Minimize J(z) with respect to z ∈ R nz and w = (w 1,0 , . . . , w 1,N , . . . , w M,0 , . . . , w

M,N ) ⊤ ∈ R (N +1) M I=1 nI subject to: h(z) = 0, w I,k ≥ 0, I = 1, . . . , M, k = 0, . . . , N, ḠI,k (z) ⊤ w I,k = 0, I = 1, . . . , M, k = 0, . . . , N, ḡI,k (z) ⊤ w I,k ≤ -ε, I = 1, . . . , M, k = 0, . . . , N, z ∈ Z := {z ∈ R nz | z ℓ ≤ z ≤ z u }.
Herein, z ℓ ≤ z u define box constraints for z, where the settings ±∞ are permitted if a component of z is not restricted from above or below.

The above nonlinear optimization problem is solved by a sequential quadratic programming (SQP) method [START_REF] Gerdts | Optimal Control of Ordinary Differential Equations and Differential-Algebraic Equations[END_REF][START_REF] Gill | SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization[END_REF]. As in [START_REF] Schittkowski | On the Convergence of a Sequential Quadratic Programming Method with an Augmented Lagrangean Line Search Function[END_REF] we use an Armijo type line-search procedure for the augmented Lagrangian function in our implementation. However, (DOCP) contains a lot of constraints: at each time step t k , k = 0, . . . , N , and for every pair of polyhedra (P (i) , Q (j) ), four anti-collision constraints are defined (compare (2)-( 3)). To overcome this difficulty, we add an active set strategy based on the following observation: the anti-collision constraints are superfluous when the robot is far from the obstacle or moves in the opposite direction. The establishment of the active set strategy is the purpose of the next section.

Backface Culling Active Set Strategy

Backface culling comes from computer graphics and consists of working only with visible objects, see [START_REF] Vaněček | Back-Face Culling Applied to Collision Detection of Polyhedra[END_REF]. We apply here the same concept to define our active set strategy and develop four criteria to determine which objects are visible.

In this section P , resp. Q, is a polyhedron belonging to the approximation of the robot, resp. obstacle. The first criterion is similar to the broad phase in [START_REF] Cohen | I-collide: An Interactive and Exact Collision Detection System for Large-Scaled Environments[END_REF][START_REF] Mirtich | V-Clip: Fast and Robust Polyhedral Collision Detection[END_REF] and consists of checking if P is far from Q. If this is the case, no collision can occur and the anti-collision constraints are superfluous.

The distance between P and Q is roughly computed by defining the axisaligned bounding box of each polyhedron. Let (x P i , y P i , z P i ), i = 1, . . . , s P , denote the vertices of P where s P is the number of vertices, and define The vertices of P evolve in time since they belong to the robot. Hence the box B P has to be determined at each grid point t k , k = 0, . . . , N . Because the obstacle does not move, the box B Q is computed only once.

x P m =
Let us assume now that Q is close enough to P and consider the situation depicted in Figure 2: P is moving downwards, v c indicates the velocity of the center of gravity and Q is generated by the faces e 1 and e 2 of Q. According to Proposition 1 Q does not collide with P if and only if

∃ w > 0, such that A C 1,2 ⊤ w = 0 and b d 1,2 ⊤ w < 0, (4)
where C 1,2 is the matrix composed of the first two rows of C and d 1,2 is the vector composed of the first two components of d. Suppose now that w exists. By setting w = ( w, 0, 0, 0), we obtain:

A C ⊤ w = 0 and b d ⊤ w < 0.
Then, Proposition 1 implies that P and Q do not collide. In summary, if no collision occurs between Q and P , then Q and P do not collide. The dimension of w is always smaller than that of w, because the polyhedra are supposed to be compact. Consequently, the problem of finding w is always smaller and there is an advantage in replacing the anti-collision constraints by [START_REF] Gerdts | Direct Shooting Method for the Numerical Solution of Higher-Index DAE Optimal Control Problems[END_REF]. In (4) only the faces visible to P are taken into consideration. The next criteria concern the determination of the visible faces of Q relative to P . The faces of Q which are located behind P , are invisible to P . P is always looking towards its velocity, v c . Hence, all objects located in the lower halfspace H generated by v c and S R , the vertex of P located furthest in the opposite direction of v c , are behind P . Then, the second criterion reads

Criterion 2 A face e of Q is invisible to P if e ⊂ H = {y ∈ R n |v ⊤ c (y -S R ) < 0}.
An example is given in Figure 2-(a) where the faces e 3 and e 4 satisfy Criterion 2. The case where all faces are located behind P means that P is moving in the opposite direction to Q. In this situation no collision can occur and we have: Not all remaining faces of Q are visible to P . Some of them can be hidden by other faces of Q. This is the case of the face e 5 in Figure 2-(b). The vertex S, which is the closest point of P to e 5 , cannot see e 5 because e 5 is hidden by e 1 .

The face e 5 will be visible to S as soon as S is no more located in the halfspace A limit case exists with Criterion 4 when P is included in Q. In that case all faces of Q are invisible to P according to Criterion 4. But in fact all these faces must be considered in the anti-collision constraints. Hence, Criterion 4 must not be applied in this particular case.

Criterion 4 can also be applied to detect which faces of P are visible to Q. Then the anti-collision constraints defined for the pair (P, Q) can be reduced as it was done in [START_REF] Gerdts | Direct Shooting Method for the Numerical Solution of Higher-Index DAE Optimal Control Problems[END_REF].

In this section criteria to determine the visible faces of Q were developed, provided Q is visible. All criteria were depending on the position of P which is given by the state variable q. In the next section we show how the backface culling is included in the SQP algorithm to solve (DOCP).

Algorithm and numerical examples

Let us recall the index transformation that associates to each pair (i, j) the new index I via the formula: I = (i -1)n Q + j and define the set of indices

K := {(I, k) | the polyhedron Q (j) is visible to P (i) at t k }.
K is determined by applying Criteria 1, 2 and 3. Let us also recall that w I belongs to R pi+qj . The first p i components of w I are associated to the faces of At the step (4), if a pair of indices (I, k) newly appears in K (ℓ+1) (i.e. if (I, k) ∈ K (ℓ+1) \ K (ℓ) ), then the variable w (ℓ+1) I,k must be initialized. We choose to take w (ℓ+1) I,k as the solution of

min w ḡI,k (z (ℓ+1) ) ⊤ w such that ḠI,k (z (ℓ+1) ) ⊤ w = 0, w c > 0, if c / ∈ J (ℓ+1) I,k
, and

w c = 0, if c ∈ J (ℓ+1) I,k
.

Thereby, we are assured to satisfy the state constraints as close as possible.

In our first example we consider an obstacle and a robot composed by a socket and 3 links. A load is fixed at the end of the last link. A complete description of the example is given in [START_REF] Gerdts | Path Planning and Collision Avoidance for Robots[END_REF]. In this example, the collision avoidance needs to be applied only between the load and the obstacle. The obstacle is always close to the load. Consequently, the number of state contraints is not reduced with the backface culling. For this example, we take 21 control grid points and ε = 10 -5 . In Figure 3 the visible faces of the obstacle are in white. We can observe that only 3 faces of the obstacle are visible. The computational time is 52 s. If we do not use the backface culling, the computational time is equal to 3 min 44. So, with about half of the unknowns, the code runs about four times faster. The second example is in 2 dimensions and uses all criteria of the backface culling. The robot P is a square and 4 obstacles, Q 1 to Q 4 , are present in the workspace. We take for this example 42 control grid points and ε = 10 -2 . Snapshots of the motion of the robot, which must reach the point F , are given in Figure 4. The visible obstacles are in white and their visible faces in gray. The computational time is equal to 36 min 50 when no backface culling is used. With the backface culling strategy the CPU time is 27 s. Hence, Criteria 1 to 4 induce a large decrease in the computational time.

Fig. 1 .

 1 Fig. 1. (a) The polyhedron P is moving downwards. The faces of Q are denoted by e1, . . . , e5. (b) The set Q is generated by the faces of Q visible to P .

Criterion 3 Remark 1 .

 31 If all faces of Q are invisible to P according to Criterion 2, then the anti-collision constraints are superfluous. If the velocity v c is zero, Criteria 2 and 3 are not applied.

H 5 =Fig. 2 .

 52 Fig. 2. (a) The faces e3 and e4 of Q are located behind P . (b) The face e5 of Q is invisible to P .

at t1 at t4 at t9 at t14 at t17 Fig. 3 .

 t173 Fig. 3. Snapshots of the motion of the robot avoiding an obstacle. The visible faces of the obstacle are in white.

Fig. 4 .

 4 Fig. 4. Snapshots of the motion of the robot P moving to F and avoiding four obstacles. The visible obstacles are in white and their visible faces in gray.

Acknowledgments. This work has been supported by the DFG Research Center MATHEON Mathematics for key technologies.

P (i) and the next q j components are related to the faces of Q (j) . Then, let us define the following set of indices for each (I, k) ∈ K J I,k := {c ∈ {1, . . . , p i } | the face c of

This set contains the index of the faces of the pair (P (i) , Q (j) ) which are invisible at t k . The invisibility of a face is determined using Criteria 2 and 4.

Backface culling involves considering the anti-collision constraints whose pair of indices (I, k) belongs to K and write these constraints according to [START_REF] Gerdts | Direct Shooting Method for the Numerical Solution of Higher-Index DAE Optimal Control Problems[END_REF]. The algorithm to solve (DOCP ) is the SQP method presented below in which the backface culling is added as an active set strategy. This means that at each iteration we update the set K and then build the quadratic problem by considering only the constraints whose pair of indices belongs to K.

Backface Culling Active Set Strategy (0) Choose ε > 0, z (0) ∈ Z and w (0) ≥ 0.

Determine the sets of indices K (0) and J

I,k for all (I, k) ∈ K (0) . Set B 0 := I, the identity matrix and ℓ := 0.

(1) If (z (ℓ) , w (ℓ) ) is a KKT point of the optimization problem, STOP.

(2) Compute a KKT point of the following linear-quadratic optimization problem: Minimize

I,k , (I, k) ∈ K (ℓ) .

Note: The constraints d w I,k ,c = 0 are only included for notational simplicity.

In practice these variables are actually eliminated.

(4) Update the sets of indices K (ℓ+1) and J (ℓ+1) I,k for (I, k) ∈ K (ℓ+1) according to Criteria 1 to 4 which depend on z (ℓ+1) . Update B ℓ+1 according to BFGS update formulas, set ℓ := ℓ + 1 and go to (1).