Kurt Helmes
email: helmes@wiwi.hu-berlin.de

Richard H Stockbridge
email: stockbri@uwm.edu

Chao Zhu
email: zhu@uwm.edu

Impulse control of standard Brownian motion: Long-term average criterion

Keywords: impulse control, long-term average criterion, infinite dimensional linear programming, expected occupation and impulse measures

This paper examines the impulse control of a standard Brownian motion under a long-term average criterion. In contrast with the dynamic programming approach, this paper first imbeds the stochastic control problem into an infinite-dimensional linear program over a space of measures and then reduces the problem to a simpler nonlinear optimization that has a familiar interpretation. One is able to easily identify the optimal cost and a family of optimal impulse control policies.

Introduction

When one seeks to control a stochastic process and every intervention incurs a strictly positive cost, one must select a sequence of separate intervention times and amounts. The resulting stochastic problem is therefore an impulse control problem in which the decision maker seeks to either maximize a reward or minimize a cost. This paper examines the impulse control of the prototypical process Brownian motion under a long-term average cost criterion; a companion paper studies the impulse control of Brownian motion under a discounted criterion. The aim of the paper is to illustrate a solution approach which first imbeds the stochastic control problem into an infinite-dimensional linear program over a space of measures and then reduces the linear program to a simpler nonlinear optimization. This approach provides a new method for determining an optimal impulse control policy.

Let W be a standard Brownian motion process with natural filtration {F t }. An impulse control policy consists of a pair of sequences (τ, Y ) := {(τ k , Y k ) : k ∈ N} in which τ k is the {F t }-stopping time of the kth impulse and the F τ kmeasurable variable Y k gives the kth impulse size. The sequence {τ k : k ∈ N} is required to be non-decreasing, a natural assumption in that intervention k + 1 must occur no earlier than intervention k. For a policy (τ, Y ), the impulsecontrolled Brownian motion process is given by

X(t) = x 0 + W (t) + ∞ k=1 I {τ k ≤t} Y k .
The goal of the decision maker is to keep the process close to zero and to minimize the long-term average cost incurred by the impulse policy. Define the running/deviation cost rate function by c 0 (x) = x 2 and set the impulse costs to be c 1 (y, z) = k 1 + k 2 |y -z| for (y, z) ∈ R 2 , in which k 1 > 0 and k 2 ≥ 0 and y denotes the pre-intervention location which will typically be far away from zero, while z denotes the post-intervention location of the process X, which should be close to zero. Note, in particular, that there is a strictly positive cost for every impulse, even one in which Y k = 0 which does not affect the value of X. Let (τ, Y ) be an impulse control policy. The quantity to be minimized is

J(τ, Y ) := lim sup t→∞ t -1 E t 0 c 0 (X(s)) ds + ∞ k=1 I {τ k ≤t} c 1 (X(τ k -), X(τ k )) . (1)
Clearly any policy (τ, Y ) for which J(τ, Y ) = ∞ is undesirable so to be admissible, we require (τ, Y ) to have a finite cost; the nonnegativity of c 0 and strict positivity of c 1 indicates that every policy will have nonnegative long-term average cost. The collection of all admissible impulse policies is denoted by A.

Similar type of problems have been extensively investigated in the literature. An incomplete list includes the now classical works on general stochastic impulse problems [START_REF] Bensoussan | Impulse control and quasi-variational inequalities[END_REF][START_REF] Harrison | Impulse control of brownian motion[END_REF][START_REF] Lions | Quasi-variational inequalities and ergodic impulse control[END_REF][START_REF] Robin | On some impulse control problems with long run average cost[END_REF] as well as their applications in various areas such as portfolio optimization, inventory control, risk management, control of a dam and exchange rate intervention [START_REF] Chen | Risk aversion in inventory management[END_REF][START_REF] Eastham | Optimal impulse control of portfolios[END_REF][START_REF] Melas | An ergodic impulse control model with applications[END_REF][START_REF] Ormeci | Impulse control of brownian motion: the constrained average cost case[END_REF][START_REF] Piunovskiy | Examples in Markov Decision Processes[END_REF]. In particular, [START_REF] Ormeci | Impulse control of brownian motion: the constrained average cost case[END_REF] explains the adoption of a Brownian motion model.

Unlike the aforementioned references, in which the primary tool is the dynamic programming principle and its associated quasi-variational inequalities, this paper aims to illustrate the utility of a different methodology, namely, the linear programming approach. In such an approach, we embed the stochastic impulse control problem into an infinite-dimensional linear program over appropriate measures. Further, the linear program, with the aid of an auxiliary linear program, is transformed into a nonlinear optimization problem. Then both the value of the impulse control problem and an optimal impulse control policy are easily determined. The linear programming approach toward stochastic control problem can be dated back to [START_REF] Manne | Linear Programming and Sequential Decisions[END_REF] for discrete time and a finite state space and to [START_REF] Stockbridge | Time-average control of martingale problems: Existence of a stationary solution[END_REF][START_REF] Stockbridge | Time-average control of martingale problems: A linear programming formulation[END_REF] for regular stochastic control problems in continuous time with general state space. It has been further developed in [START_REF] Helmes | Construction of the value function and optimal rules in optimal stopping of one-dimensional diffusions[END_REF][START_REF] Kurtz | Existence of Markov controls and characterization of optimal Markov controls[END_REF][START_REF] Kurtz | Stationary solutions and forward equations for controlled and singular martingale problems[END_REF] for optimal stopping and singular control problems. This paper aims to expand the utility of such a methodology to impulse control problems as well.

We make four important observations about impulse policies. Firstly, the "no-intervention policy" in which the process is X(t) = x 0 + W (t) for all t ≥ 0 incurs an infinite long-term average cost so is inadmissible. Secondly, let (τ, Y ) be an impulse policy and define τ ∞ := lim k→∞ τ k . Should τ ∞ < ∞ on a set of positive probability, then the fixed cost k 1 > 0 per intervention also results in an infinite long-term average cost. Thus for every admissible policy τ k → ∞ a.s. as k → ∞. Thirdly, let (τ, Y ) be a policy for which there is some k such that τ k = τ k+1 on a set of positive probability. Again due to the presence of the fixed intervention cost k 1 , the total cost up to time τ k+1 will be at least k 1 P(τ k = τ k+1 ) smaller by combining these interventions into a single intervention on this set. Hence we may restrict policies to those for which τ k < τ k+1 almost surely for each k.

The final observation is similar. Suppose (τ, Y ) is a policy such that on a set G of positive probability τ k < ∞ and |X(τ k )| > |X(τ k -)| for some k. Consider a modification of this impulse policy and resulting process X which simply fails to implement this impulse on G. Define the stopping time σ = inf{t > τ k :

|X(t)| ≤ | X(t)|}.
Notice that the running costs accrued by X over [τ k , σ) are smaller than those accrued by X. Finally, at time σ, introduce an intervention on the set G which moves the X process so that X(σ) = X(σ). This intevention will incur a cost which is no greater than the cost for the process X at time τ k . As a result, we may restrict the impulse control policies to those for which no impulse increases the distance of the process from the origin (an intuitively obvious observation).

Restricted problem and measure formulation

The initial analysis considers a restricted collection of impulse policies.

Condition 1 Let A 1 denote the set of policies such that for (τ, Y ) ∈ A 1 the resulting process X remains bounded; that is for some M < ∞, |X(t)| ≤ M for all t ≥ 0.

Intuitively, Condition 1 is not much of a restriction since unbounded processes occur by allowing the Brownian motion to diffuse which incurs an expensive running cost. This restriction, however, is needed so that a transversality condition is satisfied and a stochastic integral is a martingale. Following the initial solution, the general class of impulse policies will be analyzed.

We capture the expected behavior of the process and impulses with measures. Arbitrarily fix (τ, Y ) ∈ A 1 and let M be as given in Condition 1. For each t > 0, define the average expected occupation and average expected impulse measures µ (t) 0 and µ

(t) 1 , respectively, such that for each G, G 1 , G 2 ∈ B(R), µ (t) 0 (G) = t -1 E t 0 I G (X(s)) ds , and
µ (t) 1 (G 1 × G 2 ) = t -1 E ∞ k=1 I {τ k ≤t} I G1×G2 (X(τ k -), X(τ k )) . (2) It is immediate that µ (t) 0
is a probability measure for each t > 0 and since (τ, Y ) ∈ A 1 ⊂ A, the finiteness of J(τ, Y ) implies that the collection of measures {µ (t) 1 : t > 0} is uniformly bounded above. We also note that the second component being measured by µ (t) 1 is the post-jump location so µ (t) 1 is a measure on the product space of (pre-jump, post-jump) pairs. In light of Condition 1, the support of µ

J(τ, Y ) = lim sup t→∞ t -1 c 0 (x) µ (t) 0 (dx) + c 1 (y, z) µ (t) 1 (dy × dz) .
Let {t j : j ∈ N} be a sequence with t j → ∞ as j → ∞ such that

t -1 j c 0 (x) µ (tj ) 0 (dx) + c 1 (y, z) µ (tj ) 1 (dy × dz) → J(τ, Y ).
For i = 0, 1, the relative compactness of {µ (tj ) i

: j ∈ N} implies that there exist weak limits µ 0 and µ 1 . Note µ 0 is a probability measure whereas µ 1 is a finite measure. Since c 0 and c 1 are bounded and continuous on [-M, M ] and R, respectively,

J(τ, Y ) = c 0 (x) µ 0 (dx) + c 1 (y, z) µ 1 (dy × dz). (3)
It is now helpful to characterize the value of functions of the process.

For f ∈ D = C 2 (R), f (X(t)) = f (x 0 ) + t 0 f (X(s)) dW (s) + t 0 1 2 f (X(s)) ds + ∞ k=1 [f (X(τ k )) -f (X(τ k -))] I {τ k ≤t} . (4)
The generator A of the Brownian motion process is Af (x) = 1 2 f (x); define the jump operator B by Bf (y, z) = f (z) -f (y). First taking expectations, then dividing by t and letting t → ∞ in (4) results in lim sup

t→∞ t -1 E t 0 Af (X(s)) ds + ∞ k=1 I {τ k ≤t} Bf (X(τ k -), X(τ k )) = 0; (5)
note that the boundedness of X(t) along with f ∈ C 2 (R) implies both the transversality condition lim t→∞ t -1 E[f (X(t))] = 0 holds and that the stochastic integral exists and has mean 0. (The same argument applies by taking the limit inferior in (4), so in fact, left-hand side of ( 5) is a limit.) The fact that f and f are continuous and bounded on [-M, M ] means that

Af (x) µ (t k ) 0 → Af (x) µ 0 (dx)
and Bf (y, z) µ

(t k ) 1 (dy × dz) → B(y, z) µ 1 (dy × dz)
and hence (5) can be written in terms of these measures as

Af (x) µ 0 (dx) + Bf (y, z) µ 1 (dy × dz) = 0. ( 6
)
The restricted impulse control problem is therefore imbedded in the linear program of minimizing (3) over pairs of measures (µ 0 , µ 1 ) satisfying the constraints (6) for every f ∈ D. Since there may be pairs (µ 0 , µ 1 ) which do not correspond to any (τ, Y ) ∈ A 1 , it follows that the value of the linear program is a lower bound on the minimal long-run average cost of the impulse control problem. Observe also that by further restricting the collection of functions for which the constraint is required to be satisfied, the corresponding "auxiliary" linear program may have even more feasible pairs and hence will provide an even lower bound on the value of the impulse control problem. These observations are summarized in the following theorem.

Theorem 2 Let V denote the optimal value of the long-term average impulse control problem, V lp denote the optimal value of the linear program which seeks to minimize (3) over measures satisfying [START_REF] Kurtz | Existence of Markov controls and characterization of optimal Markov controls[END_REF] and V aux be the optimal value of an auxiliary linear program which limits (6) to a smaller collection of functions f ∈ D 1 for some D 1 ⊂ D. Then V aux ≤ V lp ≤ V.

Partial solution: First auxiliary linear program

It would be helpful to reduce the complexity of the linear program by reducing the number of constraints. We first consider the constraints [START_REF] Kurtz | Existence of Markov controls and characterization of optimal Markov controls[END_REF]. The intuition is rather straightforward. Consider a function φ for which Aφ ≡ -1. Then since µ 0 is a probability measure for each feasible pair (µ 0 , µ 1 ), the identity (6) becomes {|z|≤|y|} Bφ(y, z) µ 1 (dy × dz) = 1.

(

The general solution to the equation Aφ ≡ -1 is φ(x) = -x 2 + ax + b, in which a, b are constants. We shall select a solution φ so that Bφ(x, y) = E x [τ y ], where |x| < |y| and τ y := inf{t ≥ 0 :

|X(t)| = |x + W (t)| = |y|}.
To this end, we notice that the function u(x) := -x 2 + y 2 solves the boundary value problem

Au(x) = -1, x ∈ (-|y|, |y|), u(-|y|) = u(|y|) = 0. ( 8
)
Therefore the optional sampling theorem implies that for any x ∈ (-|y|, |y|) and t ≥ 0, we have

E x [u(X(t ∧ τ y ))] = u(x) + E x t∧τy 0 Au(X(s)) ds = u(x) -E x [t ∧ τ y ].
Since X(t ∧ τ y ) is bounded and τ y < ∞ a.s., utilizing the boundary conditions in (8), letting t → ∞ in the above equation yields E

x [τ y ] = u(x) = -x 2 + y 2 .
Therefore by selecting φ(x) = -x 2 , x ∈ R, the term Bφ(x, y) in [START_REF] Kurtz | Stationary solutions and forward equations for controlled and singular martingale problems[END_REF] gives the expected time it takes the Brownian motion process starting at x to hit the set {-|y|, |y|}.

In a similar manner, by considering the boundary value problem

Au(x) = -c 0 (x), x ∈ (-|y|, |y|), u(-|y|) = u(|y|) = 0, (9)
it follows that the function g 0 (x) := -1 6 x 4 satisfies

Bg 0 (x, y) = g 0 (y) -g 0 (x) = E x τy 0 c 0 (X(s))ds , |x| < |y|. (10)
The following proposition establishes the required identity.

Proposition 3 Let (τ, Y ) ∈ A 1 and let X denote the resulting impulse controlled process. Recall {t j : j ∈ N} is a set of times such that

t -1 j c 0 dµ (tj ) 0 + c 1 dµ (tj ) 1 → J(τ, Y ).
Let (µ 0 , µ 1 ) be a weak limit of (µ

(tj ) 0 , µ (tj )
1 ) as j → ∞. Then

c 0 (x) µ 0 (dx) = Bg 0 (y, z) µ 1 (dy × dz). (11)
Proof. Without loss of generality, assume that µ (tj ) 0

⇒ µ 0 and similarly µ (tj ) 1 ⇒ µ 1 . Using g 0 in (4) and taking expectations yields for each t j

E x0 [g 0 (X(t j ))] = g 0 (x 0 ) + E x0 tj 0 Ag 0 (X(s)) ds + ∞ k-0 I {τ k ≤tj } Bg 0 (X(τ k -), X(τ k )) .
Since for (τ, Y ) ∈ A 1 , X(t) remains bounded, dividing by t j , using the definitions of µ [c 1 (y, z) + Bg 0 (y, z)] µ 1 (dy × dz)

S.t. Bφ(y, z) µ 1 (dy × dz) = 1. ( 12
)
Notice the support of each feasible µ 1 is in the set {(y, z) ∈ R 2 : |z| ≤ |y|}. Thus Bφ(y, z) = y 2 -z 2 ≥ 0 and the constraint of ( 12) implies that Bφ(y, z) = y 2 -z 2 is a probability density for every feasible measure µ 1 of ( 12). We emphasize that the auxiliary linear program includes the cost of any impulse policy (τ, Y ) ∈ A 1 . These policies are not required to be of feedback type or even stationary. It is only required that the resulting controlled process remain bounded.

Observe that the constraint does not impose any mass restrictions on the set {(y, z) : Bψ(y, z) = 0}. However, the goal is to minimize the objective function and since the impulse cost function c 1 > k 1 > 0 any mass placed on this set will only increase the cost. We may therefore restrict the optimization to those measures µ 1 having support in {(y, z) : Bφ(y, z) > 0} = {(y, z) : |z| < |y|}. As a result the objective function can be rewritten as

[c 1 (y, z) + Bg 0 (y, z)] µ 1 (dy × dz) = c 1 (y, z) + Bg 0 (y, z) Bφ(y, z) • Bφ(y, z) µ 1 (dy × dz)
and so the problem reduces to the minimization of

F (y, z) = c 1 (y, z) + Bg 0 (y, z) Bφ(y, z) (13)
over {(y, z) : |z| ≤ |y|}. Observe that F (-y, -z) = F (y, z).

Remark 4

The minimization of F in (13) has a very familiar interpretation. Let y and z be such that 0 ≤ |z| < |y|. Let τ y = inf{t ≥ 0 : |X(t)| = |y|}. Recall from (10) that Bg 0 (y, z) represents the expected running cost for the cycle [0, τ y ). Now consider the impulse policy in which impulses occur only when the process X hits either y or -y and then jumps to z or -z, respectively. The symmetry of the fixed cost function c 1 means that c 1 (y, z) = c 1 (-y, -z) and this cost is assessed at the end of the cycle. Note also Bφ(y, z) gives the expected time it takes the Brownian motion starting from -z or z to reach -y or y. Hence the function F represents the ratio of the expected cost per cycle over the expected cycle length taken for such impulse control policies. The significance of this reformulation is that the linear programming imbedding allows arbitrary impulse policies in the class A 1 yet the resulting nonlinear optimization corresponds to minimizing the cost over a subclass of these policies.

We now determine an optimal impulse control policy.

Theorem 5 There exist values y * > z * > 0 such that

F (y * , z * ) = F (-y * , -z * ) = inf (y,z)∈R F (y, z).
Proof. First since Bφ(y, z) = y 2 -z 2 and Bg 0 (y, z) = 1 6 (y 4 -z 4 ), the function F (y, z) = k1+k2|y-z| y 2 -z 2 + 1 6 (y 2 +z 2 ). Observe that as y → ∞ or z → ∞ or |y-z| → 0, F (y, z) → ∞. Hence F achieves its minimal value at some point (y * , z * ) which by optimality requires y * and z * to have the same sign.

Theorem 6 An optimizing pair (y * , z * ) having positive components is the level set of the function h(x) = k2 2x + x 2 3 at the level F (y * , z * ). Proof. Consider the function F on the set {(y, z) : 0 < z < y} so |y -z| = y -z. the first-order optimality conditions are

[y 2 -z 2 ](k 2 + (2/3)y 3 ) -[k 1 + k 2 (y -z) + (1/6)(y 4 -z 4 )](2y) [y 2 -z 2 ] 2 = 0, ( 14
) [y 2 -z 2 ](-k 2 -(2/3)z 3 ) + [k 1 + k 2 (y -z) + (1/6)(y 4 -z 4 )](2z) [y 2 -z 2 ] 2 = 0, ( 15
)
which are satisfied for pairs (y, z) such that

F (y, z) = h(y) = h(z), (16)
where the function h is defined in the statement of the theorem. We have h (x) = 4x 3 -3k2 6x 2

. When x > 0, one observes that h strictly decreases from +∞ until it reaches a minimum at 3 3k 2 /4 after which it strictly increases to +∞. For x < 0, h is strictly decreasing from +∞ to -∞. The level sets of h consist of either a single value with x < 0 when the level lies below the minimum of h over the positive reals or three points with x < 0 < z < y when the level is above this minimum value.

It is still necessary to connect the optimal value of the linear program [START_REF] Piunovskiy | Examples in Markov Decision Processes[END_REF] to the optimal value of the long-term average impulse control problem. Examine the relative magnitudes of the roots given in the proof of Theorem 6; by definition 0 < z < y. Observe the root x < 0 is such that |x| > y and hence x < -y < -z.

Theorem 7 The optimal long-term average cost for the restricted impulse control problem is F (y * , z * ) and an optimal impulse control policy (τ * , Y * ) is defined by:

τ * 1 = 0, Y * 1 = z * -x 0 , τ * k = inf{t > τ * k-1 : X(t-) = ±y * }, Y * k = sgn(X(τ * k -))z * -X(τ * k -), k ≥ 2. ( 17
)
Remark 8 Due to the nature of the long-term average criterion, many other optimal policies exist. In fact, any impulse control policy (τ, Y ) ∈ A 1 can be used for a finite length of time so long as after some point the process is in the interval [-y * , y * ] and the policy of jumping to z * when the process hits y * and jumping to -z * at the time of hitting -y * is adopted.

Proof. Consider the impulse control policy (17) and observe that (τ * , Y * ) ∈ A 1 . Also for each t > 0 by (2), µ

1 has its support on the set {(x 0 , z * ), (-y * , -z * ), (y * , z * )} but the limiting measure µ 1 only has mass on {(-y * , -z * ), (y * , z * )}. As a result, J(τ * , Y * ) = F (y * , z * ) and hence by Theorem 2 is optimal.

General solution

The solution obtained in the previous section does not include impulse control policies which allow the process X to become unbounded in either direction. It is therefore necessary to show that such policies cannot provide a lower cost.

Theorem 9 The optimal long-term average cost over the collection of all admissible impulse control policies is F (y * , z * ) and any policy which eventually only impulses when the process reaches ±y * , with impulses respectively to ±z * , is optimal.

Proof. Let (τ, Y ) be an arbitrary admissible impulse control policy in A and let X denote the resulting process. We shall prove that F * is a lower bound on J(τ, Y ) and since F * is achieved by the impulse control policy (17), this will establish the result.

We use a localization argument. For each n ∈ N, define the stopping time σ n = inf{t ≥ 0 : |X(t)| ≥ n} and to simplify notation, let F * = F (y * , z * ). We consider the function f (x) = F * φ(x) -g 0 (x) = Taking expectations on both sides, and rearranging the terms, it follows that where the last inequality follows from the observation that f (x) = 1 6 x 4 -F * x 2 ≥ K > -∞ for some constant K. Letting n → ∞, we know σ n → ∞ almost surely so the monotone convergence theorem yields

F * E x0 [t ∧ σ n ] ≤ f (x 0 ) -E x0 [f (X(t ∧ σ n ))] + E x0
F * t ≤ f (x 0 ) -K + E x0

  in the compact set R := {(y, z) : |z| ≤ |y| ≤ M }. Similarly, each µ (t) 0 has support in the compact interval [-M, M ]. As a result, these collections are tight and hence relatively compact. Now notice the objective function (1) can be expressed as

in ( 2 )

 2 and letting j → ∞ establishes the result. We are now ready to define the first auxiliary linear program. Restrict the constraint to the single function φ(x) = -x 2 and use Proposition 3 to rewrite the objective function. The resulting linear program is

1 6 x 4 -I k=1 I

 4k=1 F * x 2 , x ∈ R. Due to the choices of φ and g 0 , we have Af (x) = c 0 (x) -F * , and for any |z| < |y|,Bf (y, z) = F * (φ(z) -φ(y)) -(g 0 (z) -g 0 (y)) ≤ c 1 (y, z) + g 0 (z) -g ( y) φ(z) -φ(y) • (φ(z) -φ(y)) -(g 0 (z) -g 0 (y)) = c 1 (y, z).Then applying Itô's formula yieldsf (X(t ∧ σ n )) = f (x 0 ) + {τ k ≤t∧σn} Bf (X(τ k -), X(τ k )) ≤ f (x 0 ) + t∧σn 0 [c 0 (X(s)) -F * ] ds + t∧σn 0 f (X(s)) dW (s) + ∞ {τ k ≤t∧σn} c 1 (X(τ k -), X(τ k )).

Ic 1 (

 1 {τ k ≤t∧σn} c 1 (X(τ k -), X(τ k )) ≤ f (x 0 ) -K + E X(τ k -), X(τ k ))I {τ k ≤t∧σn} ,

t 0 c

 0 0 (X(s)) ds+ ∞ k=0 c 1 (X(τ k -), X(τ k ))I {τ k ≤t} .Dividing by t and letting t → ∞, we obtainF * ≤ lim sup t→∞ t -1 E x0 t 0 c 0 (X(s)) ds + ∞ k=1 c 1 (X(τ k -), X(τ k ))I {τ k ≤t} = J(τ, Y ).

This research was supported in part by National Science Foundation under grant DMS-1108782 and by grant award 246271 from the Simons Foundation.


Kurt Helmes
email: helmes@wiwi.hu-berlin.de

Richard H Stockbridge
email: stockbri@uwm.edu

Chao Zhu
email: zhu@uwm.edu

Impulse control of standard Brownian motion: Discounted criterion

Keywords: impulse control, discounted criterion, infinite dimensional linear programming, expected occupation measures

This paper examines the impulse control of a standard Brownian motion under a discounted criterion. In contrast with the dynamic programming approach, this paper first imbeds the stochastic control problem into an infinite-dimensional linear program over a space of measures and derives a simpler nonlinear optimization problem that has a familiar interpretation. Optimal solutions are obtained for initial positions in a restricted range. Duality theory in linear programming is then used to establish optimality for arbitrary initial positions.

Introduction

When one seeks to control a stochastic process and every intervention incurs a strictly positive cost, one must select a sequence of separate intervention times and amounts. The resulting stochastic problem is therefore an impulse control problem in which the decision maker seeks to either maximize a reward or minimize a cost. This paper continues the examination of the impulse control of Brownian motion. It considers a discounted cost criterion while a companion paper [START_REF] Helmes | Impulse control of standard Brownian motion: Long-term Average Criterion[END_REF] studies the long-term average criterion. The aim of the paper is to illustrate a solution approach which first imbeds the stochastic control problem into an infinite-dimensional linear program over a space of measures and then reduces the linear program to a simpler nonlinear optimization. Contrasting with the long-term average paper, the dependence of the value function on the initial position of the process requires the use of duality in linear programming to obtain a complete solution.

Impulse control problems have been extensively studied using a quasi-variational approach; now classical works include [START_REF] Bensoussan | Impulse control and quasi-variational inequalities[END_REF][START_REF] Harrison | Impulse control of brownian motion[END_REF] while the recent paper [START_REF] Dai | Optimal control of Brownian inventory models with convex inventory cost, Part 2: Discount-optimal controls[END_REF] examines a Brownian inventory model. This paper extends a linear programming approach used on optimal stopping problems [START_REF] Helmes | Construction of the value function and optimal rules in optimal stopping of one-dimensional diffusions[END_REF]. See [START_REF] Helmes | Impulse control of standard Brownian motion: Long-term Average Criterion[END_REF] for additional references.

Let W be a standard Brownian motion process with natural filtration {F t }. An impulse control policy consists of a pair of sequences (τ, Y ) := {(τ k , Y k ) : k ∈ N} in which τ k is the {F t }-stopping time of the kth impulse and the F τ kmeasurable variable Y k gives the kth impulse size. The sequence {τ k : k ∈ N} is required to be non-decreasing, a natural assumption in that intervention k + 1 must occur no earlier than intervention k. For a policy (τ, Y ), the impulsecontrolled Brownian motion process is given by

X(t) = x 0 + W (t) + ∞ k=1 I {τ k ≤t} Y k .
The goal is to control the (discounted) second moment of X subject to (discounted) fixed and proportional costs for interventions. Let (τ, Y ) be an impulse control policy. Define c 0 (x) = x 2 . Let k 1 > 0 denote the fixed costs incurred for each intervention and let k 2 ≥ 0 be a cost proportional to the size of the intervention. Define the impulse cost function c 1 (y, z) = k 1 + k 2 |z -y|, in which y denotes the pre-jump location of X (typically far from 0) and z denotes the post-jump location of X which is thought to be close to 0. Let α > 0 denote the discount rate. The objective function is

J(τ, Y ; x 0 ) = E x0 ∞ 0 e -αs c 0 (X(s)) ds + ∞ k=1 I {τ k <∞} e -ατ k c 1 (X(τ k -), X(τ k )) . (1)
The controller must balance the desire to keep the process X near 0 so as to have a small second moment against the desire to limit the number and/or sizes of interventions so as to have a small impulse cost. Since the goal is to minimize the objective function, impulse control policies having J(τ, Y ; x 0 ) = ∞ are undesirable. We therefore restrict attention to the impulse policies for which J(τ, Y ; x 0 ) is finite. Denote this class of admissible controls by A.

We make five important observations about impulse policies. Firstly, "0impulses" which do not change the state only increase the cost so can be excluded from consideration. Secondly, the symmetry of the dynamics and costs means that any impulse (τ k , Y k ) which would cause sgn(X(τ k )) = -sgn(X(τ k -)) on a set of positive probability will have no greater cost (smaller cost when k 2 > 0) by replacing the impulse with one for which X(τ k ) = sgn(X(τ k -))|X(τ k )|. Thus we can also restrict analysis to those policies for which all impulses keep the process on the same side of 0. Next, any policy (τ, Y ) with lim k→∞ τ k =: τ ∞ < ∞ on a set of positive probability will have infinite cost so for every admissible policy τ k → ∞ a.s. as k → ∞. Next let (τ, Y ) be a policy for which there is some k such that τ k = τ k+1 on a set of positive probability. Again due to the presence of the fixed intervention cost k 1 , the total cost up to time τ k+1 will be at least k 1 E[e -ατ k I(τ k = τ k+1 )] smaller by combining these interventions into a single intervention on this set. Hence we may restrict policies to those for which τ k < τ k+1 a.s. for each k.

The final observation is similar. Suppose (τ, Y ) is a policy such that on a set G of positive probability τ k < ∞ and |X(τ k )| > |X(τ k -)| for some k. Consider a modification of this impulse policy and resulting process X which simply fails to implement this impulse on G. Define the stopping time σ = inf{t > τ k : |X(t)| ≤ | X(t)|}. Notice that the running costs accrued by X over [τ k , σ) are smaller than those accrued by X. Finally, at time σ, introduce an intervention on the set G which moves the X process so that X(σ) = X(σ). This intervention will incur a cost which is smaller than the cost for the process X at time τ k . As a result, we may restrict the impulse control policies to those for which every impulse decreases the distance of the process from the origin.

Restricted problem and measure formulation

The solution of the impulse control problem is obtained by first considering a subclass of the admissible impulse control pairs. Condition 1 Let A 1 ⊂ A be those policies (τ, Y ) such that the resulting process X is bounded; that is, for (τ, Y ) ∈ A 1 , there exists some M < ∞ such that |X(t)| ≤ M for all t ≥ 0.

Note that for each M > 0, any impulse control which has the process jump closer to 0 whenever |X(t-)| = M is in the class A 1 so this collection is non-empty. The bound is not required to be uniform for all (τ, Y ) ∈ A 1 . The restricted impulse control problem is one of minimizing J(τ, Y ; x 0 ) over all policies (τ, Y ) ∈ A 1 .

We capture the expected behavior of the process and impulses with discounted measures. Let (τ, Y ) ∈ A 1 be given and consider f ∈ C 2 (R). Then upon letting t → ∞ after taking expectations, the general Dynkin's formula results in

f (x 0 ) = E x0 ∞ 0 e -αs [αf (X(s)) -(1/2)f (X(s))] ds + E x0 ∞ k=0 I {τ k <∞} e -ατ k [f (X(τ k -)) -f (X(τ k ))] , (2)
in which the transversality condition lim t→∞ E x0 [e -αt f (X(t))] = 0 follows from the boundedness of X. Note the generator of the Brownian motion process is

Af (x) = (1/2)f (x). To simplify notation, define Bf (y, z) = f (y) -f (z).
Define the discounted expected occupation measure µ 0 and the discounted impulse measure

µ 1 such that for each G, G 1 , G 2 ⊂ R, µ 0 (G) = E x0 ∞ 0 e -αs I G (X(s)) ds µ 1 (G 1 × G 2 ) = E x0 ∞ k=0 I {τ k <∞} e -ατ k I G1×G2 (X(τ k -), X(τ k )) . (3)
Notice that the total mass of µ 0 is 1/α while µ 1 is a finite measure since J(τ, Y ; x 0 ) is finite. Rewriting the objective function and Dynkin's formula in terms of these measures imbeds the impulse control problem in the linear program

     Min. c 0 dµ 0 + c 1 dµ 1 S.t. (αf -Af ) dµ 0 + Bf dµ 1 = f (x 0 ), ∀f ∈ C 2 . (4)
We now wish to introduce an auxiliary linear program derived from ( 4) which only has the µ 1 measure as its variable and has fewer constraints. Define φ(x) = e - √ 2αx and ψ(x) = e √ 2αx . Notice that φ is a strictly decreasing solution while ψ is a strictly increasing solution of the homogeneous equation αf -Af = 0. For each (τ, Y ) ∈ A 1 , the resulting process X is bounded so we can use both φ and ψ in (2). This results in the two constraints

Bφ(y, z) µ 1 (dy × dz) = φ(x 0 ) and Bψ(y, z) µ 1 (dy × dz) = ψ(x 0 ) (5)
which only constrain the measure µ 1 . Note that the monotonicity and positivity of both φ and ψ require the support of µ 1 to be such that the two integrals in (5) are positive. We can also take advantage of the symmetry inherent in the problem. Define p 0 (x) = cosh( √ 2αx). Then averaging the two constraints (5

) yields Bp 0 (y, z) p 0 (x 0 ) µ 1 (dy × dz) = 1. Using g 0 (x) = (αx 2 + 1)/α 2 in (2)
, where again the boundedness of X implies that the transversality condition is satisfied, yields

E x0 ∞ 0 e -αs c 0 (X(s)) ds = αx 2 0 + 1 α 2 -E x0 ∞ k=0 I {τ k <∞} e -ατ k Bg 0 (X(τ k -), X(τ k )) . (6)
Let [c 1 -Bg 0 ] denote the sum of the two functions c 1 and Bg 0 . Using (6) in (1) establishes that

J(τ, Y ; x 0 ) = α x 2 0 +1 α 2 + E x0 ∞ k=0 I {τ k <∞} e -ατ k [c 1 -Bg 0 ](X(τ k -), X(τ k ))
and hence that

J(τ, Y ; x 0 ) = α x 2 0 +1 α 2 + [c 1 -Bg 0 ](y, z) µ 1 (dy × dz) (7)
so the objective function value only depends on the measure µ 1 . Since the objective function for each (τ, Y ) ∈ A 1 has the affine term g 0 (x 0 ), it may be ignored for the purposes of optimization but it must be included to obtain the correct value for the objective function. Now form the auxiliary linear program

     Min. [c 1 -Bg 0 ](y, z)] µ 1 (dy × dz) S.t. Bp 0 (y, z) p 0 (x 0 ) µ 1 (dy × dz) = 1. (8) Let V 1 (x 0
) denote the value of the impulse control problem over policies in A 1 , V lp denote the value of (4) and V aux denote the value of (8). The following proposition is immediate.

Proposition 2 V aux (x 0 ) ≤ V lp (x 0 ) ≤ V 1 (x 0 ).
Remark 3 Our analysis will also involve other auxiliary linear programs as well. One will replace the single constraint in (8) with the pair of constraints (5) while another will limit the constraints in (4) to a single function. Each auxiliary program will provide a lower bound on V lp (x 0 ) and hence on V 1 (x 0 ).

Nonlinear optimization and partial solution

Recall, the admissible impulse policies can be (and are) limited to those for which impulses move X closer to the origin. As a result, the integrand Bp 0 > 0 and the constraint of (8) implies that the feasible measures µ 1 of (8) are those for which Bp 0 /p 0 (x 0 ) is a probability density. For a feasible µ 1 , let μ1 be the probability measure Bp0 p0(x0) µ 1 . Thus we can write the objective function as

[c 1 -Bg 0 ] dµ 1 = c 1 -Bg 0 Bp 0 dμ 1 p 0 (x 0 ).
Since the goal is to minimize the cost, a lower bound is given by the minimal value of F scaled by the constant p 0 (x 0 ), where

F (y, z) := c 1 (y, z) -Bg 0 (y, z) Bp 0 (y, z) .
Moreover, should the infimum be attained at some pair (y * , z * ), then the probability measure μ1 (•) putting unit point mass on (y * , z * ) would achieve the lower bound and identify an optimal µ 1 measure for the auxiliary linear program. To solve the stochastic problem, one would need to connect the measure µ 1 back to an admissible impulse control policy in the class A 1 in such a way that the resulting µ 1 measure would be given by (3).

Remark 4

The objective function p 0 (x 0 )F has a natural interpretation. First observe that Bp 0 (y, z) = cosh( √ 2αy) -cosh( √ 2αz) so

p 0 (x 0 )F (y, z) = [c 1 (y, z) -Bg 0 (y, z)] • cosh( √ 2αx 0 ) cosh( √ 2αy) • ∞ n=0 cosh( √ 2αz) cosh( √ 2αy) n .
It can be shown that the first fraction gives the expected discount for the time it takes X to reach {±y} when starting at x 0 . The ratio cosh( √ 2αz) cosh( √ 2αy) then gives the expected discount for the time it takes X to again reach {±y} but this time starting at ±z so the sum represents the expected discounting for infinitely many cycles. By symmetry, the initial term gives the cost for impulsing from ±y to ±z along with the second moment. The minimization therefore optimizes the expected cost over a particular class of impulse policies. We emphasize that the linear program imbedding is not restricted to these policies.

Proposition 5 There exists pairs (y * , z * ) and (-y * , -z * ) such that

F (y * , z * ) = F (-y * , -z * ) = inf (y.z):|z|≤|y| F (y, z). (9)
Moreover, the minimizing pair (y * , z * ) having nonnegative components is unique.

Proof. First observe

F (y, z) = k 1 + k 2 |y -z| + (z 2 -y 2 )/α cosh( √ 2αy) -cosh( √ 2αz)
so there exists some pairs (y, z) for which F (y, z) < 0 since the difference of the quadratic terms is negative and will dominate the constant and linear terms in the numerator. A straightforward asymptotic analysis show that F (y, z) is asymptotically nonnegative when y → ∞, z → ∞ or |y -z| → 0. Therefore F achieves its minimum at some point (y * , z * ). Notice that F is symmetric about 0 in that F (-y, -z) = F (y, z) so it is sufficient to analyze F on the domain 0 ≤ z ≤ y. The first-order optimality conditions on F are

0 = ∂F ∂y (y * , z * ) = (k 2 -2y * /α)[cosh( √ 2α y * ) -cosh( √ 2α z * )] [cosh( √ 2α y * ) -cosh( √ 2α z * )] 2 - √ 2α [k 1 + k 2 (y * -z * ) + (z 2 * -y 2 * )/α] sinh( √ 2α y * ) [cosh( √ 2α y * ) -cosh( √ 2α z * )] 2 , 0 = ∂F ∂z (y * , z * ) = (-k 2 + 2z * /α)[cosh( √ 2α y * ) -cosh( √ 2α z * )] [cosh( √ 2α y * ) -cosh( √ 2α z * )] 2 + √ 2α [k 1 + k 2 (y * -z * ) + (z 2 * -y 2 * )/α] sinh( √ 2α z * ) [cosh( √ 2α y * ) -cosh( √ 2α z * )] 2 .
The minimizing pair (y * , z * ) will be interior to the region since ∂F ∂z (y * , 0) = -k2 cosh( √ 2α y * )-1 < 0. Simple algebra now leads to the following systems of nonlinear equations for (y * , z * ):

k 2 α -2y * = α √ 2α sinh( √ 2α y * ) • F (y * , z * ), k 2 α -2z * = α √ 2α sinh( √ 2α z * ) • F (y * , z * ). ( 10
)
The fact that the minimal value of F is negative implies y * > z * > k 2 α/2. Solving for -α √ 2α F (y * , z * ) in each equation shows that at an optimal pair (y * , z * ),

-α √ 2α F (y * , z * ) = 2z * -k 2 α sinh( √ 2α z * ) = 2y * -k 2 α sinh( √ 2α y * ) .
A straightforward analysis of the function h

(x) = [2x -k 2 α]/ sinh( √ 2α x) on the domain [k 2 α/2, ∞)
shows that the level sets of h consist of two-point sets and so on the region 0 ≤ z ≤ y, the pair (y * , z * ) is unique. Now that the lower bound given in ( 9) is determined, it is important to connect an optimizing µ * 1 with an admissible impulse control policy (τ, Y ) ∈ A 1 . The existence of two minimizing pairs (y * , z * ) and (-y * , -z * ) allows many auxiliary-LP-feasible measures µ 1 to place point masses at these two points and still achieve the lower bound. This observation leads to a solution to the restricted stochastic impulse control problem.

Theorem 6 Let (y * , z * ) be the pair having positive components that minimizes F as identified in Proposition 5. Consider initial positions -y * ≤ x 0 ≤ y * . Define the impulse control policy (τ * , Y * ) as follows:

τ * 1 = inf{t ≥ 0 : X(t-) = ±y * } and Y * 1 = sgn(X(τ * 1 -)) • z * -X(τ * 1 -) and for k = 2, 3, 4, . . ., define τ * k = inf{t > τ k-1 : X(t-) = ±y * } and Y * k = sgn(X(τ * k -)) • z * -X(τ * k -). Then (τ * , Y *
) is an optimal impulse control pair for the restricted stochastic impulse control problem and the corresponding optimal value is

V 1 (x 0 ) = αx 2 0 +1 α 2 + F (y * , z * ) • cosh( √ 2α x 0 ). ( 11
)
Proof. The measure µ * 1 defined from (τ * , Y * ) using ( 3) is concentrated on the two points (-y * , -z * ) and (y * , z * ). Since the process resulting from the admissible impulse control pair (τ * , Y * ) remains bounded, conditions (5) can be used to obtain the masses:

µ * 1 (-y * , -z * ) = φ(x 0 )[ψ(y * ) -ψ(z * )] -ψ(x 0 )[φ(y * ) -φ(z * )] [φ(-y * ) -φ(-z * )][ψ(y * ) -ψ(z * )] -[ψ(-y * ) -ψ(-z * )][φ(y * ) -φ(z * )] , µ * 1 (y * , z * ) = ψ(x 0 )[φ(-y * ) -φ(-z * )] -φ(x 0 )[ψ(-y * ) -ψ(-z * )] [φ(-y * ) -φ(-z * )][ψ(y * ) -ψ(z * )] -[ψ(-y * ) -ψ(-z * )][φ(y * ) -φ(z * )] . Recall φ(x) = e - √ 2α x and ψ(x) = e √ 2α
x so φ(-x) = ψ(x) and ψ(-x) = φ(x). As a result these expressions simplify to

µ * 1 (-y * , -z * ) = φ(x 0 )[ψ(y * ) -ψ(z * )] -ψ(x 0 )[φ(y * ) -φ(z * )] [ψ(y * ) -ψ(z * )] 2 -[φ(y * ) -φ(z * )] 2 , µ * 1 (y * , z * ) = ψ(x 0 )[ψ(y * ) -ψ(z * )] -φ(x 0 )[φ(y * ) -φ(z * )] [ψ(y * ) -ψ(z * )] 2 -[φ(y * ) -φ(z * )] 2 .
It is now straightforward to verify that J(τ * , Y * ; x 0 ) equals the value in (11).

Full solution

Theorem 6 solves the problem for initial positions x 0 with |x 0 | ≤ y * . The issue is now one of determining the optimal value and an optimal impulse control pair when |x 0 | > y * . From an intuitive point of view, |x 0 | < y * has an optimal control which waits until the state process first hits ±y * before having an impulse so one might expect an impulse to occur immediately when |x 0 | ≥ y * . Since two impulses at the same instant are no better than one, one would anticipate that the after-jump location might be z ∈ (-y * , y * ). The cost of an immediate jump from x 0 to z followed by using an optimal impulse control is

g(z) := αx 2 0 + 1 α 2 + k 1 + k 2 (x 0 -z) + z 2 -x 2 0 α + V 1 (z) = αz 2 + 1 α 2 + k 1 + k 2 (x 0 -z) + V 1 (z).
Solving g (z) = 0 to find a minimizer results in

0 = -k 2 + 2z/α + √ 2α F (y * , z * ) sinh( √ 2α z),
which is the first order condition (10) for which both y * and z * are solutions. An impulse to y * would be followed by an immediate jump to z * and incur two fixed costs whereas a single jump directly to z * would cost less. This line of reasoning indicates that a single jump to z * could be an optimal initial impulse.

The goal is to verify that this intuitive reasoning is correct. Define

V (y) =        k 1 + k 2 (|y| -z * ) + V 1 (-z * ), y ≤ -y * , V 1 (y), -y * ≤ y ≤ y * , k 1 + k 2 (y -z * ) + V 1 (z * ), y ≥ y * .
For |y| > y * , the function V is the cost associated with the process starting at initial position y, having an instantaneous jump from y to sgn(y) z * and then using the optimal impulse control policy of Theorem 6 thereafter. The following lemma is fairly straightforward so its proof is left to the reader.

Lemma 7 V ∈ C 1 (R) ∩ C 2 (R\{±y * }).
The function V therefore has sufficient regularity to use in [START_REF] Dai | Optimal control of Brownian inventory models with convex inventory cost, Part 2: Discount-optimal controls[END_REF]. We now consider the new auxiliary linear program

           Min. R\{±y * } c 0 (x) µ 0 (dx) + R 2 c 1 (y, z) µ 1 (dy × dz) S.t. R\{±y * } [α V (x) -A V (x)] µ 0 (dx) + R 2 [ V (y) -V (z)] µ 1 (dy × dz) = V (x 0 ) (12)
and its dual (having sole variable w)

       Max. V (x 0 ) • w S.t. α V (x) -A V (x) • w ≤ c 0 (x), x = ±y * , V (y) -V (z) • w ≤ c 1 (y, z), ∀y, z ∈ R. (13)
Observe that each linear program has feasible points with costs that are finite. A straightforward weak duality argument therefore shows that each value of (13) corresponding to a feasible variable w is no greater than any value of (12) for a feasible pair of measures and hence the value of ( 13) is a lower bound on the value of the restricted impulse control problem. Since V (x 0 ) > 0, one seeks as large a positive value as possible for w.

Theorem 8 The optimal value of (13) is V (x 0 ) which is achieved when w * = 1.

Proof. By symmetry, it is sufficient to examine x, y, z ≥ 0. Notice that for 0 ≤ x < y * , α V (x) -A V (x) = x 2 = c 0 (x) ≥ 0 and hence the dual variable w cannot exceed 1. The question is whether w = 1 is feasible for (13) so examine the rest of the constraints with w = 1. For x > y * , A V (x) = 0 so the first constraint of (13) requires

0 ≤ x 2 -α (k 1 + k 2 (x -z * ) + V 1 (z * )) = x 2 -αV 1 (y * ).
Since the right-hand expression is an increasing function for x ∈ [k 2 α/2, ∞), it suffices to verify its nonnegativity with x = y * :

0 ≤ y 2 * -αV 1 (y * ) = y 2 * -α αy 2 * + 1 α 2 + F (y * , z * ) cosh( √ 2α y * ) = - 1 α + [2y * -k 2 α] cosh( √ 2α y * ) √ 2α sinh( √ 2α y * )
in which (10) is used to obtain the last expression. This inequality can be rewritten as tanh(

√ 2α y * ) √ 2α ≤ y * -k 2 α/2. (14)
Since (y * , z * ) is a minimizing pair of the function F , (14) holds and the first family of constraints of (13) is satisfied with w = 1. Consider now the second family of constraints with w = 1. There are several cases to examine. When 0 ≤ y ≤ z, monotonicity of V on this range shows the condition is trivially satisfied. Next, for 0 ≤ z ≤ y ≤ y * , the constraint can be rewritten as

V 1 (y) ≤ k 1 + k 2 (y -z) + V 1 (z).
The right-hand expression gives the cost of an immediate jump from y to z followed by an optimal impulse control policy thereafter whereas the left-hand side gives the optimal cost. Hence this inequality is satisfied. Now consider y * ≤ z < y and observe that V (y) -V (z) = k 2 (y -z) < k 1 + k 2 (y -z). Finally, for 0 ≤ z < y * < y and again using the definition of V , the second set of constraints in (13) is equivalent to

k 1 + k 2 (y -z * ) + V 1 (z * ) ≤ k 1 + k 2 (y -z) + V 1 (z) or equivalently k 2 (y -y * ) + V 1 (y * ) = k 2 (y -y * ) + [k 1 + k 2 (y * -z * ) + V 1 (z * )] ≤ k 2 (y -y * ) + [k 1 + k 2 (y * -z) + V 1 (z)].
This last inequality is true by the optimality of both the pair (y * , z * ) and the function V 1 on [-y * , y * ] since the bracketed quantity on the right-hand side gives the cost associated with an initial impulse to z from y * along with optimal impulse control policy starting from z. Thus the second family of constraints in (13) hold when w = 1.

We now have the following result.

Theorem 9 Let (y * , z * ) be the optimizing pair for F having positive components. Define the impulse control policy (τ * , Y * ) as follows;

τ * 1 = inf{t ≥ 0 : |X(t-)| ≥ y * } and Y * 1 = sgn(X(τ * 1 -)) • z * -X(τ * 1 -)
and for k = 2, 3, 4, . . ., define

τ * k = inf{t > τ k-1 : X(t-) = ±y * } and Y * k = sgn(X(τ * k -)) • z * -X(τ * k -).
Then (τ * , Y * ) is an optimal impulse control pair for the restricted stochastic impulse control problem and the corresponding optimal value is V (x 0 ).

Proof. The particular choice of (τ * , Y * ) implies

V (x 0 ) ≤ V lp (x 0 ) ≤ V 1 (x 0 ) ≤ J(τ * , Y * ) = V (x 0 ).

Solution for general admissible impulse controls

The solution of Section 2.2 is restricted to those impulse control policies under which the process X remains bounded. It is necessary to show that no lower cost can be obtained by any policy which allows the process to be unbounded.

Theorem 10 The impulse control policy (τ * , Y * ) of Theorem 9 is optimal in the class of all admissible policies and V (x 0 ) is the optimal value.

Proof. This argument establishes that V (x 0 ) is a lower bound on J(τ, Y ; x 0 ) for every admissible impulse control policy. Theorem 9 then gives the existence of an optimal policy whose cost equals the lower bound.

Choose (τ, Y ) ∈ A and let X be the resulting controlled process. Suppose there exists some K > 0 such that lim inf t→∞ E x0 [e -αt V (X(t)] ≥ K. Note that lim inf Using this estimate in (1) shows J(τ, Y ; x 0 ) = ∞. Now suppose J(τ, Y ; x 0 ) < ∞ so lim inf t→∞ E x0 [e -αt V (X(t))] = 0. Then there exists a sequence {t j : j ∈ N} such that lim j→∞ E x0 [e -αtj V (X(t j ))] = 0. Note that | V | ≤ k 2 so t 0 e -αs V (X(s)) dW (s), t ≥ 0, is a martingale. Thus the dual constraints, in conjunction with the finiteness of the expected cost, implies that Dynkin's formula holds when t = t j for each j. Hence -E x0 e -αtj V (X(t j )) Letting j → ∞, an application of the monotone convergence theorem on the first expectation and the convergence to 0 of second expectation establishes that V (x 0 ) is a lower bound on the expected cost J(τ, Y ; x 0 ).

t→∞E

  x0 e -αt V (X(t)) = lim inf t→∞ E x0 e -αt V (X(t))I {|X(t)|≥y * } so the linearity of V on {x : |x| ≥ y * } implies that E x0 [|X(t)|I {|X(t)|≥y * } ] is asymptotically bounded below by Ke αt as t → ∞. Hence by Jensen's inequality for > 0 and t large,

E

  x0 X 2 (t) ≥ E x0 [ |X(t)|I {|X(t)|≥y * } ] 2 ≥ K 2 e 2αt -.

V

  (x 0 ) = E x0 tj 0 e -αs [α V (X(s)) -A V (X(s))] ds -E x0 e -αtj V (X(t j )) + E x0 ∞ k=0 I {τ k ≤tj } e -ατ k B V (X(τ k -), X(τ k )) ≤ E x0 tj 0 e -αs c 0 (X(s)) ds + ∞ k=0 I {τ k ≤t k } e -ατ k c 1 (X(τ k -), X(τ k ))

This research was supported in part by National Science Foundation under grant DMS-1108782 and by grant award 246271 from the Simons Foundation.


Peter Philip
email: philip@math.lmu.de

PDE-driven shape optimization: Numerical investigation of different descent directions and projections using penalization and regularization

Keywords: shape optimization, optimal control, fixed domain method, elliptic partial differential equation, numerical simulation

We consider shape optimization problems with elliptic partial differential state equations. Using regularization and penalization, unknown shapes are encoded via shape functions, turning the shape optimization into optimal control problems for the unknown functions. The method is designed to allow topological changes in a natural way. Based on convergence and differentiability results, numerical algorithms are formulated, using different descent directions and projections. The algorithms are assessed in a series of numerical experiments, applied to an elliptic PDE arising from an oil industry application with two unknown shapes, one giving the region where the PDE is solved, and the other determining the PDE's coefficients.

Introduction

We study an elliptic shape optimization problem motivated by the oil industry application studied in [START_REF] Woo | A direct tracking method for a grounded conductor inside a pipeline from capacitance measurements[END_REF], where one aims at monitoring the interior of a pipeline. The cross section through the pipeline is modeled by a set D ⊆ R 2 , consisting of a liquid region Ω (such that D \ Ω represents air) that is part oil (region O ⊆ Ω) and part water (Ω \ O). The shape optimization needs to reconstruct the parts O and Ω from given measurements y d of a quantity (e.g. voltage) taken in a region E adjacent to the boundary of D. Thus, we are led to the following problem, previously formulated in [START_REF] Neittaanmäki | Fixed domain approaches in shape optimization problems with Dirichlet boundary conditions[END_REF]:

min Ω, O 1 2 E |y -y d | 2 dx + 1 2 E |∇y -∇y d | 2 dx , (1a)
subject to In particular, depending on ξ, (1c) can mean homogeneous or nonhomogeneous Dirichlet conditions. The employed method (previously used, e.g., in [START_REF] Mäkinen | On a fixed domain approach for a shape optimization problem[END_REF][START_REF] Neittaanmäki | Fixed domain approaches in shape optimization problems with Dirichlet boundary conditions[END_REF][START_REF] Philip | A penalization and regularization technique in shape optimization problems[END_REF]) is based on the introduction of shape functions g and p, defined on D and encoding the unknown sets Ω and O, respectively, and on a technique for the approximation and regularization of characteristic functions. Assuming g, p : D → R to be continuous, the corresponding sets are obtained via

Ω a1 χO + a2(1 -χO) ∇y • ∇v dx + Ω b1 χO + b2(1 -χO) y v dx = Ω f v dx , ∀ v ∈ H 1 0 (Ω), (1b)
y -ξ ∈ H 1 0 (Ω), (1c)
Ω g = int{x ∈ D : g(x) ≥ 0}, O p = int{x ∈ D : p(x) ≥ 0} (2)
(Ω g and O p then are open Caratheodory sets, not necessarily connected). Enforcing the constraints E ⊆ O ⊆ Ω translates into the set of admissible pairs

U ad := (g, p) ∈ C(D) × C(D) : g ≥ p on D and p ≥ 0 on E . (3)
Let H : R → R denote the Heaviside function. Then H(g), H(p) : D → R are the characteristic functions of Ω g and O p , respectively. We use the differentiable regularization of the Heaviside function given by

H ε (r) :=          1 for r ≥ 0, ε(r + ε) 2 -2r(r + ε) 2 ε 3 for -ε < r < 0, 0 for r ≤ -ε, (4)
and obtain the following regularized fixed domain approximation of (1) (with the abovementioned constraints E ⊆ O ⊆ Ω):

min (g,p)∈U ad 1 2 E |y ε -y d | 2 dx + 1 2 E |∇y ε -∇y d | 2 dx , ( 5a
) D a 1 H ε (p) + a 2 (1 -H ε (p)) ∇y ε • ∇v + b 1 H ε (p) + b 2 (1 -H ε (p)) y ε v dx + 1 ε D 1 -H ε (g) y ε v dx = D f v dx , ∀ v ∈ H 1 0 (D), (5b)
y ε -ξ ∈ H 1 0 (D). ( 5c
)
For ε > 0 small, the penalty term with the 1/ε in (5b) forces the state y ε to be close to 0 outside Ω g (for precise, rigorous versions of this statement see [5, Th. 2], [8, Th. 2.2, Th. 3.1]). It is noted that, even though the above-described method encodes Ω g and O p , in fact, as level sets of the functions g, p, respectively, our method is essentially different from the well-known level set method of [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], since no time dependence of the functions g, p and no time evolution of the corresponding open sets Ω g , O p is assumed. In particular, we do not need to solve any Hamilton-Jacobi equations in the process. We now consider a triangular finite element partition of D, D =

T h ∈T h T h ,
h > 0, assuming the grid in D, restricted to E, provides a finite element mesh in E as well. Let V h , V h denote the corresponding finite element spaces in D constructed with piecewise affine continuous functions (with 0 trace on ∂D for elements of V h ). Defining

U h ad := (g, p) ∈ V h × V h : g ≥ p on D and p ≥ 0 on E , (6)
the discretized form of (5) reads min

(g h ,p h )∈U h ad j(g h , p h ) := 1 2 E |y ε,h -y d,h | 2 dx + 1 2 E |∇y ε,h -∇y d,h | 2 dx , (7a) D a 1 H ε (p h ) + a 2 (1 -H ε (p h )) ∇y ε,h • ∇v h dx + D b 1 H ε (p h ) + b 2 (1 -H ε (p h )) y ε,h v h dx + 1 ε D 1 -H ε (g h ) y ε,h v h dx = D f h v h dx , ∀ v h ∈ V h ⊆ H 1 0 (D), (7b)
y ε,h -ξ h ∈ V h ⊆ H 1 0 (D), (7c)
where (7b) constitutes the equation for the discretized state

y ε,h ∈ V h ⊆ H 1 (D), f h ∈ V h and ξ h ∈ V h are
) ∈ V h × V h in the direction (w h , u h ) ∈ V h × V h as 1 ε D H ε (g h )w h y ε,h q ε,h dx - D (b 1 -b 2 )H ε (p h )u h y ε,h q ε,h dx - D (a 1 -a 2 )H ε (p h )u h ∇y ε,h ∇q ε,h dx , (8)
where

q ε,h ∈ V h ⊆ H 1 0 (D) is the solution to the adjoint equation D a 1 H ε (p h ) + a 2 (1 -H ε (p h )) ∇q ε,h • ∇v h dx + D b 1 H ε (p h ) + b 2 (1 -H ε (p h )) q ε,h v h dx + 1 ε D 1 -H ε (g h ) q ε,h v h dx = E (y ε,h -y d,h )v h dx + σ E (∇y ε,h -∇y d,h ) • ∇v h dx , ∀ v h ∈ V h , (9)
and the direction of steepest descent (w d,0 , u d,0 ) is given by

w d,0 := -(1/ε) H ε (g h ) y ε,h q ε,h , u d,0 := H ε (p h ) (a 1 -a 2 ) ∇y ε,h • ∇q ε,h + H ε (p h ) (b 1 -b 2 ) y ε,h q ε,h . (10)
While ( 10) is difficult to use in practise as H ε (g h ) and H ε (p h ) are typically nonzero only in a small neighborhood of ∂Ω g h and ∂O p h , respectively, multiplication by nonnegative coefficients yields the following alternative descent directions (w d,1 , u d,1 ) and (w d,2 , u d,2 ), without such support restrictions:

w d,1 := -y ε,h q ε,h , u d,1 := (a 1 -a 2 ) ∇y ε,h • ∇q ε,h + (b 1 -b 2 ) y ε,h q ε,h , (11)
w d,2 := w d,1 χ S , u d,2 := u d,1 χ S , (12)
where χ S denotes the characteristic function of

S := {x ∈ D : w d,1 (x) ≥ u d,1 (x)} ∪ {x ∈ E : w d,1 (x) ≥ 0 and u d,1 (x) ≥ 0}.
(13) Using ( 12) has the advantage of maintaining the conditions g ≥ p on D and p ≥ 0 on E.

For the numerical results presented below, we employ four variants of an algorithm of gradient with projection type making use of (approximations of) the descent directions [START_REF] Shewchuk | Delaunay refinement algorithms for triangular mesh generation[END_REF] and [START_REF] Woo | A direct tracking method for a grounded conductor inside a pipeline from capacitance measurements[END_REF]. The two variants based on (11) will be called A1a and A1b, whereas the variants based on (12) will be called A2a and A2b. Moreover, variants A1a and A2a will use the admissible set U h ad of ( 6), whereas A1b and A2b will use the modification

U h ad,b := (g, p) ∈ U h ad : |∇g|, |∇p| ≤ 1 , (14)
enforcing uniformly bounded gradients for the shape functions, a condition suggested by the results of [START_REF] Philip | A penalization and regularization technique in shape optimization problems[END_REF]Sec. 5]. Variant A1a was previously considered in [START_REF] Philip | A penalization and regularization technique in shape optimization problems[END_REF]; the remaining three variants are new. The four algorithms are formulated below in Sec. 2.1, with a description of their implementation in Sec. 2.2. Numerical experiments comparing the performance of the four variants are then presented in Sec. 3.

Numerical algorithms

Formulation

In preparation for the numerical experiments of Sec. 3, we formulate the employed algorithms. As indicated at the end of the Introduction, we use four variants of the algorithm previously published in [START_REF] Philip | A penalization and regularization technique in shape optimization problems[END_REF], built on the earlier version of [START_REF] Neittaanmäki | Fixed domain approaches in shape optimization problems with Dirichlet boundary conditions[END_REF]. As mentioned above, we denote the four variants by A1a, A1b, A2a, A2b, where A1a is precisely the algorithm used in [START_REF] Philip | A penalization and regularization technique in shape optimization problems[END_REF]. A1a and A1b use the descent direction [START_REF] Shewchuk | Delaunay refinement algorithms for triangular mesh generation[END_REF] for line searches followed by a projection step, whereas A2a and A2b use [START_REF] Woo | A direct tracking method for a grounded conductor inside a pipeline from capacitance measurements[END_REF], which has the advantage of remaining within U h ad during the line search, avoiding the projection. Variants A1b and A2b project into the smaller space U h ad,b of (14) after each line search. The algorithms consist of the following Steps (1) -( 7):

(1): Set n := 0 and choose initial shape functions (g h,0 , p h,0 ) ∈ U h ad .

(2): Compute the solution to the state equation y n := θ ε,h (g h,n , p h,n ), where θ ε,h : V h × V h → V h denotes the control-to-state operator corresponding to (7b), (7c); and compute the solution to the corresponding adjoint equation q n := θε,h (y n ), where θε,h : V h → V h , y ε,h → q ε,h , denotes the solution operator corresponding to [START_REF] Philip | Shape optimization via control of a shape function on a fixed domain: Theory and numerical results[END_REF].

(3): Compute the descent direction (w n d , u n d ), where w n d = w d,1 (y n , q n ) and u n d = u d,1 (y n , q n ) according to [START_REF] Shewchuk | Delaunay refinement algorithms for triangular mesh generation[END_REF] for A1a and A1b, whereas w n d = w d,2 (y n , q n ) and u n d = u d,2 (y n , q n ) according to [START_REF] Woo | A direct tracking method for a grounded conductor inside a pipeline from capacitance measurements[END_REF] For all the numerical examples discussed below, we stopped the iteration and returned (g h,fin , p h,fin ) := (g h,n+1 , p h,n+1 ) if |j(g h,n , p h,n ) -j(g h,n+1 , p h,n+1 )| < 10 -5 AND g h,n -g h,n+1 2 < 10 -3 AND p h,n -p h,n+1 2 < 10 -3 , where |j(g h,n , p h,n ) -j(g h,n+1 , p h,n+1 )|/|j(g h,n+1 , p h,n+1 )| is used for |j(g h,n , p h,n )j(g h,n+1 , p h,n+1 )| if |j(g h,n+1 , p h,n+1 )| > 1 and analogous for g h,n and p h,n .

Implementation

The state equations as well as the adjoint equations that need to be solved numerically during the above algorithms are discretized linear elliptic PDE with Dirichlet boundary conditions. The numerical solution is obtained via a finite volume scheme [START_REF] Philip | Analysis, optimal control, and simulation of conductive-radiative heat transfer[END_REF]Sec. 4]. More precisely, the software WIAS-HiTNIHS 1 , originally designed for the solution of more general PDE occurring when modeling conductive-radiative heat transfer and electromagnetic heating [START_REF] Geiser | Numerical simulation of temperature fields during the sublimation growth of SiC single crystals, using WIAS-HiTNIHS[END_REF], has been 1 High Temperature Numerical Induction Heating Simulator.

adapted for use in the present context. WIAS-HiTNIHS is based on the program package pdelib [START_REF] Fuhrmann | pdelib: An open modular tool box for the numerical solution of partial differential equations. Design patterns[END_REF], it employs the grid generator Triangle [START_REF] Shewchuk | Delaunay refinement algorithms for triangular mesh generation[END_REF] to produce constrained Delaunay triangulations of the domains, and it uses the sparse matrix solver GSPAR [START_REF] Grund | Direct linear solvers for vector and parallel computers[END_REF] to solve the linear system arising from the finite volume scheme.

The numerical scheme yields discrete y n and q n (cf.

Step (2) of the above algorithms), defined at each vertex of the triangular discrete grid, interpolated piecewise affine, i.e. affinely to each triangle of the discrete grid. In consequence, the shape functions g h,n and p h,n are piecewise affine as well. Where integrals of these piecewise affine functions need to be computed (e.g. in Step (7) of the algorithms), they are computed exactly. A golden section search [10, Sect. 10.2] is used to numerically carry out the minimization (15). Note that the minimization (15) is typically nonconvex and the golden section search will, in general, only provide a local min λ n .

For some numerical examples, the stated initial shape functions (g h,0 , p h,0 ) are merely piecewise continuous (cf. the Introduction and [START_REF] Philip | Shape optimization via control of a shape function on a fixed domain: Theory and numerical results[END_REF]) and, thus, not in U h ad . However, the stated (g h,0 , p h,0 ) are only used to determine the values g h (x h i ), p h (x h i ), at the nodes x h i of the triangulation T h , and the resulting affinely interpolated functions are in U h ad . Moreover, in Step (3) of the algorithms, approximations of the descent directions are used, as for the gradients nodewise averages are computed, that are then affinely interpolated, and the conditions of (13) are enforced nodewise and affinely interpolated. In principle, in might occur that the approximated direction is no longer a descent direction, but such a case was not observed during our numerical experiments.

Numerical experiments

Numerical experiments with precomputed optimum

The numerical computations of the present section employ the circular fixed domain D := (x 1 , x 2 ) :

x 2 1 + x 2 2 < 1 ⊆ R 2 with fixed subdomain E := {(x 1 , x 2 ) ∈ D : |x 1 | > 3 4 , |x 2 | < 1 2 } ⊆ D (note E has two connected components).
We use a fixed triangular grid provided by Triangle [START_REF] Shewchuk | Delaunay refinement algorithms for triangular mesh generation[END_REF], consisting of 24458 triangles. The used regularization parameter is ε = 10 -5 (cf. [START_REF] Philip | A penalization and regularization technique in shape optimization problems[END_REF][START_REF] Philip | Shape optimization via control of a shape function on a fixed domain: Theory and numerical results[END_REF]). The settings for the remaining given quantities are a 1 := 1, a 2 := 10, b 1 := 1, b 2 := 10, f (x 1 , x 2 ) := 5, ξ(x 1 , x 2 ) := 2. The cost functional j as in (7a) depends on the given function y d,h . For the first set of numerical results, we precompute y d,h := y ε,h numerically as the solution to the state equation (7b), (7c), using The computed y d,h with the corresponding Ω g and O p is depicted in Fig. 1. Using the precomputed y d,h has the advantage that we actually know y d,h together with g h , p h as in ( 16) provides an absolute minimum in the following numerical examples, employing the cost functional j of (7a) with the precomputed y d,h from above. A series of four numerical experiments was conducted, all using the initial shape functions g h,0 (x 1 , x 2 ) := 1, p h,0 (x 1 , x 2 ) := 1 (see Fig. 2). We refer to the experiments as 1:A1a, 1:A1b, 1:A2a, and 1:A2b, depending on which variant of the the algorithm of Sec. 2.1 was used. The results for 1:A2a and 1:A2b are shown in Fig. 3. The final state and shapes for 1:A1a and 1:A1b were very similar to those of 1:A2a, with slightly higher final costs (0.69 and 0.28, respectively). All variants reduce the cost significantly, all resulting local minima being different and different from the absolute min. Variant A2a gives the best result, whereas A2b results in the highest final cost, where one also observes a symmetry breaking due to the discrete grid. Actually, for A2b, after the first line search, the cost is 0.29 with shapes resembing the final shapes of the other variants, but the projection of Step (6) can subsequently result in a cost increase, which occurs in this example.

g h (x 1 , x 2 ) :=      -1 if (x 1 , x 2 ) / ∈ E and (x 1 , x 2 ) -(-1, 0) 2 < 0.4, -1 if (x 1 , x 2 ) / ∈ E and (x 1 , x 2 ) -(1, 0) 2 < 0.4, 1 otherwise, (16a)
p h (x 1 , x 2 ) := 1 in E, -1 in D \ E. (16b) . y

Numerical experiments without precomputed optimum

In contrast to the experiments of the previous section, we now consider a setting, where we are no longer in the situation of a known precomputed optimum. For the following numerical results, the fixed domain D is still the unit disk as in Sec. 3.1. However, the fixed subdomain E is now at the bottom of D, defined by

E := {(x 1 , x 2 ) ∈ D : x 2 < -0.7} ⊆ D.

Conclusions

In a series of numerical experiments, we have studied four variants of an algorithm of gradient with projection type for shape optimization problems driven by elliptic PDE. The variants used different descent directions and different sets of admissible shape functions. Except in one situation, all variants were effective in finding local minima of significantly reduced costs. However, it did depend on both the equation and on the initial condition, which variant showed the best performance. Thus, further research seems warranted to further evaluate and improve the different variants.

  where the sets E ⊆ O ⊆ Ω ⊆ D ⊆ R 2 all are bounded and open, D also connected; χ O denotes the characteristic function of O; a 1 , a 2 , b 1 , b 2 > 0, f ∈ L 2 (D), y d , ξ ∈ H 1 (D) all given.

  ) := π h,b ( gh,n , ph,n ), where π h,b denotes the projection π h,b : U h ad → U h ad,b , obtained by dividing gh,n and ph,n by α, defined as the max of the max-norms of |∇g| and |∇p|, in case α > 1. (7): RETURN (g h,fin , p h,fin ) := (g h,n+1 , p h,n+1 ) if the change of g, p and/or the change of j(g, p) are below some prescribed tolerance parameter. Otherwise: Increment n, i.e. n := n + 1 and GO TO (2).

  given suitable approximations of f and ξ, respectively, g h , p h ∈ V h are discretized shape functions corresponding to discretizations of Ω g and O p , respectively, and y d,h is a suitable given continuous and piecewise affine approximation of y d .

	Similar to [5, Prop. 2] and [8, Cor. 5.3], one obtains the directional derivative
	of the cost functional (g, p) → j(g, p) with j as in (7a) at (g h , p h

  for A2a and A2b. (4): Set gh,n := g h,n +λ n w n d and ph,n := p h,n +λ n u n d , where λ n ≥ 0 is determined via line search, i.e. as a solution to the minimization problem For A2a and A2b, set ( gh,n , ph,n ) := (g h,n , ph,n ) (no projection is necessary to obtain ( gh,n , ph,n ) ∈ U h ad ); for A1a and A1b, set ( gh,n , ph,n ) := π h (g h,n , ph,n ), where π h denotes the projection π h: V h × V h → U had , obtained by first setting gh,n (x h For A1a and A2a, set (g h,n+1 , p h,n+1 ) := ( gh,n , ph,n ) (no second projection necessary); for A1b and A2b, set (g h,n+1 , p h,n+1

	min λ≥0	j(g h,n + λ w n d , p h,n + λ u n d ).	(15)
	(5):

i ) := max{0, gh,n (x h i )} and ph,n (x h i ) := max{0, ph,n (x h i )} for each node x h i of the triangulation T h such that x h i ∈ E, and second setting ph,n (x h i ) := min{p h,n (x h i ), gh,n (x h i )} for every node x h i of the triangulation T h . (6):

  Precomputed y d,h used in all experiments of Sec. 3.1 (left, isolevels spaced at 0.2), obtained as the solution to the state equation (7b), (7c); with the corresponding Ωg (middle) and Op (right). Initial state, shapes used in all experiments of Sec. 3.1. Left: State isolevels spaced at 0.2. Middle: Shape Ωg. Right: Shape Op. Cost: j(g h,0 , p h,0 ) = 19.0.

						g				.	p				.
					1.00					1.00					1.00
					0.50					0.50					0.50
					0.00					0.00					0.00
					-0.50					-0.50					-0.50
					-1.00					-1.00					-1.00
	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00
	Fig. 1. . y	g				.	p				.
					1.00					1.00					1.00
					0.50					0.50					0.50
					0.00					0.00					0.00
					-0.50					-0.50					-0.50
					-1.00					-1.00					-1.00
	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00
	Fig. 2.

  The numerical computations employ a Intermediate and final state, shapes for shape optimizations of Sec. 3.2, i.e. for Experiments 2:A1a (1st, 2nd row), 2:A1b (3rd, 4th row), and 2:A2b (5th, 6th row). Left: State isolevels spaced at 0.2. Middle: Shapes Ωg. Right: Shapes Op. Costs at shown intermediate states are 2.64 for 2:A1a, 1.74 for 2:A1b, 1.70 for 2:A2b. Final costs j(g h,fin , p h,fin ) are 1.72 for 2:A1a, 1.68 for 2:A1b, 1.69 for 2:A2b. Number of line searches for intermediate and for final state: 8 and 31 for 2:A1a, 2 and 29 for 2:A1b, 3 and 38 for 2:A2b.

	y				.	g				.	p				.
					1.00					1.00					1.00
					0.50					0.50					0.50
					-0.00					-0.00					-0.00
					-0.50					-0.50					-0.50
					-1.00					-1.00					-1.00
	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00
	y				.	g				.	p				.
					1.00					1.00					1.00
					0.50					0.50					0.50
					-0.00					-0.00					-0.00
					-0.50					-0.50					-0.50
					-1.00					-1.00					-1.00
	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00
	y				.	g				.	p				.
					1.00					1.00					1.00
					0.50					0.50					0.50
					-0.00					-0.00					-0.00
					-0.50					-0.50					-0.50
					-1.00					-1.00					-1.00
	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00
	y				.	g				.	p				.
					1.00					1.00					1.00
					0.50					0.50					0.50
					-0.00					-0.00					-0.00
					-0.50					-0.50					-0.50
					-1.00					-1.00					-1.00
	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00
	y				.	g				.	p				.
					1.00					1.00					1.00
					0.50					0.50					0.50
					-0.00					-0.00					-0.00
					-0.50					-0.50					-0.50
					-1.00					-1.00					-1.00
	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00
	y				.	g				.	p				.
					1.00					1.00					1.00
					0.50					0.50					0.50
					-0.00					-0.00					-0.00
					-0.50					-0.50					-0.50
					-1.00					-1.00					-1.00
	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00	-1.00	-0.50	0.00	0.50	1.00
	Fig. 5.

fixed triangular grid provided by Triangle [START_REF] Shewchuk | Delaunay refinement algorithms for triangular mesh generation[END_REF], consisting of 24623 triangles. The parameter settings are as in Sec. 3.1, except for f (x 1 , x 2 ) := 10(x 2 1 + x 2 2 ) + 5. The cost functional is as in (7a) with y d,h (x 1 , x 2 ) := x 1 +x 2 . A series of four numerical experiments was conducted, all using the initial shape functions g h,0 (x 1 , x 2 ) :=

(see Fig. 4).

We refer to the experiments as 2:A1a, 2:A1b, 2:A2a, and 2:A2b, depending on which variant of the the algorithm of Sec. 2.1 was used. Results are shown in Fig. 5, except for 2:A2a, which converged after 10 line searches to a local min almost identical to the initial condition. All other variants reduce the cost significantly, all resulting local minima being different. Here, the lowest final cost is achieved for variant A1b. One notices significant changes in shapes (including topology changes) during the optimizations, where very different shapes can result in nearly identical costs. As in Sec. 3.1, symmetry breaking can occur due to the discrete grid.


Tom Streubel
email: streubel@math.hu-berlin.de

Andreas Griewank
email: griewank@math.hu-berlin.de

Manuel Radons
email: radons@math.hu-berlin.de

Jens-Uwe Bernt
email: berntj@math.hu-berlin.de

Representation and analysis of piecewise linear functions in abs-normal form

Keywords: piecewise linearization (PL), algorithmic differentiation (AD), equation solving, semi-smooth Newton, smooth dominance, complementary piecewise linear system (CLP), linear complementarity

It follows from the well known min/max representation given by Scholtes in his recent Springer book, that all piecewise linear continuous functions y = F (x) : R n → R m can be written in a so-called abs-normal form. This means in particular, that all nonsmoothness is encapsulated in s absolute value functions that are applied to intermediate switching variables zi for i = 1, . . . , s. The relation between the vectors x, z, and y is described by four matrices Y, L, J, and Z, such that

This form can be generated by ADOL-C or other automatic differentation tools. Here L is a strictly lower triangular matrix, and therefore zi can be computed successively from previous results. We show that in the square case n = m the system of equations F (x) = 0 can be rewritten in terms of the variable vector z as a linear complementarity problem (LCP). The transformation itself and the properties of the LCP depend on the Schur complement S = L -ZJ -1 Y .

Introduction

Via algorithmic differentiation it is possible to calculate directional derivatives from evaluation procedures of vector valued functions simultaneously with their evaluation at a base point x 0 . These evaluations are exact within the limitations of machine precision. An evaluation procedure is a composition of so called elementary functions, which are aggregated as a library in their symbolic form and thus make up the atomic constituents of complex functions. Basically the selection of elementary functions for the library is arbitrary, as long as they comply with assumption (ED) (elementary differentiability, in [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF]), meaning that they are at least once Lipschitz-continuously differentiable. In the literature (see e.g. [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF][START_REF] Naumann | The art of differentiating computer programs: an introduction to algorithmic differentiation[END_REF]) the following collection is suggested as the quasi-standard for a library: Φ = {+, -, * , /, sin, cos, tan, cot, exp, log, . . . } Common software packages such as ADOL-C provide tools for the algorithmic differentiation of functions composed from the contents of this collection.

But many practical problems and most algorithms are not smooth everywhere and thus cannot be modelled via a library that consists solely of a set of functions that comply with (ED). More specifically one is likely to encounter standard functions of computer arithmetic, that are not globally differentiable, e.g. abs, max and min. Since max(x, y, z) ≡ max(max(x, y), z), max(x, y) ≡ 0.5 * (x + y + abs(x -y)) max and min can be expressed in terms of the absolute value function. As shown, this reformulation of max and min provides us with a very practical handle on the representation of piecewise linearity, since Scholtes proved in [START_REF] Scholtes | Introduction to piecewise differentiable equations[END_REF], that any scalar-valued, real piecewise linear function f : R n → R can be expressed as a finite nesting of max and min comparisons of linear functions. Here and throughout we use linear in the sense of affine, i.e. allow a constant increment.

Generally any one dimensional piecewise linear function f : R → R can be expressed in terms of absolute values. For a given set of points {(x i , y i ) : i = 0, . . . , n}, where x 0 < x 1 < • • • < x n , two outer slopes s 0 , s n+1 and n inner slopes s i = (y i -y i-1 )/(x i -x i-1 ), we obtain the formula

y = 1 2 y 0 + s 0 (x -x 0 ) + n i=0 (s i+1 -s i ) abs(x -x i ) + y n + s n+1 (x -x n )
where the two linear functions at the beginning and the end can be combined to [y 0 -s 0 x 0 + y n -s n+1 x n + (s 0 + s n+1 )x]/2. This might be helpful for the purpose of implementation. For example with a < b ∈ R we obtain the cut-off function

f (x) = max(a, min(x, b)) = 0.5 * [a + abs(x -a) -abs(x -b) + b]
Similar to linear models of smooth functions, piecewise linearizations can be used to approximate piecewise smooth functions [START_REF] Scholtes | Introduction to piecewise differentiable equations[END_REF]. The aim is to extend the principles and techniques of classic algorithmic differentiation in such a way, that these piecewise linear models can be evaluated with the same efficiency, stability and simplicity of data structures as in the linear case. Since the absolute value function is already piecewise linear, it can be modelled by itself. By proposition 3.1 from [START_REF] Griewank | On stable piecewise linearization and generalized algorithmic differentiation[END_REF] we have for the procedure (introduced in the next chapter) that the error of the piecewise linear approximation is of second order and varies Lipschitz continuously w.r.t. the developing point.

Piecewise linearization and abs-normal form

Example 1. Formula, graph and sequential code instruction of an evaluation procedure:

F (x 1 , x 2 ) = x 1 + |x 1 -x 2 | + |x 1 -|x 2 || x 2 x1 v1 |v1| v4 |v4| v0 |v0| y1 x2 y2 v6 v -1 = x 1 v 0 = x 2 v 1 = v -1 -v 0 v 2 = abs(v 1 ) v 3 = abs(v 0 ) v 4 = v -1 -v 3 v 5 = abs(v 4 ) v 6 = v 5 + v 2 y 1 = v -1 + v 6 y 2 = v 0
The v i are called intermediate values. The indices are in a dependency relation j ≺ i, if there is an edge from v j to v i . In general the values of a sequential code instruction of an evaluation procedure are denoted as a tuple

[v 1-n , v 1-(n-1) , . . . , v 0 , v 1 , v 2 , . . . , v i , . . . , v l ] where v j-n = x j for j = 1, . . . , n v i = ϕ i (v j ) j≺i for i = 1, . . . , l and ϕ ∈ Φ abs = Φ ∪ {abs}
The values of the piecewise linearization can be evaluated simultaneuosly as increments of the function value by the following set of propagation rules [START_REF] Griewank | On stable piecewise linearization and generalized algorithmic differentiation[END_REF] that implicitly defines a second code instruction.

Procedure 1

[∆v 1-n , ∆v 1-(n-1) , . . . , ∆v 0 , ∆v 1 , ∆v 2 , . . . , ∆v i , . . . , ∆v l ] where for j = 1, . . . , n : ∆v j-n = ∆x j and for i = 1, . . . , l :

∆v i = ∆v j ± ∆v k when v i = v j ± v k ∆v i = ∆v j * v k + v j * ∆v k when v i = v j * v k ∆v i = c ij * ∆v j when v i = ϕ(v j )
where ϕ ∈ Φ \ {±, * , abs} and c i,j = ϕ i (v j ) is the local partial derivative

∆v i = abs(v j + ∆v j ) -abs(v j ) when v i = abs(v j )
Then the k-th component of the piecewise linearization is determined by:

y k = v l-m+k + ∆v l-m+k
The overall costs are at most four times of those of a function evaluation [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF]. So far we have a method for a small number of evaluations at some basepoints. But for the purposes of integration, solving ODEs, optimization and solving piecewise linear equation systems (see [START_REF] Griewank | On stable piecewise linearization and generalized algorithmic differentiation[END_REF] and [START_REF] Griewank | T Solving piecewise linear equations in abs-normal form[END_REF]) we need a suitable data structure for a large number of evaluations of a single piecewise linearization. A general nonlinear concept of Barton and Khan from [START_REF] Khan | Evaluating an element of the Clarke generalized Jacobian of a piecewise differentiable function[END_REF] combined with taping technology leads to the abs-normal form:

Definition 1. For Z ∈ R s×n , L ∈ R s×s , J ∈ R m×n , Y ∈ R m×s matrices, where L is of strictly lower triangular form and vectors c ∈ R s , b ∈ R m , the system z y = c b + Z L J Y x |z| (1)
is called abs-normal form. The modulus operation |z| has to be understood componentwise here. An abs-normal form is called simply switched if L = 0.

The components of z can be evaluated successively, since L is a strictly lower triangular matrix. The control flow in the evaluation of the abs-normal form is conveniently characterised by the signature vectors and matrices

σ x ≡ σ z ≡ sign(z) ∈ {-1, 0, 1} s , Σ z = diag(σ z ) ∈ {-1, 0, 1} s×s
In particular we will use throughout the identity |z| = Σ z z. Using this relation we can eliminate z for any given x ∈ R n and obtain the explicit representation

F (x) = y = piecewise constant b + Y Σ σ (I -LΣ σ ) -1 c +J σ • x (2)
where

J σ = J + Y Σ σ (I -LΣ σ ) -1 Z (3)
On the other hand every piecewise linear function in max-min expression can be represented in abs-normal form. Thus the abs-normal form is an equivalent characterization of piecewise linear mappings, which is stable w.r.t to perturbations1 . Each signature vector σ ∈ {-1, 0, 1} s uniquely characterises the polyhedron

P σ = {x ∈ R n | σ x = σ}
The collection of these mutually disjoint and relatively open polyhedra forms a so called polyhedral decomposition or skeleton P of R n . The restriction of F to the closure of any P σ ∈ P is linear. Each P = P σ has a nonempty interior if and only if it is open, in which case we will also refer to as open. By continuity all σ that have no zero components are open, but the converse need not be the case. It can be shown, that the J σ given in (3) are limiting Jacobians in the following sense exactly if σ is open.

For general Lipschitz continuous F it follows from Rademacher's Theorem that it has a Frechet derivative F (x) at all points in a set D F , whose complement has the measure zero. The set of limiting Jacobians at any x 0 ∈ R n is defined as

∂ L F (x 0 ) = lim x→x0 x∈D F F (x) = ∅
and the set of generalized Jacobians in the sense of Clarke as

∂F (x 0 ) = conv(∂ L F (x 0 ))
The definition of ∂ L F (x 0 ) looks quite nonconstructive and in fact there is no general methodology for evaluating limiting Jacobians since the rules for propagating generalized derivatives are only inclusions. Given the abs-normal form one can compute limiting Jacobians that are also generalized Jacobians of the underlying nonlinear functions by a technique called polynomial escape [START_REF] Griewank | On stable piecewise linearization and generalized algorithmic differentiation[END_REF] and [START_REF] Khan | Evaluating an element of the Clarke generalized Jacobian of a piecewise differentiable function[END_REF]. The computational complexity is similar to that of the foward mode in the smooth case. Especially for generalized gradients where m = 1 an adaption of the much cheaper reverse mode is under development.

Throughout the remainder of this paper, we will only consider piecewise linear F in abs-normal form that are square in that m = n. Furthermore we assume w.l.o.g. that the so called smooth part J is nonsingular. If this is not a priori true one can shift terms by using the identity x = abs(x+abs(x))-abs(x). The Schur complement of J within the abs-normal form is given by S = L -ZJ -1 Y . By using the Sherman-Morrison-Woodbury formula we can characterise the nonsingularity of the generalized Jacobian J σ as follows det(J σ ) = det(J) det(I -SΣ σ ), for σ = σ x ∈ {-1, 0, 1} s

Note that the upper half of the abs-normal form, which maps x onto z, need not be surjective. Hence the mapping is maybe partially switched in that some signature vectors σ ∈ {-1, 0, 1} s do not arise as σ x for any x. In other words some P σ might be empty. On the other hand if the linear map Zx is surjective, then the abs-normal form must be totally switched in that all 3 n sign combinations of σ with corresponding nonempty P σ do arise. The following so called complementary piecewise linear mappings are always totally switched, since z ∈ R s becomes independent and ranges over all of R s .

Complementary piecewise linear systems and their relation to LCPs

In contrast the nonsingularity of the smooth part J allows the elimination of x for any given z and y.

y = b + Jx + Y |z| ⇐⇒ x = J -1 (y -b) -J -1 Y |z|
In view of solving F (x) = 0 we can set y = 0 or absorb it into b. Then substitution of x into the upper half yields the complementary piecewise linear mapping

H(z) ≡ (I -SΣ z )z -ĉ, where ĉ ≡ c -ZJ -1 b
The function H : R s → R s is still piecewise linear and has the abs-normal form z

H(z) = 0 -ĉ + I 0 I -S z |z| (5)
whose Schur complement is again S. Since the new L vanishes, the complementary piecewise linear map is always simply switched. Moreover the polyhedral decomposition consists entirely of 2 n open orthants and their faces. As shown in [START_REF] Griewank | T Solving piecewise linear equations in abs-normal form[END_REF] 4) the coherent orientation of F . Whereas the converse need not be true, i.e. F may be coherently oriented but H not. The problem of solving H(z) = 0, for some z ∈ R n can be recast as a linear complementarity problem (LCP). It turns out to have the P -matrix property if and only if H is coherently oriented [START_REF] Rump | Theorems of Perron-Frobenius type for matrices without sign restrictions[END_REF]. The reformulation requires:

M = [I + S] -1 [I -S] ⇐⇒ [I + S] 1 2 (I + M ) = 1 2 ([I + S] + [I -S]) = I ⇐⇒ S = 2(I + M ) -1 -I = (I + M ) -1 (I -M )
Now consider two vectors 0 ≤ u, w ≥ 0, such that z = u -w and u w = 0. Then by the upper half of an abs-normal form for

F u -w = c + Zx + L(u + w) = c + ZJ -1 (y -b) -ZJ -1 Y (u + w) + L(u + w) = ĉ + S(u + w) ⇐⇒ (I -S)u = ĉ + (I + S)w ⇐⇒ w = u -(I + S) -1 ĉ ⇐⇒ w = M u + q
where M ≡ (I + S) -1 (I -S) and q ≡ -(I + S) -1 ĉ. Because of the substitution of x, the solutions of this standard LCP w = q + M u, are solutions of the complementary piecewise linear system H. Any standard LCP w = q + M u, where u, w ≥ 0 and u w = 0, can be rewritten as a complementary piecewise linear equation system as

z = (I + M ) -1 (I -M )|z| -2(I + M ) -1 q
where u = 1 2 (|z| + z) and w = 1 2 (|z| -z). This was proven by Bokhoven in his thesis [START_REF] Van Bokhoven | Piecewise-linear modelling and analysis[END_REF]. To transform the complementary piecewise linear system into an LCP or vice versa one has to compute the Möbius transform of S or M , respectively. This requires in either case at least implicitly a matrix inversion and several multiplications. Therefore we consider methods for directly solving the original and complementary piecewise linear system possibly even avoiding the explicit computation of S = L -ZJ -1 Y .

Solving piecewise linear equation systems

The principal task is to find solutions x ∈ R n , such that F (x) = 0 with piecewise linear F : R n → R n . A possible nonzero right hand side can be absorbed into the vector b as described above.

There are several methods developed and discussed in detail in [START_REF] Griewank | T Solving piecewise linear equations in abs-normal form[END_REF] and [START_REF] Griewank | On stable piecewise linearization and generalized algorithmic differentiation[END_REF]. Some of them solve F (x) = 0 directly, whereas others solve the complementary piecewise linear equation System H(z) = 0. Note that there is a one-to-one solution correspondence between both representations [START_REF] Griewank | T Solving piecewise linear equations in abs-normal form[END_REF]. Now, let us give an overview of some of these methods.

Full-step Newton variants

All continuous piecewise linear functions are known to be semi smooth. Hence the result in [START_REF] Qi | A nonsmooth version of Newton's method[END_REF] ensures local convergence of the full-step iteration

x + = x -J -1 F (x), for J ∈ ∂F (x)
to a solution x * , provided that all generalized Jacobians J ∈ ∂F (x * ) are nonsingular. However this condition need not be satisfied even if F is coherently oriented. Coherent orientation in some vicinity of x * means that all limiting Jacobians J ∈ ∂ L F (x * ) are nonsingular, so that the stronger result from [START_REF] Qi | Nonsmooth equations and smoothing Newton methods[END_REF], where the J are restricted to be limiting Jacobians, is applicable.

It should be noted that both results apply here in a trivial fashion, since convergence in one step must occur from all points x 0 belonging to polyhedra P σ , whose closure contains x * . Of course finding such an initial point x 0 requires to resolve all combinatorial issues in advance.

Hence we are more interested in global convergence results. We can guarantee full step convergence for the restricted generalized Newton method in finitely many steps towards the unique solution, if either of the contractivity conditions

I -J -1 σ J σ < 1, for all σ, σ open or I -J σ J -1 σ < 1, for all σ, σ open
is satisfied w.r.t. to some induced matrix norm. The proof can be found in [START_REF] Griewank | T Solving piecewise linear equations in abs-normal form[END_REF].

Either condition is rather strong and implies bijectivity. In terms of the absnormal form they are implied by the conditions

ρ ≡ Z J -1 Y < 1 -L and ρ (1 -ρ -L )(1 -L ) < 1 2
As we have already noted suitable J σ can be computed from the abs-normal form at reasonable expense.

Naturally the generalized Newton method with or without restriction to limiting Jacobians can also be applied to the complementary piecewise linear system, yielding

z + = z -(I -SΣ z ) -1 H(z) = (I -SΣ z ) -1 ĉ
However, here the local convergence condition that all limiting Jacobians be nonsingular is no weaker than the requirement that all generalized Jacobians be nonsingular. Sufficient for global full-step convergence are either of the following independent conditions

S p < 1 3 or ρ(|S|) < 1 2
where ρ denotes the spectral radius and |S| the componentwise modulus.

If the second condition is satisfied, the calculation can be organized such that the whole solution process requires only 1 3 s 3 operations, just like a Gaussian elimination in the smooth linear case.

Piecewise Newton

Rather than taking full steps based on a local linearization one may restrict steps to stay within the closure of one polyhedron P σ . This requires some pivoting and active set managament familiar from Lemke type algorithms for LCPs. For a comparitive study of the two approaches see the dissertation of T. Munson [START_REF] Munson | Algorithms and environments for complementarity[END_REF].

In [START_REF] Griewank | On stable piecewise linearization and generalized algorithmic differentiation[END_REF] it was observed that coherent orientation implies, that the fibres

[x 0 ] ≡ {x ∈ R n : F (x) = λF (x 0 ), 0 < λ ∈ R}
are bifurcation-free piecewise linear paths for almost all x 0 ∈ R n . Then their closure contains a solution. Even in the case of singular fibres, there are strategies to reduce the residual towards a solution. An implementation is currently under development.

Modulus algorithm

Checking F for surjectivity or openess is NP-hard, because there may be 2 n possible determinants det(J σ ), for σ = σ x . An easier verifiable property is smooth dominance.

Definition 2. F : R n → R n in abs-normal form is called smooth dominant, if for some nonsingular diagonal matrix D and a

p ∈ [1, ∞] DSD -1 p < 1
Smooth dominant abs-normal forms are always injective [START_REF] Griewank | T Solving piecewise linear equations in abs-normal form[END_REF]. Nevertheless there are many practical problems which satisfy this condition.

In [START_REF] Brugnano | Iterative solution of piecewise linear systems[END_REF] Brugnano and Casulli consider unilateral constraints solve max(0, x) + T x = -e/2

where T ∈ R n×n is an irreducible, symmetric, positive semidefinite matrix and

x, e ∈ R n vectors. This class of problems is piecewise linear and its abs-normal forms are smooth dominant. Electrical engineers considered piecewise linear function as models of electrical circuits since the 50's of the last century. For example Bokhoven discussed those models in his dissertation [START_REF] Van Bokhoven | Piecewise-linear modelling and analysis[END_REF] and introduced the iteration z + = S|z| -ĉ whose convergence follows from smooth dominance, by the Banach fix point theorem. In our experience the modulus iteration is robust, but rather slow.

Alternating block Seidel iteration

Another fixed point iteration which has the potential of being significantly faster, is the following block Seidel scheme from [START_REF] Griewank | T Solving piecewise linear equations in abs-normal form[END_REF]. Solving alternatingly the upper half for z and the lower half for x, we obtain z + = h z (h x (z)), where

h z : R n → R s h z (x) = (I -LΣ x ) -1 (c + Zx) h x : R s → R n h x (z) = -J -1 b -J -1 Y Σ z z
The convergence of this method to the unique solution is ensured [START_REF] Griewank | T Solving piecewise linear equations in abs-normal form[END_REF], if

S p ≤ L p + ZJ -1 Y p < 1
for some suitable p where positive diagonal scaling may be applied.

Conclusion and outlook

We gave a short introduction to basic techniques of automatic differentiation and methods for the modelling of piecewise smooth functions via piecewise linearization with a second order error. We also discussed the solvabillity of the resulting equation systems in abs-normal form, by finitely convergent Newton variants or linearly convergent fix point solvers. Currently we are working on hybrid algorithms to obtain stable global and fast local convergence. They will then be used in the inner loop of a piecewise smooth equation solver by successive piecewise linearization. A related task to equation solving are the (un)constrained optimization of piecewise smooth objectives and the numerical integration of initial value problems with Lipschitzian right hand sides. Common utillities for manipulating abs-normal forms are developed as the linear algebra package PLAN-C, which uses abs-normal forms as objects.

Fig. 1 :

 1 Fig. 1: example of a piecewise linear function and its corresponding polyhedral decomposition.

Lemma 1 .

 1 Let M, S ∈ R s×s arbitrary, s.t. (I + S)M = (I -S), then 1. det(I + S) = 0 ⇐⇒ det(I + M ) = 0 2. S = (I + M ) -1 (I -M ) if det(I + M ) = 0 Proof.

  this implies that H is bijective if and only if it is an open map. For general PL functions and in particular the underlying F we only have the chain of implications[START_REF] Scholtes | Introduction to piecewise differentiable equations[END_REF]

F is injective =⇒ F is open =⇒ F is surjective

Furthermore Scholtes has proven in

[START_REF] Scholtes | Introduction to piecewise differentiable equations[END_REF]

that piecewise linear maps are open maps if and only if the determinants of all limiting Jacobians have the same sign (are w.l.o.g. positive). The limiting Jacobians of H are exactly the shifted identities I -SΣ for any Σ = diag(σ) with σ ∈ {-1, 1} s . Consequently, coherent orientation of H occurs if and only if all det(I -SΣ) are positive, which implies by (

Perturbations of the data Z, L, J, Y, c and b preserves the property of being a continuous, piecewise linear abs-normal form, provided L stays strictly lower triangular.