Michael Burger

Matthias Gerdts

Simone Göttlich

Michael Herty

Dynamic programming approach for discrete-valued time discrete optimal control problems with dwell time constraints

Keywords: dynamic programming, mixed-integer optimization, dwell time constraints

The article discusses a numerical approach to solve optimal control problems in discrete time that involve continuous and discrete controls. Special attention is drawn to the modeling and treatment of dwell time constraints. For the solution of the optimal control problem in discrete time, a dynamic programming approach is employed. A numerical example is included that illustrates the impact of dwell time constraints in mixed integer optimal control.

Introduction

Mixed-integer optimal control is a field of increasing importance as practical applications often include discrete decisions in addition to continuous-valued control variables. Examples of such problems can be found, e.g., in [START_REF] Göttlich | Modeling and optimizing traffic light settings on road networks[END_REF][START_REF] Howlett | Optimal Strategies for the Control of a Train[END_REF][START_REF] Gerdts | Solving Mixed-Integer Optimal Control Problems by Branch&Bound: A Case Study from Automobile Test-Driving with Gear Shift[END_REF][START_REF] Gerdts | A variable time transformation method for mixed-integer optimal control problems[END_REF][START_REF] Siburian | Numerical Methods for Robust, Singular and Discrete Valued Optimal Control Problems[END_REF][START_REF] Sager | Numerical methods for optimal control with binary control functions applied to a Lotka-Volterra type fishing problem[END_REF][START_REF] Kirches | Time-optimal control of automobile test drives with gear shifts[END_REF][START_REF] Sager | Fast solution of periodic optimal control problems in automobile test-driving with gear shifts[END_REF]. One way to approach mixed-integer optimal control problems is by solving necessary optimality conditions provided by the well-known maximum principle. These necessary conditions are valid even for discrete control sets. In [START_REF] Khmelnitsky | A combinatorial, graph-based solution method for a class of continuous-time optimal control problems[END_REF] a graph-based solution method exploiting the maximum principle was developed, which is however limited to single-state problems. The drawback of maximum principle based methods is that a very good initial guess of the switching structure is needed which is often not available for practical applications. Direct discretization methods based on variable time transformations as in [START_REF] Lee | Control parametrization enhancing technique for optimal discrete-valued control problems[END_REF][START_REF] Lee | Control parameterization enhancing technique for time optimal control problems[END_REF][START_REF] Lee | An Optimal Control Approach to Nonlinear Mixed Integer Programming Problems[END_REF][START_REF] Teo | The control parameterization enhancing transform for constrained optimal control problems[END_REF][START_REF] Siburian | Numerical Methods for Robust, Singular and Discrete Valued Optimal Control Problems[END_REF][START_REF] Gerdts | Solving Mixed-Integer Optimal Control Problems by Branch&Bound: A Case Study from Automobile Test-Driving with Gear Shift[END_REF][START_REF] Gerdts | A variable time transformation method for mixed-integer optimal control problems[END_REF], or direct discretization methods based on relaxations and sum-up-rounding strategies as in [START_REF] Sager | Numerical methods for mixed-integer optimal control problems[END_REF][START_REF] Sager | Reformulations and algorithms for the optimization of switching decisions in nonlinear optimal control[END_REF][START_REF] Sager | Direct methods with maximal lower bound for mixed-integer optimal control problems[END_REF] have shown their ability to solve difficult real-world examples. Extensions towards mixed control-state constraints, which depend on the discrete variable and lead to vanishing constraints, can be found in [START_REF] Kirches | Fast numerical methods for mixed-integer nonlinear model-predictive control[END_REF]Chapter 5] and [START_REF] Palagachev | Mathematical Programs with Blocks of Vanishing Constraints Arising in Discretized Mixed-Integer Optimal Control Problems[END_REF]. The dynamic programming approach in [START_REF] Gerdts | Optimal control of ODEs and DAEs[END_REF]Chapter 7] allows to consider switching costs to avoid frequent switches of the discrete-valued controls. This paper particularly addresses so-called dwell time constraints in mixed-integer optimal control problems in discrete time. Dwell times are an important aspect in many engineering applications, especially in the context of switched systems. After a switch has occured, the physical process typically requires a certain time period to recover and to return to normal operation. During this dwell time the process is limited in its operation and follows specified dynamics. An example for such a behavior is a truck with gear shifts. As the clutch is engaged, the motor torque does not arrive at the wheels for a short time period. Only after the clutch is released again, the motor torque is distributed at the wheels. Throughout the paper, we focus on optimal control problems in discrete time, which exist in their own rights but often they are obtained as discretizations of optimal control problems in continuous time. The purpose of the paper is twofold. Firstly, a model taking into account dwell time is suggested. Herein, we will use delays in the dynamics. Secondly, it is shown that a dynamic programming principle applies and can be used for numerical computations, if the state dimension is low. The paper is organized as follows. Section 2 defines the mixed-integer optimal control problem in discrete time with dwell time constraints. Section 3 discusses a dynamic programming approach to solve the problems. Section 4 discusses an illustrative numerical example.

Modeling dwell time constraints in mixed-integer optimal control in discrete time

Let a fixed grid G N = {t 0 < t 1 < . . . < t N } with N ∈ N be given. For i = 0, . . . , N and x ∈ R nx let ∅ = X (t i ) ⊆ R nx and ∅ = U(t i , x) ⊆ R nu be closed connected sets and V = {v 1 , . . . , v M } a discrete set of vectors v j ∈ R nv , j = 1, . . . , M . On the grid G N , grid functions are introduced. The grid function x : G N -→ R nx is called state and it is restricted by the state constraints x(t i ) ∈ X (t i ), i = 0, . . . , N . The grid function u : G N -→ R nu is called real-valued control and it is restricted by the control constraints u(t i ) ∈ U, i = 0, . . . , N . The grid function v : G N → R nv is called discrete-valued control and it is restricted by the constraints v(t i ) ∈ V, i = 0, . . . , N . We say the discrete-valued control switches at time point

t i ∈ G N with i ∈ {0, . . . , N -1}, if v(t i ) = v(t i+1
). Switches in the grid function v can be measured with the help of the discrete variation dv : G N → R nv defined by

dv(t i ) = v(t i+1 ) -v(t i ), for i = 0, . . . , N -1, v(t N ) -v(t N -1 ), for i = N.
A jump of v at t i occurs if and only if dv(t i ) = 0. In order to define dwell time constraints let L ∈ N denote a number of steps and ū ∈ U a fixed control vector. Definition 1. A dwell time constraint of time horizon length L ∈ N and control ū applies, if and only if the following conditions are met:

(a) v switches at t i for some i ∈ {0, . . . , N -1} and (b) u(t i+ ) = ū for = 1, . . . , L.

In other words, a dwell time constraint is a logical implication of type

v(t i ) = v(t i+1 ) =⇒ u(t i+ ) = ū for = 1, . . . , L,
i.e. if the discrete control switches, then control u is fixed to ū for a defined future time horizon of length L.

Remark 1. For notational simplicity we fixed the full control vector u in Definition 1 (b) to a specified value ū for all time points t i+ , = 1, . . . , L. This can be generalized such that only some components of u are fixed to specified values which may even vary with the index . This more general setting would cause a more technical notation, but otherwise can be treated with the same techniques as below.

Let the state evolve according to the discrete dynamics

x(t i+1 ) = f (t i , x(t i ), u(t i ), v(t i )), i = 0, 1, . . . , N -1.
Imposing a dwell time constraint of length L and value ū implies that the future state depends not just on the current control value but on the history of the control v, i.e.

u(t i ) ∈ Ū(v(t i ), v(t i-1 ), . . . , v(t i-L )) := U, if v(t i-) = v(t i--1 ) for = 0, 1, . . . , L -1, {ū}, otherwise.
Herein and throughout the paper we use the convention

v(t -) := v(t 0 ) for = 1, . . . , L. (1)
We arrived at the following mixed-integer optimal control problem in discrete time with a dwell time constraint of length L with value ū:

Problem 1 (DMIOCP). Minimize F (x, u, v) := ϕ(x(t N )) + N -1 i=0 f 0 (t i , x(t i ), u(t i ), v(t i ))
subject to the constraints

x(t i+1 ) = f (t i , x(t i ), u(t i ), v(t i )), i = 0, 1, . . . , N -1, x(t 0 ) = x 0 , x(t i ) ∈ X (t i ), i = 0, 1, . . . , N, u(t i ) ∈ Ū(v(t i ), v(t i-1 ), . . . , v(t i-L )), i = 0, 1, . . . , N -1, v(t i ) ∈ V, i = 0, 1, . . . , N -1.
Note that the control v(t N ) does not have any influence in the problem. Consequently, only switches at the time points t 0 , . . . , t N -1 are taken into account.

In the sequel we will make use of the convention v N = v N -1 whenever useful.

A dynamic programming approach

DMIOCP is embedded into a family of perturbed problems as follows. To this end we use the notation v L (k) := (v k-1 , . . . , v k-L ) for L ∈ N and some index k with values v k-∈ V for = 1, . . . , L. We denote with V L the cartesian product

V × • • • × V (L times). Problem 2 (DMIOCP(t k , x k , v L (k) )). For a given index t k ∈ G N , x k ∈ R nx , and v L (k) ∈ V L , minimize F k (x, u, v) := ϕ(x(t N )) + N -1 i=k f 0 (t i , x(t i ), u(t i ), v(t i )) subject to x(t i+1 ) = f (t i , x(t i ), u(t i ), v(t i )), i = k, . . . , N -1, x(t k ) = x k , x(t i ) ∈ X (t i ), i = k, . . . , N, u(t i ) ∈ Ū(v(t i ), v(t i-1 ), . . . , v(t i-L )), i = k, . . . , N -1, v(t i ) ∈ V, i = k, . . . , N -1, v(t k-) = v k-, = 1, . . . , L.                (2)
Definition 2. The function W :

G N × R nx × V L -→ R that assigns to every (t k , x k , v L (k) ) ∈ G N × R nx × V L the optimal value of DMIOCP(t k , x k , v L (k)
) is called optimal value function of DMIOCP, i.e.

W (t k , x k , v L (k) ) := inf x,u,v with (2) F k (x, u, v), if DMIOCP(t k , x k , v L (k) ) is feasible, ∞,
otherwise.

The next theorem establishes Bellman's optimality principle.

Theorem 1. Let (x, û, v) be optimal grid functions for DMIOCP and let

G k N := G N \{t 0 , . . . , t k-1 }. Then, the restrictions on G k N given by x| G k N , û| G k N , and v| G k N are optimal for DMIOCP(t k , x(t k ), vL (k) ) with vL (k) = (v(t k-1 ), . . . , v(t k-L
)) for all k = 0, . . . , N subject to the convention (1).

Proof. Assume that the restrictions x|

G k N , û| G k N , and v| G k N are not optimal for DMIOCP(t k , x(t k ), vL (k)
) for some k ∈ {0, . . . , N -1}. Then there exist feasible grid functions x :

G k N -→ R nx , ũ : G k N -→ R nu , and ṽ : G k N -→ V for DMIOCP(t k , x(t k ), vL (k) ) with F k (x, ũ, ṽ) < F k (x, û, v), (3)
x(t k ) = x(t k ), and ṽ(t k-) = v(t k-) for = 1, . . . , L. Hence, the trajectories

x : G N -→ R nx , u : G N -→ R nu , and v : G N -→ V with (x(t i ), u(t i ), v(t i )) := (x(t i ), û(t i ), v(t i )) for i = 0, 1, . . . , k -1, (x(t i ), ũ(t i ), ṽ(t i )) for i = k, k + 1, . . . , N,
are feasible for DMIOCP and satisfy

F (x, u, v) = ϕ(x(t N )) + N -1 i=0 f 0 (t i , x(t i ), u(t i ), v(t i )) = ϕ(x(t N )) + k-1 i=0 f 0 (t i , x(t i ), û(t i ), v(t i )) + N -1 i=k f 0 (t i , x(t i ), ũ(t i ), ṽ(t i )) < ϕ(x(t N )) + k-1 i=0 f 0 (t i , x(t i ), û(t i ), v(t i )) + N -1 i=k f 0 (t i , x(t i ), û(t i ), v(t i )) = ϕ(x(t N )) + N -1 i=0 f 0 (t i , x(t i ), û(t i ), v(t i )) = F (x, û, v),
where

is exploited. This contradicts the optimality of x(•), û(•), and v(•).

From the definition of the optimal value function one immediately obtains

W (t N , x N , v L (N ) ) = ϕ(x N ), if x N ∈ X (t N ), ∞, otherwise. (4)
Bellman's optimality principle allows to derive the following result.

Theorem 2. For all

(t k , x k , v L (k) ) ∈ G N × X (t k ) × V L and k = N -1, .
. . , 0, the optimal value function in Definition 2 satisfies the recursion

W (t k , x k , v L (k) ) = inf u∈ Ū (v,v k-1 ,...,v k-L ),v∈V f 0 (t k , x k , u, v) +W (t k+1 , f (t k , x k , u, v), v L (k+1) ) where v L (k) = (v k-1 , . . . , v k-L ) and v L (k+1) = (v, v k-1 , . . . , v k-L+1
) and convention (1) is used. Herein, W at t = t N is given by ( 4) and the convention

W (t k , x k , v L (k) ) = ∞ is used whenever (x k , v L (k) ) ∈ X (t k ) × V L . Proof. Let (t k , x k , v L (k) ) ∈ G N × X (t k ) × V L and k ∈ {0, . . . , N -1} be given. If (x k , v L (k) ) ∈ X (t k ) × V L , then W (t k , x k , v L (k) ) = ∞ by definition. If f (t k , x k , u, v) ∈ X (t k+1 ) for all u ∈ Ū(v, v k-1 , . . . , v k-L ) and all v ∈ V, then DMIOCP(t k , x k , v L (k) ) and DMIOCP(t k+1 , f (t k , x k , u, v), v L (k+1) ) are infeasi- ble for every u ∈ Ū(v, v k-1 , . . . , v k-L ) and every v ∈ V and hence W (t k , x k , v L (k) ) = ∞. For arbitrary u ∈ Ū(v, v k-1 , . . . , v k-L ) and v ∈ V with f (t k , x k , u, v) ∈ X (t k+1
) the definition of the optimal value function yields

W (t k , x k , v L (k) ) ≤ f 0 (t k , x k , u, v) + W (t k+1 , f (t k , x k , u, v), v L (k+1) ).
Taking the infimum over all u ∈ Ū(v, v k-1 , . . . , v k-L ) and v ∈ V yields the first part of the assertion. Now let ε > 0 and feasible grid functions x, ũ, ṽ with x(t k ) = x k and

F k (x, ũ, ṽ) ≤ W (t k , x k , ṽL (k) ) + ε
with ṽL (k) = (ṽ(t k-1 ), . . . , ṽ(t k-L )) be given. Then,

W (t k , x k , ṽL (k) ) ≥ F k (x, ũ, ṽ) -ε ≥ f 0 (t k , x k , ũ(t k ), ṽ(t k )) +W (t k+1 , f (t k , x k , ũ(t k ), ṽ(t k )), ṽL (k+1) ) -ε ≥ inf u∈ Ū (v,ṽ(t k-1 ),...,ṽ(t k-L )),v∈V f 0 (t k , x k , u, v) +W (t k+1 , f (t k , x k , u, v), v L (k+1) ) -ε,
where ṽL (k+1) = (ṽ(t k ), . . . , ṽ(t k-L+1 )) and v L (k+1) = (v, ṽ(t k-1 ), . . . , ṽ(t k-L+1 )). As ε > 0 was arbitrary, the assertion follows.

Theorem 2 allows to deduce the following dynamic programming algorithm. Note that owing to the presence of v L (k) in DMIOCP, the recursive formula for the optimal value function is not of standard type. The vector v L (k) contains the history of the discrete-valued control, which influences the decision at stage k. The vector v L (k) can be interpreted as additional states in the argument of W and leads to an additional runtime factor of M L when compared to standard dynamic programming. In general this is a prohibitively large factor, but often M and L are moderate numbers. For instance, a binary control (M = 2) and L = 3 leads to a factor of 8.

Algorithm: Dynamic Programming

Init: Let G N = {t 0 < t 1 < . . . < t N } be given. Set W (t N , x N , v L (N ) ) = ϕ(x N ), if x N ∈ X (t N ), ∞, otherwise for all x N ∈ R nx and all v L (N ) = (v N -1 , . . . , v N -L ) ∈ V L .
Phase 1 (Backward solution):

For k = N -1, . . . , 0 compute W (t k , x k , v L (k) ) = inf u∈ Ū (v,v k-1 ,...,v k-L ),v∈V f 0 (t k , x k , u, v) +W (t k+1 , f (t k , x k , u, v), v L (k+1) ) (5) with v L (k+1) = (v, v k-1 , . . . , v k-L+1 ) for all x k ∈ R nx and all v L (k) ∈ V L .
Phase 2 (Forward solution): Set x(t 0 ) = x 0 . For k = 0, . . . , N -1 find

(û(t k ), v(t k )) = arg min u∈ Ū (v,v(t k-1 ),...,v(t k-L )),v∈V f 0 (t k , x(t k ), u, v) +W (t k+1 , f (t k , x(t k ), u, v), vL (k+1) ) (6)
with vL

(k+1) = (v, v(t k-1 ), . . . , v(t k-L+1 )) and set x(t k+1 ) = f (t k , x(t k ), û(t k ), v(t k )).

Practical issues in the dynamic programming algorithm

The implementation of the dynamic programming approach works on a compact state space

Ω = {x ∈ R nx | x ≤ x ≤ x u } ⊂ R nx with lower and upper bounds x , x u ∈ R nx and -∞ < x < x u < ∞.
The bounds should be chosen such that all realistically relevant trajectories are contained in this set. The state space Ω is discretized with an equidistant partition

Ω Nx =    (x 1 , . . . , x nx ) ∈ Ω x j = x ,j + i xu,j -x ,j Nx , j = 1, . . . , n x , i ∈ {0, . . . , N x }    with N x ∈ N.
During the backward solution and forward solution phases of the dynamic programming algorithm values of the value function at points x = f (t k , x k , u, v) not in Ω Nx are approximated by polynomial interpolation of the value function on a cell of Ω Nx that contains x.

The minimization in ( 5) and ( 6) is usually done by a sufficiently dense discretization of the control set Ū and complete enumeration of all possible values. The v-component is not crucial as V is supposed to be a finite discrete set. For a fixed v the minimization w.r.t. the component u could be carried out using methods from nonsmooth optimization, e.g. Bundle methods or subgradient methods. However, such methods may result in a local minimum and therefore the previously mentioned enumeration technique is preferred. The general drawback of the dynamic programming approach is the curse of dimensions. Since the value function W depends on the state x and the control history v L (k) of length L, the approach is computationally feasible merely for low state dimensions and moderate values of L and M , where M is the number of elements in V.

Remark 2. Please note that a more efficient coding of the switching history is possible. To this end, it would be sufficient to replace v L (k) in W by v k-1 and an additional state with values in {0, . . . , L}, which encodes how long the last switch is ago (with 0 meaning longer than L).

Numerical results

We consider a car of mass m driving on a road with a given slope profile γ(•) and aim to optimize the gear shift control on a given time horizon [0, T ]. The car model is adapted from [START_REF] Gerdts | Solving Mixed-Integer Optimal Control Problems by Branch&Bound: A Case Study from Automobile Test-Driving with Gear Shift[END_REF] with the state (x, v) (position and velocity) and controls φ ∈ U := [0, 1] (gas pedal position) and µ ∈ V := {1, . . . , 5} (gear shift). The optimal control problem in discrete time with step-size h = T /N , N ∈ N, reads as follows:

Minimize -αx(t N ) + h N -1 k=0 f 0 (v(t k ), φ(t k ), µ(t k ))
subject to the constraints x(0) = 1470, v(0) = 23, and for k = 0, 1, . . . , N -1,

x(t k+1 ) = x(t k ) + hv(t k ), v(t k+1 ) = v(t k ) + h m i g (µ(t k ))i t M mot (φ(t k ), v(t k ), µ(t k )) R -F R (v(t k )) cos(γ(x(t k ))) -F L (v(t k )) -mg sin(γ(x(t k ))) , µ(t k ) ∈ V, φ(t k ) ∈ U, v(t k ) ∈ [22, 30].
Herein, α ≥ 0 is a weight factor, i g is the gear transmission coefficient, i t the motor torque transmission, R the wheel radius, F R the friction force of the rolling wheel, F L (v) = 1 2 c w ρAv 2 the drag, and

f 0 (v, φ, µ) = β 1 ω mot (v, µ) + β 2 |M mot (φ, v, µ)| + β 3 ω mot (v, µ) |M mot (φ, v, µ)|
models the fuel costs, where ω mot denotes the rotary frequency of the motor and M mot the motor torque, for details on ω mot , M mot , and F R please refer to [START_REF] Gerdts | Solving Mixed-Integer Optimal Control Problems by Branch&Bound: A Case Study from Automobile Test-Driving with Gear Shift[END_REF].

Figure 1 shows the results of the dynamic programming algorithm for both, the problem with dwell time constraint (with L = 3 and ū = φ = 0, i.e. φ is set to zero, if µ switches) and the problem without dwell time constraint. The parameters α = 1/10, β = 1/120, β 2 = 10/258, β 3 = 1/10320, c w = 0.3, ρ = 1.249512, A = 2, R = 0.302, i t = 3.91, g = 9.81, m = 2000, T = 10, N = 25, N x = (100, 45), N u = 100, i g (1) = 3.91, i g (2) = 2.002, i g (3) = 1.33, i g (4) = 1, i g (5) = 0.805 were used. The optimal value for the problem with dwell time constraint is -79.7656, while the optimal value without dwell time constraint is -88.0430. In the presence of the dwell time constraint it can be nicely seen that φ is set to zero after µ switches at approximately t = 4. The CPU time on a PC with 2.3 GHz is 1 minute 47.146 seconds without dwell time constraint and 14 minutes 38.9 seconds with dwell time constraint.

Conclusions

The paper discusses a dynamic programming principle that allows to consider dwell times in a time discrete optimal control problem with discrete-valued con- trols. Dwell times have to be taken into account in many applications where a switch in a discrete-valued control fixes the dynamics of the process for a specified time period before the process can be controlled again. A particular example is the optimization of gear shifts in a car, where a shift in gears implies that the motor torque is reduced for a given time period. A dynamic programming approach is suggested to solve such problems and a numerical example shows the applicability of the method. However, the curse of dimension associated to dynamic programming approaches limits the method to low dimensional systems with few states. Hence, it is desirable to develop further methods for this problem class and extend it to optimal control problems in continuous time. To this end, a formulation of the dwell time constraint in a continuous time setting becomes necessary, which is subject to future research.

Fig. 1 .

 1 Fig. 1. Optimal gear sequence (top left), gas pedal position (top right), velocity (bottom left), and slope profile (bottom right).


Martin Burger
email: martin.burger@wwu.de

Konstantinos Papafitsoros
email: kostas.papafitsoros@wias-berlin.de

Evangelos Papoutsellis
email: evangelos.papoutsellis@univ-orleans.fr

Carola-Bibiane Schönlieb

Infimal Convolution Regularisation Functionals of BV and L p Spaces. The Case p = ∞

Keywords: Total Variation, Infimal Convolution, L ∞ norm, Denoising, Staircasing

In this paper we analyse an infimal convolution type regularisation functional called TVL ∞ , based on the total variation (TV) and the L ∞ norm of the gradient. The functional belongs to a more general family of TVL p functionals (1 < p ≤ ∞) introduced in [5]. There, the case 1 < p < ∞ is examined while here we focus on the p = ∞ case. We show via analytical and numerical results that the minimisation of the TVL ∞ functional promotes piecewise affine structures in the reconstructed images similar to the state of the art total generalised variation (TGV) but improving upon preservation of hat-like structures. We also propose a spatially adapted version of our model that produces results comparable to TGV and allows space for further improvement.

Introduction

In the variational setting for imaging, given image data f ∈ L s (Ω), Ω ⊆ R 2 , ones aims to reconstruct an image u by minimising a functional of the type

min u∈X 1 s f -T u s L s (Ω) + Ψ(u), (1.1)
over a suitable function space X. Here T denotes a linear and bounded operator that encodes the transformation or degradation that the original image has gone through. Random noise is also usually present in the degraded image, the statistics of which determine the norm in the first term of (1.1), the fidelity term.

The presence of Ψ, the regulariser, makes the minimisation (1.1) a well-posed problem and its choice is crucial for the overall quality of the reconstruction. A classical regulariser in imaging is the total variation functional weighted with a parameter α > 0, αTV [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], where TV(u) := sup

Ω u divφ dx : φ ∈ C 1 c (Ω, R 2 ), φ ∞ ≤ 1 . (1.2)
While TV is able to preserve edges in the reconstructed image, it also promotes piecewise constant structures leading to undesirable staircasing artefacts. Several regularisers that incorporate higher order derivatives have been introduced in order to resolve this issue. The most prominent one has been the second order total generalised variation (TGV) [START_REF] Bredies | Total generalized variation[END_REF] which can be interpreted as a special type of infimal convolution of first and second order derivatives. Its definition reads TGV 2 α,β (u) := min

w∈BD(Ω) α Du -w M + β Ew M . (1.3)
Here α, β are positive parameters, BD(Ω) is the space of functions of bounded deformation, Ew is the distributional symmetrised gradient and • M denotes the Radon norm of a finite Radon measure, i.e.,

µ M = sup Ω φ dµ : φ ∈ C ∞ c (Ω, R 2 ), φ ∞ ≤ 1 .
TGV regularisation typically produces piecewise smooth reconstructions eliminating the staircasing effect. A plausible question is whether results of similar quality can be achieved using simpler, first order regularisers. For instance, it is known that Huber TV can reduce the staircasing up to an extent [START_REF] Huber | Robust regression: Asymptotics, conjectures and monte carlo[END_REF].

In [START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF], a family of first order infimal convolution type regularisation functionals is introduced, that reads TVL p α,β (u) := min

w∈L p (Ω) α Du -w M + β w L p (Ω) , (1.4)
where 1 < p ≤ ∞. While in [START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF], basic properties of (1.4) are shown for the general case 1 < p ≤ ∞, see Proposition 1, the main focus remains on the finite p case. There, the TVL p regulariser is successfully applied to image denoising and decomposition, reducing significantly the staircasing effect and producing piecewise smooth results that are very similar to the solutions obtained by TGV. Exact solutions of the L 2 fidelity denoising problem are also computed there for simple one dimensional data.

Contribution of the Present Work

The purpose of the present paper is to examine more thoroughly the case p = ∞, i.e., TVL ∞ α,β (u) := min

w∈L ∞ (Ω) α Du -w M + β w L ∞ (Ω) , (1.5)
and the use of the TVL ∞ functional in L 2 fidelity denoising min

u 1 2 f -u 2 L 2 (Ω) + TVL ∞ α,β (u). (1.6)
We study thoroughly the one dimensional version of (1.6), by computing exact solutions for data f a piecewise constant and a piecewise affine step function.

We show that the solutions are piecewise affine and we depict some similarities and differences to TGV solutions. The functional TVL ∞ α,β is further tested for Gaussian denoising. We show that TVL ∞ α,β , unlike TGV, is able to recover hatlike structures, a property that is already present in the TVL p α,β regulariser for large values of p, see [START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF], and it is enhanced here. After explaining some limitations of our model, we propose an extension where the parameter β is spatially varying, i.e., β = β(x), and discuss a rule for selecting its values. The resulting denoised images are comparable to the TGV reconstructions and indeed the model has the potential to produce much better results.

2 Properties of the TVL ∞ α,β Functional

The following properties of the TVL ∞ α,β functional are shown in [START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF]. We refer the reader to [START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF][START_REF] Papoutsellis | First-order gradient regularisation methods for image restoration. Reconstruction of tomographic images with thin structures and denoising piecewise affine images[END_REF] for the corresponding proofs and to [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF] for an introduction to the space of functions of bounded variation BV(Ω).

Proposition 1 ([5]) Let α, β > 0, d ≥ 1, let Ω ⊆ R d be
an open, connected domain with Lipschitz boundary and define for u

∈ L 1 (Ω) TVL ∞ α,β (u) := min w∈L ∞ (Ω) α Du -w M + β w L ∞ (Ω) . (2.1)
Then we have the following:

(i) TVL ∞ α,β (u) < ∞ if and only if u ∈ BV(Ω). (ii) If u ∈ BV(Ω) then the minimum in (2.1) is attained. (iii) TVL ∞ α,β (u
) can equivalently be defined as

TVL ∞ α,β (u) = sup Ω u divφ dx : φ ∈ C 1 c (Ω, R 2 ), φ ∞ ≤ α, φ L 1 (Ω) ≤ β ,
and TVL ∞ α,β is lower semicontinuous w.r.t. the strong L 1 convergence. (iv) There exist constants

0 < C 1 < C 2 < ∞ such that C 2 TV(u) ≤ TVL ∞ α,β (u) ≤ C 1 TV(u), for all u ∈ BV(Ω). (v) If f ∈ L 2 (Ω), then the minimisation problem min u∈BV(Ω) 1 s f -u 2 L 2 (Ω) + TVL ∞ α,β (u),
has a unique solution.

3 The One Dimensional L 2 -TVL ∞ α,β Denoising Problem

In order to get an intuition about the underlying regularising mechanism of the TVL ∞ α,β regulariser, we study here the one dimensional L 2 denoising problem min

u∈BV(Ω) w∈L ∞ (Ω) 1 2 f -u 2 L 2 (Ω) + α Du -w M + β w L ∞ (Ω) . (3.1)
In particular, we present exact solutions for simple data functions. In order to do so, we use the following theorem:

Theorem 2 (Optimality conditions) Let f ∈ L 2 (Ω). A pair (u, w) ∈ BV(Ω)× L ∞ (Ω) is a solution of (3.1) if and only if there exists a unique φ ∈ H 1 0 (Ω) such that φ = u -f, (3.2)
φ ∈ αSgn(Du -w), (3.3)
φ ∈ ψ ∈ L 1 (Ω) : ψ L 1 (Ω) ≤ β , if w = 0, ψ ∈ L 1 (Ω) : ψ, w = β w L ∞ (Ω) , ψ L 1 (Ω) ≤ β , if w = 0. (3.4)
Recall that for a finite Radon measure µ, Sgn(µ) is defined as

Sgn(µ) = φ ∈ L ∞ (Ω) ∩ L ∞ (Ω, µ) : φ L ∞ (Ω) ≤ 1, φ = dµ d|µ| , |µ| -a.e. . As it is shown in Sgn(µ) ∩ C 0 (Ω) = ∂ • M (µ) ∩ C 0 (Ω).
Proof. The proof of Theorem 2 is based on Fenchel-Rockafellar duality theory and follows closely the corresponding proof of the finite p case. We thus omit it and we refer the reader to [START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF][START_REF] Papoutsellis | First-order gradient regularisation methods for image restoration. Reconstruction of tomographic images with thin structures and denoising piecewise affine images[END_REF] for further details, see also [START_REF] Bredies | Properties of L 1 -TGV 2 : The onedimensional case[END_REF][START_REF] Papafitsoros | A study of the one dimensional total generalised variation regularisation problem[END_REF] for the analogue optimality conditions for the one dimensional L 1 -TGV and L 2 -TGV problems.

The following proposition states that the solution u of (3.1) is essentially piecewise affine.

Proposition 3 (Affine structures) Let (u, w) be an optimal solution pair for (3.1) and φ be the corresponding dual function.

Then |w| = w L ∞ (Ω) a.e. in the set {φ = 0}. Moreover, |u | = w L ∞ (Ω) whenever u > f or u < f .
Proof. Suppose that there exists a U ⊆ {φ = 0} of positive measure such that |w(x)| < w L ∞ (Ω) for every x ∈ U . Then

Ω φw dx ≤ Ω\U |φ||w| dx + U |φ||w| dx < w L ∞ (Ω) Ω\U |φ| dx + U |φ| dx = w L ∞ (Ω) φ L 1 (Ω) = β w L ∞ (Ω) ,
where we used the fact that φ L 1 (Ω) ≤ β from (3.4). However this contradicts the fact that φ, w = β w L ∞ (Ω) also from (3.4). Note also that from (3.2) we have that {u > f } ∪ {u < f } ⊆ {φ = 0} up to null sets. Thus, the last statement of the proposition follows from the fact that whenever u > f or u < f then u = w there. This last fact can be shown exactly as in the corresponding TGV problems, see [START_REF] Bredies | Properties of L 1 -TGV 2 : The onedimensional case[END_REF]Prop. 4.2].

Piecewise affinity is typically a characteristic of higher order regularisation models, e.g. TGV. Indeed, as the next proposition shows, TGV and TVL ∞ regularisation coincide in some simple special cases.

Proposition 4

The one dimensional functionals TGV 2 α,β and TVL ∞ α,2β coincide in the class of those BV functions u, for which an optimal w in both definitions of TGV and TVL ∞ is odd and monotone.

M + 2β w L ∞ (Ω) = TVL ∞ α,2β (u).
Exact solutions: We present exact solutions for the minimisation problem (3.1), for two simple functions f, g : (-L, L) → R as data, where f (x) = hX (0,L) (x) and g(x) = f (x) + λx, with λ, h > 0. Here X C (x) = 1 for x ∈ C and 0 otherwise. With the help of the optimality conditions (3.2)-(3.4) we are able to compute all possible solutions of (3.1) for data f and g and for all values of α and β. These solutions are depicted in Figure 1. Every coloured region corresponds to a different type of solution. Note that there are regions where the solutions coincide with the corresponding solutions of TV minimisation, see the blue and red regions in Figure 1a and the blue, purple and red regions in Figure 1b. This is not surprising since as it is shown in [START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF] for all dimensions, TVL ∞ α,β = αTV, whenever β/α ≥ |Ω| 1/q with 1/p + 1/q = 1 and p ∈ (1, ∞]. Notice also the presence of affine structures in all solutions, predicted by Proposition 3. For demonstration purposes, we present the computation of the exact solution that corresponds to the yellow region of Figure 1b and refer to [START_REF] Papoutsellis | First-order gradient regularisation methods for image restoration. Reconstruction of tomographic images with thin structures and denoising piecewise affine images[END_REF] for the rest.

Since we require a piecewise affine solution, from symmetry and (3.2) we have that φ

(x) = (c 1 -λ) x 2 2 -c 2 |x| + c 3 .
Since we require u to have a discontinuity at 0, (3.3) implies φ(0) = a while from the fact that φ ∈ H 1 0 (Ω) and from (3.4) we must have φ(-L) = 0 and φ, w = β w L ∞ (Ω) . These conditions give We also have c 1 = u = w and thus we require c 1 > 0. Since we have a jump at x = 0 we also require g(0) < u(0) < h, i.e., 0 < c 2 < h 2 and u(-L) > g(-L) i.e., φ (-L) > 0. These last inequalities are translated to the following conditions

c 1 = 6(αL -β) L 3 + λ, c 2 = 4αL -3β L 2 , c 3 = α.
β < αL + λL 3 6 , β > 4αL 3 - hL 2 6 , β > 2αL 3 , β < 4αL 3 ,
which define the yellow area in Figure 1b. We can easily compute u now:

u(x) = 6(αL-β) L 3 + λ x + h -4αL-3β L 2 , x ∈ (0, L), 6(αL-β) L 3 + λ x + 4αL-3β L 2 , x ∈ (-L, 0).
Observe that when β = αL, apart from the discontinuity, we can also recover the slope of the data g = λ, something that neither TV nor TGV regularisation can give for this example, see [8, Section 5.2].

Numerical Experiments

In this section we present our numerical experiments for the discretised version of L 2 -TVL ∞ α,β denoising. We solve (3.1) using the split Bregman algorithm, see [START_REF] Papoutsellis | First-order gradient regularisation methods for image restoration. Reconstruction of tomographic images with thin structures and denoising piecewise affine images[END_REF]Chapter 4] for more details.

First we present some one dimensional examples that verify numerically our analytical results. Figure 2a shows the TVL ∞ result for the function g(x) = hX (0,L) (x)+λx where α, β belong to the yellow region of Figure 1b. Note that the numerical and the analytical results coincide. We have also computed the TGV solution where the parameters are selected so that f -u TGV 2 = f -u TVL ∞ 2 . Figure 2b shows a numerical verification of Proposition 4. There, the TVL ∞ parameters are α = 2 and β = 4, while the TGV ones are α = 2 and β = 2. Both solutions coincide since they satisfy the symmetry properties of the proposition.

We proceed now to two dimensional examples, starting from Figure 3. There we used a synthetic image corrupted by additive Gaussian noise of σ = 0.01, cf. Figures 3a-3b. We observe that TVL ∞ denoises the image in a staircasing-free way in Figure 3c. In order to do so however, one has to use large values of α and β something that leads to a loss of contrast. This can easily be treated by solving the Bregman iteration version of L 2 -TVL ∞ α,β minimisation, that is

u k+1 = argmin u,w 1 2 f -v k -u 2 2 + α ∇u -w 1 + β w ∞ , v k+1 = v k + f -u k+1 . (4.1)
Bregman iteration has been widely used to deal with the loss of contrast in these type of regularisation methods, see [START_REF] Benning | Higher-order TV methods -Enhancement via Bregman iteration[END_REF][START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF] among others. For fair comparison we also employ the Bregman iteration version of TV and TGV regularisations.

The Bregman iteration version of TVL ∞ regularisation produces visually a very similar result to the Bregman iteration version of TGV, even though it has a slightly smaller SSIM value, cf. Figures 3g-3h. However, TVL ∞ is able to reconstruct better the sharp spike at the middle of the figure, cf. Figures 3d, 3e and3j.

The reason for being able to obtain good reconstruction results with this particular example is due to the fact that the modulus of the gradient is essentially constant apart from the jump points. This is favoured by the TVL ∞ regularisation which promotes constant gradients as it is proved rigorously in dimension one in Proposition 3. We expect that a similar analytic result holds in higher dimensions and we leave that for future work. However, this is restrictive when gradients of different magnitude exist, see Figure 4. There we see that in order to get a staircasing-free result with TVL ∞ we also get a loss of geometric information, Figure 4b, as the model tries to fit an image with constant gradient, see the middle row profiles in Figure 4c. While improved results can be achieved with the Bregman iteration version, Figure 4d, the result is not yet fully satisfactory as an affine staircasing is now present in the image, Figure 4e.

Spatially adapted TVL ∞ : One way to allow for different gradient values in the reconstruction, or in other words allow the variable w to take different values, is to treat β as a spatially varying parameter, i.e., β = β(x). This leads to the spatially adapted version of TVL ∞ : TVL

∞ s.a. (u) = min w∈L ∞ (Ω) α Du -w M + βw L ∞ (Ω) . (4.2)
The idea is to choose β large in areas where the gradient is expected to be small and vice versa, see Figure 5 for a simple illustration. In this example the Fig. 6: Best reconstruction of the "Ladybug" in terms of SSIM using TV, TVG and spatially adapted TVL ∞ regularisation. The β for the latter is computed both from the filtered (Figure 6g) and the ground truth image (Figure 6h) slope inside the inner square is roughly twice the slope outside. We can achieve a perfect reconstruction inside the square by setting β = 3500, with artefacts outside, see Figure 5a and a perfect reconstruction outside by setting the value of β twice as large, i.e., β = 7000, Figure 5b. In that case, artefacts appear inside the square. By setting a spatially varying β with a ratio β out /β in 2 and using the Bregman iteration version, we achieve an almost perfect result, visually very similar to the TGV one, Figures 5c-5d. This example suggests that ideally β should be inversely proportional to the gradient of the ground truth. Since this information is not available in practice we use a pre-filtered version of the noisy image and we set

β(x) = c |∇f σ (x)| + . (4.3)
Here c is a positive constant to be tuned, > 0 is a small parameter and f σ denotes a smoothing of the data f with a Gaussian kernel. We have applied the spatially adapted TVL ∞ (non-Bregman) with the rule (4.3) in the natural image "Ladybug" in Figure 6. There we pre-smooth the noisy data with a discrete Gaussian kernel of σ = 2, Figure 6d, and then apply TVL ∞ s.a. with the rule (4.3), Figure 6g. Comparing to the best TGV result in Figure 6f, the SSIM value is slightly smaller but there is a better recovery of the image details (objective). Let us note here that we do not claim that our rule for choosing β is the optimal one. For demonstration purposes, we show a reconstruction where we have computed β using the gradient of ground truth u g.t. , Figure 6c, as β(x) = c/(|∇u g.t. (x)|+ ), with excellent results, Figure 6h. This is of course impractical, since the gradient of the ground truth is typically not available but it shows that there is plenty of room for improvement regarding the choice of β. One could also think of reconstruction tasks where a good quality version of the gradient of the image is available, along with a noisy version of the image itself. Since the purpose of the present paper is to demonstrate the capabilities of the TVL ∞ regulariser, we leave that for future work.

Fig. 1 :

 1 Fig. 1: Exact solutions for the L 2 -TVL ∞ α,β one dimensional denoising problem Proof. Note first that for every odd and monotone bounded function w we have Dw M = 2 w L ∞ (Ω) and denote this set of functions by A ⊆ BV(Ω). For a BV function u as in the statement of the proposition we have TGV 2 α,β (u) = argmin w∈A

4 Fig. 2 :

 42 Fig. 2: One dimensional numerical examples

Fig. 3 :

 3 Fig. 3: TVL ∞ based denoising and comparison with the corresponding TV and TGV results. All the parameters have been optimised for best SSIM.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Illustration of the fact that TVL ∞ regularisation favours gradients of fixed modulus

  Gradient of the smoothed data, σ = 2, 13x13 pixels window (e) TV: α = 0.06, SSIM = 0.8608 (f) TGV: α = 0.068, β = 0.046, SSIM = 0.8874 (g) TVL ∞ s.a. : α = 0.07 and β computed from filtered version with c = 30, ε = 10 -4 , SSIM = 0.8729 (h) TVL ∞ s.a. : α = 0.5 and β computed from ground truth with c = 50, ε = 10 -4 , SSIM = 0.9300

Acknowledgments. The authors acknowledge support of the Royal Society International Exchange Award Nr. IE110314. This work is further supported by the King Abdullah University for Science and Technology (KAUST) Award No. KUK-I1-007-43, the EPSRC first grant Nr. EP/J009539/1 and the EPSRC grant Nr. EP/M00483X/1. MB acknowledges further support by ERC via Grant EU FP 7-ERC Consolidator Grant 615216 LifeInverse. KP acknowledges the financial support of EPSRC and the Alexander von Humboldt Foundation while in UK and Germany respectively. EP acknowledges support by Jesus College, Cambridge and Embiricos Trust.


Octavian Catrina

Eugen Borcoci

Piotr Krawiec

Two-Phase Multi-criteria Server Selection for Lightweight Video Distribution Systems

Keywords: Content networks, Video streaming, Server selection, Multicriteria decision algorithms

Video streaming services need a server selection algorithm that allocates efficiently network and server resources. Solving this optimization problem requires information about current resources. A video streaming system that relies entirely on the service provider for this task needs an expensive monitoring infrastructure. In this paper, we consider a two-phase approach that reduces the monitoring requirements by involving the clients in the selection process: the provider recommends several servers based on limited information about the system's resources, and the clients make the final decision, using information obtained by interacting with these servers. We implemented these selection methods in a simulator and compared their performance. The results show that the two-phase selection is effective, improving substantially the performance of lightweight service providers, with limited monitoring capabilities.

Introduction

An essential task of a video streaming service provider is to select a suitable content server. This task can be formalized as a multi-criteria optimization problem, that takes into account various static and dynamic attributes of the system's components. The goal is to allocate efficiently the network and server resources, while providing a suitable Quality of Experience (QoE) to the service users.

These problems are NP-complete, in general, but practical server selection algorithms are available. As these algorithms need information about network and server resources, a service provider has to deploy a complex monitoring infrastructure in order to use them. Lightweight service providers, without monitoring capabilities, can only make suboptimal decisions.

In this paper, we consider a more flexible approach, that splits the selection process between the service provider and the clients. The provider recommends a short list of content servers, using limited information (e.g., the lengths of the paths between servers and clients). The client refines the selection using additional information, obtained by interacting directly with these servers. The provider and the clients can use various combinations of algorithms, corresponding to different capabilities to obtain information about the system, and different tradeoffs between performance and complexity. In particular, this approach can improve the performance for lightweight service providers.

The paper has the following structure. Section 2 introduces a generic video streaming model and Section 3 describes the server selection methods. We implemented these methods in a simulator, presented in Section 4. The results of the simulations are discussed in Section 5, followed by a summary in Section 6.

Video Streaming Framework

We introduce in the following the generic model of online video streaming used in this paper. Figure 1 shows the main actors, the interactions between them, and the infrastructure used to deliver the service. We separate the management of the service, assigned to a service provider (SP) entity, from the delivery of the video content, performed by a collection of content servers (CS). An essential function of the SP is to select a content server when a client initiates a session.

We assume over-the-top delivery, using video streaming over HTTP or a similar technology [?]. The video content is stored on a collection of servers deployed in the Internet, each video file being available from several servers. These servers usually operate in a Content Delivery Network (CDN) [?]. We expect efficient CDN operation, with server placement and content replication adapted to the video streaming traffic. The video streaming system can include a CDN, or use the services of one or more CDN providers [?].

A video streaming session starts with the following interactions (Figure 1a):

-The client sends a request to the SP, specifying the desired video content.

-The SP runs a server selection algorithm and replies with a list of servers that could deliver the requested video content (or rejects the request). Table 1 lists the information used by the SP in the server selection process. -In a single-phase selection process, the client simply asks one of the servers indicated by the provider to deliver the desired video content.

Table 1. Information used by service providers for content server selection.

Static or quasi-static information Dynamic information Lists of content servers and video files.

Available capacity of the content servers. Mapping of video content to servers and Available capacity of the paths between of clients and servers to network nodes. servers and clients. Length of server-client network paths.

Admission control for servers and paths.

-In a two-phase selection process, the client chooses a server from the provider's list by running its own selection algorithm, with additional information.

For example, if the service is implemented using the DASH standard [?], the provider can answer with a Media Presentation Description (MPD) file, containing information about the encoding of the video content (including available resolutions and data rates) and the URLs of the servers that can deliver it.

We assume that the main bottlenecks of the transport network are the links between the autonomous systems (AS) of the network service providers. Therefore, the network is modeled as a collection of interconnected nodes, corresponding to an AS-level, multi-tier network topology, similar to the Internet topology (Figure 1b). The clients and the content servers are attached to these network nodes. The Border Gateway Protocol (BGP), responsible for inter-domain routing in the Internet, determines a single best path between ASes, based on network service provider policies and path length. Load balancing on multiple paths would help to allocate more efficiently the network resources to the video streaming sessions, and could be achieved using application-level overlays. We consider, therefore, this additional capability for some of the server selection algorithms.

Content Server Selection

Multi-Criteria Decision Algorithm (MCDA). Selection of a server and a path in content delivery systems can be formulated as a multi-criteria optimization problem. These problems are NP-complete in general, and are solved in practice using heuristics methods. Practical and efficient methods have been proposed for server and path selection based on Multi-Criteria Decision Algorithms (MCDA) [?]. The variant of MCDA used in this paper is described below:

-The algorithm determines a set of candidate solutions, S = {S i } i∈1,...,n , from information about the request and the delivery system. A candidate solution is a vector of decision variables, S i = (v i,j ) j∈1,...,m , representing characteristics of a server and a path that could deliver the requested content. -For each candidate solution S i and variable v i,j , the algorithm calculates the component achievement function R i,j = ri,j -vi,j ri,j -ai,j , where v i,j is the value of the variable, while r i,j and a i,j are the reservation level and aspiration level for this variable, respectively (the algorithm can be extended with more complex functions). We use two decision variables: the available (unused) capacity of the server and the available capacity (unused bandwidth) of the path. For these variables, the reservation level is a lower bound for suitable solutions, and the aspiration level is an upper bound beyond which the preference of the solutions is similar. We set a i,j to the maximum capacity of the server and of the path, respectively, and r i,j to the amount of server capacity and path capacity consumed by a session, respectively.

-The algorithm calculates the rank of each solution, R i = min j {R i,j }, and then determines the index of the best solution, argmin i {R i }.

Single-phase selection. The single-phase selection relies entirely on the service provider to find suitable servers for the clients' sessions. In its simplest variant, the provider answers a client's request by indicating a single server [?]. We consider a more general variant, in which the provider answers with an ordered list of servers. The client can start the session by connecting to one of these servers and use the others as back-up [?]. By delivering a list of servers, the provider enables the client to improve the initial server selection and/or use certain adaptive streaming techniques [?]). The selection proceeds as follows:

-The service provider maintains a resource database with the information about network and content servers that is necessary for server selection. The database contains static information and some of the dynamic information listed in Table 1, depending on the algorithm being used. -When a service request arrives, the provider determines an ordered list of suitable servers and delivers it to the client. The selection is based on information received in the service request (e.g., client and video ids) and information available in the resource database. If the provider does not find any suitable server the request is rejected.

Candidate solutions consist of a server that has a copy of the requested video file and a path from that server to the client. The selection algorithm takes as input a set of candidate solutions that (optionally) satisfy additional requirements. The provider can use a range of server selection algorithms, that offer different tradeoffs between the complexity of the system and the optimality of the solution (e.g., different requirements for the information collected by the provider about network and servers). The algorithms used in our simulation experiments are listed in Table 2. -If the request is accepted, the client starts the video session by connecting to the first server in the list. In case of failure, the client can try the next servers, without having to ask the provider to recommend another server.

This approach simplifies the functionality of the client. On the other hand, the service provider has, essentially, two main options: algorithms that use only static information and offer poor performance and algorithms that can achieve much better performance, but require an expensive monitoring infrastructure.

Two-phase selection. We consider now a more flexible, two-phase approach, that involves the clients in the selection process:

-The provider's resource database contains only the information about network and servers that is necessary for a preliminary selection: static and quasi-static information (Table 1) and, possibly, dynamic information that is easier to collect (e.g., server load). -When a request arrives, the service provider makes a preliminary server selection and delivers to the client a short, ordered list of recommended servers.

The procedure is essentially the same as for single-phase selection, using algorithms that match the limited information collected by the provider: e.g., random servers or closest servers.

-If the request is accepted, the client performs the second phase of the selection, choosing the server that will be used by the video streaming session.

The input of the client's selection algorithm is the list of servers received from the provider and additional information obtained by the client: 1. Local information: Client capabilities, processor and link load.

2. Current state of the recommended servers: Reports from the servers, indicating if they can handle the request and/or the available capacity. 3. Current state of the paths between the servers and the client: Throughput of the connection with each server. The client can use MCDA based on available server and path capacity, or choose the least loaded server, depending on the information it can obtain.

Admission control. The video streaming sessions have to ensure continuous playout. This requires timely data delivery, and hence a certain lower bound for the end-to-end data rate. An important issue, therefore, is to add some form of admission control to the server selection process. This can be achieved by removing the solutions with fully loaded servers or paths from the set of candidate solutions given as input to the selection algorithms.

Admission control is simpler for content servers. The available capacity can be estimated by a local monitoring application and reported by the servers to the provider or the clients. Moreover, a fully loaded server can reject additional requests. This offers basic admission control, independent of the selection process. However, it is more difficult to determine if the path from a server to a client can handle additional sessions. The service providers need for this purpose a network monitoring infrastructure that measures the throughput between network nodes. The clients can measure the throughput of the connections with the recommended servers. These end-to-end measurements are more relevant, but increase the delay and overhead of session establishment. Also, they may fail to detect initial congestion and avoid the QoE degradation of current sessions.

Simulation Software

We analyzed the performance of the server selection methods described in Section 3 by simulation. The simulation software developed for this purpose is written in C++ and uses the discrete event simulator OMNeT++ [?] as basic simulation engine. Figure 2 shows the main components of the simulation software. The class ResourceManager provides a resource database with information about the entire system infrastructure: network topology (nodes and links, link capacity and load), network paths, list of content servers (server capacity and load), mapping of video content to servers, mapping of clients and servers to network nodes. The simulations are driven by the RequestProcess, which generates the service requests, manages the video sessions, and collects statistics. When a session starts or ends, the RequestProcess notifies the Re-sourceManager, which updates accordingly the load of the server and the load of the links on the path from the server to the client.

The server selection process is implemented by the classes ServiceProvider and ServiceClients, and obtains information about system resources from the Re-sourceManager. The class ServiceProvider handles the clients' service requests and implements the server selection algorithms run by the service provider for single-phase and two-phase selection (Table 2). The class ServiceClients implements the server selection algorithms run by the clients in two-phase selection. The network model uses random topologies generated using the tool aSHIIP [?] and the network nodes are the autonomous systems (AS) of the network service providers (Figure 1b). The topologies have a multi-tier structure similar to the Internet. We assume that this structure corresponds to the usual business relations between network service providers: peering relations between nodes in the same tier, and customer-provider relations between nodes in different tiers, with the customer in the lower tier. The paths used in the simulations are similar to the BGP paths, according to the relations between network providers [?].

The clients and the content servers are attached to network nodes. Servers are deployed all over the network, in well-connected nodes, close to clients, taking into account the number of links of the node and their bandwidth.

The resources consumed by the video streaming sessions are measured using a metric for server capacity and a metric for link capacity. To simplify the analysis, we assume that a session consumes 1 unit of server capacity and 1 unit of link capacity. The system is provisioned by assigning a fixed maximum capacity to servers and links, which is available for video streaming. The servers have the same capacity and are deployed in clusters; better connected nodes have larger clusters. The network is provisioned by assigning capacity to links according to its hierarchical structure, with higher capacity in upper tiers.

The main performance metric used in the analysis of the server selection process is the success rate, defined as the ratio between the number of successful service requests and the total number of service requests. A request is successful if the system has sufficient server and network resources to deliver the requested video content to the client, with the desired quality, for its entire duration.

Failures can occur in the server selection phase or during the video streaming session. If admission control is available, a request is rejected when the selection process does not find a server and a path with sufficient unused capacity; the load of the system remains unchanged. Otherwise, a new session is created, which may overload the server and/or some links on the path from server to client. The RequestProcess monitors the current sessions and updates their state. If a session has been successful so far and does not have sufficient resources any more, its state changes from success to fail.

Performance Evaluation

We discuss in this section the results of simulations with single-phase and twophase selection, for different selection algorithms. The experiments use a network with 1000 nodes, structured into 4 tiers, and (up to) 3 paths for any pair of nodes. The network is provisioned so that the number of concurrent sessions can reach the total capacity of the servers, 61500 sessions, for a suitable server selection. Service requests are generated at random time intervals with exponential distribution. We measure the session success rate as a function of the request rate.

Analysis of single-phase selection. Figure 3 shows the success rate for singlephase selection and the configurations in Table 3. Service providers with network monitoring capability can use MCDA and admission control for servers and paths. The success rate of MCDA is equal to 1 for request rates up to 150 requests/sec and about 60000 concurrent sessions, close to the total server capacity (the mean session duration is 400 seconds). The performance drops drastically if the providers have only static information and select random servers or closest servers, without admission control. The degradation is due to sessions rejected by fully loaded servers and, especially, congestion of network links.

Congested sessions may be aborted due to poor QoE. The freed resources can then be used to start other sessions, which may be successful (e.g., their paths do not include congested links). Figure 3 shows the simulations results when the congested sessions are aborted and when they are allowed to continue. Important performance differences appear for randomly chosen servers and/or paths, because the allocation of network resources is less efficient. Intuitively, the performance should improve when the load is distributed on multiple network paths. We observe, however, that using the shortest path is better than the random choice of a path (e.g., RandShrt versus RandRand). In a multi-tier network with valley-free paths, the shortest path may avoid the core of the network, while the other paths are more likely to cross it (e.g., the paths between SC2 and CG1 in Figure 1). Thus, the random choice of the path concentrates the traffic in the network core. Therefore, if multiple paths are available, the choice should take into account the path load (e.g., using MCDA).

For similar reasons, the performance improves when the provider selects the servers with the shortest paths to the client (NearSrv), instead of choosing them randomly (e.g., RandRand, RandShrt): the sessions use less network resources and the paths are more likely to avoid the core. However, the performance of NearSrv selection is very sensitive to server placement.

Analysis of two-phase selection. In the two-phase selection process, the provider and the clients can use various combinations of algorithms. Our main goal is to see the effects of this approach on the performance of the video streaming services offered by lightweight providers. We assume, therefore, that the provider has global, static information about the entire system, and the clients can obtain dynamic information for a particular request and several candidate servers. The two-phase selection clearly improves the success rate with respect to single-phase selection when the provider chooses the recommended servers using NearSrv or RandShrt and the clients use MCDA. NearSrv can achieve better performance than RandShrt, for a suitable placement of the servers in the network, but the performance varies significantly depending on server placement.

However, the two-phase selection cannot reach the performance of the singlephase selection using MCDA (with admission control): the provider runs MCDA with the entire set of candidate solutions, while the client has a short list of servers, chosen without information about system load. The size of the list is a trade-off between system performance and the delay and overhead of the interactions between client and servers, that collect additional information for the final selection. The results suggest that a short list of 3 servers is sufficient.

The performance improves significantly when the clients add admission control to MCDA, becoming close to the performance of a provider that is able to run MCDA itself. The improvement is more important when the provider uses RandShrt, which is more vulnerable to network congestion.

Conclusions

A service provider that offers over-the-top video streaming needs an expensive monitoring infrastructure for efficient server selection. Lightweight providers, without such capabilities, can only use simple selection methods, based on static information. Simulation results show a large performance gap between these two approaches. We explored, therefore, a two-phase selection process, that aims at improving the performance for lightweight providers, by involving the clients.

We studied the performance of two-phase selection by simulation. The results show that the approach is effective, improving substantially the success rate. In the ideal case, when the clients can reliably apply admission control, the performance of the two-phase selection becomes close to the performance of a provider that runs MCDA using a network monitoring infrastructure.

Future work will focus on implementing the two-phase selection, using, e.g., DASH [?]. The critical issue is to enable the clients to obtain the information needed for the final selection with minimum delay and overhead, and sufficient accuracy. Ideally, this information could be collected during the initial buffering phase. The DASH standard offers several features that could be exploited for this purpose, e.g., chunks of video content can be downloaded from several servers.

Fig. 1 .

 1 Fig. 1. Video streaming model: (a) Actors and interactions. (b) Infrastructure.

Fig. 2 .

 2 Fig. 2. High-level architecture of the simulation software.

Fig. 3 .

 3 Fig. 3. Single-phase selection: (a) with session abort; (b) without session abort.

Fig. 4 .

 4 Fig. 4. Two-phase versus single-phase selection for: (a) NearSrv; (b) RandShrt.

Figure 4

 4 Figure 4 shows the simulation results for the following configurations: -TpNearSrvMcda/TpRandShrtMcda: Provider: NearSrv/RandShrt. Client: MCDA for available server and path capacity. Congested sessions are aborted. -TpNearSrvMcdaClac/TpRandShrtMcdaClac: Similar to the previous configrations, with client admission control for servers and paths. -The provider sends a list of 3 recommended servers, out of 8 candidates. -The results for single-phase selection using NearSrv, RandShrt, and MCDA are added to Figure 4 to facilitate the performance comparison.

Table 2 .

 2 Server selection algorithms used by the service provider.

	Algorithm Input	Output (list of servers).
	Random	Set of feasible solutions.	Randomly chosen servers; shortest or
	servers		random path.
	Closest	Set of feasible solutions. Length of	Servers with the shortest paths to the
	servers	the paths from servers to client.	client.
	MCDA	Set of feasible solutions. Available ca-	Servers in the solutions with the
		pacity of the servers and the paths.	highest MCDA rank.

Table 3 .

 3 Configurations used for server selection by the service provider.

	Name	Server and path selection algorithm	Admission control (AC)
	RandRand Random servers and paths	Server reject
	RandShrt Random servers, shortest paths	Server reject
	NearSrv	Servers with shortest paths to client	Server reject
	MCDA	MCDA for available server and path capacity AC for servers and paths

Acknowledgments. This work was partially supported by the Research Project DISEDAN, No.3-CHIST-ERA C3N, 2014-2015.


Domingo Alberto Tarzia
email: dtarzia@austral.edu.ar

DOUBLE CONVERGENCE OF A FAMILY OF DISCRETE DISTRIBUTED MIXED ELLIPTIC OPTIMAL CONTROL PROBLEMS WITH A PARAMETER

Keywords: Double convergence, Distributed optimal control problems, Elliptic variational equalities, Mixed boundary conditions, Numerical analysis, Finite element method, Fixed points, Optimality conditions, Error estimations

The convergence of a family of continuous distributed mixed elliptic optimal control problems ( P α ), governed by elliptic variational equalities, when the parameter α → ∞ was studied in Gariboldi -Tarzia, Appl. Math. Optim., 47 (2003), 213-230 and it has been proved that it is convergent to a distributed mixed elliptic optimal control problem ( P ). We consider the discrete approxi- mations ( h P α ) and ( h P ) of the optimal control problems ( P α ) and ( P ) respec- tively, for each > 0 h

and α > 0 . We study the convergence of the discrete distributed optimal control problems ( h P α ) and ( h P ) when → 0 h

, α → ∞ and α → +∞ ( , ) (0, ) h obtaining a complete commutative diagram, including the diagonal convergence, which relates the continuous and discrete distributed mixed elliptic optimal control problems ( ) ( ) ( ) , , h h P P P α α

and ( P ) by taking the corresponding limits. The convergent corresponds to the optimal control, and the system and adjoint system states in adequate functional spaces.

Introduction

The purpose of this paper is to do the numerical analysis, by using the finite element method, of the convergence of the continuous distributed mixed optimal control problems with respect to a parameter (the heat transfer coefficient) given in [START_REF] Gariboldi | Convergence of distributed optimal controls on the internal energy in mixed elliptic problems when the heat transfer coefficient goes to infinity[END_REF][START_REF] Gariboldi | A new proof of the convergence of the distributed optimal controls on the internal energy in mixed elliptic problems[END_REF] obtaining a double convergence when the parameter of the finite element method goes to zero and the heat transfer coefficient goes to infinity.

We consider a bounded domain n Ω ⊂ R whose regular boundary Γ = ∂Ω = Γ ∪Γ consists of the union of two disjoint portions 1 Γ and 2 Γ with meas 1 ( ) 0 Γ > . We consider the following elliptic partial differential problems with mixed boundary conditions, given by:

1 2 in ; on ; on u u g u b q n ∂ -Δ = Ω = Γ - = Γ ∂ , ( 1
)
1 2 in ; ( ) on ; on u u u g u b q n n α ∂ ∂ -Δ = Ω - = - Γ - = Γ ∂ ∂ ( 2
)
where g is the internal energy in Ω , . 0 b Const = > is the temperature on 1 Γ for the system (1) and the temperature of the external neighborhood on 1 Γ for the system (2)

respectively, q is the heat flux on 2 Γ and α > 0 is the heat transfer coefficient on 1 Γ . The systems ( 1) and ( 2) can represent the steady-state two-phase Stefan problem for adequate data [START_REF] Tabacman | Sufficient and/or necessary condition for the heat transfer coefficient on 1 Γ and the heat flux on 2 Γ to obtain a steady-state two-phase Stefan problem[END_REF][START_REF] Tarzia | An inequality for the constant heat flux to obtain a steady-state two-phase Stefan problem[END_REF]. We consider the following continuous distributed optimal control problem ( ) P and a family of continuous distributed optimal control problems ( ) P α for each parameter α > 0 , defined in [10], where the control variable is the internal energy g in Ω, that is: Find the continuous distributed optimal controls 2 ( ) J g J g P J g J g

op g H L ∈ = Ω and
α α α α ∈ ∈ = = (3)
where the quadratic cost functional 0 , : [2, 18, 26]:

J J H α + → R are defined by
2 2 2 2 1 1 ) ( ) , ) ( ) 2 2 2 2 g d g d H H H H M M a J g u z g b J g u z g α α = - + = - + (4)
with 0 M > and d z H

∈ given, g u K ∈ and g u V α ∈ are the state of the systems de- fined by the mixed ellliptic differential problems (1) and ( 2) respectively whose elliptic variational equality are given by [START_REF] Kinderlehrer | An introduction to variational inequalities and their applications[END_REF]:

( ) ( ) 2 0 : , , , g g u K a u v g v qvd v V γ Γ ∈ = - ∀∈ ∫ (5) ( ) ( ) 2 1
: , , ,

g g u V a u v g v qvd bvd v V γ α γ Γ Γ ∈ = - + ∀∈ ∫ ∫ (6)
and their adjoint system states g p V ∈ and g p V α ∈ are defined by the following elliptic variational equalities:

( ) ( ) ( ) ( ) 0 ) : , , , ; ) : , , , g o g g d g g g d a p V a p v u z v v V b p V a p v u z v v V α α α α ∈ = - ∀∈ ∈ = - ∀∈ (7)
with the spaces and bilinear forms defined by:

{ } 1 2 2 0 1 0 2 ( ), , / 0 , , ( ), ( ) V H V v V v K b V H L Q L = Ω = ∈ Γ = = + = Ω = Γ (8) ( ) ( ) ( ) 1 , . , , , , ( , ) a u v u vdx a u v a u v uvd u v uvdx α α γ Ω Γ Ω = ∇ ∇ = + = ∫ ∫ ∫ (9)
where the bilinear, continuous and symmetric forms a and α a are coercive on 0 V and V respectively, that is [16]:

2 0 0 such that ( , ), V v a v v v V λ λ ∃ > ≤ ∀ ∈ (10) 2 1 min(1, ) 0 such that ( , ), V v a v v v V α α α λ λ α λ ∃ = > ≤ ∀ ∈ (11)
and λ > 1 0 is the coercive constant for the bilinear form 1 a [16, 21].

The unique continuous distributed optimal energies op g and op g α have been characterized in [START_REF] Gariboldi | Convergence of distributed optimal controls on the internal energy in mixed elliptic problems when the heat transfer coefficient goes to infinity[END_REF] as a fixed point on H for a suitable operators W and W α over their optimal adjoint system states

0 op g p V ∈ and op g p V α α
∈ defined by: ( ) ( )

1 1 , : such that ) , ) g g W W H H a W g p b W g p M M α α α → = - = - . ( 12
)
The limit of the optimal control problem ( P α ) when α → ∞ was studied in [START_REF] Gariboldi | Convergence of distributed optimal controls on the internal energy in mixed elliptic problems when the heat transfer coefficient goes to infinity[END_REF] and it was proven that: lim 0, lim 0, lim 0

op op op op op op g g g g H V V u u p p g g α α α α α α α α →∞ →∞ →∞ - = - = - = (13)
for a large constant > 0 M by using the characterization of the optimal controls as fixed points through operators (12a) and (12b); this restrictive hypothesis on data was eliminated in [START_REF] Gariboldi | A new proof of the convergence of the distributed optimal controls on the internal energy in mixed elliptic problems[END_REF] by using the variational formulations. We can summary the conditions [START_REF] Hintermüller | Moreau-Yosida regularization in state constrained ellliptic control problems: Error estimates and parameter adjustement[END_REF] saying that the distributed optimal control problems ( P α ) converges to the distributed optimal control problem (P) when α → +∞ . Now, we consider the finite element method and a polygonal domain n Ω ⊂ R with a regular triangulation with Lagrange triangles of type 1, constituted by affineequivalent finite element of class 0 C being h the parameter of the finite element approximation which goes to zero [START_REF] Brenner | The mathematical theory of finite element methods[END_REF][START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. Then, we discretize the elliptic variational equalities for the system states ( 6) and ( 5), the adjoint system states (7a) and (7b), and the cost functional (4a,b) respectively. In general, the solution of a mixed elliptic boundary problem belongs to

( ) r H Ω with 3 1 ( 0 ) 2 r ε ε < ≤ - >
but there exist some examples which solutions belong to

( ) r H Ω with 2 r
≤ [START_REF] Azzam | On solutions of elliptic equations satisfying mixed boundary conditions[END_REF][START_REF] Lanzani | The mixed problem in p L for some twodimensional Lipschitz domain[END_REF][START_REF] Shamir | Regularization of mixed second order elliptic problems[END_REF]. Note that mixed boundary conditions play an important role in various applications, e.g. heat conduction and electric potential problems [START_REF] Haller-Dintelmann | Hölder continuity and optimal control for nonsmooth elliptic problems[END_REF].

The goal of this paper is to study the numerical analysis, by using the finite element method, of the convergence results [START_REF] Hintermüller | Moreau-Yosida regularization in state constrained ellliptic control problems: Error estimates and parameter adjustement[END_REF] corresponding to the continuous distrib-uted elliptic optimal control problems ( ) and ( ) P P α when α → +∞ . The main result of this paper can be characterized by the following result: Theorem 1

We have the following complete commutative diagram which relates the continuous distributed mixed optimal control problems ( ) P α and ( ) P , with the discrete dis- tributed mixed optimal control problems ( ) are respectively the optimal control, the system and the adjoint system states of the discrete distributed mixed optimal control problem ( ) h P α for each 0 h > and 0 α > , and the double convergence is the diagonal one. The study of the limit 0 h → of the discrete solutions of optimal control problems can be considered as a classical limit, see [4-6, 8, 9, 13-15, 19, 23, 24, 27, 28] but the limit α → +∞ , for each 0 h > , and the double limit α → +∞ ( , ) (0, ) h can be considered as a new ones.

The paper is organized as follows. In Section II we define the discrete elliptic variational equalities for the state systems hg u and α h g u , we define the discrete dis- tributed cost functional h J and α h J , we define the discrete distributed optimal control problems ( ) h P and ( ) h P α , and the discrete elliptic variational equalities for the ad- joint state systems hg p and α h g p for each 0 h > and 0 α > , and we obtain properties for the discrete optimal control problems ( ) h P and ( )

h P α .
In Section III we study the classical convergences of the discrete distributed optimal control problems ( ) h P to ( ) P , and ( ) h P α to ( ) P α when 0 h → (for each 0 α > ) and the estimations for the discrete cost functional h J and α h J . In Section IV we study the new convergence of the discrete distributed optimal control problems ( ) and ( ) P by taking the limits 0 h → and α → +∞ . In Section V we study the new double conver- gence of the discrete distributed optimal control problems ( )

h P α to ( ) P when α → +∞ ( , ) (0, ) h
and we obtain the diagonal convergence in the previous commutative diagram.

Discretization by Finite Element Method and Properties

We consider the finite element method and a polygonal domain n Ω ⊂ R with a regular triangulation with Lagrange triangles of type 1, constituted by affineequivalent finite element of class 0 C being h the parameter of the finite element approximation which goes to zero [START_REF] Brenner | The mathematical theory of finite element methods[END_REF][START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. We can take h equal to the longest side of the triangles h T τ ∈ and we can approximate the sets 0 , and V V K by: ( )

( ) { } { } 0 1 0 1 0 / , , / 0; h h h h h h h h h h V v C v T P T T V v V v K b V τ = ∈ Ω ∈ ∀ ∈ = ∈ Γ = = + (14)
where 1 P is the set of the polymonials of degree less than or equal to 1. Let 0 : ( )

h h C V π Ω →
be the corresponding linear interpolation operator. Then there exists a constant 0 0 c > (independent of the parameter h ) such that [3]:

( ) ( ) ( ) 1 0 0 ) ; ) ; , 1 2 r r r h h r r H V a v v c h v b v v c h v v H r π π - - ≤ - ≤ ∀ ∈ Ω < ≤ . ( 15
)
We define the discrete cost functional 0 , :

h h J J H α + → R by the following expressions: ( ) ( ) 2 2 2 2 1 1 ) , ) 2 2 2 2 h h g d h h g d H H H H M M a J g u z g b J g u z g α α = - + = - + (16)
where hg u and h g u α are the discrete system states defined as the solution of the fol- lowing discrete elliptic variational equalities [START_REF] Kinderlehrer | An introduction to variational inequalities and their applications[END_REF][START_REF] Tarzia | Numerical analysis of a mixed elliptic problem with flux and convective boundary conditions to obtain a discrete solution of non-constant sign[END_REF]:

( ) ( ) 2 0 : , , , hg h hg h h h h h u K a u v g v qv d v V γ Γ ∈ = - ∀ ∈ ∫ , (17)
( ) ( ) 2 1
: , , ,

h g h h g h h h h h h u V a u v g v qv d bv d v V γ α γ Γ Γ ∈ = - + ∀ ∈ ∫ ∫ . ( 18
)
The corresponding discrete distributed optimal control problems consists in finding

, op op h h g g H α ∈ such that: ( ) ( ) ( ) ( ) ) Problem ( ) :
, ) Problem ( ) :

op op h h h h h h h h g H g H a P J g M i nJ g b P J g M i nJ g α α α α ∈ ∈ = = (19)
and their corresponding discrete adjoint states hg p and h g p α are defined respectively as the solution of the following discrete elliptic variational equalities:

( ) ( ) 0 0 : , , , hg h hg h hg d h h h p V a p v u z v v V ∈ = - ∀ ∈ (20) ( ) ( ) : , , , h g h h g h h g d h h h p V a p v u z v v V α α α α ∈ = - ∀ ∈ (21)
Remark 1. We note that the discrete (in the n-dimensional space) distributed optimal control problem ( h P ) and ( h P α ) are still infinite dimensional control problems since the control space is not discretized.

Lemma 2 (i) There exist unique solutions

hg h u K ∈ and 0 hg h p V ∈ , and h g h u V α ∈ and h g h p V α ∈
of the elliptic variational equalities ( 17) and ( 20), [START_REF] Lions | Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles[END_REF], and ( 21) respec-

tively g H ∀ ∈ , q Q ∀ ∈ , 1 0 b on > Γ.
(ii) The operators Proof. We use the Lax-Milgram Theorem, the variational equalities ( 17), ( 18), ( 20) and ( 21), the coerciveness [START_REF] Gariboldi | Convergence of distributed optimal controls on the internal energy in mixed elliptic problems when the heat transfer coefficient goes to infinity[END_REF] and [START_REF] Gariboldi | A new proof of the convergence of the distributed optimal controls on the internal energy in mixed elliptic problems[END_REF] and following [START_REF] Gariboldi | Convergence of distributed optimal controls on the internal energy in mixed elliptic problems when the heat transfer coefficient goes to infinity[END_REF][START_REF] Lions | Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles[END_REF][START_REF] Tarzia | A commutative diagram among discrete and continuous boundary optimal control problems[END_REF].

hg g H u V ∈ → ∈ ,

W Theorem 3 (i) The discrete cost functional h

J and h J α are H -elliptic and strictly convexe applications, that is

[ ] ( ) 1, 2 , 0,1 g g H t ∀ ∈ ∀ ∈ : ( ) ( ) ( ) ( ) ( ) ( ) 2 2 1 1 2 2 1 1 1- 1 2 h h h H t t t J g tJ g J tg t g M g g - + - + - ≥ - (22)
( ) ( ) ( ) ( ) ( ) ( ) 2 2 1 1 2 2 1 1 1- 1 2 h h h H t t t J g tJ g J tg t g M g g α α α - + - + - ≥ - (23)
(ii) There exist a unique optimal controls (iii) h J and h J α are Gâteaux differenciable applications and their derivatives are given by the following expressions:

( ) ( ) ) , ) , , 0 h h g h h g a J g Mg p b J g Mg p g H h α α ʹ ʹ = + = + ∀ ∈ ∀ > (24)
(iv) The optimality condition for the optimization problems (19a) and (19b) are given by: ( ) ( )

1 1 a) 0 ; b) 0 op op h op op h op op h h h h g h h h h g J g g p J g g p M M α α α α α ʹ ʹ = ⇔ = - = ⇔ = - (25)
(v) h J ʹ and h J α ʹ are Lipschitzians and strictly monotone operators.

Proof. We use the definitions (16a,b), the elliptic variational equalities ( 17) and ( 18) and the coerciveness ( 10) and [START_REF] Gariboldi | A new proof of the convergence of the distributed optimal controls on the internal energy in mixed elliptic problems[END_REF], following [START_REF] Gariboldi | Convergence of distributed optimal controls on the internal energy in mixed elliptic problems when the heat transfer coefficient goes to infinity[END_REF][START_REF] Lions | Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles[END_REF][START_REF] Tarzia | A commutative diagram among discrete and continuous boundary optimal control problems[END_REF].

W

We define the operators:

( ) ( ) 1 1 , : such that a) , ) h h h h g h h g W W H H W g p b W g p M M α α α → = - = - . ( 26
)

Theorem 4

We have that: (i) h W and h W α are Lipschitzian operators, and h

W ( h W α ) is a contraction opera-
tor if and only if M is large, that is:

2 2 1 1 a) M > , ) M > b α λ λ . ( 27
)

(ii) If M verifies the inequalities (27a,b) then the discrete distributional optimal control

op h g H ∈ ( op h g H α ∈
) is obtained as the unique fixed point of h W ( h W α ), i.e.:

( ) ( )

1 1 , op h op op op h op op op op h h g h h h h h g h h h g p W g g g p W g g M M α α α α α α = - ⇔ = = - ⇔ = . ( 28
)
Proof. We use the definitions (26a,b), and the properties (25a,b) and Lemma 2. , ,

op op h h op op op op r r r h h g g h g g H V V g g ch u u ch p p ch - - - - ≤ - ≤ - ≤ (29) 1 1 1
, ,

op op h h op op op op r r r h h g g h g g H V V g g ch u u ch p p ch α α α α α α α α α α - - - - ≤ - ≤ - ≤ ( 30
)
where c's are constants independents of h.

Proof. It is useful to use the restriction 1 α > by splitting a α by [21, 24,25] ( ) ( )

1 1 , , ( 1 )
a u v a u v uvd α α γ Γ = + -∫ (31)
but then it can be replaced by 0 α α ≥ for any 0 0 α > . We follow a similar method to the one developed in [START_REF] Tarzia | A commutative diagram among discrete and continuous boundary optimal control problems[END_REF] for Neumann boundary optimal control problems by using the elliptic variational equalities ( 17), ( 18), ( 20) and ( 21), the thesis holds.

( ) ( ) ( ) ( ) 2 2 2( 1) 2 ( 1)
, 2 2

h o p o p o p h o p o p o p op op r r h h H H M M g g J g J g Ch g g J g J g Ch α α α α α α - - - ≤ - ≤ - ≤ - ≤ (32) ( ) ( ) ( ) ( ) 2 2 2( 1) 2( 1)
; 2 2

h o p o p h h o p o p h op op op op r r h h h h H H M M g g J g J g Ch g g J g J g Ch α α α α α α - - - ≤ - ≤ - ≤ - ≤ (33) ( ) ( ) ( ) ( ) 1 1
,

op op op op r r h h h J g J g Ch J g J g Ch - - - ≤ - ≤ (34) ( ) ( ) ( ) ( ) 1 1
,

op op op op r r h h h J g J g Ch J g J g Ch α α α α α α - - - ≤ - ≤ (35)
where C's are constants independents of h and α .

Proof. Estimations (32) and (33) follow from the estimations (29), and the equalities (similar relationship for and J J α ):

( ) ( )

- = - = - = ∀ > . ( 39
)
Proof. We omit this proof because we prefer to prove the next one with more details.

Double Convergence of the Discrete Distributed Optimal

Control Problem ( ) h P α to ( )

P when α → +∞ ( , ) (0, ) h
For the discrete distributed optimal control problem ( ) h P α we will now con- sider the double limit α → +∞ ( , ) (0, ) h

.

Theorem 10

We have the following limits:

( , ) (0, ) ( , ) (0, ) ( , ) (0, ) lim lim lim 0

h h o p o p op op op op h g g h g g h h h h H V V u u p p g g α α α α α α α α → +∞ → +∞ → +∞ - = - = - = (40)
Proof. From now on we consider that c's represent positive constants independents simultaneously of 0 h > and 0 α > (see (31)). We show a sketch of the proof by obtaining the following estimations (for

0 h ∀ > and 1 α ∀ > ):
( ) ( )

1 2 0 1 0 2 0 3 , , 1
h h h V V u c u c u b d c α α α γ Γ ≤ ≤ - - ≤ ∫ (41) 4 5 6
, ,

op h op op h h g h H H H g c u c g c α α α ≤ ≤ ≤ (42) ( ) ( ) 1 2 7 8 9 , , 1
h h h op op op hg h g h g V V u c u c u b d c α α α α α γ Γ ≤ ≤ - - ≤ ∫ (43) ( ) 1 2 10 11 12 , , 1
h h h op op op hg h g h g V V p c p c p d c α α α α α γ Γ ≤ ≤ - ≤ ∫ . ( 44
)
For example, the constant 11 c is a positive constant independent simultaneously of 0 h > and 0 α > , and it is given by the following expression:

11 1 1 1 1 1 1 1 1 0 2 2 1 d H Q c z M M b M M M q M λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ γ λ λ λ λ λλ λ λ λ λ ⎡ ⎤ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ = + + + + + + ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ ⎡ ⎤ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎛ ⎞ + + + + + + + ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ + + + + + + + + ⎜ ⎟ ⎜ ⎟ ⎜ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ 2 1 1 1 1 1 1 M λ λ λ ⎡ ⎤ ⎡ ⎤ ⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎢ ⎥ + + + ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎟ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ ⎣ ⎦ (45)
Therefore, from the above estimations we have that:

/ in weak as ( , ) (0, ) 30), (39), (40) and [START_REF] Gariboldi | Convergence of distributed optimal controls on the internal energy in mixed elliptic problems when the heat transfer coefficient goes to infinity[END_REF][START_REF] Gariboldi | A new proof of the convergence of the distributed optimal controls on the internal energy in mixed elliptic problems[END_REF]. Remark 3. We note that this double convergence is a novelty with respect to the recent results obtained for a family of discrete Neumann boundary optimal control problems [START_REF] Tarzia | A commutative diagram among discrete and continuous boundary optimal control problems[END_REF].

op h f H g f H h α α ∃ ∈ ⎯⎯ → → +∞ (46) 1 / i n w e a k ( s t r o n g ) a s (, ) ( 0 , ) w i t h / h op h g V u V H h b α α η η α η ∃ ∈ ⎯⎯ → → +∞ Γ = (47) 1 / i n w e a k ( s t r o n g ) a s (, ) ( 0 , ) w i t h / 0 h op h g V p V H h α α ξ ξ α ξ ∃ ∈ ⎯⎯ → → +∞ Γ = (48) / i

  for each α > 0 ) such that:

Fig. 1 .h

 1 Fig. 1. Relationship among optimal control problems ( ) ( ) ( ) , , h h P P P α α

  h P α to ( ) h P when α → +∞ for each 0 > h and we obtain a commutative diagram which relates the continuous and discrete distributed mixed optimal control problems ( ) ( ) ( )

  satisfy the optimization problems (19a) and (19b) respectively.

Theorem 6

 6 We obtain the following error estimations between the continuous and discrete solutions: We suppose the continuous system states and adjoint system states have the regularities If M verifies the inequalities (27a,b) then we have the following error bonds:

W Remark 2 .Lemma 7

 27 If M verifies the inequalities (27a,b) we can obtain the convergence in Theorem 6 by using the characterization of the fixed point (28a,b), and the uniqueness of the optimal controls op give some estimations for the discrete cost functional h J α and h J . If M verifies the inequality (27a,b) and the continuous system states and adjoint system states have the regularities we have the following error bonds:

.Proof of Theorem 1 .

 1 and therefore the three double limits (40) hold when α → +∞ ( , ) (0, ) h It is a consequence of the properties (29), (

  and

					g H u ∈ →	h g α	V ∈ are Lipschitzians.
	The operators	g H ∈ →	hg p	0 V ∈ , and g	g H ∈ →	h g p α	h V ∈ are Lipschitzians and
	strictly monotone operators.

α α α

Acknowledgements

The present work has been partially sponsored by the Projects PIP No 0534 from CONICET -Univ. Austral, Rosario, Argentina, and AFOSR-SOARD Grant FA9550-14-1-0122.