Vyacheslav Maksimov
email: maksimov@imm.uran.ru

Regularized Extremal Shift in Problems of Stable Control

Keywords: parabolic equation, stable control, extremal shift

We discuss a technical approach, based on the method of regularized extremal shift (RES), intended to help solve problems of stable control of uncertain dynamical systems. Our goal is to demonstrate the essence and abilities of the RES technique; for this purpose we construct feedback controller for approximate tracking a prescribed trajectory of an inaccurately observed system described by a parabolic equation. The controller is "resource-saving" in a sense that control resource spent for approximate tracking do not exceed those needed for tracking in an "ideal" situation where the current values of the input disturbance are fully observable.

Introduction

In the present work, the problem of tracking a solution of a system with distributed parameters is discussed. The essence of this problem can be formulated in the following way. A parabolic equation is considered on a given time interval T = [t 0 , ϑ], ϑ < +∞. The solution of this equation w(•) = w(•; v(•)) depends on a time-varying control v = v(•). This solution is inaccurately measured at frequent enough time moments. It is required to organize a control process for the equation by the feedback principle in such a way that it is possible to preserve given properties of the solution. The quality of the solution constructed is estimated by the distance from a given (prescribed, standard) solution x(•). The latter is a solution of the parabolic equation generated by some input u = u(•). The problem in question is treated as the problem of constructing a control v = v(•) providing the retention of the trajectory w(•) = w(•; v(•)) nearby x(•). This is the conceptual statement of the control problem under consideration.

Problem Statement

Let H and V be real Hilbert spaces. The space V is a dense subspace of H and V ⊂ H ⊂ V * algebraically and topologically, (•, •) stands for the inner product in H, •, • stands for the duality relation between V and V * . We consider a system Σ which is described by the parabolic equation ẇ(t) + Aw(t) = Bv(t) + f (t), for a. a. t ∈ T, w(t 0 ) = w 0 .

(1)

Here A : V → V * is a linear continuous (A ∈ L(V ; V * )) and symmetrical operator satisfying (for some c * > 0 and real ω) the coercitivity condition

Aw, w + ω|w| 2 H ≥ c * |w| 2 V ∀y ∈ V, (2)
U is a Hilbert space, f ∈ L 2 (T ; H) is a given function, | • | H , | • | U and | • | V stand
for the norms in H, U and V , respectively, B : U → H is a linear continuous operator (B ∈ L(U ; H)). Let the following condition be fulfilled.

Condition 1. Operator B is invertible. Let w(t 0 ) = w 0 ∈ D(A H ), where D(A H ) = {w ∈ V : A H ∈ H}.
It is known that under such conditions, for any v(•) ∈ L 2 (T ; U ), there exists a unique solution w(•) = w(•; t 0 , w 0 , v(•)) of equation ( 1) with the following properties [START_REF] Bensoussan | Representation and Control of Infinite Dimensional Systems[END_REF]:

w(•) ∈ W (T ) = W 1,2 (T ; H) ∩ L 2 (T ; V ). Here, W 1,2 (T ; H) = {w(•) ∈ L 2 (T ; H) : ẇ(•) ∈ L 2 (T ; H)}, the derivative ẇ(•) is understood in the sense of distributions.
Assume that along with equation ( 1) we have another equation of the same form:

ẋ(t) + Ax(t) = Bu(t) + f (t) for a.a. t ∈ T (3)
with an initial state x(t 0 ) = x 0 ∈ D(A H ). This equation (in what follows, we call it reference) is subject to the action of some reference control u(•) ∈ L 2 (T ; U ). The reference control as well as the corresponding solution x(•) = x(•; t 0 , x 0 , u(•)) of equation ( 3) are a priori unknown. At discrete, frequent enough, time moments τ i ∈ ∆ = {τ i } m i=0 (τ 0 = t 0 , τ m = ϑ, τ i+1 = τ i + δ), the states w(τ i ) = w(τ i ; t 0 , w 0 , v(•)) of equation [START_REF] Bensoussan | Representation and Control of Infinite Dimensional Systems[END_REF] as well as the states x(τ i ) = x(τ i ; t 0 , x 0 , u(•)) of reference equation (3) are measured. The states w(τ i ) are measured with an error. The results of measurements are elements ξ h i ∈ H satisfying the inequalities

|w(τ i ) -ξ h i | H ≤ h, i ∈ [1 : m -1]. (4)
Here, the value h ∈ (0, 1) is the measurement accuracy. It is required to design an algorithm for forming the control v = v h (•) in equation ( 1) allowing us to track the solution x(•) of equation ( 3) by the solution w(•) of equation [START_REF] Bensoussan | Representation and Control of Infinite Dimensional Systems[END_REF]. Thus, we consider the problem consisting in constructing an algorithm, which (on the basis of current measurements of the values w(τ i ) and x(τ i )) forms in real time mode (by the feedback principle) the control v = v h (•) in the righthand part of inequality (1) such that the deviation of w(•) = w(•; t 0 , w 0 , v h (•)) from x(•) = x(•; t 0 , x 0 , u(•)) in metric of the space C(T ; H) ∩ L 2 (T ; V ) is small if the measurement accuracy h is small enough. We also want the constructed algorithm to be resource-saving. This means that the resources of the synthetic

control v = v h (•) (i.e., the value ϑ t0 |v h (τ )| 2 U dτ
) should exceed the resources of the reference control by a small value depending on the measurement accuracy h. This value tends to zero as h tends to zero. Thus, we require the validity of the inequality

ϑ t0 |v h (τ )| 2 U dτ ≤ ϑ t0 |u(τ )| 2 U dτ + ϕ(h), (5)
where ϕ(h) → 0 as h → 0.

In the case when the reference control u as well as the control v in inequality (1) are subject to instantaneous constraints (u ∈ P , v ∈ P , where P ⊂ U is a given bounded and closed set), the problem above can be solved by means of the method of extremal shift [START_REF] Krasovskii | Game-theoretical Control Problems[END_REF]. Namely, if the control v = v h (•) in the right-hand part of ( 1) is calculated by the formula

v h (t) = v(τ i , ξ h i , x(τ i )) = arg min{(x(τ i ) -ξ h i , Bv) : v ∈ P } for t ∈ [τ i , τ i+1 ), (6)
then, as it follows from [START_REF] Vaysburd | Differential pursuit game for systems with distributed parameters[END_REF], for any ε > 0 one can find numbers h 1 > 0 and

δ 1 > 0 such that the inequality sup t∈T |w(t; t 0 , w 0 , v h (•)) -x(t; t 0 , x 0 , u(•))| H ≤ ε is fulfilled if h ∈ (0, h 1 ) and δ ∈ (0, δ 1 )
. The last inequality is valid for any reference control, i.e., for any Lebesgue measurable function u(t) ∈ P for almost all t ∈ T . Here and below, we assume that ω > 0 and

w 0 ∈ D ⊂ V , where D is a bounded set, |w 0 -x 0 | H ≤ h. (7)
Thus, the method of extremal shift allows us to solve the problem of tracking the solution of the reference equation under instantaneous constraints on the controls (v, u ∈ P ). In the present paper, we assume that any function from the space L 2 (T ; U ) can be the admissible control (both reference, u(•), and real, v(•)). No additional information on the functions v(•) and u(•) is required. We construct a corresponding modification of the method of extremal shift, using, according to [START_REF] Osipov | Inverse Problems for Ordinary Differential Equations: Dynamical Solutions[END_REF][START_REF] Maksimov | Dynamical Inverse Problems of Distributed Systems[END_REF][START_REF] Osipov | Methods of Dinamical Reconstruction of Inputs of Controlled Systems[END_REF][START_REF] Maksimov | Dynamical reconstruction of inputs for constraction semigroup systems: the boudary input case[END_REF][START_REF] Maksimov | On reconstruction of unbounded controls in nonlinear dynamical systems[END_REF][START_REF] Maksimov | Dynamical state and control reconstruction for a phase field model. Dynamics of continuous, discrete and impulsive systems[END_REF], the idea of its local regularization. Along with measuring the phase states at discrete time moments (see ( 4)), we also consider the case of "continuous" measuring of the states x(t) and w(t). Namely, it is assumed that, at every time t ∈ T , the phase states of equations ( 1) and (3) are measured; as a result, we have functions ξ h (t) ∈ H with the properties

|ξ h (t) -w(t)| H ≤ h, t ∈ T. (8)
The functions ξ h (t), t ∈ T , are Lebesgue measurable.

In control theory for distributed systems, a linear quadratic control problem (LQP) is widely known. Its solution methods have been studied rather well (see, for example, [START_REF] Lasiecka | Control theory for partial differential equations: continuous and approximation theories[END_REF][START_REF] Curtain | Infinite dimensional linear systems theory[END_REF]). This problem consists in the minimization of some quadratic functional depending on a phase trajectory and control (for example, in the minimization of the deviations in L 2 -norm from a reference control and state trajectory). The problem in question, which is in essence close to LQP, has, at the same time, several distinctive features. Among them, at the first turn, it is worth while noticing the following two features. Firstly, an LQP solution, as a rule, does not guarantee that we find a control generating a trajectory that is close to the reference trajectory in uniform metric. In addition, some apriori information on the reference control is rather often required. The second distinction is in the essence of solving methods. Namely, the method suggested in the present paper is based on constructions of the well-known in the theory of guaranteed control principle of extremal shift.

Control algorithm. Case of continuous measuring of solutions

First, we consider the case of "continuous" measuring of solutions of equations ( 1) and ( 3). In this case, inequalities (8) are valid (for simplicity, we set ξ h (t 0 ) = w 0 ). The problem consists in designing a rule forming (by the feedback principle) the control v = v(t, ξ h (t), w(t)). Fix a function α = α(h) : (0, 1) → (0, 1). Let the control v α,h (t) in equation ( 1) be defined by the formula

v = ṽα,h (t) = v α,h (t) + ṽh (t), (9)
where

ṽh (t) = cB -1 (x(t) -ξ h (t)), v α,h (t) = α -1 B * (x(t) -ξ h (t)). (10)
Here, B * denotes the adjoint operator, c = const > 2ω. Thus, we obtain system (1), (3); i.e., we have the pair of equations

ẋ(t) + Ax(t) = Bu(t) + f (t), ẇα,h (t) + Aw α,h (t) = α -1 BB * (x(t) -ξ h (t)) + c(x(t) -ξ h (t)) + f (t)
with the initial condition x(t 0 ) = x 0 , w α,h (t 0 ) = w 0 . Here, we denote by w α,h (•) the solution of equation ( 1) corresponding to the function v = v α,h (•) of form [START_REF] Maksimov | Dynamical state and control reconstruction for a phase field model. Dynamics of continuous, discrete and impulsive systems[END_REF]. The second formula in ( 10) is an analog of relation [START_REF] Osipov | Methods of Dinamical Reconstruction of Inputs of Controlled Systems[END_REF]. If the constraint in the form of the set P is absent then the application of formula (6) for calculating the control v is impossible, since in this case it is required to solve the problem of minimization of the linear functional l i (u) = (ξ h i -y(τ i ), Bu) over the whole space U . It is natural to replace this problem by a new regularized problem with a smoothing functional of the form α(h)|v| 2 U , i.e., to replace problem ( 6) by the problem of finding the function v α,h (t) by the rule

v α,h (t) = arg min{α|v| 2 U -2(B * (x(t) -ξ h (t)), v) U : v ∈ U }.
Formula (10) provides the solution of the new problem. Thus, to calculate v α,h (t), we realize the regularization of the method of extremal shift by means of the method of smoothing functional, which is known in the theory of ill-posed problems.

Theorem 1. Let α = α(h) → 0. Then the following inequalities

|x(t) -w α,h (t)| 2 H + 2c t t0 |x(τ ) -w α,h (τ )| 2 V dτ ≤ d 0 (h + α(h)), t ∈ T, (11)
ϑ t0 |ṽ α,h (τ )| 2 U dτ ≤ ϑ t0 |u(τ )| 2 U dτ + d * (hα -2 (h) + h 1/2 + α 1/2 (h)) (12)
are fulfilled. Here, d 0 , d * = const > 0 are constants, which do not depend on h ∈ (0, 1).

Proof. Due to [START_REF] Samarskii | Introduction to the Theory of Difference Schemes[END_REF], it holds that |v α,h (t

)| 2 U ≤ 2b 2 α -2 (h 2 + |µ α,h (t)| 2 H ), t ∈ T , where µ α,h (t) = x(t) -w α,h (t), b = |B * | L(H;U )
is the norm of the linear operator B * ∈ L(H; U ). In this case, we have

t t0 |v α,h (τ )| 2 U dτ ≤ 2b 2 α -2 ̺ h (t) + c 1 h 2 α -2 , ̺ h (t) = t t0 |µ α,h (τ )| 2 H dτ. (13)
Due to coercivity condition (2), we obtain,

εh (t) ≤ -2(v α,h (t), B * (x(t)-ξ h (t))) U +α|v α,h (t)| 2 U +2(u(t), B * (x(t)-ξ h (t))) U - -α|u(t)| 2 U + 2bh{|u(t)| U + |v α,h (t)| U } + (2ω -c)|µ α,h (t)| 2 H + 2ch 2 , ( 14
)
where 14) and [START_REF] Samarskii | Introduction to the Theory of Difference Schemes[END_REF], it follows that

ε h (t) = |µ α,h (t)| 2 H + 2c * t t0 |µ α,h (τ )| 2 V dτ + α t t0 {|v α,h (τ )| 2 U -|u(τ )| 2 U } dτ . From (
ε h (t) ≤ ε h (t 0 ) + c 2 + t t0 2bh{|u(τ )| U + |v α,h (τ )| U }dτ + (2ω -c)̺ h (t). (15)
From ( 13), (15), and the inequality 2ω -c < 0, we derive

ε h (t) ≤ c 3 (h + h 3 α -2 ) + c 4 (hα -2 + h)̺ h (t). (16)
Therefore, from (16) we get the bound

|µ α,h (t)| 2 H ≤ c 5 (h + α + h 3 α -2 ) + c 4 (hα -2 + h)̺ h (t). ( 17
)
By the Gronwall inequality and (17), we obtain

|µ α,h (t)| 2 H ≤ c 5 (h + α + h 3 α -2 ) exp{c 4 (t -t 0 )(hα -2 + h)}. ( 18
)
Note that hα -2 ≤ const. Then

|µ α,h (t)| 2 H ≤ c 6 (h + α). (19)
From ( 16) and ( 19) we derive

ε h (t) ≤ c 3 (h + h 3 α -2 ) + c 7 (hα -2 + 1)(h + α) ≤ c 9 (h + α). (20)
Relation [START_REF] Lasiecka | Control theory for partial differential equations: continuous and approximation theories[END_REF] follows from [START_REF] Samarskii | Introduction to the Theory of Difference Schemes[END_REF]. Let us verify [START_REF] Curtain | Infinite dimensional linear systems theory[END_REF]. By virtue of inequality (19), from ( 16) we obtain

ε h (t) ≤ c 8 {h + h 3 α -2 + (hα -2 + h)(h + α)} ≤ c 9 {h + h 2 α -2 + hα -1 }. (21)
Using ( 21), we get for t ∈ T

t t0 |v α,h (τ )| 2 U dτ ≤ t t0 |u(τ )| 2 U dτ + c 10 hα -2 , ( 22
)
|ṽ h (t)| U ≤ c 11 (h + α) 1/2 . ( 23
)
Relation [START_REF] Curtain | Infinite dimensional linear systems theory[END_REF] follows from ( 22) and ( 23). The theorem is proved.

Control algorithm. Case of discrete measuring of solutions

Let us describe the algorithm for solving the problem in the case of discrete measuring of phase states. In this case, we assume that relations (4) are fulfilled. Let l(•) :

W 1,2 (T ; H) ∩ L 2 (T ; V ) → R + , l(y(•)) = |y(•)| C(T ;H) + | ẏ(•)| L2(T ;H) + |y(•)| L2(T ;V )
. In a standard way, we establish the validity of the following lemma.

Let a family of partitions

∆ h = {τ h,i } m h i=0 , τ h,0 = t 0 , τ h,m h = ϑ, τ h,i+1 = τ h,i + δ(h)
and a function α(h) : (0, 1) → (0, 1) be fixed. First, before the moment t 0 , the value h and the partition ∆ h of the interval T are chosen and fixed. The work of the algorithm is decomposed into m -1 (m = m h ) identical steps. At the ith step, which is carried out on the time interval δ i = [τ i , τ i+1 ), τ i = τ h,i , the following sequence of actions is fulfilled. First, at the moment τ i , the element

v h i = α -1 B * (x(τ i ) -ξ h i ) (24)
is calculated. Then, the control defined by the formula

v = ṽh (t) = v h i + cB -1 (x(τ i ) -ξ h i ), t ∈ δ i , (25)
is fed onto the input of equation [START_REF] Bensoussan | Representation and Control of Infinite Dimensional Systems[END_REF], where c = const > 2ω. Under the action of this control, instead of the state

w h (τ i ) = w h (τ i ; τ i-1 , w h (τ i-1 ), v h i-1 ), the state w h (τ i+1 ) = w h (τ i+1 ; τ i , w h (τ i ), v h i ) is realized.
The work of the algorithm stops at the time moment ϑ.

Let the family of partitions ∆ h of the time interval T and the function α(h) have the following properties:

hδ -1 (h) ≤ C 1 , δ(h)α -2 (h) → 0, hα -1 (h) → 0, (26)
α(h) → 0, δ(h) → 0 as h → 0+.
Here C 1 = const > 0 is a constant, which does not depend on h.

Theorem 2. Uniformly with respect to h ∈ (0, 1), the inequalities

λ h (t) ≡ |x(t)-w h (t)| 2 H +2c t t0 |x(τ )-w h (τ )| 2 V dτ ≤ d 1 (h+α+δ) ∀t ∈ T, (27)
ϑ t0 |ṽ h (τ )| 2 U dτ ≤ ϑ t0 |u(τ )| 2 U dτ + d 2 (hα -1 + δα -2 ) + d 3 (h + α + δ) 1/2 (28) are true. Here, d 1 , d 3 (d 1 -d 3 = const > 0)
are constants, which do not depend on h, α = α(h), and δ = δ(h).

Proof. First, we verify inequality (27). Using the invertibility of the operator B as well as coercitivity condition (2), we obtain for a.a. t ∈ δ i the inequality 0.5

d|µ h (t)| 2 H dt + c * |µ h (t)| 2 V -ω|µ h (t)| 2 H ≤ (B(u(t)-v h (t)) -c(x(τ i )-ξ h i ), µ h (t)) U ,
where µ h (t) = x(t) -w h (t) for t ∈ T , v h (t) = v h i for t ∈ δ i . From the inequality

c(ξ h i -x(τ i ), µ h (t)) ≤ -0.5c|µ h (t)| 2 H +4ch 2 +8c(t-τ i ) t τi {| ẋ(τ )| 2 H +| ẇh (τ )| 2 H } dτ,
we have for a.a.

t ∈ δ i (B(u(t) -v h (t)) -c(x(τ i ) -ξ h i ), µ h (t)) U ≤ (B(u(t) -v h (t)), x(τ i ) -ξ h i ) U + ̺ i (t, h) + χ i (t, h) -0.5c|µ h (t)| 2 H . Here, χ i (t, h) = 4ch 2 + 8c(t -τ i ) t τi {| ẋ(τ )| 2 H + | ẇh (τ )| 2 H } dτ, ̺ i (t, h) = b(|u(t)| U + |v h (t)| U )(h + t τi {| ẇh (τ )| H + | ẋ(τ )| H } dτ ).
For a.a. t ∈ δ i , we deduce that

εh (t) ≤ -2(v h (t), B * (x(τ i ) -ξ h i )) U + α|v h (t)| 2 U + + 2(u(t), B * (x(τ i ) -ξ h i )) U -α|u(t)| 2 U + 2̺ i (t, h) + 2χ i (t, h) + (2ω -c)|µ h (t)| 2 H ,
where

ε h (t) = |µ h (t)| 2 H + 2c * t t0 |µ h (τ )| 2 V dτ + α t t0 {|v h (τ )| 2 U -|u(τ )| 2 U } dτ .
Therefore, by virtue of the rule of forming the control ṽh (•) (see ( 24) and ( 25)), we conclude that, for a.a. t ∈ δ i ,

ε h (t) ≤ ε h (τ i ) + c 1 h 2 + c 2 δ t τi {|u(τ )| 2 U + |v h (τ )| 2 U } dτ + (29) + c 3 δ t τi {| ẇh (τ )| 2 H + | ẋ(τ )| 2 H }dτ + (2ω -c) t τi |µ h (τ )| 2 H dτ.
Summing the right-hand and left-hand parts of (29) over i and taking into account Lemma 1, we obtain for t ∈ T

ε h (t) ≤ ε h (t 0 ) + c 4 h 2 δ -1 + c 6 δ + c 7 γ h,δ (t). (30)
Here, γ h,δ (t) = δ 2 i(t) j=0 |v h j | 2 U . Using (4) and the rule of forming v h i (see ( 24)), we get

|v h i | 2 U ≤ 2b 2 (̺ h i + h 2 )α -2 ≤ c 8 (̺ h i + h 2 )α -2 , (31)
where

̺ h i = |x(τ i ) -w h (τ i )| 2 H .
Due to [START_REF] Maksimov | Dynamical reconstruction of inputs for constraction semigroup systems: the boudary input case[END_REF], we have ε h (t 0 ) ≤ h 2 . Therefore, we derive from (30) the estimate

λ h (t) ≤ c 9 (δ + h 2 δ -1 + α + γ h,δ (t)). ( 32
)
Note that ̺ h i ≤ λ h i , where λ h j = λ h (τ j ). Therefore, for t ∈ [τ i , τ i+1 ], due to (31), the inequality

γ h,δ (t) ≤ c 8 δ 2 i(t) j=0 (λ h j + h 2 )α -2 (33)
is valid. Consequently, (32) implies the inequality

λ h i ≤ c 10 (δ + h 2 δ -1 + α) + c 11 δh 2 α -2 + c 12 δ 2 α -2 i j=0 λ h j . (34)
By the discrete Gronwall inequality [START_REF] Samarskii | Introduction to the Theory of Difference Schemes[END_REF], (34), and the inequalities hδ -1 (h) ≤ C 1 , δα -2 (h) ≤ C 2 as h → 0 (see ( 26)), we have

λ h i ≤ c 14 (h + δ + α), i ∈ [0 : m].
This and (33) imply γ h,δ (t) ≤ c 15 (h + δ + α) t ∈ T . Moreover, from the last inequality and (32), we have λ h (t) ≤ c 16 (δ +h 2 δ -1 +α+γ h,δ (t)) ≤ c 17 (h+δ +α). Relation (27) follows from the last inequality. The proof of (28) is similar to the proof of ( 12). The theorem is proved. It follows from Theorems 1 and 2 that the algorithms presented above are resource-saving.

Example

The second algorithm was tested. The parabolic equation w t (t, η) -∂ 2 w(t, η)/∂η 2 = v(t, η), η ∈ [0, 1] (

with the boundary w(t, 0) = w(t, 1), t ∈ T = [0, 2] and initial w(0, η) = 0, η ∈ [0, 1] conditions was considered. The reference equation (see ( 3)) was of the form x t (t, η) -∂ 2 x(t, η)/∂η 2 = u(t, η), η ∈ [0, 1] (36)

x(t, 0) = x(t, 1) = 0, t ∈ T, x(0, η) = 0, η ∈ [0, 1].

Equations ( 35) and (36) were solved by the grid method [START_REF] Samarskii | Introduction to the Theory of Difference Schemes[END_REF]. The grid {η j } n j=0 , η 0 = 0, η n = 1 with the step γ N = 1/n was taken on the interval [0, 1]. The control v = v h (t, η) in the right-hand part of (35) was calculated by formula (25) taking the form v h (t, η j ) = (α -1 + c)(x(τ i , η j ) -ξ h i (η j )), t ∈ [τ i , τ i+1 ), j ∈ [0 : n].

During the experiment, we assumed that ξ h i (η j ) = w(τ i , η j ) + h. In figs. 1-4, the cross-sections of the trajectories x (dashed line) and w (solid line) by the hyperplane η = 0, 4 are presented, as well as the variations of the values p(t) = |x(τ i , η j ) -w h (τ i , η j )| = 0, 01926, in the first case 0, 00972, in the second case.

Lemma 1 .

 1 There exists a number K = K(ω, D, c * , |B| L(U ;H) ) such that the inequality l(x(•; t 0 , x, u(•))) ≤ K(1 + |u(•)| L2(T ;U ) ) is fulfilled uniformly for any x ∈ D and u(•) ∈ L 2 (T ; U ).

Fig. 3 .Fig. 4 .

 34 Fig. 1. Fig. 2.

Research (No. 11-01-12112-ofi-m), by the Program on Basic Research of the Presidium of the Russian Acad. Sci. (No. 12-P-1-1012) and by the Program for support of leading scientific schools of Russia (No. 6512.2012.1).


Daniel Wachsmuth

Gerd Wachsmuth

Necessary conditions for convergence rates of regularizations of optimal control problems

Keywords: optimal control problem, inequality constraints, Tikhonov regularization, source condition

We investigate the Tikhonov regularization of control constrained optimal control problems. We use a specialized source condition in combination with a condition on the active sets. In the case of high convergence rates, these conditions are necessary and sufficient.

Introduction

In this article, we investigate regularization schemes for the following class of optimization problems:

Minimize 1 2
Su -z 2 Y + β u L 1 (Ω) such that u ∈ L 2 (Ω) and u a ≤ u ≤ u b a.e. on Ω.

(P)

Here, Ω is a measurable subset of R n , n ≥ 1, Y is a Hilbert space, S : L 2 (Ω) → Y a bounded linear operator, and the function z ∈ Y is given. The parameter β is assumed to be non-negative. The control constraints u a , u b ∈ L ∞ (Ω) satisfy u a ≤ 0 ≤ u b . This model problem can be interpreted as an optimal control problem as well as an inverse problem. In the point of view of inverse problems, the unknown u has to be constructed in order to reproduce given measurements z. The inequality constraints on u reflect certain a-priori knowledge about the solution u † of the linear ill-posed equation Su = z. If the problem at hand is seen as an optimal control problem, then u is the control, Su the state of the system, which has to be close to a desired state z, the inequality constraints restrict the feasible set and may hinder the state Su to reach the target z. If the parameter β is positive, then the resulting optimal control will be sparse, that is, its support is a possibly small subset of Ω. The resulting optimization problem (P) is nevertheless ill-posed if S is not continuously invertible. Due to the control constraints, problem (P) still possesses a solution, which is even unique if S is injective. However, the solution may be unstable with respect to perturbations in the problem data, for instance in the given state z. Here small perturbations due to measurement errors may lead to large changes in the solution. Consequently, any numerical approximation of (P) is challenging to solve and numerical approximations of solutions may converge arbitrarily slow. Let us note, that a positive value of β does not make the problem well-posed. This is due to the fact, that L 1 (Ω) is not a dual space and hence bounded sets in L 1 (Ω) are not compact w.r.t. the weak(-star) topology, see also the discussions in [7,[START_REF] Wachsmuth | On the regularization of optimization problems with inequality constraints[END_REF]. In order to overcome this difficulty, we apply common ideas from inverse problem theory. We will study a regularization of the type

Minimize 1 2 Su -z 2 Y + β u L 1 (Ω) + α 2 u 2 L 2 (Ω)
such that u ∈ L 2 (Ω) and u a ≤ u ≤ u b a.e. on Ω,

(P α )
where α > 0 is given. Clearly, the problems (P α ) are uniquely solvable for α > 0. Now, the question arises, whether their solutions u α converge (weakly or strongly) to a solution u 0 of (P) for α → 0. Moreover, in the case of convergence, one is interested in proving convergence rates of u α -u 0 L 2 (Ω) and Su α -Su 0 Y under suitable assumptions.

In this work, we will prove necessary conditions for convergence rates. In some parts, the necessary conditions are similar to sufficient conditions found in earlier works [7,[START_REF] Wachsmuth | On the regularization of optimization problems with inequality constraints[END_REF]. Moreover, the result of Theorem 3 leads to a weakened sufficient condition for convergence rates.

Standing assumptions and notation

Let us fix the standing assumptions on the problem (P). We assume that S : L 2 (Ω) → Y is linear and continuous. In many applications this operator S is compact. Furthermore, we assume that the Hilbert space adjoint operator S ⋆ maps into L ∞ (Ω), i.e., S ⋆ ∈ L(Y, L ∞ (Ω)).

U ad := {u ∈ L 2 (Ω) : u a ≤ u ≤ u b a.e. on Ω}.
The problem (P α ) is uniquely solvable for α > 0. We will denote its solution by u α , with the corresponding state y α := Su α and adjoint state p α := S ⋆ (z -y α ).

There is a unique solution of (P) with minimal L 2 (Ω) norm, see [START_REF] Wachsmuth | On the regularization of optimization problems with inequality constraints[END_REF]Thm. 2.3,Lem. 2.7]. This solution and the associated state and adjoint state will be denoted by u 0 , y 0 and p 0 , respectively. Note that the weak convergence

u α ⇀ u ⋆ in L 2 (Ω),
where u ⋆ is a solution of (P) already implies u ⋆ = u 0 , see [START_REF] Wachsmuth | On the regularization of optimization problems with inequality constraints[END_REF]Rem. 3.3].

Optimality conditions

As both problems (P) and (P α ) are convex, their solutions can be characterized by the following necessary and sufficient optimality conditions:

Theorem 1 ([7, Lemma 2.2]). Let α ≥ 0 be given, and let u α be a solution of (P α ) (or (P) in the case α = 0).

Then, there exists a subgradient

λ α ∈ ∂ u α L 1 (Ω) , such that the variational inequality (α u α -p α + β λ α , u -u α ) ≥ 0 ∀u ∈ U ad , (1)
is satisfied, where p α = S ⋆ (z -Su α ) is the associated adjoint state.

Here, (•, •) refers to the scalar product in L 2 (Ω).

Standard arguments (see [6, Section 2.8]) lead to a pointwise a.e. interpretation of the variational inequality, which in turn implies the following relation between u α and p α in the case α > 0:

u α (x) =                u a (x) if p α (x) < α u a (x) -β 1 α (p α (x) + β) if α u a (x) -β ≤ p α (x) ≤ -β 0 if |p α (x)| < β 1 α (p α (x) -β) if β ≤ p α (x) ≤ α u b (x) + β u b (x) if α u b (x) + β < p α (x)
a.e. on Ω.

(2)

In the case α = 0, we have

u 0 (x)                = u a (x) if p 0 (x) < -β ∈ [u a (x), 0] if p 0 (x) = -β = 0 if |p 0 (x)| < β ∈ [0, u b (x)] if p 0 (x) = β = u b (x) if β < p 0 (x)
a.e. on Ω.

Note that if β = 0, one obtains u 0 (x) ∈ [u a (x), u b (x)] where p 0 (x) = 0 in (3). This implies that u 0 (x) is uniquely determined by p 0 (x) on the set, where it holds |p 0 (x)| = β.

Sufficient conditions for convergence rates

Let us first recall the sufficient conditions for convergence rates as obtained in [START_REF] Wachsmuth | On the regularization of optimization problems with inequality constraints[END_REF]. We will work with the following assumption. There we denote by proj

(x) =      proj [ua(x),0] (S ⋆ w)(x) if β > 0, p 0 (x) ≤ -β 2 , proj [0,u b (x)] (S ⋆ w)(x) if β > 0, p 0 (x) ≥ β 2 , proj [ua(x),u b (x)] (S ⋆ w)(x) if β = 0. (4)
2. (structure of active set) A = Ω \ I and for all ǫ > 0

meas {x ∈ A : 0 < |p 0 (x)| -β < ǫ} ≤ c ǫ κ if w = 0, meas {x ∈ A : 0 < |p 0 (x)| -β < ǫ} ≤ c ǫ κ if w = 0. (5)
Some remarks are in order. The first part of the assumption is analogous to source conditions in inverse problems: we assume that on the set I ⊂ Ω the solution u 0 is the restriction to I of a certain pointwise projection of an element in the range of S ⋆ . This part of the condition is different from other conditions in the literature: in our earlier work [START_REF] Wachsmuth | On the regularization of optimization problems with inequality constraints[END_REF] we used the assumption u 0 (x) = proj [ua(x),u b (x)] (S ⋆ w)(x) on I. However, in the light of the derivation of necessary conditions it turns out that such a condition can be weakened without losing anything with respect to convergence rates. In works on inverse problems [START_REF] Engl | Regularization of inverse problems[END_REF][START_REF] Neubauer | Tikhonov-regularization of ill-posed linear operator equations on closed convex sets[END_REF], the source condition u 0 = proj U ad (S ⋆ w) is used, which is retained as the special case I = Ω in Assumption 2.

The assumption (5) (without the second alternative) on the active sets was already employed to obtain regularization error estimates [7,8], error estimates for finite-element discretizations of (P) [START_REF] Deckelnick | A note on the approximation of elliptic control problems with bang-bang controls[END_REF], as well as stability results of bang-bang controls [START_REF] Felgenhauer | On stability of bang-bang type controls[END_REF]. Note that in the case β = 0, both conditions in (5) are equivalent. However, if β > 0 and w = 0 (in particular, if I has measure zero), the second alternative provides a weaker condition than the first one. Hence, condition [START_REF] Neubauer | Tikhonov-regularization of ill-posed linear operator equations on closed convex sets[END_REF] is weaker than the condition used in our earlier work [START_REF] Wachsmuth | On the regularization of optimization problems with inequality constraints[END_REF].

Theorem 3. Let Assumption 2 be satisfied. Let d be defined as

d =      1 2-κ if κ ≤ 1, 1 if κ > 1 and w = 0, κ+1 2 if κ > 1 and w = 0.
Then there is α max > 0 and a constant c > 0, such that

y 0 -y α Y ≤ c α d p 0 -p α L ∞ (Ω) ≤ c α d u 0 -u α L 2 (Ω) ≤ c α d-1/2
holds for all α ∈ (0, α max ].

Under the assumptions of the theorem, one can prove also convergence rates for u α -u 0 L 1 (A) [START_REF] Wachsmuth | On the regularization of optimization problems with inequality constraints[END_REF].

Proof. The proof is analogous to the proof of [START_REF] Wachsmuth | On the regularization of optimization problems with inequality constraints[END_REF]Thm. 3.14]. We have to take into account the modification of the source condition (4) in the case β > 0 and the modification of (5) in the case w = 0. By [8, Lemma 2.12], we have

y 0 -y α 2 Y + α u 0 -u α 2 L 2 (Ω) ≤ α (u 0 , u 0 -u α ). (6)
Since U ad is bounded, we obtain p 0 -p α L ∞ (Ω) ≤ c α 1/2 for some c > 0 independent of α.

Let now α be small enough such that p 0 -p α L ∞ (Ω) < β/2. This implies that p 0 and p α have the same sign on the set {x ∈ I : |p 0 (x)| ≥ β/2}. Consequently, u 0 and u α have the same sign on this set, too. Moreover, on the set {x ∈ I : |p 0 (x)| < β/2} it holds |p α | < β, and hence u α = 0 = u 0 on this set. This yields

(χ I u 0 , u 0 -u α ) ≤ (χ I S ⋆ w, u 0 -u α )
for α > 0 small enough. Note that in case of w = 0, the right-hand side in the previous estimate vanishes and it remains to estimate (χ A u 0 , u 0 -u α ). Taking into account that u 0 (x) = 0 whenever |p 0 (x)| < β, the weakend estimate ( 5) is sufficient in this case. Arguing as in the proof of [8, Thm. 3.14] proves the claim. ⊓ ⊔

3 Necessary conditions for convergence rates

u 0 (x) =      proj [ua(x),0] (S ⋆ w)(x) if β > 0, p 0 (x) = -β, proj [0,u b (x)] (S ⋆ w)(x) if β > 0, p 0 (x) = +β, proj [ua(x),u b (x)] (S ⋆ w)(x) if β = 0, p 0 (x) = 0. If moreover y α -y 0 Y = o(α), then u 0 = 0 on {x ∈ Ω : |p 0 (x)| = β}, i.e. w = 0.
This result shows that the source condition ( 4) is necessary on the set {x ∈ Ω :

|p 0 (x)| = β}.
Proof. Let us prove the claim in the case β > 0. The result in the case β = 0 can be proved with obvious modifications. Let us take a test function u ∈ U ad defined as

u(x)                = u a (x) if p 0 (x) < -β, ∈ [u a (x), 0] if p 0 (x) = -β, = 0 if |p 0 (x)| < β, ∈ [0, u b (x)] if p 0 (x) = β, = u b (x) if p 0 (x) > β.
Due to the relation

λ 0 = proj [-1,1] 1 β p 0 , (7)
which is a consequence of the necessary optimality condition, see [START_REF] Casas | Optimality conditions and error analysis of semilinear elliptic control problems with L 1 cost functional[END_REF] for a proof, we obtain λ 0 = ±1 where p 0 = ±β. Hence it holds

(-p 0 , u -u 0 ) = (-β λ 0 , u -u 0 ) = β u 0 L 1 (Ω) -β u L 1 (Ω) (8)
for u as above. Since λ 0 ∈ ∂ u 0 L 1 (Ω) , we obtain

(λ 0 , u α -u 0 ) ≤ u α L 1 (Ω) -u 0 L 1 (Ω) . (9)
Using the optimality of u α and the relation -p α = -p 0 + S ⋆ S(u α -u 0 ) we get

(-p 0 + S ⋆ S(u α -u 0 ) + α u α , u -u α ) + β u L 1 (Ω) -β u α L 1 (Ω) ≥ 0.
Adding (-p 0 + β λ 0 , u α -u 0 ) ≥ 0 to the left-hand side yields

(S ⋆ S(u α -u 0 ) + α u α , u -u α ) + (-p 0 , u -u 0 ) + (β λ 0 , u α -u 0 ) +β u L 1 (Ω) -β u α L 1 (Ω) ≥ 0.
Using ( 8) and ( 9) we obtain

(S ⋆ S(u α -u 0 ) + α u α , u -u α ) ≥ 0.
Due to the assumptions of the theorem, the functions 1 α (S(u α -u 0 )) = 1 α (y α -y 0 ) are uniformly bounded for α ց 0. As a consequence, α ց 0 implies (S ⋆ ẏ0 + u 0 , u -u 0 ) ≥ 0 for any weak subsequential limit ẏ0 of 1 α (y α -y 0 ). Due to the construction of the test function u, we obtain

u 0 = proj [ua,0] (S ⋆ ẏ0 ) where p 0 = -β, proj [0,u b ] (S ⋆ ẏ0 ) where p 0 = +β. If y α -y 0 Y = o(α) then 1 α (y α -y 0 ) → 0 strongly in Y for α → 0, hence ẏ0 = 0, and u 0 = 0 on the set {|p 0 | = β}. ⊓ ⊔
As can be seen from the proof, the element that realizes the source condition can be interpreted as the (weak) directional derivate of α → y α at α = 0. The result of the theorem resembles known results of necessity of the source condition in linear inverse problems, see e.g. [START_REF] Engl | Regularization of inverse problems[END_REF][START_REF] Neubauer | Tikhonov-regularization of ill-posed linear operator equations on closed convex sets[END_REF].

Necessity of the condition (5) on the active set

In this section, we want to prove the necessity of (5) in the case of high convergence rates d > 1. In this case, we have w = 0, see Theorem 4. It remains to show that the second condition in ( 5) is necessary to obtain convergence rates d > 1. Hence, we derive a bound on

µ(ǫ) := |{x ∈ Ω : 0 < |p 0 (x)| -β < ǫ}| ,
which is the measure of a subset of

A := {x ∈ Ω : β < |p 0 (x)|}.
For α > 0 let ũα denote the unique solution of

min u∈U ad -(u, p 0 ) + α 2 u 2 L 2 (Ω) + β u L 1 (Ω) . (P aux α )
Analogous to (2), we have the representation

ũα (x) =                u a (x) if p 0 (x) < α u a (x) -β 1 α (p 0 (x) + β) if α u a (x) -β ≤ p 0 (x) ≤ -β 0 if |p 0 (x)| < β 1 α (p 0 (x) -β) if β ≤ p 0 (x) ≤ α u b (x) + β u b (x) if α u b (x) + β < p 0 (x)
a.e. on Ω. (10)

Let us first prove a relation between the convergence rates of u 0 -ũα L 2 (A) and µ(ǫ) for α → 0 and ǫ → 0, respectively. Proof. Due to the pointwise representations of ũα and u 0 in (10) and (3), respectively, it holds

u 0 -ũα 2 L 2 (A) = {β<p0<αu b +β} (u b -α -1 (p 0 -β)) 2 + {αua-β<p0<-β} (u a -α -1 (p 0 + β)) 2 .
Due to the assumption on the control constraints we have

{β<p0<αu b +β} (u b -α -1 (p 0 -β)) 2 ≥ {β<p0<ασ/2+β} (u b -α -1 (p 0 -β)) 2 ≥ {β<p0<ασ/2+β} (σ/2) 2 ≥ (σ/2) 2 {x ∈ Ω : 0 < p 0 (x) -β < ασ/2} .
Similarly, we obtain {αua-β<p0<-β}

(u a -α -1 (p 0 +β)) 2 ≥ (σ/2) 2 {x ∈ Ω : 0 < -p 0 (x)-β < ασ/2} . This implies u 0 -ũα 2 L 2 (A) ≥ (σ/2) 2 µ(ασ/2). Hence if u 0 -ũα L 2 (A) = O(α d ) holds, then µ(ασ/2) ≤ O(α 2d ),
for α → 0, which proves the claim.

⊓ ⊔

Using the same arguments, we can prove the following result.

Corollary 6. Let the requirements of Lemma 5 be satisfied. Let p ∈ [1, ∞) be given. Then it holds

σ 2 p µ σ 2 α ≤ u 0 -ũα p L p (A) ≤ M p µ(M α) with M = max( u a L ∞ (Ω) , u b L ∞ (Ω) ).
Lemma 7. Let ũα be defined as above. Then it holds

α ũα -u α 2 L 2 (Ω) + y 0 -y α 2 Y ≤ (p 0 -p α , ũα -u 0 ).
Proof. Since u α and ũα solve (P α ) and (P aux α ), respectively, we have

(αu α -p α + βλ α , ũα -u α ) ≥ 0, (αũ α -p 0 + β λα , u α -ũα ) ≥ 0,
with some λα ∈ ∂ ũα L 1 (Ω) . Due to the monotonicity of the subdifferential we have (λ α -λα , u α -ũα ) ≥ 0. This gives

α ũα -u α 2 L 2 (Ω) ≤ (p 0 -p α , ũα -u α ).
The identity

(p 0 -p α , ũα -u α ) = (p 0 -p α , ũα -u 0 + u 0 -u α ) = (p 0 -p α , ũα -u 0 ) -y 0 -y α 2 Y
finishes the proof. ⊓ ⊔ Theorem 8. Let us assume that there is σ > 0 such that u a ≤ -σ < 0 < σ ≤ u b a.e. on Ω. The we have the following implication: If

u 0 -u α L 2 (Ω) = O(α d-1/2 ), y 0 -y α Y = O(α d ) for α → 0 holds with d > 1, then it follows µ(ǫ) ≤ O(ǫ 2d-1 ) for ǫ → 0.
Proof. Let us begin with In our opinion, this condition is too strong and has to be relaxed in order to obtain a characterization for convergence rates d ≤ 1.

u 0 -ũα 2 L 2 (Ω) ≤ 2( u 0 -u α 2 L 2 (Ω) + u α -ũα 2 L 2 (Ω) ) ≤ O(α 2d-1 ) + α -1 (p 0 -p α , ũα -u 0 ) ≤ O(α 2d-1 ) + O(α d-1 ) u 0 -ũα L 2 (Ω) , which gives u 0 -ũα L 2 (Ω) = O(α d-

Necessary conditions for exact reconstruction with α > 0

Let us now investigate the case of exact reconstruction. That is, the solutions of the regularized problem u α coincide with the (minimal L 2 -norm) solution u 0 of the original problem.

Lemma 9. Let us assume that u α * = u 0 a.e. on Ω for some α * > 0. Then u α = u 0 a.e. on Ω for all α ∈ (0, α * ).

Proof. The claim follows from known monotonicity results: The mapping α → u α L 2 is monotonically decreasing, while α → 1 2 y α -y d 2 Y + β u α L 1 is monotonically increasing from (0, +∞) to R, see e.g. [START_REF] Wachsmuth | On the regularization of optimization problems with inequality constraints[END_REF]Lemma 2.8].

⊓ ⊔

Assumption 2

 2 [a,b] (v) the projection of the real number v onto the interval [a, b]. Let u 0 be a solution of (P). Let us assume that there exist a measurable set I ⊂ Ω, a function w ∈ Y , and positive constants κ, c such that it holds: 1. (source condition) I ⊃ {x ∈ Ω : |p 0 (x)| = β}, and for almost all x ∈ I u 0

3. 1 4 ) 4 .

 144 Necessity of the source condition (Theorem Let us suppose that y α -y 0 Y = O(α) with y 0 = Su 0 . Then there exists w ∈ Y such that

Lemma 5 .

 5 Let us assume that there is σ > 0 such that u a ≤ -σ < 0 < σ ≤ u b a.e. on Ω. Then it holds: If u 0 -ũα L 2 (A) = O(α d ), d > 0, for α → 0, then µ(ǫ) = O(ǫ 2d ) for ǫ → 0.

1 ). 2 L 2 2 L 2

 12222 Hence by Lemma 5, we obtain µ(ǫ) = O(ǫ 2d-2 ). Let us note that the convergence rates imply u 0 (x) = 0 if |p 0 (x)| = β by Theorem 4. Moreover, we have u 0 = ũα = 0 on {x ∈ Ω : |p 0 (x)| ≤ β} by (3) and (10). This implies u 0 = ũα = 0 on the set {x ∈ Ω : |p 0 (x)| ≤ β} = Ω \ A, cf. (10). Using the convergence rate p 0 -p α L ∞ (Ω) = O(α d ) and Corollary 6, we findα -1 |(p 0 -p α , ũα -u 0 )| = O(α d-1 ) ũα -u 0 L 1 (Ω) = O(α d-1 ) ũα -u 0 L 1 (A) ≤ O(α d-1 )µ(M α).Since by the above considerations we already got µ(ǫ) = O(ǫ 2d-2 ) this givesu 0 -ũα (Ω) = O(α 2d-1 ) + O(α 3(d-1) ).Repeating this processk times until k(d -1) ≥ 2d -1 yields u 0 -ũα (Ω) = O(α 2d-1 ),which finishes the proof.⊓ ⊔Together with Theorem 4, this result shows that the requirements of Theorem 3 for convergence rates d > 1 are sharp.It is an open question, whether the requirement (5) on the active set is also necessary for convergence rates d ≤ 1.

Theorem 10. Let us assume that there is σ > 0 such that u a ≤ -σ < 0 < σ ≤ u b a.e. on Ω. Then the exact recovery u α * = u 0 a.e. on Ω for some α * > 0 is equivalent to

for some ǫ > 0.

Proof. Let us assume u α * = u 0 for some α * > 0. Lemma 9 and Theorem 4 imply

, where ũα * is defined by (10). Hence, Corollary 6 implies µ(σ α * /2) = 0.

To prove the converse, let (11) be satisfied for some ǫ > 0. Using ( 6) we obtain

where

Arguing similarly as in [8, Corollary 3.13], we have |A α | = 0, and hence u 0 -u α L 2 (A) = 0 holds for α > 0 small enough.

⊓ ⊔

In many applications, the adjoint state p 0 belongs to C(Ω). In this case, the result of Theorem 10 shows that an exact reconstruction is only possible if |p 0 (x)| = β for all x ∈ Ω. This in turn implies either u 0 ≡ u a or u 0 ≡ 0 or u 0 ≡ u b on every connected component of Ω.


Eugenio Mijangos
email: eugenio.mijangos@ehu.es

An Algorithm for Two-Stage Stochastic Quadratic Problems

Keywords: Stochastic Programming, Mixed-Integer Quadratic Problems, Branch-and-Fix Coordination, Perspective Cuts

An algorithm for solving quadratic, two-stage stochastic problems is developed. The algorithm is based on the framework of the Branch and Fix Coordination (BFC) method. These problems have continuous and binary variables in the first stage and only continuous variables in the second one. The objective function is quadratic and the constraints are linear. The nonanticipativity constraints are fulfilled by means of the twin node family strategy. On the basis of the BFC method for two-stage stochastic linear problems with binary variables in the first stage, an algorithm to solve these stochastic quadratic problems is designed. In order to gain computational efficiency, we use scenario clusters and propose to use either outer linear approximations or (if possible) perspective cuts. This algorithm is implemented in C++ with the help of the Cplex library to solve the quadratic subproblems. Numerical results are reported.

Introduction

Two-stage stochastic mixed integer programs are among the most interesting problems since the complexity generated by the integrality of variables and the high dimensionality. Stochastic parameters can exist anywhere in the problem. In order to model the uncertainty a finite set of scenarios, Ω, is used, where each ω ∈ Ω has an associated probability of ocurrence p ω . In a two-stage program decisions on the first and second stage variables must be taken. First-stage variables are chosen before knowing the realization of the uncertain parameters. After having decided on first stage and having known each realization of the uncertain parameters, the second stage decision must be taken. The first-stage variables take the same value in each scenario, which yields nonanticipativity constraints. If we consider a finite number of scenarios, a general two-stage program can be expressed regarding the first-stage variables being equivalent to a large programming problem suggested in citeWets66 and known as deterministic equivalent model (DEM). A general two-stage problem can include binary firststage variables. The simplest form of two-stage stochastic integer problems have first-stage binary and second-stage continuous variables. In [START_REF] Laporte | An integer L-shaped algorithm for the capacited vehicle routing problem with stochastic demands[END_REF] a branch-and-cut method is used for those problems, which is based on the Benders decomposition method. An efficient branch-and-fix coordination (BFC) method for solving twostage programs is provided in [START_REF] Alonso-Ayuso | BFC, a branch-and-fix coordination algorithm framework for solving some types of stochasticpure and mixed 0-1 programs[END_REF], where the first-stage has only binary variables, and where the uncertainty only appears in the coefficients of the objective function and in the right-hand-side of the constraints. If the first stage involves pure binary variables, finite termination is justified by using branching over the 0-1 first-stage variables, see among others [START_REF] Sen | Decomposition with branch-and-cut approaches for twostage stochastic mixed-integer programming[END_REF] and [START_REF] Sherali | A modification of Bender's decomposition algorithm for discrete subproblems: an approach for stochastic programs with integer recourse[END_REF]. Escudero et al. [START_REF] Escudero | A general algorithm for solving two-stage stochastic mixed 0-1 first-stage problems[END_REF] study general two-stage stochastic mixed 0-1 problems, where the first stage only involves binary variables and continuous variables and the second stage continuous variables. They use a specialization of the BFC scheme and the twin-node-family (TNF) concept, which was introduced in [START_REF] Alonso-Ayuso | BFC, a branch-and-fix coordination algorithm framework for solving some types of stochasticpure and mixed 0-1 programs[END_REF]. Their scheme is specifically designed for coordinating the node branching selection and pruning and the 0-1 variable branching selection and fixing at each branch-and-fix (BF) tree. Also, they suggest to decompose the set of scenarios in clusters.

On the other hand, real problems with this structure exist and have a high dimensionality. They often need to be solved and it is important to find the procedure that will solve them with the highest efficiency. An example of this type is the Iberian Electricity Market (MIBEL), which comprises Spanish and Portuguese electricity systems, see [START_REF] Heredia | Solving electricity market quadratic problems by Branch and Fix Coordination methods[END_REF].

In this paper we consider the two-stage mixed 0-1 quadratic problem minimize c t δ + q(x, y)

subject to : l a ≤ A δ x ≤ u a , l t ≤ T   δ x y   ≤ u t , x ≥ 0, y ≤ y ≤ y, δ ∈ {0, 1} n δ ,
where δ are first-stage-binary variables, x ∈ R nx are first-stage continuous variables, y ∈ R ny are second-stage continuous variables, c is the coefficient vector for δ, and q is the quadratic function defined as follows

q(x, y) = b t x y + [x t y t ]Q x y ,
where Q is a positive-definite matrix and b and Q are partitioned as

x y ; i.e., b = b x b y and Q = Q xx Q xy Q yx Q yy .
In addition, l a and u a are the bounds for the first-stage constraints and l t and u t are the bounds for the second-stage constraints.

Let us suppose that some of the coefficient in b y , Q xy , Q yx , Q yy , l t , u t and T are uncertain. The uncertainty is given by the scenarios ω in the finite set Ω and p ω is the probability of that occurs ω ∈ Ω. Therefore, the initial problem given in a stochastic way can be written as the so-called Deterministic Equivalent Model (DEM)

minimize c t δ + ω∈Ω p ω q ω (x, y ω ) (1)
subject to :

l a ≤ A δ x ≤ u a , (2)
l ω t ≤ T ω   δ x y ω   ≤ u ω t , ω ∈ Ω, (3)
x ≥ 0, y ≤ y ω ≤ y, ω ∈ Ω, (4)
δ ∈ {0, 1} n δ , (5)
As is shown by [START_REF] Escudero | A general algorithm for solving two-stage stochastic mixed 0-1 first-stage problems[END_REF] the compact representation DEM can be written as a splitting variable representation; i.e., δ and x are respectively replaced by δ ω and x ω , for ω ∈ Ω. So, we have

(MIQP) minimize ω∈Ω p ω (c t δ ω + q ω (x ω , y ω )) subject to : l a ≤ A δ ω x ω ≤ u a , ω ∈ Ω, l ω t ≤ T ω   δ ω x ω y ω   ≤ u ω t , ω ∈ Ω, x ω ≥ 0, y ≤ y ω ≤ y, δ ω ∈ {0, 1} n δ , ω ∈ Ω, (NAC δ ) δ ω -δ ω ′ = 0, ∀ω, ω ′ ∈ Ω : ω = ω ′ , (NAC x ) x ω -x ω ′ = 0, ∀ω, ω ′ ∈ Ω : ω = ω ′ ,
where NAC δ and NAC x are the nonanticipativity constraints.

Note that the relaxation of the NACs in the model MIQP gives rise to |Ω| independent MIQP ω submodels minimize p ω (c t δ ω + q ω (x ω , y ω )), [START_REF] Heredia | Solving electricity market quadratic problems by Branch and Fix Coordination methods[END_REF] subject to :

l a ≤ A δ ω x ω ≤ u a , (7)
l ω t ≤ T ω   δ ω x ω y ω   ≤ u ω t , (8)
x ω ≥ 0, y ≤ y ω ≤ y, (9)

δ ω ∈ {0, 1} n δ , (10)
and these models are linked by the NACs, which force the equality of the firststage variables.

In this work to solve the original quadratic problem DEM a Branch-and-Fix-Coordination scheme (BFC) is used for each scenario ω ∈ Ω to fulfill the integrality condition (IC) given by [START_REF] Wets | Programming under uncertainty: the equivalent convex program[END_REF], so that the NAC δ are also satisfied when selecting branching nodes and branching variables by the Twin-Node-Families concept (TNF), which was introduced by [START_REF] Alonso-Ayuso | BFC, a branch-and-fix coordination algorithm framework for solving some types of stochasticpure and mixed 0-1 programs[END_REF]. A similar approach to that suggested in [START_REF] Escudero | A general algorithm for solving two-stage stochastic mixed 0-1 first-stage problems[END_REF] is used in this work to coordinate the selection of the branching node and branching variable for each scenario-related BF tree, such that the NAC δ are satisfied when fixing δ ω , ∀ω ∈ Ω, either to 1 or to 0. A TNF integer set is a set of integer BF nodes (i.e. they verify IC), one per BF tree, in which the NAC δ are verified.

For each TNF integer set we use two quadratic submodels. The quadratic model QP T N F obtained after fixing in DEM δ = δ ∈ {0, 1} n δ for a TNF integer set

(QP T N F ) Z T N F = c t δ+ min ω∈Ω p ω q ω (x, y ω ) subject to : l a ≤ A δ x ≤ u a , l ω t ≤ T ω   δ x y ω   ≤ u ω t , ω ∈ Ω, x ≥ 0, y ≤ y ω ≤ y, ω ∈ Ω,
It gives a feasible solution and a possible incumbent solution.

The second quadratic submodel to solve at a TNF integer set corresponds to the case where not all the δ variables have been branched on in the current TNF, but all of them hold the integrality condition. Then, the other quadratic submodel (QP f ) is obtained from problem DEM with δ = δ δ f , where δ j , for j ∈ {1, . . . , k}, are fixed to 0-1 values and the componentes of δ f j are in the interval [0, 1].

(QP f ) Z f = min c t δ + ω∈Ω p ω q ω (x, y ω ) subject to : l a ≤ A δ x ≤ u a , l ω t ≤ T ω   δ x y ω   ≤ u ω t , ω ∈ Ω, x ≥ 0, y ≤ y ω ≤ y, ω ∈ Ω, δ j = δ j fixed to 0-1, for j ∈ {1, . . . , k} δ j = δ f j ∈ [0, 1], for j ∈ {k + 1, . . . , n δ }
This model contributes strong lower bounds of the solution value of the descendent nodes from a given node, by satisfying the NAC x .

Outline of BFC

This method branches on the δ-variables, obtaining the solution of the quadratic submodels MIQP ω and coordinating the selection of the branching node and branching variable for the BF trees, such that the NAC δ constraints are fulfilled once fixed the suitable variables δ to 1 or to 0. . A sequence of lower bounds Z i is computed, where Z i = ω∈Ω z ω i and z ω i is the solution to the quadratic relaxation (QP ω ) of MIQP ω once the previous variables δ have been fixed to 0 or to 1.

If the optimal solution obtained in each node of the TNF satisfies the IC (integrality constraints) and the NAC δ , two cases can happen with respect to the NAC x . If NAC x are satisfied, the incumbent solution is updated and the TNF's branch is pruned; if the set of active nodes is empty, that solution is the optimum. Otherwise, to satisfy NAC x we solve the TNF quadratic problem obtained by fixing the δ-variables that verified IC and NAC δ ; if this problem is feasible, the incumbent solution is updated, and if the TNF cannot be pruned, we continue with the tree examination. For more details about the BFC method for two-stage stochastic problems see [START_REF] Escudero | A general algorithm for solving two-stage stochastic mixed 0-1 first-stage problems[END_REF].

Scenario clusters

When the number of scenarios is very large, in order to gain computational efficiency we can take scenario clusters; see in [START_REF] Escudero | An algorithmic framework for solving large-scale multistage stochastic mixed 0-1 problems with nonsymmetric scenario trees[END_REF] an information structuring for scenario cluster partitioning of nonsymmetric scenario trees.

Let p be the number of clusters and Ω 1 , . . . , Ω p , where Ω p ∩ Ω p ′ = ∅ for p, p ′ = 1, . . . , p, such that p = p ′ , and ∪ p p=1 Ω p = Ω. So, instead of the submodel MIQP ω for ω ∈ Ω we can consider the following submodel for the scenario cluster p = 1, . . . , p

(MIQP p ) minimize ω∈Ω p p ω (c t δ p + q ω (x p , y ω )) subject to : l a ≤ A δ p x p ≤ u a , l ω t ≤ T ω   δ p x p y ω   ≤ u ω t , ω ∈ Ω p x p ≥ 0, y ≤ y ω ≤ y, ω ∈ Ω p , δ p ∈ {0, 1} n δ
These models are linked by the NACs δ pδ p ′ = 0 and x px p ′ = 0, for all p, p ′ ∈ {1, . . . , p} such that p = p ′ . However, since the number of branches to test can be huge, the BFC method has some troubles: the number of feasible solutions can be too high, a high number of quadratic problems QP p , QP T N F , and QP f can exist to solve, and QP T N F and QP f can have very high dimensions. Hence, in order to gain computational efficiency, we propose to use either outer linear approximations or (if possible) perspective cuts to solve QP p in each TNF (i.e., MIQP p where the previous branching variables have been fixed and the rest is relaxed in [0, 1])

Outer linear approximations (OLA)

Let the problem min y∈Y g(y), where g is convex and Y is a polyhedral set. The optimal value of that problem is not smaller than that of minimize η subject to : η ≥ g(y i ) + ∇g(y i ) t (yy i ),

y ∈ Y.
Therefore, the value of η * gives us an underestimate of g(y * ) and, so, we can use it instead of g(y * ) in the comparison with the current upper bound, in order to prune (or not) the current branch.

In our problem, we use as y i the solution in the previous node.

Perspective cuts

When Q xy and Q yx are zero matrices, the quadratic function q is defined as follows q(x, y) = b t x x + b t y y + x t Q xx x + y t Q yy y. This kind of model can be found in liberalized electricity markets [START_REF] Heredia | Solving electricity market quadratic problems by Branch and Fix Coordination methods[END_REF] and [START_REF] Corchero | A new optimal electricity market bid model solved through perspective cuts[END_REF].

For each scenario ω ∈ Ω we can write the objective function of the submodel MIQP ω as follows

p ω (c t δ ω + b t x x + (b ω y ) t y ω + x t Q xx x + (y ω ) t Q ω yy y ω ) = p ω b t x x + x t Q xx x + (y ω ) t Q ω yy y ω + (b ω y ) t y ω + c t δ ω
and, if n := n δ = n y and Q yy is diagonal (as in [START_REF] Heredia | Solving electricity market quadratic problems by Branch and Fix Coordination methods[END_REF]), we can write the last bracket as

n i=1 q ω ii (y ω i ) 2 + b ω i y ω i + c i δ ω i .
For notational simplicity in this paragraph we drop the indices. The issue is then how to best represent the quadratic function f (y, δ) = qy 2 + by + cδ by means of a piecewise-linear one. There is an effective way based on ideas developed by Frangioni and Gentile [START_REF] Frangioni | Perspective cuts for a class of convex 0-1 mixed integer programs[END_REF]. The function f (y, δ) is only relevant at points (y, δ) of its (disconnected) domain D = [0, 0] ∪ [y, y] × {1} . Standard branch-and-cut approaches typically solve the continuous relaxation of the mixed problem, where δ ∈ [0, 1] instead of {0, 1}, in order to obtain lower bounds on the optimal value. This makes sense to use the convex envelope of f (y, δ) over D, that is, the convex function with the smallest (in set-inclusion sense) epigraph containing that of f (y, δ). As is showed in [START_REF] Frangioni | Perspective cuts for a class of convex 0-1 mixed integer programs[END_REF] the convex envelope is

h(y, δ) =          0, if (y, δ) = (0, 0) qy 2 δ + by + cδ, if δy ≤ y ≤ δy, for δ ∈ (0, 1] +∞, otherwise.
This function is strongly related with the perspective-function

f (y, δ) = δf (y/δ) of f (y) = qy 2 + by + c, which is convex if f (y) is convex. h(y, δ) ≥ f (y, δ) for 0 < δ ≤ 1, i.e.
h is a tighter objective function than f for the continuous relaxation. As is well-known, every convex function is the point-wise supremum of affine functions. In fact, the epigraph of h is composed of all and only triples (v, y, δ) satisfying δy ≤ y ≤ δy, 0 ≤ δ ≤ 1 and the infinite system of linear inequalities v ≥ (2q y + b)y + (cq y 2 )δ taking y ∈ [y, y]. For each y we have an inequality so-called a perspective cut (PC), which is the unique supporting hyperplane to the function passing by (0, 0) and ( y, 1).

PC formulation (PCF)

PC formulation (PCF) lies in choosing these supporting hyperplanes and using as an objective function the polyhedral function that is the point-wise maximum of the corresponding linear functions. A small set of initial PCs is chosen to solve the problem with the continuous relaxation. When δ * > 0, check whether the solution (v * , y * , δ * ) satisfies the PC for y = y * /δ * ; if not, the obtained cut can be added to PCF.

PCF starts with only two pieces, the ones corresponding with y and y; additional cuts are then dynamically generated when needed as described in the previous paragraph.

Therefore, the objective function of MIQP ω for PCF becomes

p ω b t x x + x t Q xx x + n i=1 v ω i ,
and the initial PCs added to the constraints of MIQP ω for each i ∈ {1, . . . , n}

v ω i ≥ (2q ω ii y i + b ω i )y ω i + (c i -q ω ii y 2 i )δ ω i v ω i ≥ (2q ω ii y i + b ω i )y ω i + (c i -q ω ii y 2 i )δ ω i
We can extend this formulation to scenario clusters obtaining the MIQP p submodels for p = 1, . . . , p in this way min

ω∈Ω p p ω b t x x + x t Q xx x + n i=1 v ω i s.t.: v ω i ≥ (2q ω ii y i + b ω i )y ω i + (c i -q ω ii y 2 i )δ ω i , i ∈ {1, . . . , n}, ω ∈ Ω p v ω i ≥ (2q ω ii y i + b ω i )y ω i + (c i -q ω ii y 2 i )δ ω i , i ∈ {1, . . . , n}, ω ∈ Ω p l a ≤ A δ p x p ≤ u a , l ω t ≤ T ω   δ p x p y ω   ≤ u ω t , ω ∈ Ω p
x p ≥ 0, y ω ∈ [y, y], ω ∈ Ω p , and δ p ∈ {0, 1} n δ 4 Implementation

These methods have been implemented in C++ with the help of Cplex 12.1 to solve only the quadratic subproblems QP p in each node of the BF tree, for each p ∈ {1, . . . , p}, and the QP T N F and QP f subproblems. These algorithmic alternatives have been considered:

⊲ QBFC: coordination of δ in the TNF of the BF trees for clusters p ∈ {1, . . . , p} without using neither OLAs nor PCs, i.e. solving the quadratic subproblems QP p .

⊲ QBFC-PC: coordination of δ in the TNF of the BF trees for clusters p ∈ {1, . . . , p} using PCs. ⊲ QBFC-OLA: coordination of δ in the TNF of the BF trees for clusters p ∈ {1, . . . , p} using OLAs.

For our instances the number of scenarios in each cluster is the same, |Ω p | = |Ω|/ p. Each cluster contains |Ω p | consecutive scenarios, starting from the first one and following in natural order.

Numerical tests

In order to obtain a computational comparison of the performance of the algorithmic alternatives QBFC, QBFC-PC, and QBFC-OLA some computational tests are carried out, which consist in solving two-stage stochastic problems, where the objective function is convex quadratic with linear constraints using QBFC code with those algorithmic choices. Therefore, these problems have a unique primal solution and the duality gap is zero. The tests have been performed on HP Compact with Intel Core 2 Quad Q9550 2.83GHz 4 CPU under Linux 2.6.38-8-generic-pae (x86 64).

The test problems have been randomly generated by using a C++ code developed by this author. This generator provides the scenarios set together with the associated probability of occurrence for two-stage stochastic mixed quadratic problems where Q xy = 1l, Q xy = 0 0, and Q yx = 0 0. Moreover, as can be seen in Table 1, in some problems Q xx = 0 0 and in the rest of problems Q xx = 1l. Also, "# var" means the number of continuous variables, "# bin" the number of binary variables, "# constr" the number of constraints for DEM, see (1)-( 5), and "dens" constraint matrix density %. Table 2 presents the main results of the computational experimentation for given values of the number of scenario clusters. Below the heading QBFC-1 are the times in CPU-seconds used for solving problems with 1 only scenario cluster (i.e. p = 1) and by solving the quadratic subproblem QP p for each node using Cplex; the heading QBFC-5OLA indicates the CPU-times for 5 scenario cluster (i.e. p = 5) and by solving the quadratic subproblem QP p for each node using outer linear approximations (OLA). Finally, PC means perspective cuts are used.

As can be observed in Table 2, the best efficiency is mainly obtained when the problems are solved with 5 clusters, except in the case of QBFC-5 because of the computational cost of solving a quadratic problem in each node of the BF tree for the different clusters. In addition, the outer linear approximations give us a higher efficiency than the perspective cuts.

Conclusions

An algorithm to solve two-stage stochastic quadratic problems based on the Twin Node Family concept involved in the Branch-and-Fix Coordination has been implemented in C++ with the help of Cplex library to solve only the quadratic subproblems. When the problem's structure makes it possible, the algorithm uses perspective cuts or OLA to linearize the MIQ subproblems in each BF tree. The preliminary numerical results show a bit better efficiency with OLA than with PCF.

The path started from this work has the aim of solving nonlinear (nonquadratic) stochastic problems with nonlinear constraints and for two stage or more.

Table 1 .

 1 Test problems Prob. |Ω| # var # bin # constr Qxx dens.

	P1	20 420	20	840	0 0	38
	P2	30 620	20	1240	0 0	37
	P3	40 820	20	1640	0 0	42
	P4	50 1020	20	2040	0 0	43
	P5	60 1220	20	2440	0 0	39
	P6	70 1420	20	2840	0 0	37
	P7	30 930	30	1860	0 0	2
	P8	40 1230	30	2460	0 0	2
	P9	50 1530	30	3060	0 0	1
	P10	60 1830	30	3660	0 0	1
	P11	70 2130	30	4260	0 0	1
	P12 100 3030	30	6060	0 0	1
	P13	30 930	30	1550	1l	13
	P14	40 1230	30	2050	1l	10
	P15	50 1530	30	2550	1l	10
	P16	60 1830	30	3050	1l	9
	P17	70 2130	30	3550	1l	9

⋆ This work was partially supported by the Ministry of Science and Technology of Spain through MICINN Project DPI2008-02153.


Richard H Stockbridge
email: stockbri@uwm.edu

Chao Zhu
email: zhu@uwm.edu

Harvesting in Stochastic Environments: Optimal Policies in a Relaxed Model

Keywords: Singular stochastic control, linear programming, relaxed control AMS subject classification. 93E20, 60J60

This paper examines the objective of optimally harvesting a single species in a stochastic environment. This problem has previously been analyzed in [1] using dynamic programming techniques and, due to the natural payoff structure of the price rate function (the price decreases as the population increases), no optimal harvesting policy exists. This paper establishes a relaxed formulation of the harvesting model in such a manner that existence of an optimal relaxed harvesting policy can not only be proven but also identified. The analysis imbeds the harvesting problem in an infinite-dimensional linear program over a space of occupation measures in which the initial position enters as a parameter and then analyzes an auxiliary problem having fewer constraints. In this manner upper bounds are determined for the optimal value (with the given initial position); these bounds depend on the relation of the initial population size to a specific target size. The more interesting case occurs when the initial population exceeds this target size; a new argument is required to obtain a sharp upper bound. Though the initial population size only enters as a parameter, the value is determined in a closed-form functional expression of this parameter.

Introduction

This paper examines the problem of optimally harvesting a single species that lives in a random environment. Let X be the process denoting the size of the population and Z denote the cumulative amount of the species harvested. We assume X(0-) = x 0 > 0, Z(0-) = 0, and X and Z satisfy dX(t) = b(X(t))dt + σ(X(t))dW (t) -dZ(t), [START_REF] Alvarez | Singular stochastic control in the presence of a state-dependent yield structure[END_REF] in which W (•) is a 1-dimensional standard Brownian motion that provides the random fluctuations in the population's size, and b and σ are real-valued continuous functions. We assume that b and σ are such that in the absence of harvesting the population process X takes values in R + and that ∞ is a natural boundary so that the population will not explode to ∞ in finite time. The boundary 0 may be an exit or a natural boundary point but may not be an entrance point; this indicates that the species will not spontaneously reappear following extinction. Note that X(0) may not equal X(0-) due to an instantaneous harvest Z(0) at time 0 and the process Z is restricted so that ∆Z(t) := Z(t) -Z(t-) ≤ X(t-) for all t ≥ 0. This latter condition indicates that one cannot harvest more of the species than exists. Let r > 0 denote the discount rate and f denote the marginal yield for harvesting. The objective is to select a harvesting strategy Z so as to maximize the expected discounted revenue

J(x 0 , Z) := E x0 τ 0 e -rs f (X(s-))dZ(s) , (2)
where τ = inf {t ≥ 0 : X(t) = 0} denotes the extinction time of the species.

As a result of developments in stochastic analysis and stochastic control techniques, there has been a resurgent interest in determining the optimal harvesting strategies in the presence of stochastic fluctuations (see, e.g., [START_REF] Alvarez | Singular stochastic control in the presence of a state-dependent yield structure[END_REF][START_REF] Lungu | Optimal harvesting from a population in a stochastic crowded environment[END_REF]). In particular, [START_REF] Alvarez | Singular stochastic control in the presence of a state-dependent yield structure[END_REF] examines the current problem using dynamic programming techniques and determines the value function. The paper indicates the lack of an optimal policy in the admissible class of (strict) harvesting policies by commenting that a "chattering" policy will be optimal. The problem of optimal harvesting of a single species in a random environment is also studied in [START_REF] Song | On optimal harvesting problems in random environments[END_REF] in which the model is extended to regime-switching diffusions so as to capture different dynamics such as for drought and non-drought conditions. The paper also adopts a dynamic programming solution approach to determine the value function while at the same time exhibiting ǫ-optimal harvesting policies since, as in the static environment of [START_REF] Alvarez | Singular stochastic control in the presence of a state-dependent yield structure[END_REF], no optimal harvesting policy exists. In light of the complexities of the regime-switching model, it further identifies a condition under which the value function is shown to be continuous and a viscosity solution to the variational inequality.

The focus of this paper is on developing a relaxed formulation for the harvesting problem under which an optimal harvesting control exists and on establishing optimality using a linear programming formulation instead of dynamic programming. In addition, it is sufficient to have a weak solution to (1) rather than placing Lipschitz and polynomial growth conditions on the coefficients b and σ that guarantee existence of a strong solution. Intuitively, relaxation completes the space of admissible harvesting rules by allowing measure-valued policies. A benefit of the linear programming solution methodology is the analysis concentrates on the optimal value for a single, fixed initial condition, rather than seeking the value function and thus no smoothness properties need to be established about the value as a function of the initial position.

To set the stage for the relaxed singular control formulation of the model, let D = C 2 c (R + ) and for a function g ∈ D, define the operators A and B by

Ag(x) = 1 2 σ 2 (x)g ′′ (x) + b(x)g ′ (x), and (3)
Bg(x, z) = g(x-z)-g(x) z , if z > 0, -g ′ (x), if z = 0, (4)
where x, z ∈ R + . Itô's formula then implies

g(X(t)) = g(x 0 ) + t 0 Ag(X(s)) ds + t 0
Bg(X(s), ∆Z(s)) dZ(s)

+ t 0 σ(X(s))g ′ (X(s)) dW (s), ∀g ∈ D.
It therefore follows that for any g ∈ D

g(X(t)) -g(x 0 ) - t 0 Ag(X(s)) ds - t 0 Bg(X(s), ∆Z(s)) dZ(s) (5)
is a mean 0 martingale. In fact, requiring (5) to be a martingale for a sufficiently large collection of functions g is a way to characterize the processes (X, Z) which satisfy [START_REF] Alvarez | Singular stochastic control in the presence of a state-dependent yield structure[END_REF]. We turn now to a precise formulation of the model in which the processes are relaxed solutions of a controlled martingale problem for the operators (A, B).

Formulation of the Relaxed Model

For a complete and separable metric space S, we define M (S) to be the space of Borel measurable functions on S, B(S) to be the space of bounded, measurable functions on S, C(S) to be the space of continuous functions on S, C(S) to be the space of bounded, continuous functions on S, M(S) to be the space of finite Borel measures on S, and P(S) to be the space of probability measures on S. M(S) and P(S) are topologized by weak convergence. Recall, the amount of harvesting is limited by the size of the population. Define R = {(x, z) : 0 ≤ z ≤ x, x ≥ 0}; R denotes the space on which the paired process (X, Z) evolves when considering solutions of [START_REF] Alvarez | Singular stochastic control in the presence of a state-dependent yield structure[END_REF].

The formulation of the population model in the presence of "relaxed" harvesting policies adapts the relaxed formulation for singular controls given in [START_REF] Kurtz | Stationary solutions and forward equations for controlled and singular martingale problems[END_REF] to the particulars of the harvesting problem. This adaptation sets the state space E to be R + and the control space

U = R + , with U = R ⊂ R + × R + .
Let X be an R + -valued process and Γ be an

L(R)-valued random variable. Let Γ t denote the restriction of Γ to R × [0, t]. Then (X, Γ ) is a relaxed solution of the harvesting model if there exists a filtration {F t } such that (X, Γ t ) is {F t }-progressively measurable, X(0-) = x 0 , and for every g ∈ D, g(X(t)) -g(x 0 ) - t 0 Ag(X(s)) ds - R×[0,t] Bg(x, z) Γ (dx × dz × ds) (6)
is an {F t }-martingale, in which the operators A and B are given by ( 3) and ( 4), respectively. Throughout the paper we assume that a relaxed solution (X, Γ ) exists and is strong Markov. Let A denote the set of measures Γ for which there is some X such that (X, Γ ) is a relaxed solution of the harvesting model.

We turn now to the extension of the reward criterion (2) to the relaxed framework. Specifically, f : R + → R + represents the instantaneous marginal yield accrued from harvesting. Assume f is continuous and non-increasing with respect to x. Thus f (x) ≥ f (y) whenever x ≤ y; this assumption indicates that the price when the species is plentiful is smaller than when it is rare. Moreover, we assume 0 < f (0) < ∞. Let (X, Γ ) be a solution to the harvesting model [START_REF] Lungu | Optimal harvesting from a population in a stochastic crowded environment[END_REF]. Let S = (0, ∞) denote the survival set of the species. Then the expected total discounted value from harvesting is

J(x 0 , Γ ) := E R×[0,τ ] e -rs f (x) Γ (dx × dz × ds) . ( 7
)
The goal is to maximize the expected total discounted value from harvesting over relaxed solutions (X, Γ ) of the harvesting model and to find an optimal harvesting strategy Γ * . Thus, we seek

V (x 0 ) = J(x 0 , Γ * ) := sup Γ ∈A J(x 0 , Γ ). (8)
We emphasize that the initial position x 0 is merely a parameter in the problem and that V is not to be viewed as a function with any particular properties but merely is the value of the harvesting problem when the initial population size is x 0 . We do, however, obtain the value in functional form for x 0 in two regions.

Linear Programming Formulation and Main Result

Throughout this paper, we assume the equation (A -r)u(x) = 0 has two fundamental solutions ψ and φ, where ψ is strictly increasing and φ is strictly decreasing; without loss of generality we may assume ψ(0) = 0 (see [START_REF] Alvarez | Singular stochastic control in the presence of a state-dependent yield structure[END_REF]).

The main result of this paper is summarized in the following theorem.

Theorem 1 Assume that there exists some b ≥ 0 such that

(i) f (x) ψ ′ (x) ≤ f ( b) ψ ′ ( b) , ∀x ≥ 0, (ii) the function f /ψ ′ is nonincreasing on [ b, ∞), and (iii) the function f is continuously differentiable on ( b, ∞). Put b * = inf{ b ≥ 0 : b satisfies (i)-(iii)}.
Then the value is given by

V (x 0 ) = f (b * ) ψ ′ (b * ) ψ(x 0 ∧ b * ) + x0∨b * b * f (y)dy (9)
and an optimal relaxed harvesting policy is given by We begin the task of reformulating the harvesting problem with the following observation. Let τ be any {F t }-stopping time. The optional sampling theorem along with the requirement that (6) be a mean 0 martingale for each g ∈ D implies

Γ * (dx × dz × dt) = I (b * ,∞) (x 0 )λ [b * , x0] (dx)δ {0} (dz)δ {0} (dt) + Γ b * (dx × dz × dt), (10
e -r(t∧τ ) g(X(t ∧ τ )) -g(x 0 ) - t∧τ 0 e -rs [A -r]g(X(s)) ds - R×[0,t∧τ ] e -rs Bg(x, z) Γ (dx × dz × ds)
is also a martingale. Recall g ∈ D means g has compact support and hence is bounded. So taking expectations and letting t → ∞ yields

g(x 0 ) = E e -rτ I {τ <∞} g(X(τ )) -E τ 0 e -rs [A -r]g(X(s)) ds (11) -E R×[0,τ ] e -rs Bg(x, z) Γ (dx × dz × ds) .
The initial analysis takes τ = τ ; later we will need (11) for a different stopping time.

The measures involved in the infinite-dimensional linear program are expected discounted occupation measures corresponding to relaxed solutions (X, Γ ) of the harvesting model. Indeed, for any Borel measurable G 1 ⊂ S and G ⊂ R, we define

µ τ (G 1 ) = E e -rτ I G1 (X(τ ))I {τ <∞} , µ 0 (G 1 ) = E τ 0 e -rs I G1 (X(s))ds , µ 1 (G) = E R×[0,τ ] e -rs I G (x, z)Γ (dx × dz × ds) . (12)
Using these measures, the singular control problem of maximizing (7) over relaxed solutions of the harvesting problem (6) can be written in the form

         Maximize f dµ 1 , subject to gdµ τ -(A -r)gdµ 0 -Bgdµ 1 = g(x 0 ), ∀g ∈ D, µ τ , µ 0 , µ 1 ∈ M(S), µ τ (S) ≤ 1, µ 0 (S) ≤ 1 r . (13)
Since each relaxed solution (X, Γ ) defines measures µ τ , µ 0 and µ 1 by (12), the harvesting problem is embedded in (13). There might be feasible measures which do not arise in this manner. Consequently, letting V lp (x 0 ) denote the value of the LP problem (13) with initial condition X(0-) = x 0 > 0, we have V (x 0 ) ≤ V lp (x 0 ).

The Proof of Theorem 1

This section is devoted to the proof of Theorem 1 and involves two steps.

Step 1: Universal Upper Bound

The proof follows along the lines of the arguments used in [START_REF] Helmes | Thinning and harvesting in stochastic forest models[END_REF]. The general argument involves finding an upper bound for V lp (x 0 ) by reducing the number of constraints in the linear program (13). We state the results and leave the proofs to the reader.

Proposition 2 Let b * be defined as in Theorem 1. Then for every x 0 ≥ 0,

V (x 0 ) ≤ f (b * ) ψ ′ (b * ) ψ(x 0 ). ( 14
)
Notice the bound in ( 14) holds for all initial positions x 0 . The following result shows that this bound is sharp for x 0 ≤ b * . (15)

Then J(x 0 , L b * ) = J(x 0 , Γ b * ) = f (b * ) ψ ′ (b * ) ψ(x 0 ).
Since ∆L b * (s) = 0 for every s ≥ 0, an optimal strategy is to harvest just enough of the population (using the local time of X * at b * ) so that the population size "reflects" at b * .

The value function has been determined for intial population sizes x 0 that are smaller than b * . It therefore remains to prove the validity of (9) when x 0 > b * .

Step 2: Return of Stochasticity for a Refined Upper Bound

This step is the more interesting of the two and requires a new argument and also a different type of harvesting policy than appears in the literature.

When dealing with singular control problems, one usually takes the so-called reflection strategy, namely,

Z(t) = (x 0 -b * ) + + L b * (t)
, where one follows an immediate jump from x 0 to b * by using the local time process L b * at b * . Such a reflection strategy is used in [START_REF] Choulli | A diffusion model for optimal dividend distribution for a company with constraints on risk control[END_REF], [START_REF] Pham | Continuous-time stochastic control and optimization with financial applications[END_REF] and others. The corresponding income is

J(x 0 , Z) = f (x 0 )(x 0 -b * ) + f (b * ) ψ ′ (b * ) ψ(b * ).
When f is strictly decreasing, the reflection strategy is not optimal. Our purpose is to find an optimal relaxed harvesting strategy.

To develop a sharp upper bound, it is beneficial to revisit the definitions of the occupation measures in (12) so that the connection between the measures, the initial position and the harvesting strategy is more clearly displayed. Let x 0 ∈ R + and (X, Γ ) be a relaxed solution of the harvesting model. Modify the notations of the measures to indicate their dependence on x 0 and Γ by writing µ τ (G; x 0 , Γ ), µ 0 (G; x 0 , Γ ) and µ 1 (G; x 0 , Γ ).

Proposition 4 For x 0 > b * , V (x 0 ) ≤ x0 b * f (y) dy + f (b * ) ψ ′ (b * ) • ψ(b * ). ( 16
)
Proof. The proof of ( 16) is broken into two parts, with a technical lemma between the parts.

Γ τ b * (G × [0, t]) = I {τ b * <τ } Γ (G × [τ b * , τ b * + t]), G ∈ B(R), t ≥ 0.
Notice that Γ τ b * captures all harvesting using the measure Γ from time τ b * onwards. Also define the measures µ 0,τ b * and µ 1,τ b * by

µ 0,τ b * (G; x 0 , Γ ) = E x0 τ b * 0 e -rs I G (X(s)) ds , µ 1,τ b * (G; x 0 , Γ ) = E x0 R×[0,τ b * )
e -rs I G (x, z) Γ (dx × dz × ds) .

Note carefully that any harvesting at the time τ b * is excluded from the measure µ 1,τ b * . Also observe that the total mass of µ 0,τ b * equals r -1 (1

-E x0 [e -rτ b * ]).
Using the strong Markov property of (X, Γ ), for each G ∈ B(R) it follows that

E x0 R×[0,τ ] e -rs I G (x, z) Γ (dx × dz × ds) = E x0 R×[0,τ b * ) e -rs I G (x, z) Γ (dx × dz × ds) +E x0 E x0 I {τ b * <τ } R×[τ b * ,τ ] e -rs I G (x, z) Γ (dx × dz × ds) F τ b * = E x0 R×[0,τ b * ) e -rs I G (x, z) Γ (dx × dz × ds) +E x0 e -rτ b * I {τ b * <τ } E X(τ b * ) R×[0,τ ] e -rs I G (x, z) Γ τ b * (dx × dz × ds) .
As a result, for each G ∈ B(R), this identity can be written in terms of the measures as

µ 1 (G; x 0 , Γ ) = µ 1,τ b * (G; x 0 , Γ ) + E x0 e -rτ b * I {τ b * <τ } µ 1 (G; X(τ b * ), Γ τ b * ) .
Notice, in particular, that the expectation term involves the measure µ 1 evaluated at the random initial position X(τ b * ). Hence

f (y) µ 1 (dy; x 0 , Γ ) = f (y) µ 1,τ b * (dy; x 0 , Γ ) + E x0 e -rτ b * I {τ b * <τ } f (y) µ 1 (dy; X(τ b * ), Γ τ b * ) ≤ f (y) µ 1,τ b * (dy; x 0 , Γ ) + E x0 e -rτ b * I {τ b * <τ } f (b * ) ψ ′ (b * ) • ψ(b * ), (17)
(A -r)h(x) ≤ r f (b * ) ψ ′ (b * ) ψ(b * ), for x ≥ b * and (18)
-Bh(x, z) ≥ f (x), for all (x, z) ∈ R. (19)
Proof. Since by assumption the function f /ψ ′ is nonincreasing and differentiable on (b * , ∞), we have

0 ≥ d dx f (x) ψ ′ (x) = f ′ (x)ψ ′ (x) -f (x)ψ ′′ (x) (ψ ′ (x)) 2 , x > b * .
But ψ is strictly increasing and so ψ ′ (x) > 0. Hence it follows that

f ′ (x)ψ ′ (x) - f (x)ψ ′′ (x) ≤ 0, or equivalently f ′ (x) ≤ f (x) ψ ′ (x) ψ ′′ (x), for x > b * . It then follows that for each x > b * (A -r)h(x) ≤ 1 2 σ 2 (x) f (x) ψ ′ (x) ψ ′′ (x) + b(x)f (x) -r f (x) ψ ′ (x) (ψ(x) -ψ(b * )) = f (x) ψ ′ (x) 1 2 σ 2 (x)ψ ′′ (x) + b(x)ψ ′ (x) -rψ(x) + r f (x) ψ ′ (x) ψ(b * ) = r f (x) ψ ′ (x) ψ(b * ) ≤ r f (b * ) ψ ′ (b * ) ψ(b * ).
Turning to a consideration of (19), since f is nonincreasing, for any 0 ≤ x 1 < x 2 , we have

f (x 2 )[x 2 -x 1 ] ≤ x2 x1 f (y)dy = h(x 2 ) -h(x 1 ).
Hence it follows that for (x, z) ∈ R, we have

-Bh(x, z) = h ′ (x) = f (x), if z = 0 h(x)-h(x-z) z ≥ f (x), if z > 0.
The relation ( 19) is therefore established.

Part 2:

The goal is of this part of the proof is to estimate f (y) µ 1,τ b * (dy × dz) of ( 17). Using Itô's formula, one obtains for each t > 0, We have derived an upper bound for the value V (x 0 ) in Proposition 4. The following proposition exhibits an optimal relaxed harvesting policy. The proof is left to the reader. We observe that the manner in which this optimal harvesting policy differs from the typical "reflection" strategy occurs at the initial time. Whereas the reflection strategy has the process X instantaneously jump from x 0 to b * , the optimal relaxed harvesting policy obtains this relocation in an instantaneous but continuous manner.

-E R×[0,t∧τ b * ) e -rs Bh(x, z) Γ (dx × dz × ds) = h(x 0 ) -E e -r(t∧τ b * ) h(X((t ∧ τ b * )-)) + E t∧τ b * 0 e -rs [A -r]h(X (
Finally we note that the combination of Propositions 2 and 6 establishes Theorem 1. Moreover, the optimal relaxed harvesting policy in (10) unifies the two cases.

  ) where λ [b * , x0] (•) denotes Lebesgue measure on [b * , x 0 ] and Γ b * is defined in Proposition 6.

Proposition 3 0 I

 30 For x 0 ≤ b * , let L b * denote the local time process of X at b * . Define the random measure Γ b * for Borel measurable G ⊂ R and t ≥ 0 by Γ b * (G × [0, t]) = t G (X(s-), ∆L b * (s)) dL b * (s).

Part 1 :

 1 Define the stopping time τ b * = inf{t ≥ 0 : X(t) ≤ b * } to be the first time the process X takes value at most b * and note that τ b * ≤ τ . For the harvesting measure Γ , define Γ τ b * by

  in which the inequality follows from Step 1. This concludes Part 1 of the proof. Part 2 concentrates on estimating the first term of the right-hand side of (17); a technical lemma is required. Lemma 5 Assume the conditions in Theorem 1. Define the function h by h(x) := x b * f (y)dy for x ≥ 0. Then the following estimates hold:

  s)) ds , in which the deliberate choice of the half-open interval [0, t ∧ τ b * ) in the integral with respect to Γ leads to the use of X((t ∧ τ b * )-) for the location of the process just before any harvest occurs at time τ b * . This is extremely important since h(X((t ∧ τ b * )-)) ≥ 0 and hence the right-hand side is not decreased by dropping the first expectation. Letting t → ∞ yields-E R×[0,τ b * ) e -rs Bh(x, z) Γ (dx × dz × ds) ≤ h(x 0 ) + E τ b * 0 e -rs [A -r]h(X(s)) ds .Using the estimates (18) and (19) and the definition of the measures µ 1,τ b * and µ 0,τ b * , we obtainR f (y) µ 1,τ b * (dy × dz; x 0 , Γ ) ≤ -R Bh(y, z ) µ 1,τ b * (dy × dz; x 0 , Γ ) ≤ h(x 0 ) + r • f (b * )ψ(b * ) ψ ′ (b * ) µ 0,τ b * (dx; x 0 , Γ ) = h(x 0 ) + 1 -E x0 e -rτ b * f (b * ) ψ ′ (b * ) • ψ(b * ), (20) in which the last equality follows from the mass of µ 0,τ b * . Combining (17) and (20) produces the desired relation f (y) µ 1 (dy × dz; x 0 , Γ ) ≤ x0 b * f (y) dy + f (b * ) ψ ′ (b * ) • ψ(b * ).

Proposition 6

 6 Let λ [b * ,x0] (•) denote Lebesgue measure on [b * , x 0 ]. Also let L b * denote the local time process of Proposition 3 with x 0 taken to be b * and denote by Γ b * the random measure defined in (15). Finally, define the relaxed harvesting strategy byΓ * (dx × dz × dt) = λ [b * ,x0] (dx)δ {0} (dz)δ {0} (dt) + Γ b * (dx × dz × dt). Then V (x 0 ) = J(x 0 , Γ * ) = x0 b * f (y)dy + f (b * ) ψ ′ (b * ) ψ(b * ).(21)

Acknowledgements. This research was supported in part by the U.S. National Security Agency under Grant Agreement Number H98230-09-1-0002, the National Science Foundation under DMS-1108782, a grant from the UWM Research Growth Initiative, and City University of Hong Kong (SRG) 7002677. The United States Government is authorized to reproduce and distribute reprints notwithstanding any copyright notation herein.